Maxima/Numerical methods

"If you are doing purely numerical computation and are concerned about speed, use a compiled numerical programming language. Maxima is intended for use if you have symbolic mathematical symbols, too, and while it works for numbers, most components of the system are on the lookout for non-numeric inputs, and this checking takes time.   It is possible to speed up certain kinds of numeric computations in Maxima by using compile and mode_declare  together. " RJF

Newton method
Functions
 * newton for equation with function of one variable
 * mnewton is an implementation of Newton's method for solving nonlinear equations in one or more variables.

newton
load : load("newton"); (%o1)                                                                                                  /home/a/maxima/share/numeric/newton.mac (%i2) load("newton1"); (%o2)                                                                                                  /home/a/maxima/share/numeric/newton1.mac

Code :

/* Maxima CAS code from /maxima/share/numeric/newton1.mac input : exp =  function of one variable, x var = variable x0 = initial value of variable eps =

The search begins with x = x_0 and proceeds until abs(expr) < eps (with expr evaluated at the current value of x).

output : xn = an approximate solution of expr = 0 by Newton's method

newton(exp,var,x0,eps):= block( [xn,s,numer],   numer:true,   s:diff(exp,var),   xn:x0,   loop,     if abs(subst(xn,var,exp))<eps          then return(xn),    xn:xn-subst(xn,var,exp)/subst(xn,var,s),   go(loop) )$

mnewton
One can use Newton method for solving systems of multiple nonlinear functions. It is implemented in mnewton function. See directory :

/Maxima..../share/contrib/mnewton.mac

which uses linsolve_by_lu defined in :

/share/linearalgebra/lu.lisp

See this image for more code. (%i1) load("mnewton")$

(%i2) mnewton([x1+3*log(x1)-x2^2, 2*x1^2-x1*x2-5*x1+1],             [x1, x2], [5, 5]); (%o2) x1 = 3.756834008012769, x2 = 2.779849592817897 (%i3) mnewton([2*a^a-5],[a],[1]); (%o3)            a = 1.70927556786144 (%i4) mnewton([2*3^u-v/u-5, u+2^v-4], [u, v], [2, 2]); (%o4) u = 1.066618389595407, v = 1.552564766841786