A Procedure for Newton's Method in R^n

Version 0.8

This worksheet:

1. contains a proceudre to do n iterations of Newtons method for solving f=0 where
f:R^n->R^n is given as a liast of expressions

2. optionally show the function values at each iteration.

3. optionally show an estimate |f| |D2f| |Df^(-1)|^2 used in checking the convergence

of Newton's method near this point, where the norms are sums of squares of entries.

If the quantity (replacing the second derivative term by a second derivative norm

bound in a suitable ball) is smaller than 1/2 then convergence is guaranteed.

4. will be available at the URL http://www.mathlab.cornell.edu/local_maple/newton_conv_est.mws when completed.

> restart;

> with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> with(plots):

Warning, the name changecoords has been redefined

Change this value if you'd like a different floating point accuracy. More digits are fine.

> Digits:=10;

This error tolerance is used to decide if a matrix is too close to 0 (as measured by its det) to be considered invertible. Reduce it if you like.

> tol := 10.0^(-8);

Change this to false and re-execute the line if you don't want to see function values.

> show_f_val := true;

Change this between true and false and re-execute the line depending on whether you want to see the test for f convergence for the Newton's method near this point.

> show_conv_estimate := true;

Typical Usage: newton([x^2-y^2-16, 2*x*y-30], [x,y], [4,4],10);

to do 10 iterations of Newton's method for solving

x^2-y^2-16 = 0
2*x*y-30 = 0
using the variables [x,y] and starting with the point [4,4].
The sequence of iterations are maintained as p[0], p[1], ... in a global table p.

> newton := proc(f,vars,p0,n)
global p,d,tol, show_f_val, show_radius_estimate;
local i, j, m, u, f_val, Df, Df_val, Df_detval, Df_invval,var_subs_set, err_found,q,q_pre, rad_est, H2f, H2f_val, D2f_norm;
p:= 'p';
p[0] := eval(p0):
q:= map(evalf,p0);
m:= nops(vars);
Df := jacobian(f,vars);
var_subs_set := {seq(vars[j]=q[j],j=1..m)};
f_val := map(evalf,map2(subs,var_subs_set,f));
Df_val := map(evalf,map2(subs,var_subs_set,Df));
Df_detval := det(Df_val);
err_found := evalb(abs(Df_detval) < tol);
printf(`Iteration starts with %a.\n`,convert(q,list));
if (show_f_val) then
printf(`\t\tThe function value here is %a.\n`, f_val);
fi:
if (show_conv_estimate) then
H2f:=[seq(hessian(f[i],vars),i=1..m)];
H2f_val:= map2(subs,var_subs_set,H2f):
Df_invval := inverse(Df_val);
D2f_norm := sqrt(add(trace(evalm(H2f_val[i] &* transpose(H2f_val[i]))),i=1..m)):
rad_est:= D2f_norm* sqrt(trace(evalm(f_val&*transpose(f_val))))*(trace(evalm(Df_invval&*transpose(Df_invval))))^2;
printf(`\t\tThe convergence estimate here is %a.\n`, evalf(rad_est));
fi:
printf(`\n`);
q_pre := eval(q);
for i from 1 to n while (not err_found) do
q := map(evalf,evalm(q_pre - Df_val^(-1) &* f_val));
u := evalm(q - q_pre);
d := sqrt(1.0*dotprod(u,u));
printf(`Iteration %d gives %a, a distance %g from the previous.\n`,i,convert(q,list),d);
var_subs_set := {seq(vars[j]=q[j],j=1..m)};
f_val := map(evalf,map2(subs,var_subs_set,f));
if (show_f_val) then
printf(`\t\tThe function value here is %a.\n`, f_val);
fi:
Df_val := map(evalf,map2(subs,var_subs_set,Df));
Df_detval := det(Df_val);
if (show_conv_estimate) then
H2f_val:= map2(subs,var_subs_set,H2f):
Df_invval := inverse(Df_val);
D2f_norm := add(trace(evalm(H2f_val[j] &* transpose(H2f_val[j]))),j=1..m):
rad_est:= D2f_norm* sqrt(trace(evalm(f_val&*transpose(f_val))))*(trace(evalm(Df_invval&*transpose(Df_invval))))^2;
printf(`\t\tThe convergence estimate here is %a.\n`, rad_est);
fi:
err_found := evalb(abs(Df_detval) < tol);
q_pre := eval(q);
p[i] := convert(eval(q_pre),list);
printf(`\n`);
od:
eval(q);
end:

Enter your function as a list of expressions.

> f := [x^2-y^2-16,2*x*y-30];

> newton(f,[x,y],[4,4],10);

Iteration starts with [4., 4.].

The function value here is [-16., 2.].

The convergence estimate here is .1574659717e-1.

Iteration 1 gives [4.875000000, 2.875000000], a distance 1.425219 from the previous.

The function value here is [-.50000000, -1.96875000].

The convergence estimate here is .7919095779e-2.

Iteration 2 gives [5.001402439, 3.002378049], a distance .179451 from the previous.

The function value here is [-.24759e-3, .3220180e-1].

The convergence estimate here is .1112429410e-3.

Iteration 3 gives [5.000000023, 3.000000653], a distance .002760 from the previous.

The function value here is [-.369e-5, .666e-5].

The convergence estimate here is .2634571979e-7.

Iteration 4 gives [5.000000000, 3.000000001], a distance 6.524055e-07 from the previous.

The function value here is [-.1e-7, 0.].

The convergence estimate here is .3460207611e-10.

Iteration 5 gives [5.000000001, 3.000000001], a distance 1.000000e-09 from the previous.

The function value here is [0., .2e-7].

The convergence estimate here is .6920415222e-10.

Iteration 6 gives [5.000000000, 3.000000000], a distance 1.414214e-09 from the previous.

The function value here is [0., 0.].

The convergence estimate here is 0..

Iteration 7 gives [5.000000000, 3.000000000], a distance 0 from the previous.

The function value here is [0., 0.].

The convergence estimate here is 0..

Iteration 8 gives [5.000000000, 3.000000000], a distance 0 from the previous.

The function value here is [0., 0.].

The convergence estimate here is 0..

Iteration 9 gives [5.000000000, 3.000000000], a distance 0 from the previous.

The function value here is [0., 0.].

The convergence estimate here is 0..

Iteration 10 gives [5.000000000, 3.000000000], a distance 0 from the previous.

The function value here is [0., 0.].

The convergence estimate here is 0..

> p[3];

> pts := [seq(p[i],i=0..10)];

> plot(pts,`sequence of points in the plane`,thickness=3);

>

>

>