There is something somewhat new here. On the equations with x^2 below, both
newton2d and mnewton give answers.
Either I don't know how to use mnewton or it gives the wrong answer.
On the equations with x^3, neither give answers.
On the equations with x^4, newton2d gives an answer, but mnewton does
not.
Here are the files, newton2d.mac and newton_test_output.mac so you
can easily load them and play.
Thanks for any ideas,
-sen
(Incidentally, I plan to change newton2d as suggested by Richard
Fateman and Barton Willis, but I have not yet had the time to do so.)
/********************** begin file newton2d.mac *****************/
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) )$
/***********end file newton2d.mac *******************/
/********* begin file newton_test_output.mac ***********/
build_info();
load("newton2d");
load("mnewton");
var:[x,y];
ic:[4,2];
l2:[x^2-y^2+sin(y),x+exp(y)];
l3:[x^3-y^2+sin(y),x+exp(y)];
l4:[x^4-y^2+sin(y),x+exp(y)];
r1:newton2d(l2,var,ic,9.9999999999999995E-7);
r2:mnewton(l2,var,ic);
ev(l2,x = r1[1],y = r1[2]);
prt("OK, this looks like a root");
ev(l2,x = part(r2,1,2,2),y = part(r2,1,2,2));
prt("Apparently no root");
prt("Now for the cubic and quartic terms");
newton2d(l3,var,ic,9.9999999999999995E-7);
mnewton(l3,var,ic);
newton2d(l4,var,ic,9.9999999999999995E-7);
ev(l4,x = %[1],y = %[2]);
newton2d(l4,var,ic,1.0000000000000001E-9);
ev(l4,x = %[1],y = %[2]);
mnewton(l4,var,ic);
prt("So, mnewton does not work here");
/********* end file newton_test_output.mac ***********/
----------------------------------------------------------------------------
Notice the following outputs using newton2d and mnewton.
(%i3) var:[x,y];
(%o3) [x, y]
(%i4) ic:[4,2];
(%o4) [4, 2]
(%i5) l2:[x^2-y^2+sin(y),x+exp(y)];
2 2 y
(%o5) [sin(y) - y + x , E + x]
(%i6) l3:[x^3-y^2+sin(y),x+exp(y)];
2 3 y
(%o6) [sin(y) - y + x , E + x]
(%i7) l4:[x^4-y^2+sin(y),x+exp(y)];
2 4 y
(%o7) [sin(y) - y + x , E + x]
(%i8) r1: newton2d(l2,var,ic,9.9999999999999995E-7);
(%o8) [- 0.69637915328756, - 0.36185483763238]
(%i9) r2: mnewton(l2,var,ic);
(%o9) [[x = - 0.91071428571429, y = - 0.48809523809524]]
(%i10) ev(l2,x = r1[1],y = r1[2]);
(%o10) [- 4.5613267799504698E-6, 4.295887681249333E-6]
(%i11) ev(l3,x = part(r2,1,2,2),y = part(r2,1,2,2));
(%o11) [- 0.82346368930204, 0.12569917555302]
(%i12)
-------------------------------------------------------------------------
Now for the cubic and quartic versions
Thanks for your help.
-sen
---------------------------------------------------------------------------
| 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 Sat, 7 Apr 2007, Robert Dodier wrote:
> Sheldon, I don't see any recursive function call here or other stuff
> that would be pushing stuff onto a stack, so I don't know where the
> bind stack overflow error is coming from. It could well be some
> internal error in Maxima. If you want to post an example of a
> call to newton2d which triggers the bind stack overflow error,
> I would be interested to see it.
>
> Robert
>
> PS. Here is the code from before.
>> >/***************************/
>> >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) )$
>> >/*****************************/
>