Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie

Aug 2016

Lecture topic: Course Introduction

Error in Derivative Approximations

f[x_] = Cos[x] BesselJ[4,x]
Plot[f[x],{x,0,10},PlotStyle->{Thick},BaseStyle->20,PlotLegends->{f[x]}]

ApproxDer1[x_,h_] = ((f[x+h]-f[x])/h);
ApproxDer2[x_,h_] = ((f[x+h]-f[x-h])/(2h));
error1[x_,h_]= Abs[D[f[x],x]- ApproxDer1[x,h]];
error2[x_,h_]= Abs[D[f[x],x]- ApproxDer2[x,h]];

Approximation 1

$Post=(#/.x_/;x≤0.01:>ScientificForm[x])&;
a = 3.0;
N[D[f[x],x] /. x->a,4]
stepsize = 0.001;
Table[{ApproxDer1[a,h] ,error1[a,h]},{h,-10*stepsize,10*stepsize,stepsize}];
TableForm[%,TableHeadings->{None,{"((f (a + h) - f (a))/h)","Absolute Error"}}]
$Post=.
list1 = Table[{h,error1[a,h]},{h,-10*stepsize,10*stepsize,stepsize}];
plot1=ListPlot[list1,PlotStyle->{Blue,PointSize[.02]},AxesLabel->{"h","Absolute Error"},BaseStyle->20]

Approximation 2

$Post=(#/.x_/;x≤0.01:>ScientificForm[x])&;
a = 3.0;
N[D[f[x],x] /. x->a,4]
stepsize = 0.001;
Table[{ApproxDer2[a,h] ,error2[a,h]},{h,-10*stepsize,10*stepsize,stepsize}];
TableForm[%,TableHeadings->{None,{"((f (a + h) - f (a - x))/(2  h))","Absolute Error"}}]
$Post=.
list1 = Table[{h,error2[a,h]},{h,-10*stepsize,10*stepsize,stepsize}];
plot2=ListPlot[list1,PlotStyle->{Red,PointSize[.02]},AxesLabel->{"h","Absolute Error"},BaseStyle->20]

Show[plot1,plot2]

Roundoff Error

Plot[error1[a,h],{h,-10^(-5),10^(-5)},AxesLabel->{"h","Absolute Error"},PlotStyle->{Blue},BaseStyle->20]

Plot[error1[a,h],{h,-10^(-6),10^(-6)},AxesLabel->{"h","Absolute Error"},PlotStyle->{Blue},BaseStyle->20]

Plot[error1[a,h],{h,-10^(-11),10^(-11)},AxesLabel->{"h","Absolute Error"},PlotStyle->{Blue},BaseStyle->20]

There are ways in Mathematica to set the precision of a result. Check this out:

SetPrecision[1.2,20]
% - 12/10

SetPrecision[1.2000000000000000000,20]
%-12/10

N[12/10,20]
%-12/10

Taylor Series

Example f(x)=cos(x)

f[x_] :=Cos[x]

PT[a_,n_,x_] :=(∑_(j=0)^n)Derivative[j][f][a]((((x-a))^j)/(j!))

R[a_,ξ_,n_,x_]:=Simplify[Derivative[n+1][f][ξ]((((x-a))^(n+1))/((n+1)!))]

TableForm[Table[{n,PT[0,n,x],R[0,ξ,n,x]},{n,0,10}],TableHeadings->{None,{"n","(Subscript[P, n])(x)","(Subscript[R, n])(x)"}}]

TableForm[Table[{n,PT[Pi/2,n,x],R[Pi/2,ξ,n,x]},{n,0,3}],TableHeadings->{None,{"n","(Subscript[P, n])(x)","(Subscript[R, n])(x)"}}]

plotPT[n_]:= Plot[{f[x],PT[0,n,x]},{x,-2Pi,2Pi},AxesLabel->{"x","y"},PlotStyle->{Red,Blue},PlotRange->{{-2Pi,2Pi},{-1,1}},BaseStyle->20]

Manipulate[Show[plotPT[n]],{n,0,20,1}]

Integral (∫_0^1)cos(x)dx

The integral we want to do:

actual=SetPrecision[Integrate[f[x],{x,0,0.1}],12]

Error Estimate, n=2:

Abs[R[0,ξ[x],2,x]]

We can get a bound on sin(ζ(x)) since we know |sin(ζ(x))|≤1,

error=(1/6) Abs[(x^3)]

SetPrecision[Integrate[PT[0,2,x],{x,0,0.1}],12]
SetPrecision[Integrate[error,{x,0,0.1}],12]

Error Estimate, n=3:

Abs[R[0,ξ[x],3,x]]

We can get a bound on cos(ζ(x)), since |cos(ζ(x))|≤1,

error=(1/24) (x^4)

SetPrecision[Integrate[PT[0,3,x],{x,0,0.1}],12]
SetPrecision[Integrate[error,{x,0,0.1}],12]

Actual Error:

Integrate[f[x],{x,0,0.1}]-Integrate[PT[0,2,x],{x,0,0.1}] //N

Notice that the actual error is slightly better than our error bound.

Another Example: Integral We Can Not Do Exactly

Example (∫_0^1)(e^(-(x^4)))cos(x)dx

Here is a simple looking integral.

f[x_] = Exp[-x^4]Cos[x]
Integrate[f[x],{x,0,1}]
SetPrecision[NIntegrate[f[x],{x,0,1}],20]

How many decimals are correct in our answer? Without further analysis, we simply don't know.

Let's try and do it by using Taylor Series. First, define the general Taylor series:

PT[a_,n_,x_] :=(∑_(j=0)^n)Derivative[j][f][a]((((x-a))^j)/(j!))
R[a_,ξ_,n_,x_]:=Simplify[Derivative[n+1][f][ξ]((((x-a))^(n+1))/((n+1)!))]

Just for kicks, let's choose to expand about a=1/3. In practice you would probably expand about a=1/2.

x0=1/3;
TableForm[Table[{n,PT[x0,n,x],R[x0,ξ[x],n,x]},{n,0,1}],TableHeadings->{None,{"n","(Subscript[P, n])(x)","(Subscript[R, n])(x)"}}]

That remainder is nasty looking!

Here is a plot of how the approximation looks for various values of n:

plotPT[n_]:= Plot[{f[x],PT[x0,n,x]},{x,0,1},AxesLabel->{"x","y=f(x)"},PlotStyle->{Red,Blue},PlotRange->{{0,1},{0,1}}]

Manipulate[Show[plotPT[n]],{n,0,15,1}]

Notice the Taylor Polynomial approximates the function very well once n=15.

Here are some approximations to the integral. Notice that I do the integration symbolically (since it is an integral over a polynomial, it can be done easily), and then approximate as a decimal at the end.

TableForm[Table[{n,SetPrecision[Integrate[PT[x0,n,x],{x,0,1}] ,10]},{n,0,20}],TableHeadings->{None,{"n","(∫_0^1)(Subscript[P, n])(x) \[DifferentialD]x"}}]

This is all fine, but we need to estimate the error in these approximation. For each value of n, you will get a different error estimate using the remainder term from the Taylor polynomial.

Estimating the Error for a given n=15

Here is the remainder for n=15.

R[x0,ξ,15,x]

To get an estimate for the error in  (∫_0^1)(Subscript[P, 15])(x) dx, we need to estimate (∫_0^1)|(Subscript[R, 15])(ζ(x))| dx.

Notice that we can separate the parts that depend on ζ and x:

ξpart[ξ_]=((E^(-(ξ^4))) ((726745175521-51987401704800 (ξ^2)-777963501631680 (ξ^4)+3317347889854080 (ξ^6)+19679165616015360 (ξ^8)-25108410579732480 (ξ^10)-102949325633889280 (ξ^12)+56414063392727040 (ξ^14)+193731468867256320 (ξ^16)-52719313863344128 (ξ^18)-167954270687723520 (ξ^20)+24146461175316480 (ξ^22)+76578563140485120 (ξ^24)-5900076192890880 (ξ^26)-19784184700600320 (ξ^28)+796530701238272 (ξ^30)+3007804184985600 (ξ^32)-58809510789120 (ξ^34)-271116119572480 (ξ^36)+2198486384640 (ξ^38)+14092861440000 (ξ^40)-32212254720 (ξ^42)-386547056640 (ξ^44)+4294967296 (ξ^48)) Cos[ξ]+64 ξ (-212500688190+1268971933319 (ξ^2)+41529673784340 (ξ^4)-33898345331760 (ξ^6)-548228261380240 (ξ^8)+154273575980160 (ξ^10)+1850205427150080 (ξ^12)-233227667382528 (ξ^14)-2444079141504000 (ξ^16)+153267369615360 (ξ^18)+1545120360509440 (ξ^20)-49903463301120 (ξ^22)-520674143109120 (ξ^24)+8541273128960 (ξ^26)+98941548625920 (ξ^28)-771474063360 (ξ^30)-10803213238272 (ξ^32)+34351349760 (ξ^34)+664126095360 (ξ^36)-587202560 (ξ^38)-21139292160 (ξ^40)+268435456 (ξ^44)) Sin[ξ]))/900657498850357248000;
xpart[x_]=(((1-3 x))^16) ;

error=ξpart[ξ]*xpart[x]

We need a bound on the ζ part over the interval ζ between 0 and 1 (the integration region).

Plot[Abs[ξpart[ξ]],{ξ,0,1},AxesLabel->{"ξ","ξ part of (Subscript[R, 15])(x))"},PlotStyle->{Red},BaseStyle->20]

ξpart[1]
N[%]

Notice the maximum value occurs at the left endpoint of the interval. If it didn't, we could examine the graph more closely to get an upper bound on the ξ part.

Plot[{Abs[ξpart[ξ]],ξpart[1]},{ξ,0,1},AxesLabel->{"ξ","ξ part of (Subscript[R, 15])(x))"},PlotStyle->{Red,Blue},ImageSize->{400,Automatic}]

Whoops! Actually, the error is not bounded by the value at the endpoint! We could use calculus to try and get the value at which the maximum occurs, or we can use a sketch to justify an upper bound. Let's do the latter.

Plot[{Abs[ξpart[ξ]],7.5*(10^(-7))},{ξ,0,1},AxesLabel->{"ξ","ξ part of (Subscript[R, 15])(x))"},PlotStyle->{Red,Blue},ImageSize->{400,Automatic}]

So an upper bound for the ξ part is 7.5*(10^(-7)). Therefore, we can estimate the remainder by

Integrate[ xpart[x],{x,0,1}]
N[%]

Integrate[ 7.5*(10^(-7))*xpart[x],{x,0,1}]