images_newton/sm121_newton1.mws

NOTE on the syntax: As you'll see almost all commands must end in either a ":" (to suppress output) or a ";" and almost always an "=" sign is preceeded by a ":".

EXAMPLE 1: Let's begin by examining the function sin(x).

> linear_approx:=proc(f,x,a)
local f0,Df0,Df;
f0:=limit(f(x),x=a);
Df := unapply(diff(f(x),x),x);
Df0:=limit(Df(x),x=a);
RETURN(f0+Df0*(x-a));
end;

> linear_approx(sin,x,0);
#Caution: the wrong syntax:
#linear_approx(sin(x),x,0);
#crashes MAPLE!

> linear_approx(tan,x,0);

> linear_approx(cos,x,0);

> f:=x->(1-cos(x))/x;
linear_approx(f,x,0);

> f:=x->sqrt(x);
linear_approx(f,x,4);

> linear_approx(exp,x,0);

Let us plot the values of sin(x) and it's linear approximation about x=0:

> f:=sin:
print(`x`.` `.f(x),` `,linear_approx(f,x,0));
for i from -10 to 10 do
print(evalf(i/100),evalf(f(i/100)),evalf(subs(x=i/100,linear_approx(f,x,0))));
od;

We define f(x) (in MAPLE's "arrow" notation ->) to be sqrt(4+x) as follows.

> f:=x->sqrt(x);

> print(`x`.` `.f(x),` `,linear_approx(f,x,4));
for i from -10 to 10 do
print(evalf(4+i/100),evalf(f(4+i/100)),evalf(subs(x=4+i/100,linear_approx(f,x,4))));
od;

> L:=x->linear_approx(f,x,4);

> plot1 := plot({f(x),L(x)},x=0..10):
plot2 := plots[textplot]({[4,2.5,`y=1+x/4`],[1.2,1,`y=sqrt(x)`]},align = RIGHT):

> plots[display]({plot1,plot2});

EXERCISE for Example 1(optional): Do all the same commands but replace sqrt(x), a=4, by tan(x), a=Pi/4.

> newtons_method:=proc(f,x,x0,N)
#interates newton's method N steps
#starting at x0
local i,Df,f0,Df0,next;
f0:=f(x0);
Df := unapply(diff(f(x),x),x);
Df0:=Df(x0);
next:=evalf(x0-f0/Df0):
for i from 1 to N do
print(i,next);
next:=evalf(next-limit(f(x),x=next)/limit(Df(x),x=next)):
od:
RETURN(next);
end;

> plot_newtons_method:=proc(f,x,x0,N)
#interates Newton's method N steps
#starting at x0
local ff,i,Df,next,L,k,xx;
if x1=x0 then ERROR(`can't have x0=x1`); fi;
ff[0]:=f(x0);
Df := unapply(diff(f(x),x),x);
xx[0]:=x0;
xx[1]:=evalf(xx[0]-ff[0]/Df(xx[0])):
for i from 1 to N do
ff[i]:=f(xx[i]);
xx[i+1]:=evalf(xx[i]-ff[i]/Df(xx[i])):
od:
L:=[seq([xx[k],f(xx[k])],k=0..N)];
#print(L);
plots[pointplot](L,thickness=3,color=red,symbol=circle,axes=boxed);
end;

Example 2a: Let us try to find a root of f(x)=x^2-1=0 using Newton's method, where our first guess is 1.5.

We use 5 steps, then tabulate and plot the resulting points (x_n,f(x_n)) from the data.

> f:=x->x^2-1;
newtons_method(f,x,1.5,4);
plot_newtons_method(f,x,1.5,4);

Example 2b: Let us try to find a root of f(x)=x^2-1=0 using Newton's method, where our first guess is 10.

We use 4 steps, then tabulate and plot the resulting points (x_n,f(x_n)) from the data

> f:=x->x^2-1;
newtons_method(f,x,10,4);
plot_newtons_method(f,x,10,4);

Example 2c: Let us try to find a root of f(x)=x^2-1=0 using Newton's method, where our first guess is -10.

We use 4 steps, then tabulate and plot the resulting points (x_n,f(x_n)) from the data.

> f:=x->x^2-1;
newtons_method(f,x,-10,4);
plot_newtons_method(f,x,-10,4);

Example 2c: Let us try to find a root of f(x)=x^2-1=0 using Newton's method, where our first guess is 0.001, not a bad guessbut Newton's method, as you see from the data below, does notwork well in this case. We must use 12 steps, then we tabulate and plot the resulting points (x_n,f(x_n)).

> f:=x->x^2-1;
newtons_method(f,x,0.001,12);
plot_newtons_method(f,x,0.001,12);

Example 3: Let us find all the real roots of x^3+5x-3=0.

> f:=x->x^3+5*x-3;
the_roots:=solve(f(x)=0,x):
print(`first root `,evalf(the_roots[1]));
print(`second root `,evalf(the_roots[2]));
print(`third root `,evalf(the_roots[3]));
newtons_method(f,x,3,5);
plot_newtons_method(f,x,3,5);

EXERCISE for Examples 2,3(optional): Use Newton's method to find all the real solutions of x=sin(x).

> newtons_method_complex:=proc(f,z,z0,N)
#interates newton's method N steps
#starting at z0=x0+I*y0
local i,Df,f0,Df0,next;
f0:=f(z0);
Df := unapply(diff(f(z),z),z);
Df0:=Df(z0);
next:=evalc(z0-f0/Df0):
for i from 1 to N do
#print(i,next);
next:=evalf(next-limit(f(z),z=next)/limit(Df(z),z=next)):
od:
RETURN([[Re(next),Im(next)],next]);
end:

> newtons_method_complex1:=proc(f,z,z0,N)
#interates newton's method N steps
#starting at z0=x0+I*y0
local i,Df,f0,Df0,next;
f0:=f(z0);
Df := unapply(diff(f(z),z),z);
Df0:=Df(z0);
next:=evalc(z0-f0/Df0):
for i from 1 to N do
print(i,next);
next:=evalf(next-limit(f(z),z=next)/limit(Df(z),z=next)):
od:
RETURN([[Re(next),Im(next)],next]);
end:

Example 4a: Newton's method can try to find a (possibly complex) solution to z^3-1=0, with initial guess -1+I, where MAPLE uses the notation I=sqrt(-1).

> f:=z->z^3-1;
newtons_method_complex1(f,z,-1+I,5);

Example 4b: Newton's method can try to find a (possibly complex) solution to z^3-1=0, with initial guess -1-I.

> f:=z->z^3-1;
newtons_method_complex1(f,z,-1-I,5);

Example 4c: Newton's method can try to find a (possibly complex) solution to z^3-1=0, with initial guess -1.

> f:=z->z^3-1;
newtons_method_complex1(f,z,-1,7);

Computer art and the cube roots of 1.

> cuberoots_of_one:=proc(z0)
local x,y,colorname,f,z;
f:=z->z^3-1;
x:=newtons_method_complex(f,z,z0,10)[1][1]:
y:=newtons_method_complex(f,z,z0,10)[1][2]:
if (-.51<x and x <-.49 and .86<y and y<.87) then RETURN(`red`); fi;
if (-.51<x and x <-.49 and -.87<y and y<.86) then RETURN(`blue`); fi;
if (.9<x and x<1.1 and -0.1<y and y<.1) then RETURN(`green`); fi;
RETURN(`yellow`);
end;

> cuberoots_of_one_plot:=proc(N)
local m,n,P;
for m from 1 to 2*N do
for n from 1 to 2*N do
P[m,n]:=plots[pointplot]([(m-N-.1)/N,(n-N-.1)/N],color=cuberoots_of_one((m-N-.1)/N+I*(n-N-.1)/N),thickness=2):
od;
od:
print(`green points are Newton's method initial guesses for the root z=1 of z^3-1=0`);
print(`red points are Newton's method initial guesses for the root z=-1/2-I*sqrt(3)/2 of z^3-1=0`);
print(`blue points are Newton's method initial guesses for the root z=-1/2+I*sqrt(3)/2 of z^3-1=0`);
plots[display]([seq(seq(P[m,n],m=1..2*N),n=1..2*N)],axes=none);
end;

You can vaguely see that the complex plane is divided into regions which are either red, blue or green, depending on which root that point (as a first guess in Newton's method for z^3-1=0) leads to.

> cuberoots_of_one_plot(10);# about 13 seconds on a 300Mz pentium

>

> t0:=time():
cuberoots_of_one_plot(50);
time()-t0;#almost 8 minutes on a 300Mz

>

> secant_method:=proc(f,x,x0,x1,N)
#interates the secant method N steps
#starting at x0<>x1
local ff,i,Df,next;
if x1=x0 then ERROR(`can't have x0=x1`); fi;
ff[0]:=f(x0);ff[1]:=f(x1);
Df:=(ff[1]-ff[0])/(x1-x0);
x[1]:=x1;
print(0,x0);
print(1,x1);
x[2]:=evalf(x1-ff[1]/Df):
for i from 2 to N do
print(i,x[i]);
ff[i]:=f(x[i]);
if x[i-1]=x[i] then print(`stopped after step `,i); RETURN(x[i]); fi;
Df:=(ff[i]-ff[i-1])/(x[i]-x[i-1]);
x[i+1]:=evalf(x[i]-ff[i]/Df):
od:
RETURN(x[N]);
end:

> plot_secant_method:=proc(f,x,x0,x1,N)
#interates the secant method N steps
#starting at x0<>x1
local ff,i,Df,next,L,k,xx;
if x1=x0 then ERROR(`can't have x0=x1`); fi;
ff[0]:=f(x0);ff[1]:=f(x1);
Df:=(ff[1]-ff[0])/(x1-x0);
xx[0]:=x0;
xx[1]:=x1;
xx[2]:=evalf(xx[1]-ff[1]/Df):
for i from 2 to N do
ff[i]:=f(xx[i]);
if xx[i-1]=xx[i] then print(`stopped after step `,i); RETURN(); fi;
Df:=(ff[i]-ff[i-1])/(xx[i]-xx[i-1]);
xx[i+1]:=evalf(xx[i]-ff[i]/Df):
od:
L:=[seq([xx[k],f(xx[k])],k=0..N)];
#print(L);
plots[pointplot](L,thickness=3,color=red,symbol=circle);
end:

Example 5: Let us try to find a root of f(x)=x^2-1=0 using the secant method, where our first guesses are x0=-10, x1=-11.We use 10 steps, then tabulate and plot the resulting points (x_n,f(x_n)) from the data.

> f:=x->x^2-1;
secant_method(f,x,-10,-11,10);
plot_secant_method(f,x,-10,-11,10);

Example 6: Let us try to find a root of f(x)=x^2-2=0 using the secant method, where our first guesses are x0=10, x1=20.We use 10 steps, then tabulate and plot the resulting points (x_n,f(x_n)) from the data.

> f:=x->x^2-2;
secant_method(f,x,10,20,10);
plot_secant_method(f,x,10,20,10);

>

EXERCISE for Examples 5,6(optional): Use the secant method to find all the real solutions of x=sin(x)