I have not looked at your code, but, in case you are interested, here
is a simple 2d routine I wrote. Looking at it now, it is easy to
extend it to higher dimension.
HTH,
-sen
/* File newton2d.mac
This is a 2-d version of the Newton root finding algorithm.
expr is a 2d list of functions of a vector var.
x0 is an initial 2d vector guess
eps is the desired error.
Example: newton2d([x^2 - y^2 + sin(y), x + exp(y)],[x,y],
[2,1],1e-19);
Note: For high degree in x, y, this produces bind stack overflows.
*/
newton2d(expr,var,x0,eps):=
block([u,v,s,numer],
numer:true,
s:matrix([diff(expr[1],var[1]),diff(expr[1],var[2])],
[diff(expr[2],var[1]),diff(expr[2],var[2])]),
u: matrix(x0),
loop, if ( subst(u[1][2],var[2],subst(u[1][1],var[1],expr[1]) )^2
+ subst(u[1][2],var[2],subst(u[1][1],var[1],expr[2]))^2 < eps)
then return([u[1][1],u[1][2]]),
v: subst(u[1][2],var[2],subst(u[1][1],var[1],expr)),
u: u - v
.
subst(u[1][2],var[2],subst(u[1][1],var[1],invert(transpose(s)))),
go(loop) )$
---------------------------------------------------------------------------
| Sheldon E. Newhouse | e-mail: sen1 at math.msu.edu |
| Mathematics Department | |
| Michigan State University | telephone: 517-355-9684 |
| E. Lansing, MI 48824-1027 USA | FAX: 517-432-1562 |
---------------------------------------------------------------------------
On Thu, 5 Apr 2007, Wolfgang Lindner wrote:
> dear maxima experts,
>
> making my first steps into Maxima I want to translate A. Pinkus YaCAS-Code for 2-dim
> Newton-Iteration. Here is the original YaCAS-Code from the file Algo.Book.pdf S.18 (3.9
> Newton iteration) of the YaCAS 1.063 distribution:
>
>>> f(a,x):={Sin(a*x),a-2}
>>> matrix(a,x):=Eval({D(a)f(a,x),D(x)f(a,x)})
>>> {a,x}:={1.5,1.5}
>>> For(ii:=1,ii<100,ii++)[{a,x}:={a,x}+ N(SolveMatrix(matrix(a,x),-f(a,x)));]
>>> {a,x}
> / \
> | 2. |
> | |
> | 0.596673098e-1 |
> \ /
>
> And here is my beginner's translation to Maxima (yes, I know of mnewton.mac and
> newton1.mac) searching for a very simple version of newton. But I get stuck.
>
> numer : true;
> f(a,x) := matrix([Sin(a*x)],[a-2]);
> J(a,x) : =matrix(''diff(f(a,x),a),''diff(f(a,x),x));
>
> /* What is wrong with the following line? */
> [a,x]: [1.5,1.5];
>
> /*
> How to manage the following line?
> [a,x] is type list I think, the product is type matrix.
> So I have to coerce it someway - but how?
> */
> for i:1 thru 100 do ( [a,x] : [a,x] + J(a,x)^^(-1).(-f(a,x)) );
> [a,x];
>
> Sorry if my question is too simple, but it would be nice if someone could help me.
> The wxMaxima sheet ist attached.
>
> Thanks Wolfgang
>
>