Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Sept 2016
Lecture topic: Numerical Differentiation

Deriving Some (n+1)-Point Formulas for first derivative

Clear[n,k,x,xp,i,Lag,f]

Here is the wiggly polynomial defined. Notice I am using Mathematica's derivative abilities to calculate the derivative of the wiggly polynomial, which is trival to do since it is just a polynomial, but a pain for larger values of n since it is easy to make a simple mistake.

Lag[n_,k_,x_]:=((∏_(i=0)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))
Lagp[n_,k_,x_]:=D[Lag[n,k,x],x]

Here is the n+1-point formula that approximates the derivative:

fp[n_]:=DerivativeApprox[xp[j]]==(∑_(k=0)^n)f[xp[k]]*Lagp[n,k,xp[j]]+(derfc/((n+1)!))((∏_(k=0)^(j-1))(xp[j]-xp[k]))((∏_(k=j+1)^n)(xp[j]-xp[k]))

Note I want to track the order of the error term, so I am inlcuding the piece that depends on c.

2-Point Formulas: Equally Spaced Nodes

Define our equally spaced nodes:

xp[0]=(Subscript[x, 0])
xp[1]=(Subscript[x, 0])+h

Here is the simplification when n=1 and j=0:

fp[1]
fp[1]/.j->0

Here is the simplification when n=1 and j=1:

fp[1]
fp[1]/.j->1
fp[1]/.j->1 /.(Subscript[x, 0])->(Subscript[x, 0])-h

These are precisely the formulas we derived! We rule!

3-Point Formulas: Equally Spaced Nodes

Define our equally spaced nodes:

xp[0]=(Subscript[x, 0])
xp[1]=(Subscript[x, 0])+h
xp[2]=(Subscript[x, 0])+2h

Here is the simplification when n=2 and j=0:

fp[2]/.j->0

Here is the simplification when n=2 and j=1: (this is equation 5.7 in the text, when we set (Subscript[x, 0])+h=x)

fp[2]/.j->1

Here is the simplification when n=2 and j=2:

fp[2]/.j->2

Here are the simplifications:

fp[2]/.j->0
fp[2]/.j->1/.(Subscript[x, 0])->(Subscript[x, 0])-h

5-Point Formulas: Equally Spaced Nodes

Define our equally spaced nodes:

xp[0]=(Subscript[x, 0])
xp[1]=(Subscript[x, 0])+h
xp[2]=(Subscript[x, 0])+2h
xp[3]=(Subscript[x, 0])+3h
xp[4]=(Subscript[x, 0])+4h

Here is the simplification when n=4 and j=0:

fp[4]/.j->0

And of course, we could generate more if we liked: Here is the simplification when n=4 and j=2:

fp[4]/.j->2
fp[4]/.j->2/.(Subscript[x, 0])->(Subscript[x, 0])-2h

Notice this uses only four of the five points!

The Various Derivative Formulas

Here is a module to create all the  n-point derivative formulas.

Note that since it is defined as n-point instead of n+1-point, some of the product limits have shifted to n-1.

DerivativeFormulas[nn_]:=Module[{Lag,Lagp,fp,error,i,ii,n,k,xp},
Print[Style[nn,Bold,24],Style["-Point Derivative Formulas for f'((Subscript[x, 0])):",Bold,24]];
Do[xp[ii]=(Subscript[x, 0])+ii h,{ii,0,nn-1}];
Lag[n_,k_,x_]:=((∏_(i=0)^(k-1))((x-xp[i])/(xp[k]-xp[i]))) (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i]));
Lagp[n_,k_,x_]:=Simplify[((∏_(i=0)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))(((∑_(i=0)^(k-1))(1/(x-xp[i])))+((∑_(i=k+1)^n)(1/(x-xp[i])) ))];
fp[n_]:=(∑_(k=0)^(n-1))f[xp[k]] Lagp[n-1,k,xp[j]];
error[n_]:=((((∏_(k=0)^(j-1))(xp[j]-xp[k])) (∏_(k=j+1)^(n-1))(xp[j]-xp[k]))/(n!));results=Table[{ii,Simplify[fp[nn]/.j->ii/.(Subscript[x, 0])->(Subscript[x, 0])-ii h],error[nn]/.j->ii},{ii,0,nn-1}];Style[TableForm[results,TableHeadings->{None,{"j","Derivative Approximation","Error Goes As"}}]
,FontSize->22,FontFamily->"Times",FontColor->Black]]

DerivativeFormulas[2]

DerivativeFormulas[3]

DerivativeFormulas[4]

DerivativeFormulas[5]

You can also create other crazy derivative formulas:

DerivativeFormulas[6]

(n+1)-Point Formula for first derivative: Error Plots

Here is our favourite function:

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

Here are some derivative approximations:

ApproxDer1[x_,h_] = ((f[x+h]-f[x])/h);
ApproxDer2[x_,h_] = ((f[x+h]-f[x-h])/(2h));
ApproxDer3[x_,h_] = ((-3f[x]+4f[x+h]-f[x+2h])/(2h));
ApproxDer4[x_,h_] = ((f[x-2h]-8f[x-h]+8f[x+h]-f[x+2h])/(12h));

error1[x_,h_]= Abs[D[f[x],x]- ApproxDer1[x,h]];
error2[x_,h_]= Abs[D[f[x],x]- ApproxDer2[x,h]];
error3[x_,h_]= Abs[D[f[x],x]- ApproxDer3[x,h]];
error4[x_,h_]= Abs[D[f[x],x]- ApproxDer4[x,h]];

Here is the error:

a = 3.0;

Plot[{error1[a,h],error2[a,h],error3[a,h],error4[a,h]},{h,-0.01,0.01},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Purple,Thick}},AxesLabel->{"h","Error"},ImageSize->{500,Automatic},PlotRange->{All,{0,(10^(-4))}}]

And the error on the log scale:

LogPlot[{error1[a,h],error2[a,h],error3[a,h],error4[a,h]},{h,-0.01,0.01},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Purple,Thick}},AxesLabel->{"h","Log Error"},ImageSize->{500,Automatic}]

Notice that roundoff error has an important effect in the supposedly more accurate five point formula. This is why we do not typically proceed with higher point formulas. Zooming in shows the trouble clearly:

LogPlot[{error1[a,h],error2[a,h],error3[a,h],error4[a,h]},{h,-0.00005,0.00005},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Purple,Thick}},AxesLabel->{"h","Log Error"},ImageSize->{500,Automatic}]

This shows that the noise is due to round off error:

{error3[3.0,(10^(-6))],error4[3.0,(10^(-6))]}
SetPrecision[{error3[3,(10^(-6))],error4[3,(10^(-6))]},30]

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.