Numerical Methods in Engineering (ENGR20005)
Lecture 13 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
Differentiation
Derivatives are prevalent in many problems and a numerical approximation is often required.
Recap: Fluid flows represented by Navier-Stokes eqs.:
∂ρ + ∂ (ρuk) = 0 ∂ ∂∂t∂xk
∂t (ρui) + ∂xk [ρuiuk + ρδik − τik] = 0 ∂(ρE)+ ∂ [ρukH+qk−ui(τik−ρσik)]=0
∂t ∂xk
τik = 2μ(Sik − 1Sjjδik)
Sik = 1(∂ui + ∂uk) 2∂xk ∂xi
Re 3 qk = − κ ∂T
PrEcRe ∂xk
2
3
∂ρ + ∂ (ρuk) = 0 ∂ ∂∂t∂xk
∂t (ρui) + ∂xk [ρuiuk + ρδik − τik] = 0
∂(ρE)+ ∂ [ρukH+qk−ui(τik−ρσik)]=0 ∂t ∂xk
L13.1:
Approximating First Derivatives: Equally Spaced Data
4
Differentiation of Equally Spaced Data
In many applications, one needs to find the derivatives of a function. This section of the notes describes a procedure that one could use to find the spatial derivatives of a function.
Suppose you know the values of a function at the three points shown below
fi−1 fi fi+1 ΔΔ
xi−1
We would like to approximate the derivative at x = xi as
xi
fi′ = Δ1 (a−1fi−1 + a0 fi + a1fi+1)
fi′ is the approximation of the derivative df/dx(x = xi), at x = xi. We will assume for the
moment that all data is equally spaced i.e.
Δ = xi+1 − xi = constant .
Our task is to find a−1,a0 and a1 such that fi′ best approximates the derivative df/dx(x = xi) 5
xi+1
(13.1)
fi−1 fi fi+1 ΔΔ
′ 1 xi−1 xi xi+1 fi = Δ(a−1fi−1 + a0 fi + a1fi+1)
(13.1)
fi−1, fi and fi+1 are defined to be
fi−1 = f(xi − Δ) fi = f(xi)
fi+1 = f(xi + Δ)
So performing Taylor series expansion on Eqs. (13.2)-(13.4) gives
(13.2)
(13.3) (13.4)
fi−1=f(xi−Δ)=f(xi)−Δf′(xi)+ 1Δ2f′′(xi)− 1Δ3f′′′(xi)+ 1Δ4fiv(xi)±HOT. 2! 3! 4!
(13.5) (13.6)
(13.7)
f(xi) = f(xi)
fi+1=f(xi+Δ)=f(xi)+Δf′(xi)+ 1Δ2f′′(xi)+ 1Δ3f′′′(xi)+ 1Δ4fiv(xi)±HOT. 2! 3! 4!
6
fi′ = Δ1 (a−1fi−1 + a0 fi + a1fi+1) Substituting (13.5)-(13.7) into the RHS of Eq. (13.1) gives
(13.1)
(13.5) (13.6)
(13.7)
fi−1=f(xi−Δ)=f(xi)−Δf′(xi)+ 1Δ2f′′(xi)− 1Δ3f′′′(xi)+ 1Δ4fiv(xi)±HOT. 2! 3! 4!
f(xi) = f(xi)
fi+1=f(xi+Δ)=f(xi)+Δf′(xi)+ 1Δ2f′′(xi)+ 1Δ3f′′′(xi)+ 1Δ4fiv(xi)±HOT. 2! 3! 4!
fi′=Δ1(a−1+a0+a1)f(xi) +
+1(a 1+a1)Δ2d2f(x) Δ −12! 12! dx2 i
+1(−a−11+a11)Δ3d3f(xi) Δ 3! 3! dx3
+O(Δ3)
Error term
a−1+a0+a1=0
1 1 d2f i
−a 1+a1=0 : d3f(x) −13! 13! dx3 i
1 (−a−1+a1)Δdf(xi) Δ dx
: f(xi) −a−1+a1=1 : df(xi)
7
dx a−12+a12!=0 : dx2(x)
1
Δ −1 0 1
a +a +a =0 1 −a +a =1
Error
1 −1 11
We only have 3 unknowns so we at most only satisfy 3 equations.
Δ a−1 2 + a1 2! = 0 1 1
Δ2 −a−1 3! + a1 3! = 0
We have 3 unknowns, a−1, a0 and a1, thus we can satisfy the first 3 equations at most. If
we only choose to satisfy the first two equations, then we will have a numerical scheme that is first order accurate (the error term is O(Δ)). We will also have one free variable. Let’s set
in
+ a0 + a1 = 0 to give
+ a1 = 1
With this set of coefficients, the derivative of f(x) can be approximated as
fi′=Δ1(a−1fi−1+ fi+a1fi+1)=Δ1(−fi+fi+1)
This is called the 1st order (because the error is O(Δ)) forward (because it requires
information at xi+1) difference scheme.
8
a−1 = 0, a−1
−a−1
a0
a−1 = 0 a0 = − 1 a1 = 1
1 Δ1
Error
a +a +a =0 −1 0 1
−a + a = 1 1 −1 11
a−1 2 + a1 2! = 0
If we had chosen to set a1 = 0, then we can only satisfy the first two conditions
Δ
We only have 3 unknowns so we at most only satisfy 3 equations.
(because we only have two other unknowns).
a−1 + a0 + a1 = 0 −a−1 + a1 = 1
a−1 = − 1 to give
a0 = 1
a1 = 0
With this set of coefficients, the derivative of f(x) can be approximated as
f i ′ = Δ1 ( a − 1 f i − 1 + a 0 f i + a 1 f i + 1 )
This is called the 1st order (because the error is O(Δ)) backward (because it requires
information at xi−1) difference scheme.
= Δ1 ( − f i − 1 + f i )
9
1 Δ1
Error
If we now choose to satisfy all 3 equations above, we will get the solution for a−1, a0 and a1 given below. The formula we get will be second order accurate because we are now also
satisfying the 3rd equation above.
a − 1 = − 12 With the corresponding formula for fi′,
fi′=1(a−1fi−1+a0fi+a1fi+1) = 1 (fi+1−fi−1). Δ 2Δ
This scheme is called the 2nd order (because the error is O(Δ2)) centred (because it requires the same amount of information at both xi−1 and xi+1) difference scheme.
Δ
a−1 + a0 + a1 = 0 −a−1 + a1 = 1
a−1 1 + a1 1 = 0 2 2!
a0 = 0
a 1 = 12
10
In general the derivative of a function, f(x), can be approximated in terms the function values at neighbouring points as
f i ′ = Δ1 ∑M a − j f i − j + a 0 f i + ∑M a j f i + j j=1 j=1
(13.8)
f i ′ = Δ1 ( a − M f i − M + a − M + 1 f i − M + 1 + . . . . + a − 1 f i − 1 + a 0 f i + a 1 f i + 1 + . . . + a M − 1 f i + M − 1 + a M f i + M )
where fi′ is the approximation of the derivative df/dx(x = xi) at x = xi. We will assume for the
moment that all data is equally spaced i.e.
Δ=xi+1−xi =xi−xi−1 =constant.
In Eq. (13.8)
fi±j = f(xi ± jΔ) 11
f i ′ = Δ1 ( a − M f i − M + a − M + 1 f i − M + 1 + . . . . + a − 1 f i − 1 + a 0 f i + a 1 f i + 1 + . . . + a M − 1 f i + M − 1 + a M f i + M ) fi±j = f(xi ± jΔ)
Performing Taylor series expansion gives
(13.8)
f(xi±jΔ)=f(xi)±(jΔ)f′(xi)+ 1(jΔ)2f′′(xi)± 1(jΔ)3f′′′(xi)+ 1(jΔ)4fiv(xi)±HOT. 2! 3! 4!
Substituting into the RHS of Eq. (13.8) gives
f i ′ = Δ1 ( a − M + a − M + 1 + . . . + a − 1 + a 0 + a 1 + . . . a M − 1 + a M ) f ( x i ) +1 (a−M(−M)+…+a−1(−1)+a1+…+aMM))Δdf(xi)
Δ dx
+1 a−M(−M)2 +…+a−1(−1)2 +a112 +…+aMM2 Δ2 d2f(xi)
Δ(2! 2!2! 2!)dx
2
Δ(3! 3!3! 3!)dx +O(Δ3)
+1 a−M(−M)3 +…+a−1(−1)3 +a113 +…+aMM3 Δ3d2f(xi)
12
2
f i ′ = Δ1 ( a − M + a − M + 1 + . . . + a − 1 + a 0 + a 1 + . . . a M − 1 + a M ) f ( x i ) +1 (a−M(−M)+…+a−1(−1)+a1+…+aMM))Δdf(xi)
Δ dx
+1 a−M(−M)2 +…+a−1(−1)2 +a112 +…+aMM2 Δ2 d2f(xi)
Δ(2! 2!2! 2!)dx
2 +1 a−M(−M)3 +…+a−1(−1)3 +a113 +…+aMM3 Δ3d3f(xi)
Δ(3! 3!3! 3!)dx +O(Δ3)
3
So if we want fi′ to approximate df(xi)/dx as closely as possible, we would need to satisfy the
following equations
1 Δ1
a−M
The more equations we choose to satisfy, the higher the accuracy of the approximation. For
O(Δ). If we satisfy the first three equations, then the approximation is O(Δ2). 13
1 2
+a12!+…+aM 2! =0
dx3 example, if we only satisfy the first two equations, then the accuracy of the approximation is
Δ Δ2
+…+a−1
3! 3! 3! 3!
dx2
: d3f (xi)
M 2 a−M(−M)3 +…+a−1(−1)3 +a113 +…+aMM3 =0
( − M ) 2 2!
( − 1 ) 2 2!
: f(xi)
: df (xi)
a−M + a−M+1 + . . . + a−1 + a0 + a1 + . . . aM−1 + aM = 0
a−M(−M) + … + a−1(−1) + a1 + … + aMM = 1
dx
: d2f (xi)
Exercise
fi fi+1 fi+2 ΔΔΔΔ
fi−2 fi−1
xi−2 xi−1 xi xi+1 xi+2
A five-point stencil finite difference formula can be obtained by setting M = 2 in Eq. (13.8) to obtain
fi′ = Δ1 (a−2fi−2 +a−1fi−1 +a0fi +a1fi+1 +a2fi+2)
Show that in order to obtain a fourth order accurate finite difference formula, one must satisfy the following 5 equations
a−2+a−1+a0+a1+a2 =0 −2a−2 − a−1 + a1 + 2a2 = 1
a =0 0 2
2a−2 + 12 a−1 + 12 a1 + 2a2 = 0 4 1 1 4
which can be solve to give
a1 = − a−1 = 3
a2 =−a−2 =−12
−3a−2−6a−1+6a1+3a2 =0 2 a−2 + 1 a−1 + 1 a1 + 2 a2 = 0
1
3 24 24 3
fi′= 1 (1 fi−2− 8 fi−1+ 8 fi+1− 1 fi+2) 12Δ 12 12 12 12
14
Exercise
f f f
i i+1 i+2
fi−2 fi−1
Δ Δ Δ Δ
x x x x x i−2 i−1 i i+1 i+2
A three-point stencil, one-sided finite difference formula can be written as
f i ′ = Δ1 ( a 0 f i + a 1 f i + 1 + a 2 f i + 2 )
Show that in order to obtain a second order accurate finite difference formula, one must
satisfy the following 3 equations
a0 + a1 + a2 = 0 1 a1 + 2a2 = 1
which can be solve to give
a 0 = − 32 a1 = 1 1
a2 = − 2
2 a1 + 2a2 = 0
So we get the 2nd order forward differencing scheme.
f i ′ = Δ1 ( − 32 f i + 2 f i + 1 − 12 f i + 2 ) 15
Some popular finite difference schemes
Central difference scheme Error
fi′ = 1 (−fi−1 + fi+1) O(Δ2) 1 2Δ
f i′ = 12Δ ( fi−2 − 8fi−1 + 8fi+1 − fi+2) O(Δ)4 Forward difference scheme Error
f i ′ = Δ1 ( − f i + f i + 1 ) O ( Δ ) fi′ = 1 (−3fi + 4fi+1 − fi+2) O(Δ2)
1 2Δ
fi′ = 6Δ (−11fi + 18fi+1 − 9fi+2 + 2fi+3) O(Δ3)
One can also come up with a partially one sided scheme, for example
fi′ = 1 (−2fi−1 − 3fi + 6fi+1 − fi+2) 6Δ
16
Backward difference scheme Error
f i ′ = Δ1 ( f i − f i − 1 ) O ( Δ ) fi′ = 1 (3fi − 4fi−1 + fi+2) O(Δ2)
1 2Δ
fi′ = 6Δ (−2fi−3 + 9fi−2 − 18fi−1 + 11fi) O(Δ3)
One can also come up with a partially one sided scheme, for example
fi′ = 1 (fi−2 − 6fi−1 + 3fi + 2fi+1) 6Δ
17
Example L13.1: This example is adapted from the book “Numerical Methods for Engineers” by Chapra & Canale. Approximate the value of the derivative at x = 2 of the function
f(x) = x2 cos(x)
using the 2nd and 4th order central differencing schemes. Compare with the exact answer
df =x(2cos(x)−xsin(x)) dx
Plot log10 of the error vs log10 Δ for Δ = 0.001,0.002,0.005,0.01,0.1,0.2,0.4,1.0.
18
Lecture13M03.m
MyFunc=@(x) x.^2.*cos(x) MyDeriv=@(x)x.*(2.*cos(x)-x.*sin(x))
% Define the vector of Delta values
DeltaVec=[0.001 0.002 0.005 0.01 0.1 0.2 0.4 1.0] ;
%Defining data point
xi=2.0;
for k=1:length(DeltaVec)
Delta=DeltaVec(k) xim2=xi-2*Delta xim1=xi-Delta; xip1=xi+Delta; xip2=xi+2*Delta;
%2nd order centred
dfdx2ndCentered=(1./(2*Delta))*(-MyFunc(xim1)+MyFunc(xip1)); error2ndCentered(k)=abs(MyDeriv(xi)-dfdx2ndCentered);
%4th order centred
′ 1
fi = 2Δ(−fi−1 +fi+1)
dfdx4thCentered=(1./(12*Delta))*(MyFunc(xim2)-8*MyFunc(xim1)+8*MyFunc(xip1)-MyFunc(xip2));
error4thCentered(k)=abs(MyDeriv(xi)-dfdx4thCentered);
end
hold off plot(log10(DeltaVec),log10(error2ndCentered),’ko-‘) hold on plot(log10(DeltaVec),log10(error4thCentered),’go-‘) xlabel(‘log(\Delta)’,’FontSize’,14) ylabel(‘log(Error)’,’FontSize’,14)
fi′ =
1 (fi−2 −8fi−1 +8fi+1 −fi+2) 12Δ
2 1
19
4 1
End of Example L13.1
20
Richardson’s Extrapolation
The exact value of an integral can be decomposed into its numerical approximation and an error term.
I = I(Δ) + E(Δ)
(11.14)
Throwback Thursday #TBT
I is the exact value of the integral
I(Δ) is the integral predicted by the integration scheme for spacing Δ E(Δ) is the error.
If we make two predictions of I, one with spacing Δ1 and also with spacing Δ2, then
I = I(Δ1) + E(Δ1) = I(Δ2) + E(Δ2)
We will assume that Δ2 < Δ1. If we look at the multiple application of the
Trapezoidal rule, we know from error analysis (see Eq. (11.8)) that
E(Δ1) ≈ Δ21 E(Δ2) Δ2
if we assume that the double derivative, f (2)(ξ) term is constant. 21
(11.15)
(11.16)
I = I(Δ1) + E(Δ1) = I(Δ2) + E(Δ2)
Substitute Eq. (11.16) into Eq. (11.15) gives
I(Δ1)+E(Δ2)Δ21 ≈I(Δ2)+E(Δ2) Δ2
E(Δ2) ≈ I(Δ1) − I(Δ2) 1 − ( Δ1 )2
(11.16)
(11.15) E(Δ1) ≈ Δ1 E(Δ2) Δ2
Thro2wback Thursday #TBT
Substitute back into Eq. (11.15) gives
I = I(Δ2) + E(Δ2) ≈ I(Δ2) + I(Δ1) − I(Δ2) 1 − ( Δ1 )2
Δ2
Δ2
1
approximated values of the integral, I(Δ1) and I(Δ2). 22
−I(Δ1)
Trapezoidal rule. It gives you a much better estimate of the integral if you have two
I≈I(Δ2) 1+
Equation (11.17) is the Richardson extrapolation formula for the multiple application of the
( Δ1 )2 − 1 Δ2
( Δ1 )2 − 1 Δ2
(11.17)
1
Can we do this for Differentiation?
23
L13.2:
Approximating Higher order derivatives
24
Generally, the finite difference formula for the second derivative, f′′ ≈ d2f/dx2(x)
can be expressed as
i i fi =Δ2 ∑a−jfi−j+a0fi+∑ajfi+j .
′′1M M j=1 j=1
(13.9)
where
You can expand the right hand side using Taylor series gives
fi±j = f(xi ± jΔ)
a−M + a−M+1 + . . . + a−1 + a0 + a1 + . . . aM−1 + aM = 0
f(xi±jΔ)=f(xi)±(jΔ)f′(xi)+ 1(jΔ)2f′′(xi)± 1(jΔ)3f′′′(xi)+ 1(jΔ)4fiv(xi)±HOT. 2! 3! 4!
Substitute into Eq. (13.9) gives
a−M(−M) + ... + a−1(−1) + a1 + ... + aMM = 0 a−M(−M)2 +...+a−1(−1)2 +a112 +...+aMM2 =1
2! 2! 2! 2!
a−M(−M)3 +...+a−1(−1)3 +a113 +...+aMM3 =0 3! 3! 3! 3!
As for the case as the first derivative, the more equations you satisfy, the more accurate your finite difference approximation will be.
25
Consider the case where M = 1.
f′′=1(a f +af+af )
(13.9)
i Δ2 −1 i−1 0 i 1 i+1 fi−1, fi and fi+1 are defined to be
fi−1 = f(xi − Δ) fi = f(xi)
(13.2)
(13.3) (13.4)
fi+1 = f(xi + Δ)
Performing Taylor series expansion on the right hand side of Eqs.(13.2)-(13.4)
f(xi) = f(xi)
(13.5) (13.6)
(13.7)
fi−1=f(xi−Δ)=f(xi)−Δf′(xi)+ 1Δ2f′′(xi)− 1Δ3f′′′(xi)+ 1Δ4fiv(xi)±HOT. 2! 3! 4!
fi+1=f(xi+Δ)=f(xi)+Δf′(xi)+ 1Δ2f′′(xi)+ 1Δ3f′′′(xi)+ 1Δ4fiv(xi)±HOT. 2! 3! 4!
26
f′′=1(a f +af+af )
(13.10)
i Δ2 −1 i−1 0 i
f(xi) = f(xi)
Substituting (13.5)-(13.7) into the RHS of Eq. (13.10) gives
1 i+1
fi−1=f(xi−Δ)=f(xi)−Δf′(xi)+ 1Δ2f′′(xi)− 1Δ3f′′′(xi)+ 1Δ4fiv(xi)±HOT. 2! 3! 4!
(13.5) (13.6)
(13.7)
fi+1=f(xi+Δ)=f(xi)+Δf′(xi)+ 1Δ2f′′(xi)+ 1Δ3f′′′(xi)+ 1Δ4fiv(xi)±HOT. 2! 3! 4!
Δ1 ( a − 1 + a 0 + a 1 ) f ( x i ) 2
f ′′ = i
+1 (−a−1+a1)Δdf(xi) Δ2 dx
+ 1 (−a−1 1 + a1 1 ) Δ3 d3f (xi) Δ2 3! 3! dx3
Compare
coeffs both sides give
a−1 + a0 + a1 = 0 −a−1 + a1 = 0
Solve to get
a−1 = a1 = 1 a0 = − 2
+ 1 (a−1 1 + a1 1 ) Δ2 d2f (xi) Δ2 2! 2! dx2
+O(Δ2)
Note that this approximation is second order accurate (error is O(Δ2)) since the fact that
a−1 = a1 makes the coefficient of d3f/dx3 zero. 27
f′′=1(f −2f+f ) i Δ2 i−1 i i+1
a−1 1 + a1 1 = 1 2 2!
Some popular finite difference schemes for 2nd order derivatives Central difference scheme Error
f′′=1(f −2f+f )O(Δ2) i Δ2 i−1 i i+1
Forward difference scheme Error
f′′=1(f −2f +f)O(Δ) i Δ2 i+2 i+1 i
Backward difference scheme Error
f′′=1(f −2f +f)O(Δ) i Δ2 i−2 i−1 i
28
In order to find finite difference formula for higher order derivatives of order n, one would need to try the general formula
n1M M
fi =Δn ∑a−jfi−j+a0fi+∑ajfi+j
j=1
j=1
29
Example L13.2: This example is taken from the book “Numerical Methods for Engineers” by Chapra & Canale.
A plane is tracked by radar and data is taken every second in polar coordinates r and θ. The readings are given in the table below.
t(s) 100 101 102 103 104 105
q(rad) 0.75 0.72 0.70 0.68 0.67 0.66
r(m) 1560 1636 1694 1767 1838 1902
Calculate the velocity and acceleration at t = 103s. The velocity and acceleration polar coordinates are ·r ·θ ···2r ····θ
v⃗=re⃗+rθe⃗ and a⃗=(r−rθ)e⃗+(rθ+2rθ)e⃗.
Plot all components of velocity and acceleration for t ∈ [100,105].
30
Lecture13M02.m
fi′= 1(fi+1−fi−1) 2Δ
rdot4=(radius(5)-radius(3))/(2*Delta); thetadot4=(theta(5)-theta(3))/(2*Delta);
rdoubledot4=(radius(5)-2*radius(4)+radius(3))/(Delta.^2); thetadoubledot4=(theta(5)-2*theta(4)+theta(3))/(Delta.^2);
v ⃗ = r· e r ⃗ + r θ· e θ ⃗
ur4=rdot4; utheta4=radius(4)*thetadot4;
f′′=1(f −2f+f ) i Δ2 i−1 i i+1
ar4=rdoubledot4-radius(4)*thetadot4.^2; atheta4=radius(4)*thetadoubledot4+2*rdot4*thetadot4;
···2r ····θ a ⃗ = (r − rθ )e ⃗ + (rθ + 2rθ)e ⃗
sprintf('The radial and azimuthal velocity at t=103s is %f %f ',ur4,utheta4) sprintf('The radial and azimuthal acceleration at t=103s is %f %f ',ar4,atheta4)
%Get data to plot magnitude of velocity and acceleration
ur=zeros(size(t));
utheta=zeros(size(t));
ar=zeros(size(t));
atheta=zeros(size(t));
rdot=zeros(size(t));
thetadot=zeros(size(t));
rdoubledot=zeros(size(t));
thetadoubledot=zeros(size(t)); 31
t=[100 101 102 103 104 105]
theta=[0.75 0.72 0.7 0.68 0.67 0.66] radius=[1560 1636 1694 1767 1838 1902]
%Delta=1s
Delta=1
figure(1)
plot(t,theta,'ko-') xlabel('t(s)');ylabel('\theta') figure(2) plot(t,radius,'bo-') xlabel('t(s)');ylabel('r(m)')
%Index 4 correspond to t=103
f i ′ = Δ1 ( − f i + f i + 1 )
rdot(1)=(radius(2)-radius(1))/Delta; thetadot(1)=(theta(2)-theta(1))/Delta;
rdoubledot(1)=(radius(3)-2*radius(2)+radius(1))/(Delta.^2); thetadoubledot(1)=(theta(3)-2*theta(2)+theta(1))/(Delta.^2);
f′′=1(f −2f +f) i Δ2 i+2 i+1 i
rdot(i)=(radius(i+1)-radius(i-1))/(2*Delta); thetadot(i)=(theta(i+1)-theta(i-1))/(2*Delta);
rdot(N)=(radius(N)-radius(N-1))/Delta; thetadot(N)=(theta(N)-theta(N-1))/Delta;
Central differencing schemes
f i ′ = Δ1 ( f i − f i − 1 )
Forward differencing schemes
rdoubledot(i)=(radius(i+1)-2*radius(i)+radius(i-1))/(Delta.^2); thetadoubledot(i)=(theta(i+1)-2*theta(i)+theta(i-1))/(Delta.^2);
rdoubledot(N)=(radius(N-2)-2*radius(N-1)+radius(N))/(Delta.^2); thetadoubledot(N)=(theta(N-2)-2*theta(N-1)+theta(N))/(Delta.^2);
ur=rdot; utheta=radius.*thetadot;
ar=rdoubledot-radius.*thetadot.^2; atheta=radius.*thetadoubledot+2.*rdot.*thetadot;
v ⃗ = r· e r ⃗ + r θ· e θ ⃗
fi′= 1(fi+1−fi−1) 2Δ
f′′=1(f −2f+f ) i Δ2 i−1 i i+1
f′′=1(f −2f +f) i Δ2 i−2 i−1 i
Backward differencing schemes
···2r ····θ a ⃗ = (r − rθ )e ⃗ + (rθ + 2rθ)e ⃗
32
N=length(t) %number of data points for i=2:N-1
end
figure(3)
subplot(221)
plot(t,ur,'ko-') xlabel('t(s)');ylabel('U_r(m/s)') subplot(222)
plot(t,utheta,'ko-') xlabel('t(s)');ylabel('U_{\theta}(m/s)') subplot(223)
plot(t,ar,'ko-') xlabel('t(s)');ylabel('A_r(m/s^2)') subplot(224)
plot(t,atheta,'ko-') xlabel('t(s)');ylabel('A_{\theta}(m/s^2)')
33
figure(3)
subplot(221)
plot(t,ur,'ko-') xlabel('t(s)');ylabel('U_r(m/s)') subplot(222)
plot(t,utheta,'ko-') xlabel('t(s)');ylabel('U_{\theta}(m/s)') subplot(223)
plot(t,ar,'ko-') xlabel('t(s)');ylabel('A_r(m/s^2)') subplot(224)
plot(t,atheta,'ko-') xlabel('t(s)');ylabel('A_{\theta}(m/s^2)')
End of Example L13.2
34
L13.3:
Derivatives of Unequally Spaced Data
35
β = 0.2
β = 0.5
Flow Visualisation: Contours of span wise vorticity
Unequally spaced grid used for simulations
Differentiation Unequally Spaced Grids
In order to find a formula to differentiate a function on an unequally spaced grid, consider the Taylor series 1
where
3!
xi+1 = xi + Δ
df dx
f(xi+1) = f(xi)+f′(xi)(xi+1 − xi) + 2! f′′(xi)(xi+1 − xi)2 + 1 f′′′(xi)(xi+1 − xi)3 + H . O . T
(13.11)
f′(x) = 38
Differentiation
Re-arrange Eq. (13.11) to obtain
f(xi+1) − f(xi) (xi+1 − xi)
1 2!
f′′(xi)(xi+1 − xi)
f′(xi) = −
+ 1 f′′′(xi)(xi+1 − xi)2 + ⋯
(13.12)
(13.13)
3! Hence f′(xi) can be approximated as
f′(xi) =
f(xi+1) − f(xi) (xi+1 − xi)
This is called a Forward Difference Approximation (FDA). Higher-order terms neglected are error in approximating f′(xi). Leading term in truncation error is
EFDA =
39
1 2!
f′′(xi)(xi+1 − xi)
A Taylor series expansion can also be used to obtain an expression for f(xi−1). ( )
f(xi−1) = f(xi) − f′(xi) xi − xi−1
+ 1 f′′(xi)(xi − xi−1)2 − 1 f′′′(xi)(xi − xi−1)3⋯
(13.14)
2! 3!
Rearranging gives
f′(x)=f(xi)−f(xi−1)+1f′′(x)(x−x ) i (xi − xi−1) 2! i i i−1
− 1 f′′′(xi)(xi − xi−1)2 + ⋯ 3!
40
Differentiation
f′(xi) = + 1 f′′(xi)(xi − xi−1) 2!
− 1 f′′′(xi)(xi − xi−1)2 + ⋯ 3!
Hence we can get yet another approximation for f′(xi) s f′(xi) =
(13.15)
f(xi) − f(xi−1) (xi − xi−1)
This is called Backward Difference Approximation (BDA). Leading term in truncation error for the BDA is
EBDA= 1 f′′(xi)(xi−xi−1) 2! 41
f(xi) − f(xi−1) (xi − xi−1)
(13.16)
Differentiation
We can find a more accurate estimate of f′(xi) on an unequally spaced grid
(13.11)
Subtract Eq. (13.14) from Eq. (13.11) to get:
f(xi+1) = f(xi)+f′(xi)(xi+1 − xi) + 1 f′′(xi)(xi+1 − xi)2
1 2!
+ 3! f′′′(xi)(xi+1 − xi)3 + ⋯
f(xi−1) = f(xi)−f′(xi)(xi − xi−1) + 1 f′′(xi)(xi − xi−1)2 1 2!
− 3! f′′′(xi)(xi − xi−1)3 + ⋯ 42
(13.14)
Differentiation
Subtract Eq. (13.14) from Eq. (13.11) to obtain:
f(xi+1) − f(xi−1) xi+1 − xi−1
f′(xi) =
This formula is called the Central Difference Approximation (CDA)
(13.17)
and its leading order error is given by
E = f′′(x ) (xi+1 − xi)2 − (xi − xi−1)2 CDA i 2!(xi+1 − xi−1)
+f′′′(x ) (xi+1 − xi)3 + (xi − xi−1)3
i 3!(xi+1 − xi−1) 43
Differentiation
Assuming all the xi's are equally spaced, i.e. xi+1 − xi = Δ = const, the FDA, BDA and CDA simplifies to
f(xi+1) − f(xi) Δ
f(xi) − f(xi−1) Δ
f(xi+1) − f(xi−1) 2Δ
and it can be shown that leading error term is O(Δ) in the FDA and BDA, and O(Δ2) in the CDA.
f′(xi) = f′(xi) = f′(xi) =
44
Error of the order of the grid spacing, xi+1 − xi.
Error of the order of the grid spacing, xi − xi−1 and xi+1 − xi. However, if
xi − xi−1 = xi+1 − xi = Δ, then error is of the order of Δ2.
EFDA = 1 f′′(xi)(xi+1 − xi) 2!
E = f′′(x ) (xi+1 − xi)2 − (xi − xi−1)2 CDA i 2!(xi+1 − xi−1)
+f′′′(x ) (xi+1 − xi)3 + (xi − xi−1)3 i 3!(xi+1 − xi−1)
EBDA= 1 f′′(xi)(xi−xi−1) 2!
45
Error of the order of the grid spacing, xi − xi−1.
Actual gradient
f(xi) f (xi+1)
f (xi−1)
xi−1
xi
xi+1
x
f′(xi) =
BDA
f(xi) − f(xi−1) (xi − xi−1)
f′(xi) =
CDA
f(xi+1) − f(xi−1) xi+1 − xi−1
f′(xi) =
FDA
f(xi+1) − f(xi) (xi+1 − xi)
46
Results from the above exercise illustrates that the most accurate scheme is the
CDA because the error in the CDA approximation is O(Δ2). This mean that if we halve the grid, the error will reduce to a quarter its original value. For the
FDA and BDA approximations, halving Δ will only halve the error.
In most practical engineering calculations, it is inefficient to have equally spaced
grid. Take for example, in the situation illustrated in the next example.
47
Example L13.3: Suppose we have a function
f(x) = e−x2 sin(10x)
Use (CDA, FDA or BDA) to calculate the derivative of this function using equally and unequally spaced grid Note that all the 'action' of this function is confined to the region 0 ≤ x ≤ 3. Thus, if we want to use the difference schemes to calculate the derivatives,
48
Lecture13M03.m
f′(x)= f(xi+1)−f(xi) i (xi+1 − xi)
f′(xi) = f(xi+1) − f(xi−1) xi+1 − xi−1
df(i)=(f(i+1)-f(i-1))/(x(i+1)-x(i-1));
df(N)=(f(N)-f(N-1))/(x(N)-x(N-1))
f′(x)= f(xi)−f(xi−1) i (xi − xi−1)
clear all
xveryfine=linspace(0,5,1000);
fxveryfine=exp(-xveryfine.^2).*sin(10*xveryfine); dfactualveryfine=exp(-xveryfine.^2).*(-2*xveryfine.*sin(10*xveryfine)+10*cos(10*xveryfine))
x=GenGrid(0,5,30,1.0); fx=exp(-x.^2).*sin(10*x); df=DerivCDA(fx,x);
figure(1)
subplot(211)
hold off plot(x,fx,'bo-',xveryfine,fxveryfine,'k-') xlabel('x');ylabel('f(x)')
subplot(212)
hold off plot(x,df,'bo-',xveryfine,dfactualveryfine,'k-') xlabel('x');ylabel('df(x)/dx')
function df=DerivCDA(f,x) df=zeros(size(f)) N=length(f) df(1)=(f(2)-f(1))/(x(2)-x(1)) for i=2:N-1
end end
function x=GenGrid(a,b,N,r) x=linspace(a,b,N);
if abs(r-1) > 1.0e-8
dx0=(b-a)*(1-r)/(1-r.^(N-1)); x(1)=a;
dx=dx0;
for i=2:N
x(i)=x(i-1)+dx;
dx=r*dx;
end end
end
49
Equally spaced grid
Unequally spaced grid:
More grid points clustered closer to x = 0
50
End of Example L13.3
51
Approximating 2nd Derivative:
fi+1 xi−1 xi−1/2 xi xi+1/2 xi+1
fi−1 fi
A formula for the 2nd derivative can be obtained by evaluating f ′(x) at
points halfway between xi+1 and xi and between xi and xi−1.
d2f dx2 i
df (xi+1/2) − df (xi−1/2) dx dx
(x) =
(xi + 12(xi+1 −xi))−(xi−1 + 12(xi −xi−1))
(13.18)
Using Eq. (13.14) (CDA) to approximate the first derivatives, f′
Eq. (13.18) can be written as
d2f (xi) = −
d x 2 12 ( x i + 1 − x i − 1 )
and f′ , i+1/2
i−1/2
52
f′ i+1/2
f′ i−1/2
f(xi+1) − f(xi) xi+1 − xi
f(xi) − f(xi−1) xi − xi−1
Gives us the following formula for approximating the 2nd derivative on an unequally spaced grid.
d2f (xi) = f(xi+1)(xi − xi−1) − f(xi)(xi+1 − xi−1) + f(xi−1)(xi+1 − xi) dx2 12 (xi+1 − xi−1) (xi+1 − xi) (xi − xi−1)
In addition, if one assumes that xi+1 − xi = xi − xi−1 = Δ, then
2
(13.19)
f(x )−2f(x)+f(x ) i+1 Δi i−1
f′′(xi) =
This formula is again a Central Difference Approximation (CDA)
and its leading order error is O(Δ2). 53