Math 4401 Numerical Methods
Root Finding

Notebook Author/Course Instructor: Barry McQuarrie
Date: Aug 2016

Function

f[x_] = Cos[x] BesselJ[4,x]
Plot[f[x],{x,0,10},AxesLabel->{"x","y"},ImageSize->{500,Automatic},PlotStyle->{Blue},BaseStyle->20,PlotLegends->{f[x]}]

Bisection Method

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

Algorthim

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

f[x_] = Cos[x] BesselJ[4,x]

a=4.0;
b=6.0;
maxit=500;
tol=(10.0^(-6));

data = {};
i=1;
fa=f[a];
Print["Bisection Method"]
Print["Input:  ","a= ",a," b= ",b," tol= ",tol]
If[fa*f[b]≥0,{Print["Condition f(a)f(b)<0 is not satisfied! Check your initial points."]}]
While[i≤maxit,
p=a+((b-a)/2.0);
fp=f[p];
data = Append[data,{a,b,p,fp}];
If[Abs[fp]<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 Estimate

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

Another Example, With Further Error Analysis

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

f[x_] = x^3+4x^2-10
Plot[f[x],{x,1,2},AxesLabel->{"x","y"},BaseStyle->20]

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

data = {};
i=1;
fa=f[a];
Print["Bisection Method"]
Print["Input:  ","a= ",a," b= ",b," tol= ",tol]
If[fa*f[b]≥0,{Print["Condition f(a)f(b)<0 is not satisfied! Check your initial points."]}]
While[i≤maxit,
p=a+((b-a)/2.0);
fp=f[p];
data = Append[data,{a,b,p,fp}];
If[Abs[fp]<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=Table[{n,Abs[data[[n,4]]]},{n,1,Length[data]}];
ListLogPlot[error,PlotStyle->{PointSize[0.02],Blue},AxesLabel->{"n","Logarithm of error"},ImageSize->{500,Automatic},BaseStyle->20]

Notice that the ninth iteration was close to the root, but the bisection method does not pick up on that and the next iteration is further from the root.

Here is the Error Bound we found analytically:

analyticerror=LogPlot[((2-1)/(2^n)),{n,1,Length[data]},PlotStyle->{Red},ImageSize->{500,Automatic},BaseStyle->20]

Here is the error at each stage, |(Subscript[p, n])-p|, which isn't what we plotted above as the error (above, we plotted |f((Subscript[p, n]))|  as an estimate of the error, since as we stepped through the process of finding the root we didn't know the value of p until we were done.

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

Show[errorplot,analyticerror]

So the analytic error bound we found is always larger than the actual error.

Another Example: Just be careful!

f[x_] = Expand[Product[(x-i),{i,1,20}]]
Plot[f[x],{x,1,3},AxesLabel->{"x","y"},BaseStyle->20]
f[1]
f[3]

Notice this function has zeros at the integers, if we don’t notice that we can get poor results. Just problem below is the condition f(a)f(b)<0 is not satisfied!

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

f[x_] = Expand[Product[(x-i),{i,1,20}]]

a=1.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]
If[fa*f[b]≥0,{Print["Condition f(a)f(b)<0 is not satisfied! Check your initial points."]}]
While[i≤maxit,
p=a+((b-a)/2.0);
fp=f[p];
data = Append[data,{a,b,p,fp}];
If[Abs[fp]<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]"}}]

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

f[x_] = Product[(x-i),{i,1,20}]

a=1.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]
If[fa*f[b]≥0,{Print["Condition f(a)f(b)<0 is not satisfied! Check your initial points."]}]
While[i≤maxit,
p=a+((b-a)/2.0);
fp=f[p];
data = Append[data,{a,b,p,fp}];
If[Abs[fp]<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]"}}]

Fixed Point Iteration

Here we want to solve equations of the form g(x)=x for x.

Cobweb diagram: Convergence

fixedPointPlot[f_,x_,x0_,n_,a_,b_,opts_]:=
Module[{t0=x0,y0,points={{x0,0}}},
Do[y0=f/.x->t0;points={points,{t0,y0}};
t0=y0;points={points,{t0,y0}},{n}];
points=Partition[Flatten[points],2];
Plot[{f,x},{x,a,b},Epilog->{Red,Line[points]},PlotStyle->{Blue,Black},opts]]

fixedPointPlot[Cos[x],x,0.2,8,0,1.1,{AspectRatio->Automatic,AxesLabel->{"x","y"},BaseStyle->20}]

Cobweb diagram: Divergence

fixedPointPlot[1-x^3,x,0.74,8,0,1.1,{AspectRatio->Automatic,AxesLabel->{"x","y"},BaseStyle->20}]

Algorithm

g[x_]=Cos[x]
Plot[Abs[g'[x]],{x,0,1},AxesLabel->{"x","y"},BaseStyle->20]

p0=0.2
data =Table[{n,SetPrecision[Nest[g,p0,n],10],If[n>0,Abs[Nest[g,p0,n]-Nest[g,p0,n-1]]]},{n,0,25}];

TableForm[data,TableHeadings->{None,{"i","(Subscript[x, i])","Error"}}]

Error Estimate

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

Another Example

g[x_]=x-Cos[x]BesselJ[4,x]
Plot[Abs[g'[x]],{x,4,6},AxesLabel->{"x","y"}]

This should converge since the absolute value of the derivative is less that one.

p0=4.0;
data =Table[{n,SetPrecision[Nest[g,p0,n],10],If[n>0,Abs[Nest[g,p0,n]-Nest[g,p0,n-1]]]},{n,0,25}];

TableForm[data,TableHeadings->{None,{"i","(Subscript[x, i])","Error"}}]

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

Comparison of Fixed Point to Bisection

Calculate (4^(1/3))to 8 decimals.

Here is the Bisection Method, where we find the root of f(x)=(x^3)-4.

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

f[x_] = (x^3)-4

a=1.0;
b=2.0;
maxit=500;
tol=(10.0^(-9));

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[Abs[fp]<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,12],TableHeadings->{None,{"a","b","p","f[p]"}}]

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

Solving for (4^(1/3))is equivalent to solving the fixed point problem x = 2/(Sqrt[x]).

g[x_]=2/Sqrt[x];
Plot[Abs[g'[x]],{x,1,2},BaseStyle->20,AxesLabel->{"x","|g'(x)|"}]
p0=1.0;
data =Table[{n,SetPrecision[Nest[g,p0,n],10],If[n>0,Abs[Nest[g,p0,n]-Nest[g,p0,n-1]]]},{n,0,25}];

TableForm[data,TableHeadings->{None,{"i","(Subscript[x, i])","Error"}}]

error=Table[{n-1,Abs[data[[n,3]]]},{n,2,Length[data]}];
eplot2=ListLogPlot[error,PlotStyle->{Red,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"},BaseStyle->Large]

Show[eplot1,eplot2,ImageSize->{500,Automatic}]

Since the value of k for the fixed point is about 1/2 in the region, we should not be surprised that the number of iterations to get the desired accuracy is about the same as the bisection method.

Notice the method fails if we use a different decomposition of the g(x):

g[x_]=4/(x^2);
Plot[Abs[g'[x]],{x,1,2},BaseStyle->20,AxesLabel->{"x","|g'(x)|"}]
p0=1.0;
data =Table[{n,SetPrecision[Nest[g,p0,n],10],If[n>0,Abs[Nest[g,p0,n]-Nest[g,p0,n-1]]]},{n,0,25}];

TableForm[data,TableHeadings->{None,{"i","(Subscript[x, i])","Error"}}]

References:

Heikki Ruskeepaa, Mathematica Navigator, Graphics and Methods of Applied Mathematics (1999) San Diego,Academic Press.