/* file: mymnewton.mac e. woollett, jan.,09 Maxima ver. 5.17.1 */ disp("working version mymnewton, assumes two dimensional problem only, syntax: mymnewton( exprlist, guesslist, numiter )$ exprlist should have the form: [expr1, expr2], to find (x,y) such that simultaneously expr1=0,expr2=0, expr1, expr2 should be explicit functions of x and y, guesslist should be in the form: [xguess,yguess], numiter = number of iterations ")$ mymnewton( exprlist, guesslist, numiter) := block([numer,ratprint,fpprintprec, small : 1.0e-30 , g,x,y,gv,h,hv,b,v,ls,d ] , local(f,j), numer:true, ratprint:false, fpprintprec:8, /* g = col vec: row 1=expr 1, row 2=expr 2 depends on (x,y) */ g : matrix( [ exprlist[1] ],[ exprlist[2] ] ), display(g), gv : ev(g, x=v[1,1], y=v[2,1] ), /* v is generic col vector */ define( f(v), gv ), /* h is jacobian matrix associated with col vec g(x,y) */ h : matrix( [diff( g[1,1], x), diff( g[1,1], y)], [diff( g[2,1], x), diff(g[2,1], y) ] ), hv : ev(h, x=v[1,1], y=v[2,1] ), define( j(v), hv ), /* b is col vec containing (x,y) values */ b : matrix([ guesslist[1] ],[ guesslist[2] ]), /* start iterations */ for i:0 thru numiter do ( print(" "), print(" i = ",i," b = ",b," f(b) = ",f(b) ), if ( i = 0) then print(" starting values ") else print(" condition number = ", second(ls) ), /* check jacobian determinant */ d : ( determinant( j(b) ), float(%%) , abs(%%) ), if ( d < small ) then ( print(" abs(det(jacobian)) is ", d," program halt " ), return() ), /* using 'floatfield arg gets condition number returned */ ls : linsolve_by_lu(j(b),-f(b), 'floatfield ), /* improved (hopefully) estimates of values of (x,y) which simultaneously satisfy expr1=0 and expr2=0 */ b : b + first(ls) ) /* end do */ )$ /* end block */