Your code has no protection against infinite looping.
I think the stack overflow happens when the newton
sequence doesn't converge.
Another thing: Instead of invert(transpose(s)), I suggest you
use linsolve_by_lu:
(%i9) m : matrix([0.001,1],[1,2])$
(%i10) b : matrix([1],[2])$
(%i11) linsolve_by_lu(m,b,'floatfield);
(%o11) [matrix([0.0],[1.0]),9.021036072144288]
The condition number is about 9.02. The function linsolve_by_lu
is new to 5.12.0. With 5.11.0, you'll need something like
(%i11) lu_backsub(lu_factor(m,'floatfield),b);
(%o11) matrix([0.0],[1.0])
Barton
-----maxima-bounces at math.utexas.edu wrote: -----
>To: Richard Fateman
>From: sen1 at math.msu.edu
>Sent by: maxima-bounces at math.utexas.edu
>Date: 04/06/2007 06:47PM
>cc: Maxima at math.utexas.edu
>Subject: Re: [Maxima] bind stack overflow
>
>OK, here it is. This is just a simple newton method for two
>dimensional systems. I have not tried to study it carefully, but I
>was curious about the 'bind stack overflow in some cases, but not
>others'
>
>Here is the code
>
>/***************************/
>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) )$
>/*****************************/