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

Boundaries of Components of Mandelbrot set by Newton method. Image and Maxima CAS code

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]]