Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Aug 2016
Lecture topic: Newton Methods

Function

f[x_] = 2x Cos[2x]-(x-2)^2
Plot[f[x],{x,0,10},AxesLabel->{"x","y"},PlotStyle->Blue,ImageSize->{500,Automatic},BaseStyle->20]

Newton's Method

We want to solve the equation f(x)=0 for x.

Algorthim

Clear[p0,p1,tol, max, i,p,q]

f[x_] =2x Cos[2x]-(x-2)^2
p0=5.0;
tol=10.0^(-6);
max=300;

data={};
i=1;
Print["Newton's Method"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[i≤max ,
p=p0-f[p0]/f'[p0];
data=Append[data,{p,f[p]}];
If[Abs[p-p0]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",Abs[p-p0], " Steps= ",i],Break[]}];
p0=p;
i=i+1]
If[i>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"p","f[p]"}}]

Error

error=Table[{n,Abs[data[[n,2]]]},{n,1,Length[data]}];
nmerrplot=ListLogPlot[error,PlotStyle->{Red,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{500,Automatic},BaseStyle->20]

Secant Method

We want to solve the equation f(x)=0 for x.

Algorthim

Clear[p0,p1,tol, max, i,p,q]

f[x_] =2x Cos[2x]-(x-2)^2
p0=5.0;
p1=4.99;
tol=10.0^(-6);
max=300;

data={};
i=2;
q0=f[p0];
q1=f[p1];
Print["Secant Method"]
Print["Input:  ","p0= ",p0," p1= ",p1," tol= ",tol]
While[i≤max ,
p=p1-q1(p1-p0)/(q1-q0);
data=Append[data,{p,f[p]}];
If[Abs[p-p1]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",Abs[p-p1], " Steps= ",i],Break[]}];
p0=p1;
q0=q1;
p1=p;
q1=f[p];
i=i+1]
If[i>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"p","f[p]"}}]

Error

error=Table[{n,Abs[data[[n,2]]]},{n,1,Length[data]}];
smerrplot=ListLogPlot[error,PlotStyle->{Blue,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{500,Automatic},BaseStyle->20]

Bisection Method

We want to solve the equation f(x)=0 for x.

Algorthim

Clear[a,b,tol,fa,fp]

f[x_] =2x Cos[2x]-(x-2)^2

a=2.0;
b=3.0;
maxit=500;
tol=10.0^(-6);

data = {};
i=1;
fa=f[a];
Print["Bisection Method"]
Print["Input:  ","a= ",a," b= ",b," tol= ",tol]
While[i≤maxit,
p=a+(b-a)/2.0;
fp=f[p];
data = Append[data,{a,b,p,fp}];
If[fp==0||(b-a)/2.0<tol,{Print["Result: "," p= ",SetPrecision[p,10]," Error ",Abs[(b-a)/2.0], " Steps= ",i],Break[]}];
i=i+1;
If[Sign[fa]*Sign[fp]>0.0,{a=p,fa=fp},{b=p}]
]
If[i>maxit,Print["Method Failed to Converge after ",maxit, " steps."]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"a","b","p","f[p]"}}]

Error

error=Table[{n,Abs[data[[n,4]]]},{n,1,Length[data]}];
bmerrplot=ListLogPlot[error,PlotStyle->{Purple,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{500,Automatic},BaseStyle->20]

Comparison of Errors

Show[bmerrplot,smerrplot,nmerrplot,PlotRange->{Automatic,Automatic}]

Notice that once the secant and Newton's methods get close to the root they converge extremely quickly. Getting good initial approximations to the root is therefore important.

Refinement: Bisection Method intializes Secant Method

We want to solve the equation f(x)=0 for x. We will use the bisection method to get a rough approximation, then the secant method to refine that approximation.

f[x_] =2x Cos[2x]-(x-2)^3
Plot[f[x],{x,0,5},PlotStyle->Blue,ImageSize->{500,Automatic}]

Refinement Algorthim

Clear[a,b,tol,fa,fp]

f[x_] =2x Cos[2x]-(x-2)^3

a=3.0;
b=4.0;
maxit=500;
tol=10.0^(-2);

data = {};
i=1;
fa=f[a];
Print["Bisection Method"]
Print["Input:  ","a= ",a," b= ",b," tol= ",tol]
While[i≤maxit,
p=a+(b-a)/2.0;
fp=f[p];
data = Append[data,{p,fp}];
If[fp==0||(b-a)/2.0<tol,{Print["Result: "," p= ",SetPrecision[p,10]," Error ",Abs[(b-a)/2.0], " Steps= ",i],Break[]}];
i=i+1;
If[Sign[fa]*Sign[fp]>0.0,{a=p,fa=fp},{b=p}]
]
If[i>maxit,Print["Method Failed to Converge after ",maxit, " steps."]]

Clear[p0,p1,tol, max, i,p,q]

p0=Part[data,Length[data],1];
p1=Part[data,Length[data]-1,1];
tol=10.0^(-12);
max=300;

i=2;
q0=f[p0];
q1=f[p1];
Print["Secant Method"]
Print["Input:  ","p0= ",p0," p1= ",p1," tol= ",tol]
While[i≤max ,
p=p1-q1(p1-p0)/(q1-q0);
data=Append[data,{p,f[p]}];
If[Abs[p-p1]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",Abs[p-p1], " Steps= ",i],Break[]}];
p0=p1;
q0=q1;
p1=p;
q1=f[p];
i=i+1]
If[i>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"p","f[p]"}}]

Error

error=Table[{n,Abs[data[[n,2]]]},{n,1,Length[data]}];
raerrplot=ListLogPlot[error,PlotStyle->{Purple,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{500,Automatic},BaseStyle->20]

Bisection Algorthim Alone

Clear[a,b,tol,fa,fp]

f[x_] =2x Cos[2x]-(x-2)^3

a=3.0;
b=4.0;
maxit=500;
tol=10.0^(-12);

data = {};
i=1;
fa=f[a];
Print["Bisection Method"]
Print["Input:  ","a= ",a," b= ",b," tol= ",tol]
While[i≤maxit,
p=a+(b-a)/2.0;
fp=f[p];
data = Append[data,{a,b,p,fp}];
If[fp==0||(b-a)/2.0<tol,{Print["Result: "," p= ",SetPrecision[p,10]," Error ",Abs[(b-a)/2.0], " Steps= ",i],Break[]}];
i=i+1;
If[Sign[fa]*Sign[fp]>0.0,{a=p,fa=fp},{b=p}]
]
If[i>maxit,Print["Method Failed to Converge after ",maxit, " steps."]]

Error

error=Table[{n,Abs[data[[n,4]]]},{n,1,Length[data]}];
bmerrplot=ListLogPlot[error,PlotStyle->{Purple,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{500,Automatic},BaseStyle->20]

Secant Algorthim Alone

Clear[p0,p1,tol, max, i,p,q]

f[x_] =2x Cos[2x]-(x-2)^3

p0=3.0;
p1=2.9;
tol=10.0^(-12);
max=300;

data={};
i=2;
q0=f[p0];
q1=f[p1];
Print["Secant Method"]
Print["Input:  ","p0= ",p0," p1= ",p1," tol= ",tol]
While[i≤max ,
p=p1-q1(p1-p0)/(q1-q0);
data=Append[data,{p,f[p]}];
If[Abs[p-p1]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",Abs[p-p1], " Steps= ",i],Break[]}];
p0=p1;
q0=q1;
p1=p;
q1=f[p];
i=i+1]
If[i>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"p","f[p]"}}]

Whoops, we got the wrong root!

Newton’s Method Module

Here is Newton’s method based module that will find a root to a desired tolerance.

NewtonsMethod[f_,ps_,tol_]:=Module[{p,p0,fp,i,max=300},

fp=D[f,x];
p0=ps;
p=p0-(f/fp/.x->p0);

i=1;
While[i≤max &&Abs[p-p0]>tol,
p0=p;
p=p0-(f/fp/.x->p0);
i=i+1];
If[i≥ max,Print["Method Failed to Converge"]];
Return[SetPrecision[p,12]];
]

f[x_] =2x Cos[2x]-(x-2)^2;
NewtonsMethod[f[x],2.0,(10.0^(-4))]

FindRoot[f[x]==0,{x,2.0}]

Newton’s Method Module Multiplicity m

Here is Newton’s method based module that will find a root to a desired tolerance.

NewtonsMethod[f_,ps_,tol_]:=Module[{p,p0,fp,i,max=300,error},

fp=D[f,x];
p0=ps;
p=p0-(f/fp/.x->p0);

error={};
i=1;
While[i≤max &&Abs[p-p0]>tol,
p0=p;
p=p0-(f/fp/.x->p0);
AppendTo[error,{i,Abs[f/.x->p0]}];
i=i+1
];
Print[ListLogPlot[error,PlotStyle->{Red,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{400,Automatic},BaseStyle->20]];
If[i≥ max,Print["Method Failed to Converge"]];
Return[SetPrecision[p,12]]
]

NewtonsMethodMultiplicity[f_,ps_,tol_]:=Module[{p,p0,fp,fpp,i,max=300,error},

fp=D[f,x];
fpp=D[f,{x,2}];
p0=ps;
p=p0-(f/fp/.x->p0);

error={};
i=1;
While[i≤max &&Abs[p-p0]>tol,
p0=p;
p=p0-(f*fp/((fp^2)-f*fpp)/.x->p0);
AppendTo[error,{i,Abs[f/.x->p0]}];
i=i+1
];
Print[ListLogPlot[error,PlotStyle->{Red,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{400,Automatic},BaseStyle->20]];
If[i≥ max,Print["Method Failed to Converge"]];
Return[SetPrecision[p,12]]
]

Multiplicity 2

Notice that Newont’s method without the modification has linear convergence when the root has multiplicity 2.

f[x_]=Exp[x]-x-1;
NewtonsMethod[f[x],2.0,(10.0^(-7))]
NewtonsMethodMultiplicity[f[x],2.0,(10.0^(-7))]
FindRoot[f[x]==0,{x,2.0}]

Multiplicity 3

f[x_]=Sin[x]+(x^2)Cos[x]-(x^2)-x;
NewtonsMethod[f[x],2.0,(10.0^(-9))]
NewtonsMethodMultiplicity[f[x],2.0,(10.0^(-9))]
FindRoot[f[x]==0,{x,2.0}]

Newton's Method: Complex Roots

We want to solve the equation f(x)=0 for x using Newton's Method, and see if we can locate a quadratic root. If we initialize with a complex number, Newton’s method will find complex valued roots!

Newton's Method

Clear[p0,p1,tol, max, i,p,q]

f[x_] =16x^4-40x^3+5x^2+20x+6

p0=-0.0+0.1I;
tol=10.0^(-8);
max=300;

data={};
i=1;
Print["Newton's Method"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[i≤max ,
p=p0-f[p0]/f'[p0];
data=Append[data,{i,p,f[p]}];
If[Abs[f[p]]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",Abs[p-p0], " Steps= ",i],Break[]}];
p0=p;
i=i+1]
If[i>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"n","p","f[p]"}}]

Error in Newton's Method

error=Table[{n,Abs[data[[n,3]]]},{n,1,Length[data]}];
nmerrplot=ListLogPlot[error,PlotStyle->{Purple,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{500,Automatic},BaseStyle->20]

Further Example: Better than NSolve?

f[x_]=π^(3x^2)-5^(Cos[x])
Solve[f[x]==0,x]
NSolve[f[x]==0,x]
Plot[f[x],{x,0,1},PlotStyle->Blue,ImageSize->{500,Automatic}]

Clear[p0,p1,tol, max, i,p,q]

f[x_]=π^(3x^2)-5^(Cos[x])
p0=1.0;
tol=10.0^(-12);
max=300;

data={};
i=1;
Print["Newton's Method"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[i≤max ,
p=p0-f[p0]/f'[p0];
data=Append[data,{p,f[p]}];
If[Abs[p-p0]≤tol,{Print["Result: ","p= ",SetPrecision[p,15]," Error ",Abs[p-p0], " Steps= ",i],Break[]}];
p0=p;
i=i+1]
If[i>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,15],TableHeadings->{None,{"p","f[p]"}}]

f[0.618029948531105599585089294123463`15]

error=Table[{n,Abs[data[[n,2]]]},{n,1,Length[data]}];
nmerrplot=ListLogPlot[error,PlotStyle->{Red,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},ImageSize->{500,Automatic},BaseStyle->20]

References:

Heikki Ruskeepaa, Mathematica Navigator, Graphics and Methods of Applied Mathematics (1999) San Diego,Academic Press.
Richard Burden & J. Douglas Faires, Numerical Analysis, 7th ed. (2001) Pacific Grove, Brooks/Cole.