Blasius Equation
Blasius Equation
Blasius Equation
A. Ooi
October 14, 2005
Contents
1 Ordinary Differential Equations
1.1 Eulers method . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Backward (Implicit) Euler method . . . . . . . . . . . . . . . .
1.3 Trapezoidal or Crank-Nicolson method . . . . . . . . . . . . . .
1.4 Linearization of Crank-Nicolson method . . . . . . . . . . . . .
1.5 Runge-Kutta methods . . . . . . . . . . . . . . . . . . . . . . .
1.5.1 Second Order Runge-Kutta Method . . . . . . . . . . . .
1.5.2 4th Order Runge-Kutta Scheme (RK-4) . . . . . . . . .
1.6 Stability and error analysis . . . . . . . . . . . . . . . . . . . . .
1.7 Systems of Ordinary Differential Equations . . . . . . . . . . . .
1.8 Runge-Kutta-Fehlberg method: Runge-Kutta with error control
1.9 Boundary Value problem . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
9
10
11
13
14
16
17
24
28
32
2 Fourier series
37
2.1 Some properties of the Fourier coefficients . . . . . . . . . . . . . . . 40
2.2 Aliasing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3 Finite Difference
3.1 General Finite Difference Schemes . . . . . . . . . . . . . . . .
3.1.1 First Derivatives . . . . . . . . . . . . . . . . . . . . .
3.1.2 Some popular differencing schemes . . . . . . . . . . .
3.1.3 Higher order derivatives . . . . . . . . . . . . . . . . .
3.1.4 Summary of finite difference formula . . . . . . . . . .
3.2 Centred Difference Schemes . . . . . . . . . . . . . . . . . . .
3.2.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Solving PDEs using finite difference schemes . . . . . . . . . .
3.4 Fourier Analysis of error . . . . . . . . . . . . . . . . . . . . .
3.4.1 Fourier analysis of central differencing scheme . . . . .
3.5 Stability analysis using the modified wavenumber . . . . . . .
3.6 Dispersion-Relation-Preserving Scheme . . . . . . . . . . . . .
3.7 General Finite Difference Schemes For The Second Derivative
3.7.1 Some popular differencing schemes . . . . . . . . . . .
3.8 Multidimensional problems . . . . . . . . . . . . . . . . . . . .
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
47
47
47
48
51
52
53
55
56
60
65
68
71
72
73
75
CONTENTS
3.8.1 Steady problem . . . . . . . . . . . . .
3.8.2 Unsteady problem . . . . . . . . . . .
3.8.3 Modified wavenumber stability analysis
3.9 Test Problems . . . . . . . . . . . . . . . . . .
3.10 Euler equations . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
75
77
81
83
83
95
95
96
98
6 Collocation Method
101
6.1 Matrix operator for Fourier spectral numerical differentiation . . . . . 101
7 Some numerical examples
103
7.1 Heating of an oil droplet in air . . . . . . . . . . . . . . . . . . . . . 103
7.2 Blasius solution:
Contributed by Mr. M. Giacobello . . . . . . . . . . . . . . . . . . . 109
Chapter 1
Ordinary Differential Equations
In many engineering problems, you will need to solve differential equations that look
something like
dx
= f (t, x)
(1.1)
dt
in the domain
atb
with the initial condition
x(t = a) =
1.1
Eulers method
Eulers method is probably the simplest method used to solve Eq. (1.1). Consider
Taylors theorem
dx
(tn+1 tn )2 d2 x
(tn ) +
(n )
dt
2
dt2
where is somewhere in between tn and tn+1 .
x(tn+1 ) = x(tn ) + (tn+1 tn )
(1.2)
xn+1 = xn + hf (tn , xn )
(1.4)
(a)
2.5
x(t)
x(t)
1.5
1.5
0.5
0.5
(b)
2.5
10
10
10
(c)
2.5
(d)
2.5
x(t)
x(t)
1.5
1.5
0.5
0.5
10
Figure 1.1: Figure showing the solution to Exercise 1.2. h = 2.0 (a), h = 1.0 (b),
h = 0.5 (c), h = 0.1 (d).
exact solution, computed solution.
The solution to Exercise 1.2 is shown in Figure 1.1. Note that the numerical solution
gets closer to the exact solution for smaller values of h. By comparing, Eqs. (1.3)
and (1.4), one can conclude that if the exact solution is known at time t = tn , the
the numerical solution at time t = tn+1 will have an error which is O(h2 ). One can
say that the Eulers method has a local truncation error of order 2. This is because
if we reduce h by 2, then it can be expected that the error will be approximately
reduced by 4. However, it can be shown that (see [3] and [1] and Exercise 1.3 ) the
errors actually accumulate over time. This is illustrated in Fig. 1.2 which shows
the error associated with the Eulers method when the solution is computed using
h = 0.2, 0.1 and 0.05. Figure 1.2 (a) plots the error
Eh (tn ) = |x(tn ) xn |
(1.5)
for various values of h. The gaps between all three curves are the same which
indicates that the error is proportional to h. This fact is confirmed in Fig. 1.2(b)
which shows the ratio of the errors, Eh=0.1 /Eh=0.2 error Eh=0.05 /Eh=0.2 . From the
graphse, it is clear that the error is approximately halved when the h is halved.
Thus, the global error is shown to be O(h) i.e. Eulers method is only first order
accurate. In summary, Eulers method is the easiest method to implement but it is
not very accurate.
10
(a)
!1
10
Error
!2
10
!3
10
!4
10
!5
10
0.5
1.5
2.5
3.5
4.5
4.5
(b)
Error ratio
0.8
0.6
0.4
0.2
0.5
1.5
2.5
3.5
Figure 1.2: Figure showing the error associated to the solution to Exercise 1.2. The
, h = 0.1
,h = 0.05
. (b) shows
error is shown in (a) with h = 0.2
the ratio between the error,
Eh=0.1 /Eh=0.2 ,
Eh=0.05 /Eh=0.2
Exercise 1.3: Subtract Eq. (1.4) from Eq. (1.3) and show that
h2 00
x (n ).
(1.6)
2
where n = x(tn ) xn , the error at time t = tn . For the case where f (t, x) = 2t,
use the above equation and show that the error at time, tn = nh is O(h).
n+1 = n + h(f (tn , x(tn )) f (tn , xn )) +
Note: Even though you have only proven that the error is O(h) for the specific
case where f (t, x) = 2t, this result is more general and applies to all cases where
f (t, x) is a smooth function.
1.2
In many applications (see later) the Euler method described in section 1.1 is very
unstable. A more stable method is the backward Euler scheme given by the following
formula:
xn+1 = xn + hf (tn+1 , xn+1 )
(1.7)
The Backward Euler method is an implicit method because it has xn+1 on both
sides of the equation. For simple problems, there is generally no real difference
between the implicit Eulers method and the more conventional explicit Eulers
method because it is possible to obtain an explicit expression for xn+1 from Eq.
(1.7). As an example, for the problem discussed in Exercise 1.2, one can use Eq.
(1.7) to show that
xn+1 =
h + xn
1+h
(1.8)
Using the Eq. (1.8), one can generate the results shown in Fig. 1.3. Comparing
Figs. 1.3(a) and 1.1(a), it is clear that the Backward Euler method gives a better
solution h = 2 in this situation. This is one property of the implicit Eulers method.
For a given value of h, it usually gives a better solution than the explicit Eulers
method.
3
2.5
x(t)
x(t)
1.5
1.5
0.5
0.5
(b)
2.5
(a)
10
10
Figure 1.3: Figure showing the solution to Exercise 1.2 using the Backward Euler
method. h = 2.0 (a), h = 1.0 (b).
exact solution, computed solution.
10
(1.9)
with the initial condition x(0) = 0.0. The exact solution to Eq. (1.9) is
2t
e 1
x(t) = arctan
2et
(1.10)
We will now show how the backward Euler method can be used to obtain an
approximate solution to Eq. (1.9). Applying Eq. (1.7) to Eq. (1.9) gives
xn+1 h cos(xn+1 ) xn = 0.
(1.11)
xn is known from the previous time step. One then needs to find xn+1 . It is not
easy to find an explicit expression for xn+1 . Root finding methods (see 2nd year
computational mechanics lecture notes) must be used to find xn+1 numerically. If
we used the Newton-Raphson method, the iterative scheme could be written as
(k+1)
(k+1)
(k+1)
(k)
xn+1 = xn+1
xn+1 h cos(xn+1 ) xn
(k+1)
1 + h sin(xn+1
(k)
where xn+1 is the kth guess for the value of xn+1 that satisfies Eq. (1.11). While
(k)
(0)
one may use any value for xn+1 , it would be logical to use xn+1 = xn .
Figure 1.4 shows the numerical solution to Eq. (1.9) obtained using both the
implicit and the explicit Eulers method. For large h = 2.0, the solution obtained
using explicit Euler is oscillatory even though the exact solution does not oscillate.
The backward Euler gives a solution that is in better agreement with the exact
solution. For h = 0.5 both implicit and explicit Euler solution converges to the
exact analytical solution.
1.3
The integration can be approximated using the trapezoidal rule (see 2nd year, Computational mechanics lecture notes)
Z tn+1
h
f (t, x(t)) = (f (tn , x(tn )) + f (tn+1 , x(tn+1 ))) + O(h3 )
(1.12)
2
tn
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
x(t)
x(t)
1
0.8
11
1
0.8
(b)
0.6
0.6
(a)
0.4
0.4
0.2
0.2
10
10
Figure 1.4: Figure showing the solution to Example 1.1 for h = 2.0 (a) h = 0.5 (b).
exact solution,
explicit Euler and
backward Euler solution.
where h = tn+1 tn = h. Ignoring the error term, one will obtain the following
formula for the Crank-Nicolson (trapezoidal) method
xn+1 = xn +
h
(f (tn , xn ) + f (tn+1 , xn+1 ))
2
(1.13)
By comparing Eqs. (1.12) with (1.13) one can conclude that the error associated
with the Crank-Nicolson method is O(h3 ). However, similar to the Eulers method,
the error accumulate over time. Thus, the global error of the Crank-Nicolson method
is O(h2 ).
Exercise 1.4: Solve
dx
+ 2x = 0
dt
for 0 < t < 3 with x(t = 0) = 1 using
(a) Euler method
(b) Backward Euler method
(c) Crank Nicolson method
1.4
The main problem with implicit method is that they are difficult to solve, i.e. it
is difficult to obtain xn+1 from the resulting discretized equations. Sometimes, it is
possible to linearized implicit methods. Consider the Crank-Nicolson scheme
xn+1 = xn +
h
(f (tn , xn ) + f (tn+1 , xn+1 )) + O(h3 )
2
(1.14)
12
The difficulty comes from the term f (tn+1 , xn+1 ). Lets consider a Taylor series
expansion about of f (tn+1 , xn+1 ) about xn
2f
f
(xn , tn+1 )+(xn+1 xn )2 2 (xn , tn+1 )+HOT
x
x
(1.15)
f
(xn , tn+1 ) + O(h2 )
x
(1.16)
xn+1
h
= xn +
2
f
f (tn+1 , xn ) + (xn+1 xn ) (xn , tn+1 ) + f (tn , xn )
x
(1.17)
(1.18)
This method has got good stability characteristics (see later). The main problem
with this scheme is that you have to find the derivative of f with respect to x. This
is not always possible.
Example 1.2: Consider again the example
dx
= cos(x)
dt
Applying the Crank Nicolson method will give an implicit expression for xn+1
h
h
cos(xn+1 ) xn cos(xn ) = 0
(1.19)
2
2
Since xn is known from the previous time step, it is possible to use root finding
methods to find a solution to Eq. (1.19). On the other hand, applying the linearized
Crank-Nicolson method (Eq. (1.18)) gives an explicit expression for xn+1
xn+1
xn+1 = xn +
h cos(xn )
1 + (h/2) sin(xn )
(1.20)
13
As you can probably imagine, writing a computer program to implement Eq. (1.19)
is a lot simpler than writing a computer program to implement Eq. (1.20). The
result for both the implicit and linearized Crank-Nicolson methods are shown in Fig.
1.5. As you can see, both methods give very similar results.
2
1.8
1.6
1.4
x(t)
1.2
1
0.8
0.6
0.4
0.2
0
10
t
Figure 1.5: Comparing the linearized and implicit Crank-Nicolson methods for h =
1.
linearized Crank-Nicolson and
implicit Crank-Nicolson method.
1.5
Runge-Kutta methods
These methods are probably the most popular methods in solving initial value problems. However, many variations of the Runge-Kutta methods exist. We will derive
one set below and see why the various formulae are not unique. Lets start with the
problem we would like to solve
dx
= f (t, x)
dt
In general, the Runge-Kutta schemes can be written as
xn+1 = xn + (xn , tn , h)h
(1.21)
14
= a1 k1 + a2 k2 + a3 k3 + a4 k4 + aN kN
(1.22)
where
k1
k2
k3
k4
..
.
..
.
..
.
kN
=
=
=
=
..
.
..
.
..
.
=
f (tn , xn )
f (tn + p1 h, xn + q11 k1 h)
f (tn + p2 h, xn + q21 k1 h + q22 k2 h)
f (tn + p3 h, xn + q31 k1 h + q32 k2 h + q33 k3 h)
..
.
..
.
..
.
f (tn + pN 1 h, xi + qN 1,1 k1 h + qN 1,2 k2 h + + qN 1,N 1 kN 1 h)
For N = 1, we get the first order Runge-Kutta scheme. This is just the same as the
Euler integration scheme presented earlier.
1.5.1
xn+1 = xn + (a1 k1 + a2 k2 ) h
xn+1 = xn + a1 f (tn , xn )h + a2 f (tn + p1 h, xn + q11 k1 h)h
(1.23)
Thus in order to find a numerical scheme, we need to find values for the following
constants, a1 , a2 , p1 and q11 . The last term in the above equation is a 2 variable
function and the Taylor series expansion (to the lineariazed approximation) for a
two variable function is given by
f
f
+
+ HOT
t
x
Using this relationship and ignoring the higher order terms (HOT),
f (t + h, x + ) = f (t, x) + h
f
f
(tn , xn ) + q11 k1 h (tn , xn )
t
x
(1.24)
15
k1 = f (tn , xn )
Hence Eq. (1.24) can be written as
f
f
+ q11 f (tn , xn )h
t
x
(1.25)
f 2
f
h + a2 q11 f (tn , xn ) h2
t
x
(1.26)
d2 x
h2
dx
(tn ) h + 2 (tn )
+ HOT
dt
dt
2!
where HOT stands for higher order terms. Lets assume that they are small. So the
above equation becomes
xn+1 = xn +
dx
d2 x h2
h+ 2
dt
dt 2!
(1.27)
df (tn , xn ) h2
dt
2
f
f
dt +
dx
t
x
df
f
f dx
=
+
dt
t
x dt
df =
(1.28)
Hence
xn+1 = xn + f (tn , xn )h +
1 f 2 1 f
h +
f (tn , xn )h2
2 t
2 x
(1.29)
Comparing Eqs. (1.29) with Eq. (1.26) will give you the following three equations
16
a1 + a2 = 1
1
a2 p 1 =
2
1
a2 q11 =
2
(1.30)
We have four unknowns (a1 , a2 , p1 and q11 ) but only the above three equations.
So there cannot be a unique solution. You have an infinite number of solutions. One
possible solution is to set
1
a2 =
2
then
1
2
= 1
= 1
a1 =
p1
q11
(1.31)
Hence, one possible second order Runge-Kutta time stepping scheme is
1
1
xn+1 = xn +
k1 + k2 h
2
2
(1.32)
where
k1 = f (tn , xn )
k2 = f (tn + h, xn + hk1 )
(1.33)
1.5.2
By far, the most popular numerical method for solving ODE is the 4th order RungeKutta scheme
1
1
1
xn+1 = xn +
k1 + (k2 + k3 ) + k4 h
(1.34)
6
3
6
where
17
k1 = f (tn , xn )
h
h
k2 = f tn + , xn + k1
2
2
h
h
k3 = f tn + , xn + k2
2
2
k4 = f (tn + h, xn + hk3 )
1.6
When discussing stability of the numerical methods, one usually considers the model
problem
dx
= x
(1.35)
dt
is a constant which can be a complex number. In most engineering problems, the
real part of is usually negative. This means that the solution that you are after
will typically decay with time. Applying the Euler method (with timestep t = h)
to Eq. (1.35) gives
xn+1 = xn + hxn
(1.36)
18
is usually called the amplification factor. If || < 1 then the error will decay with
time. The opposite is true if || > 1. Hence, in order to ensure the stability of the
numerical method, || < 1.
||2 = (1 + hR )2 + (hI )2 < 1
This is just a circle of radius 1 centred on (-1,0). This plot is called the stability
diagram and is shown in Fig. 1.6
Ih
Rh
If is real and negative, then in order for the numerical method to be stable,
h
2
||
(1.38)
(1.39)
19
where I is the imaginary unit and is a real number. The model problem for
stability, Eq. (1.35) becomes
dx
= Ix
dt
For illustrative purposes, lets use x(t = 0) = 1. The exact solution to this problem
is
x(t) = x0 eIt
(1.40)
where x0 = x(t = 0). If the Euler method is used to solve this equation, we know
that the amplitude will grow with time as the value of is not within the stability
region. Lets now analyse the phase error. The amplification factor, for this
problem can be written as
= 1 + Ih
= ||eI
(1.41)
where
= arctan(h)
Lets now compare the exact solution at time t = nh
x(t = nh) = x0 eInh
with the approximated solution given by the Eulers equation.
xn = n x0 = An eIn x0
Dividing the two equations and a little bit of algebra will give you
xn = x(nh)An eIn(h)
So the approximated solution that you will get will be amplified by An and its
phase will be shifted by n( h). Thus A is usually called the amplitude error and
( h) is called the phase error.
At any time step, h, the error associated with the phase is h
h = h arctan(h)
(h)3 (h)5 (h)7
= h h
+
+ ......
3
5
7
(h)3 (h)5 (h)7
=
+
+ ......
3
5
7
(1.42)
20
So the phase error has a leading order term of (h)3 which is very small. In contrast,
the amplitude error grows like An .
We can also apply the stability analysis to the model equation using the implicit
Euler method
xn+1 = xn + xn+1 h
Rearranging the above equation gives
xn+1 =
xn
1 h
hence
xn = n x0
where
1
1 h
Exercise 1.6: Show that the amplification factor, for the implicit Euler scheme
can be written as
= Aei
where
1
A= p
(1 R h)2 + (I h)2
and
= arctan
I h
1 R h
Since we are only interested in problems which have R < 0, A is always smaller than
unity. Thus the implicit Euler scheme is always stable (as long as R < 0), no matter
what the value of h. This means that the implicit Euler scheme is unconditionally
stable, which is a property of most implicit numerical scheme.
21
Exercise 1.7: Compare the implicit and explicit Euler scheme and show that for
the model problem
dx
= Ix,
dt
the two schemes have the different amplitude error but the phase error is the
same. Show that the magnitude of the solution given by the explicit Euler scheme
is growing and the magnitude of the solution given by the implicit Euler scheme
decreases with time. The solution to this exercise with h = 0.1 is shown in Fig.
1.7.
5
4
3
2
x(t)
1
0
!1
!2
!3
!4
!5
10
12
14
16
18
20
exact solution,
22
Exercise 1.8:
(a) For the model problem
dx
= x,
dt
perform a stability analysis on the Crank-Nicolson method and show that it
is unconditionally stable, i.e. the method is stable for all values of h.
(b) Show that the linearized Crank-Nicolson and the implicit Crank-Nicolson
methods have the same stability characteristics.
(c) Also, show that when using the Crank-Nicolson to solve the model problem
dx
= Ix,
dt
there is no amplitude error associated with the numerical method. Show that
the phase error has a leading order term that looks like
PE =
(h)3
+ ...
12
Exercise 1.9: Perform stability analysis for the second order Runge-Kutta
method on the model equation
dx
= x
dt
and show that the region of stability can be obtained by solving the following
equation
2 h2
ei = 0
(1.43)
2
Show also that for the 4th order Runge-Kutta method, the stability region is
obtained by solving
1 + h +
2 h2 3 h3 4 h4
+
+
+ 1 ei = 0
2
6
24
The solutions to Eqs. (1.43) and (1.44) are shown in Fig. 1.8.
h +
Example 1.3:
For the model problem
(1.44)
23
!I h
!1
!2
!3
!3
!2
!1
!R h
dx
= I,
dt
The solution given by the numerical schemes can be written as
xn+1 = xn
where
|Euler | =
1 + (h)2
(1.45)
(h)4
4
(1.46)
(h)6 (h)8
+
72
576
(1.47)
|RK2 | =
r
|RK4 | =
1+
24
for the Eulers, second order Runge-Kutta (RK2) and fourth order Runge-Kutta
(RK4) schemes respectively. A plot of the magnitude of for the various schemes
is shown in Fig. 1.9. It is seen that the magnitude of is very close to 1 for both
the RK4 and RK2 schemes while it diverges away from 1 very quickly for the Euler
scheme. Thus, the magnitude error for the RK4 and RK2 schemes (especially the
RK4) is very small.
To get an idea how small the error is, if we take h = 0.5, then
Euler =
1 + 0.52 = 1.118
r
RK2 =
1+
0.54
= 1.0078
4
0.56 0.58
+
= 0.9998948784
72
576
These numbers might seem small but consider this. If we carry out simulations with
h=0.5, then after 100 time steps, the RK4 scheme produce a numerical solution
whose magnitude is 98.9% of the true solution, the RK2 scheme will give you a
solution that is 217% of the true solution, and the Eulers method will be out by
7000000% (yes 7 million percent) !
It is interesting to point out here that since RK4 is given by Eq. (1.47), the
stability curve for the 4th order Runge-Kutta scheme (shown in Fig. 1.8) intercepts
the (I h axis when
RK4 =
I h = h
r
576
=
72
8
=
= 2.8284
1.7
In many engineering problems, it is essential to solve, not just one, but a set of
ordinary differential equations. This problem can be expressed in matrix-vector
notation as
d
{x} = [A] {x}
dt
(1.48)
25
1.05
1.04
1.03
1.02
"
1.01
1
0.99
0.98
0.97
0.96
0.95
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
!h
Figure 1.9: Figure showing the magnitude of for Euler (
RK4 (
)
), RK2 (
) and
where the curly brackets ({}) denotes a vector and the square brackets [] represents
a matrix. {x} is the vector of dependent variables and [A] is usually a matrix of
constant coefficients. Equation (1.48) can be solved using the methods described
above. Applying the explicit Euler method
{x}n+1 {x}n
= [A] {x}n
h
{x}n+1 = {x}n + h [A] {x}n
At every time step, one would need to perform a matrix-vector multiplication in
order to obtain the values of the vector {x} at the next time step.
If an implicit scheme is required to solve the problem, then one would need
to solve an system of linear algebraic equations at every time step. For example,
applying the trapezoidal/Crank Nicolson scheme to Eq. (1.48) gives
26
{x}n+1 {x}n
1
1
=
[A] {x}n + [A] {x}n+1
h
2
2
h
h
n+1
n+1
n
{x}
[A] {x}
= {x} + [A] {x}n
2
2
h
h
I [A] {x}n+1 = I + [A] {x}n
2
2
(1.49)
Equation (1.49) can now be solved using methods for solving a system of differential (e.g. LU decomposition, Gauss-Seidel etc.)
Exercise 1.10: Solve the following second order ordinary differential equation
d2 x
+ 2x = 0
dt2
(1.50)
x (t = 0) = p
dx
(t = 0) = 0.
dt
Outline the numerical algorithm for the Eulers and Trapezoidal method.
Multi-dimensional RK-4
In many situations, you might be required to solve a multi-dimensional system using
RK-4. The RK-4 can easily be extended for a multi-dimensional system. An example
of using RK-4 to solve a system with 4 dependent variables is given below.
In the example below, it will be assumed that the system you are interested in
is autonomous, i.e. the right hand side is not a function of time. Suppose you are
asked to solve a system of ordinary differential equations that looks like
da1
dt
da2
dt
da3
dt
da4
dt
= f1 (a1 , a2 , a3 , a4 )
= f2 (a1 , a2 , a3 , a4 )
= f3 (a1 , a2 , a3 , a4 )
= f4 (a1 , a2 , a3 , a4 )
an+1
1
an+1
=
2
=
an+1
3
=
an+1
4
1
+
k11 +
6
1
n
a2 +
k12 +
6
1
n
a3 +
k13 +
6
1
n
a4 +
k14 +
6
an1
1
(k21 + k31 ) +
3
1
(k22 + k32 ) +
3
1
(k23 + k33 ) +
3
1
(k24 + k34 ) +
3
1
k41 h
6
1
k42 h
6
1
k43 h
6
1
k44 h
6
27
(1.51)
(1.52)
(1.53)
(1.54)
k11
k12
k13
k14
k21 = f1 a1 +
k22 = f2 a1 +
k23 = f3 a1 +
k24 = f4 a1 +
=
=
=
=
f1 (a1 , a2 , a3 , a4 )
f2 (a1 , a2 , a3 , a4 )
f3 (a1 , a2 , a3 , a4 )
f4 (a1 , a2 , a3 , a4 )
h
k11 , a2 +
2
h
k11 , a2 +
2
h
k11 , a2 +
2
h
k11 , a2 +
2
h
k12 , a3 +
2
h
k12 , a3 +
2
h
k12 , a3 +
2
h
k12 , a3 +
2
h
k13 , a4 +
2
h
k13 , a4 +
2
h
k13 , a4 +
2
h
k13 , a4 +
2
h
k14
2
h
k14
2
h
k14
2
h
k14
2
28
k31 = f1 a1 +
k32 = f2 a1 +
k33 = f3 a1 +
k34 = f4 a1 +
h
k21 , a2 +
2
h
k21 , a2 +
2
h
k21 , a2 +
2
h
k21 , a2 +
2
h
k22 , a3 +
2
h
k22 , a3 +
2
h
k22 , a3 +
2
h
k22 , a3 +
2
h
k23 , a4 +
2
h
k23 , a4 +
2
h
k23 , a4 +
2
h
k23 , a4 +
2
h
k24
2
h
k24
2
h
k24
2
h
k24
2
and
k41
k42
k43
k44
=
=
=
=
1.8
In this method, the step size, h is adjusted every time level in order to keep the
error of the numerical scheme under control. In order to see how this can be done,
recall that the N th order Runge-Kutta scheme can be written as,
xn+1 = xn + h(tn , xn , h)
(1.55)
(1.56)
n , xn , h)
xn+1 = xn + h(t
(1.57)
(1.58)
Lets assume that at time level, tn , there is hardly any error in the numerical solutions, i.e. xn xn x(tn ). Rearranging Eq. (1.56), one obtains an expression for
the local truncation error,
x(tn+1 ) x(tn )
(tn , x(tn ))
h
x(tn+1 ) xn
(tn , x(tn ))
=
h
x(tn+1 ) (xn + h(tn , x(tn )))
=
.
h
n+1 (h) =
Note that n+1 (h) is O(hN ). Since the Runge-Kutta scheme requires that the numerical value of = , we can continue the above as
x(tn+1 ) (xn + h(tn , xn ))
h
1
=
(x(tn+1 ) xn+1 )
h
n+1 (h) =
(1.59)
Similarly, starting from Eq. (1.58), one can obtain an expression for the local truncation
1
(x(tn+1 ) xn+1 )
h
which is O(hN +1 ). Going back to Eq. (1.59) and using Eq. (1.60)
n+1 (h) =
(1.60)
1
(x(tn+1 ) xn+1 )
h
1
([x(tn+1 ) xn+1 ] + [
xn+1 xn+1 ])
=
h
1
= n+1 (h) + [
xn+1 xn+1 ]
h
n+1 (h) =
Recall that n+1 (h) is O(hN ) and n+1 (h) is O(hN +1 ). Thus, if h is small, n+1 (h)
can be simply approximated as
n+1 (h)
1
(
xn+1 xn+1 )
h
(1.61)
30
Lets now see how this information can be used to control the local truncation error.
Since n+1 (h) is O(hN ), one can write
n+1 (h) = KhN
(1.62)
If we increase or decrease the time step h by a factor of r, then the local truncation
N
error would be n+1 (rh) = K(rh)N = rN KhN = rN n+1 (h) rh (
xn+1 xn+1 ),
using Eq. (1.61). Thus if we want to bound the local truncation error to a small
value , then
r
h
xn+1 xn+1
1/N
h
2
1/N
(1.63)
(1.64)
xn+1 xn+1
One popular method to implement the above algorithm is called the RungeKutta-Fehlberg method. In this method xn+1 and xn+1 are approximated as
25
1408
2197
1
xn+1 = xn + h
k1 +
k3 +
k4 k5
(1.65)
216
2565
4104
5
xn+1 = xn + h
16
6656
28561
9
2
k1 +
k3 +
k4 k5 + k6
135
12825
56430
50
55
(1.66)
where
k1 = f (tn , xn )
h
h
k2 = f tn + , xn + k1
4
4
3h
3h
9h
k3 = f tn + , xn + k1 + k2
8
32
32
12h
1932h
7200h
7296h
k4 = f tn +
, xn +
k1
k2 +
k3
13
2197
2197
2197
439h
3680h
845h
k5 = f tn + h, xn +
k1 8hk2 +
k3
k4
216
513
4104
h
8h
3544h
1859h
11h
k6 = f tn + , xn k1 + 2hk2
k3 +
k4
k5
2
27
2565
4104
40
It can be shown that the global error associated with xn+1 is O(h4 ) and the global
error associated with xn+1 is O(h5 ). So N = 4 and r is calculated as
r = 0.84
h
xn+1 xn+1
1/4
(1.67)
32
if t(n)>=tmax
FLAG=0;
elseif t(n)+h>tmax
h=tmax-t(n);
elseif h<hmin
%time step is too small. exiting program....
FLAG=0;
end
end
1.9
The methods introduced above can only be used to solve initial value problems.
Meaning that all the information you are given is at time t = 0, say, and you are
asked to predict the solution at a later point in time, say t = tf . What if you are
given some information at t = tf ? This kind of problems are called boundary value
problems. Can the methods above be used to solve boundary value problems? It
turns out that it is possible to use the methods introduced in the preceding chapters
to solve boundary value problems. These family of methods are called shooting
methods and in general they guess the information at t = 0 in order to give the
required conditions at t = tf . The method will be introduced by considering a
system of 3 ODEs. Generalization to a system of N ODEs is straightforward.
Consider the system of 3 ordinary differential equations.
da
= fa (a, b, c)
dt
db
= fb (a, b, c)
dt
dc
= fc (a, b, c)
dt
33
(1.68)
(1.69)
You are asked to find a(t), b(t) and c(t). It is important to note that c(0) is not
defined and for the purpose of the following discussion, let c(0) = . The idea behind
the shooting methods is that because c(0) is not defined, we are free to choose any
value for that will give us c(tf ) = P . But we do not know before hand what
value of will give you c(tf ) = P . So we need iterate through different values of
until the computed value of c(tf ) = P . As you can imagine, this is a very inefficient
method. There are more intelligent methods in choosing the value of that will
give you c(tf ) = P . This is described below.
Using any numerical scheme for ODEs, Eq. (1.68) can be solved with the following initial conditions
a(0) = 0, b(0) = 0, c(0) =
(1.70)
Lets say we take N steps of step size h to approach t = tf , then the approximate
value of c(tf ) can be denoted as cN . The value of cN is dependent on the value of
. Since we are only guessing the value of , it is very likely that cN P 6= 0. For
the following analysis, it will be convenient to define a function
g() = cN () P.
(1.71)
In order for the numerical solution to satisfy the original boundary conditions defined
in Eq. (1.69), we must ensure that
g() = 0.
Thus this becomes a root finding problem, i.e we have to iteratively find the value
of such that g() = 0. Each value will give you a numerical solution. Only
the numerical solution computed with the value of that ensures that g() = 0 is
the correct solution to the original problem. Since the problem has been recast as a
root finding problem, the Secant formula is usually used to provide a better guess
value of ,
i+1 = i
(i i1 )g(i )
g(i ) g(i1 )
(1.72)
34
i
g(i )
0
-6.7733
1
-5.9639
8.3682 3.2575
5.7654 0.6666
5.0957 -0.1424
5.2136 0.0085
5.2069 0.0001
The numerical solution computed using different values of i is shown in Fig. 1.10.
Note that the i are chosen according to Eq. (1.72) and it forces the numerical
35
10
x(t)
!3
!4
!1
!0
!2
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
t
Figure 1.10: Figure showing the solution to Example 1.4. The numerical solution
computed using 0 , 1 , 2 , 3 , 4 are indicated on the diagram.
36
Chapter 2
Fourier series
Any periodic function can be, f (x) can expressed as a complex Fourier series with
coefficients ck . To prove this, recall that in Engineering Analysis, any periodic
function can be expressed as a sin and cosine series
f (x) = a0 +
(2.1)
l=1
where
2
,
L
L is the length of the domain and k0 is the fundamental wavenumber. The term
k0 =
lk0 = kl
is the fundamental harmonics. The coefficients, al and bl in Eq. (2.1) can be found
by evaluating the following integrals
a0
al
bl
Z
1 L
=
f (x) dx
L 0
Z
2 L
=
f (x) cos (lk0 x) dx
L 0
Z
2 L
=
f (x) sin (lk0 x) dx
L 0
Note that
al = al
bl = bl
(2.2)
37
38
Recall that sin and cos functions can be expressed in terms of complex exponentials
sin =
eI eI
2i
(2.3)
cos =
eI + eI
2
(2.4)
Ilk0 x
eIlk0 x + eIlk0 x
e
eIlk0 x
f (x) = a0 +
al
+ bl
2
2I
l=1
X
al
bl
al
bl
Ilk0 x
Ilk0 x
= a0 +
e
+
+e
2
2I
2
2I
l=1
X
al Ibl
al + Ibl
= a0 +
eIlk0 x +
eIlk0 x
2
2
l=1
Define cl to be
c 0 = a0
al Ibl
cl =
2
al + Ibl
al Ibl
=
cl =
2
2
Substituting the definition for c into the expression for f (x) gives
f (x) = c0 +
cl eIlk0 x +
l=1
cl eIlk0 x
l=1
cl eIlk0 x +
l=0
cl eI(l)k0 x
l=1
Observe that
X
l=1
Thus
I(l)k0 x
cl e
X
l=1
cl eIlk0 x
39
f (x) =
=
cl eIlk0 x +
l=0
cl eIlk0 x
l=1
cl eIlk0 x
(2.5)
l=
(2.6)
f (x) =
cl eIlk0 x
(2.7)
l=
Exercise 2.1: Use the definition of al and bl above and show that
al Ibl
2
Z
1 L
=
f (x) eIlk0 x dx
L 0
cl =
(2.8)
cl =
al Ibl
2 Z
Z
2I L
1 2 L
f (x) cos (lk0 x) dx
f (x) sin (lk0 x) dx
2 L 0
L 0
Z
Z
1 L
I L
f (x) cos (lk0 x) dx
f (x) sin (lk0 x) dx
L 0
L 0
Ilk0 x
Z
1 L
e
+ eIlk0 x
f (x)
dx
L 0
2
Ilk0 x
Z
i L
e
eIlk0 x
f (x)
dx
L 0
2I
Z
1 L
f (x) eIlk0 x dx
L 0
We usually cannot sum to infinity - problem with Eq. (2.7) In general you cannot
do the integral specified in (2.8) exactly. The Nyquist critical wavenumber i.e. the
largest wavenumber that is resolvable by an equally sampled data in physical space
is
N
k0
2
(2.9)
40
1
f (xj ) =
N
2
X
l= N
2
cl e+Ilk0 xj
(2.10)
+1
Ik0 (lm)xj
j=0
=
if
l m = pN
otherwise
N
0
where
p = 0, 1, 2, ...
(2.11)
If we multiply both sides of Eq. (2.10) by eIk0 mxj and summing from j = 0 to
N 1 gives
cl =
N
1
X
f (xj ) eIlk0 xj
(2.12)
j=0
2.1
NP
1
j=0
NP
1
cl =
j=0
If f (x) is a real function, then both c0 and cN/2 are real numbers
If f (x) is a real function, then cN l = cl .
41
Proof:
cl =
NP
1
f (xi ) eIlk0 xj ,
j=0
NP
1
cl =
f (xi ) eIlk0 xj
j=0
cN l =
N
1
X
j=0
N
1
X
j=0
where
xj = j
L
= j
N
(2.13)
Lj
2
eiN k0 xj = eIN ( L )( N )
= eI2j
= cos (2j) + I sin (2j) = 1
| {z } | {z }
=1
=0
Hence,
cN l = cl
(2.14)
1
f (xj ) =
N
2
X
cl e+Ilk0 xj
l= N
+1
2
N 1
1 X Ilk0 xj
=
cl e
N l=0
42
1
f (xj ) =
N
2
X
cl e+Ik0 lxj
l= N
+1
2
1
2
X
1 X Ilk0 xj
=
cl e
+
cl eIlk0 xj
N l=0
+1
l= N
2
N
N
1
2
2
X
X
1
cl eIlk0 xj +
cl eIlk0 xj
=
N l=0
l=1
But we know that cl = cN l . So
f (xj ) =
f (xj ) =
1
cl eIlk0 xj +
N l=0
1
N
N
1
X
2
X
1
2
X
cN l eI(N l)k0 xj
l=1
cl eIlk0 xj
l=0
Note that
eI(N l)k0 xj = eIN k0 xj eIlk0 xj
= eIN k0 (Lj/N ) eIlk0 xj
= eIlk0 xj
(2.15)
Thus the discrete Fourier transform pair can be written as
N 1
1 X Ilk0 xj
f (xj ) =
cl e
N l=0
cl =
N
1
X
f (xj ) eIlk0 xj
(2.16)
(2.17)
j=0
43
Proof:
1
f (x) =
N
N/2
X
cl eIkl x
l=N/2+1
1
f (x + t) =
N
1
=
N
N/2
X
cl eIkl (x+t)
l=N/2+1
N/2
X
cl eIkl t eIkl x
l=N/2+1
N 1
1 X Ikl x
=
bl e
N k=0
f (xj ) =
2
X
cl e+Ikl xj
l= N
+1
2
N
2
2 X
= p0 +
(pl cos(kl xj ) ql sin(kl xj ))
N l=0
where kl = lk0 and pl and ql are the real and imaginary parts of cl respectively.
If we only consider situations where f is defined in the domain
0 < x < 2
i.e. L = 2,
k0 = 2/L
= 2/2
= 1
(2.18)
With this assumption, the discrete Fourier series
44
N
2
X
1
f (xj ) =
N
cl =
cl e+Ilxj
(2.19)
l= N
+1
2
N
1
X
f (xj ) eIlxj
(2.20)
j=0
(2.22)
where kl = lk0 .
If we set U = 1, the solutions at time t = 0 and t = 20 are shown in Fig.
2.2
Aliasing
For any sample interval , there is a special wave number kc called the Nyquist
critical wave number which is defined to be
2
2
kc =
2/N
N
=
2
kc =
kc is important for the following reasons. Any energy in wave number that exist
beyond kc will be reflected to k < kc with kc as the point of reflection.
2.2. ALIASING
45
0.8
f(x,t)
0.6
0.4
0.2
0.2
10
10
x
15
20
25
30
Figure 2.1:
initial condition at time t = 0.
Solution at t = 20 for = 0,
solution at t = 20 for = 0.01,
solution at t = 20 for = 0.1.
46
Chapter 3
Finite Difference
3.1
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.
3.1.1
First Derivatives
We would like to express the derivative of a function, f (x), in terms the function
values at neighbouring points. The general explicit formula to for doing this can be
written as
1
0
fi =
M
X
aj fij + a0 fi +
j=1
M
X
!
aj fi+j
(3.1)
j=1
48
fi =
+
+
+
+
1
(aM + aM +1 + ... + a1 + a0 + a1 + ...aM 1 + aM )f (xi )
1
df
(aM (M ) + ... + a1 (1) + a1 + ... + aM M )) (xi )
dx
d2 f
(M )2
(1)2
12
M2
1
aM
+ ... + a1
+ a1 + ... + aM
2 2 (xi )
2!
2!
2!
2!
dx
3
3
3
3
d2 f
(M )
(1)
1
M
1
aM
+ ... + a1
+ a1 + ... + aM
3 2 (xi )
3!
3!
3!
3!
dx
4
O( )
0
= 0
= 1
= 0
= 0
.
= ..
The more equations we choose to satisfy, the higher the accuracy of the approximation. For example, if we only satisfy the first two equations, then the accuracy
of the approximation is O(). If we satisfy the first three equations, then the
approximation is O(2 ).
3.1.2
Lets now see what happens when we choose M = 1. The finite diffrencing can be
written as
1
(a1 f1 + a0 f0 + a1 f1 )
0
The equations that we need to satisfy in order for fi to approximate df /dx(xi ) can
be written as
0
fi =
a1 + a0 + a1
a1 + a1
1
1
a1 + a1
2
2!
1
1
a1 + a1
3!
3!
= 0
= 1
(3.2)
(3.3)
= 0
(3.4)
= 0
(3.5)
49
a1 = 0
a0 = 1
a1 = 1
With this set of coefficients, the derivative of f can be approximated as
0
fi =
1
(fi+1 fi )
(3.6)
a1 = 0
a0 = 1
a1 = 1
With this set of coefficients, the derivative of f can be approximated as
0
fi =
1
(fi fi1 )
(3.7)
a1 =
1
2
a0 = 0
1
a1 =
2
with the corresponding formula being the approximation for f ,
0
fi =
1
(fi+1 fi1 ) .
2
(3.8)
50
Show that in order to obtain a fourth order accurate finite difference formula, one
must satisfy the following 5 equations
0
fi =
a2 + a1 + a0 + a1 + a2
2a2 a1 + a1 + 2a2
1
1
2a2 + a1 + a1 + 2a2
2
2
4
1
1
4
a2 a1 + a1 + a2
3
6
6
3
2
1
1
2
a2 + a1 + a1 + a2
3
24
24
3
Solve the 5 equations above and show that
a0 = 0
2
a1 = a1 =
3
a2 = a2 =
1
12
= 0
= 1
= 0
= 0
= 0
51
Show that in order to obtain a second order accurate finite difference formula, one
must satisfy the following 3 equations
0
fi =
a0 + a1 + a2 = 0
a1 + 2a2 = 1
1
a1 + 2a2 = 0
2
Solve the 3 equations above and show that
a0 =
3
2
a1 = 2
a2 =
3.1.3
1
2
Generally, the finite difference formula for the second derivative, fi d2 f /dx2 (xi )
is expressed as
!
M
M
X
1 X
00
fi = 2
aj fij + a0 fi +
aj fi+j .
(3.9)
j=1
j=1
Expanding the right hand side in terms of Taylor series gives
aM + aM +1 + ... + a1 + a0 + a1 + ...aM 1 + aM
aM (M ) + ... + a1 (1) + a1 + ... + aM M
(M )2
(1)2
12
M2
aM
+ ... + a1
+ a1 + ... + aM
2!
2!
2!
2!
(M )3
(1)3
13
M3
aM
+ ... + a1
+ a1 + ... + aM
3!
3!
3!
3!
..
.
= 0
= 0
= 1
= 0
.
= ..
As for the case as the first derivative, the more equations you satisfy, the more
accurate your finite difference approximation will be.
52
a1 + a0 + a1 = 0
a1 + a1 = 0
1
1
a1 + a1
= 1
2
2!
Solving the above equations will give you
a1 = a1 = 1
a0 = 2
fi =
1
fin = n
3.1.4
M
X
j=1
aj fij + a0 fi +
M
X
!
aj fi+j
(3.11)
j=1
In the interior points, it is common to use central difference scheme. This is because
for a given stencil, the central difference scheme always gives a higher order accuracy.
Some commonly used central difference scheme
fi =
1
(fi2 8fi1 + 8fi+1 fi+2 ) O()4
12
At the end points, it is common to use one-sided difference scheme. Forward difference scheme is
53
1
0
fi =
(3fi + 4fi+1 fi+2 ) O(2 )
2
0
fi =
1
(11fi + 18fi+1 9fi+2 + 2fi+4 ) O(3 )
6x
One can also come up with a partially one sided scheme, for example
0
fi =
1
(2fi1 3fi + 6fi+1 fi+2 )
6x
(3.12)
1
0
(3fi 4fi1 + fi+2 ) O(2 )
fi =
2
0
fi =
1
(2fi3 + 9fi2 18fi1 + 11fi ) O(3 )
6x
One can also come up with a partially one sided scheme, for example
0
fi =
3.2
1
(fi2 6fi1 + 3fi + 2fi+1 )
6x
(3.13)
As you have seen in the exercises above, the most accurate scheme arise when the
coefficients of the finite difference scheme (the aj s and aj s) are antisymmetrical
(i.e. aj = aj ). These finite difference schemes are called the centred difference
scheme and one would start with the following formula
M
df
1 X
0
aj (fi+j fij )
(xi ) fi =
dx
j=1
(3.14)
where fi+j is the numerical approximation to f (xi + j). It is known that one can
expand f (x + j) as a Taylor series
54
f (xi + m) = f (xi )
1 0
+
f (xi )(m)
1!
1 00
f (xi )(m)2
+
2!
1 000
+
f (xi )(m)3
3!
1 iv
f (xi )(m)4
+
4!
1 v
+
f (xi )(m)5
5!
1 vi
+
f (xi )(m)6
6!
+ ...
f (xi )(m)
1!
1 00
+
f (xi )(m)2
2!
1 000
f (xi )(m)3
3!
1 iv
+
f (xi )(m)4
4!
1 v
f (xi )(m)5
5!
1 vi
+
f (xi )(m)6
6!
...
55
2 0
df
(xi ) =
f (xi )(1a1 + 2a2 + 3a3 + . . . + M aM )()
dx
1!
2 00
+
f (xi )(0 + 0 + 0 + . . . + 0)()2
2!
2 000
f (xi )(1a1 + 8a3 + 27a3 + . . . + M 3 aM )()3
+
3!
2 iv
+
f (xi )(0 + 0 + 0 + . . . + 0)()4
4!
2 v
+
f (xi )(1a1 + 32a2 + 243a3 + . . . + M 5 aM )()5
5!
2 vi
+
f (xi )(0 + 0 + 0 + . . . + 0)()6
6!
+ ...
Thus, if we want a formula to calculate the derivative at xi , then all we have to do
is to solve the following system of equations
1a1 + 2a2 + 3a3 + . . . + M aM
1a1 + 8a2 + 27a3 + . . . + M 3 aM
1a1 + 32a2 + 243a3 + . . . + M 5 aM
..
.
= 1
= 0
= 0
.
= ..
The number of equations we need to use will depend on the number of coefficients
you are using in your formula. The higher the value of M , the more accurate the
answer. The error will be of the order of ()2M
3.2.1
Example
56
3.3
Lets now consider using the finite difference technique to solve to solve the convection equation
f
f
= U
t
x
(3.15)
2
Lets first discretize the domain into N + 1 data points. The grid is your series
of x values, x0 , x1 ..... xN , and the corresponding value of f is denoted as f0 , f1 , ....
, fN . For this example, we will use 2nd order central spatial discretization at the
interior nodes. So at the interior points,
df1
dt
df2
dt
df3
dt
..
.
f2 f0
2
f3 f1
U
2
f4 f2
U
2
..
.
fN fN 2
U
2
= U
=
=
=
dfN 1
=
dt
Since, the function f is only defined inside the domain, we do not have the value of
the function at x = x1 and x = xN +1 . Thus we cannot use the central difference
formula at the end points. It is common to use a lower order scheme at the end
points. For this example, we will use 1st order one-sided scheme at the end points to
approximate the spatial derivatives. Note that we only need to find the x derivative
at x = xN since we have been given that f0 = 0 as the boundary condition to our
problem.
dfN
fN fN 1
= U
dt
The above set of N ordinary differential equations can be put into the form
d
{x} = [A] {x} + {c} .
dt
For example, if N = 5, the above system will look like
f
0
1/2
0
0
0
f
1/2
0
1/2
0
0
d 2
U
f3
1/2
0
1/2 0
= 0
dt
f4
0
0
1/2 0 1/2
f5
0
0
0
1 1
f1
f2
f3
f4
f5
f0 /2
0
0
+
57
From previous lectures, we have already seen that the eigenvalues of [A] will determine the stability of the time-integration scheme. The eigenvalues for the above
example is
U
(0.0712 0.8331I)
U
= (0.2500 0.4330I)
U
5 = (0.3576)
1,2 =
3,4
where I is the imaginary unit. Remember that from our stability analysis, the
stability of a system is determined by plotting h where h is the temporal step size.
Thus the stability of our system is dependent on the parameter
Uh
(3.16)
This parameter is called the CFL number. For stability, we would like the CFL to
be as small as possible. Note that the error of our spatial discretization is dependent
on . However, if we make small, then CFD becomes big and the system becomes
unstable !
The stability plots for N = 5 is shown in Fig. 3.1. The eigenvalues are plotted
with CFL=1.0 on both the Euler and Runge-Kutta stability diagrams. It is easily
seen that the Euler scheme is unstable and the 4th order Runge-Kutta scheme would
be stable.
CFL =
58
0.8
2
0.6
0.4
1
I h
I h
0.2
0.2
1
0.4
0.6
2
0.8
1
2
1.8
1.6
1.4
1.2
1
R h
0.8
0.6
0.4
0.2
3
3
2.5
1.5
1
R h
0.5
0.5
Figure 3.1: Convection equation stability plots. Euler scheme (left) and RungeKutta scheme (right)
0.5
0.5
(b)
(a)
0.3
0.3
0.2
0.2
f(x)
0.4
f(x)
0.4
0.1
0.1
0.0
0.0
-0.1
-0.1
-0.2
100
200
300
400
500
600
-0.2
100
200
300
400
500
600
Figure 3.2: Solution of the convection equation at time t = 400. (a) computed with
h = 1.0 and (b) computed with h = 0.001. The exact solution is shown as a solid
line and the computed solution is shown as a dashed line.
59
0.5
0.4
0.3
f(x)
0.2
0.1
0.0
-0.1
-0.2
100
200
300
400
500
600
Figure 3.3: Numerical solution at t = 400 computed with = 0.1, h = 0.01. The
numerical solution is plotted with dashed line and the exact solution is shown as a
solid line. They lie on top of each other.
60
0.6
0.5
0.6
(a)
(b)
0.5
0.3
0.3
f(x)
0.4
f(x)
0.4
0.2
0.2
0.1
0.1
0.1
380
385
390
395
400
x
405
410
415
420
0.1
380
385
390
395
400
x
405
410
415
420
Figure 3.4:
Exact solution, numerical solution. (a) is for = 1.0 and (b) is
computed with = 0.5.
h = 0.01 and small value of = 0.1. The result at t = 400 is shown in Fig. 3.3
together with the exact solution. One cannot distinguish between the two plots.
If we use the 4th order central scheme at the interior points and Kennedy and
Carpenters one-sided scheme at the end points with the classical 4th order Runge
Kutta time step scheme to solve the convection equation given above, with h = 0.1,
the results are shown in Fig 3.4. The wiggles present in Fig 3.4 (a) is due to
dispersion error, waves with different wave numbers travel at different speeds. Some
get left behind. If we reduce , then we get the exact solution.
3.4
The above analysis only gives you an idea how accurate the calculation of derivative
is with respect to the grid size. This does not give all information regarding the
characteristics of the numerical scheme. Many physical processes exhibit wave-like
motions. Hence, a Fourier analysis would provide additional information about the
resolution characteristics of a particular numerical scheme.
Consider a domain, length L. Lets assume that the depedent variable, f (x) is
periodic. Thus, f (x) can be represented as a Fourier series
f (x) =
N/2
X
l=N/2+1
where
cl eIkl x
61
kl = k 0 l
2
=
l
L
kl [0, ]. L is the length of the domain of interest. Direct differentiation of f (x)
gives
df
(x) =
dx
N/2
X
cl Ikl eikl x
(3.17)
l=N/2+1
Note that the Fourier coefficients of the derivative of a periodic function is just the
Fourier coefficients of the function multiplied by Ikl .
In order to analyse the error introduced by the discretisation scheme, note that
N
X
f (x + j) =
cl eIkl (x+j)
l=N +1
N
X
f (x j) =
cl eIkl (xj)
l=N +1
N/2
M
X
1 X
0
fn =
aj
cl eIkl (xn j) + a0
j=1
l=N/2+1
N/2
cl
l=N/2+1
M
X
aj eIkl j + a0 +
j=1
N/2
cl eIkl xn +
j=1
l=N/2+1
M
X
M
X
aj
N/2
l=N/2+1
!
aj eIkl j eIkl xn
j=1
fn =
N/2
X
cl Ikl eIkl xn
(3.18)
l=N/2+1
Ikl =
M
X
j=1
M
X
j=1
aj eIkl j + a0 +
M
X
aj eIkl j
j=1
(aj + aj ) cos(kl j) + a0 + I
M
X
j=1
(aj aj ) sin(kl j)
62
M
X
(aj aj ) sin(kl j) I
j=1
M
X
!
(aj + aj ) cos(kl j) + a0
j=1
If we split kl into its real and imaginary parts, l and -l respectively, then
l =
M
X
l =
(aj aj ) sin(kl j)
j=1
M
X
!
(aj + aj ) cos(kl j) + a0
j=1
Example 3.1:
Recall that the modified wavenumber can be written as kl = l Il For a
second order scheme, one will have a1 = 1/2, a0 = 0 and a1 = 1/2.
l = sin(kl )
l = 0
l =
N/2
X
l=N/2+1
Ikl xn
cl e
= U
N/2
X
cl Ikl eIkl xn
l=N/2+1
Comparing the coefficient of eIkl xn on both sides of the equations will give
dcl
= U Ikl cl
dt
If you solve this equation using analytical techniques
63
cl = cl (0)eIU l t eU l t
(3.19)
N/2
X
Ikl xn
N/2
X
= U
cl e
l=N/2+1
Ikl xn
cl Ikl e
l=N/2+1
N/2
N/2
X
cl kl2 eIkl xn
l=N/2+1
cl IU kl kl2 eIkl xn
l=N/2+1
Comparing the coefficient of eIkl xn on both sides of the equations will give
dcl
= U Ikl + kl2 cl
dt
Solving the above equation using analytical techniques will give you
2
cl = cl (0)eIU kl t ekl t
(3.20)
Comparing Eqs. (3.19) with (3.20) will give you the numerical dissipation
= U
l
kl2
(3.21)
l
kl
(3.22)
Example 3.2:
Lets look at the case where M = 1. Thus the finite difference scheme can be
writen as
1
(a1 fi1 + a0 fi + a1 fi+1 )
For the central difference scheme which is 2nd order accurate, then
0
fi =
a1 =
a0 = 0
1
a1 =
2
1
2
64
(3.23)
For the backward difference scheme which is 1st order accurate, then
a1 = 1
a0 = 1
a1 = 0
thus
l = sin(kl )
l = 1 cos(kl )
(3.24)
Plots of l and l for the central, forward and backward differencing schemes
are shown in Fig. 3.5. Note that in Fig. 3.5 (b), it is shown that the l for the
forward differencing scheme is always negative. This indicates that if you solve
f
f
= U ,
t
x
using the forward differencing scheme, then it will blow up. However, if you solve the
above equation using the backward differencing scheme, then l is always positive
and looking at Eqs. (3.21) and (3.19) would suggest that the solution will decay
with time. Figure 3.6 shows the numerical solution together with the exact solution.
It is clear that the magnitude of the numerical solution is less than than the exact
solution which indicates a diffusion in the numerical scheme.
65
N
U /U
11 0.756827
51 0.989506
Since U /U is always less than 1, then the numerical wave speed will always be
slower than the true wave speed. When we use more number of grid points, then
U /U will approach 1 and the dispersion error will decrease to zero.
3.4.1
1 X
df
(xi ) =
aj
dx
j=1
1
=
1
=
M
X
N
X
l=N +1
N
X
aj
j=1
l=N +1
M
X
N
X
j=1
!
Ik (xi +j)
Ik
(x
j)
cl e l
e l i
aj
!
cl eIkl xi eIkl (j) eIkl (j)
!
cl eIkl xi 2I sin(kl )
l=N +1
N
1 X
=
cl I
l=N +1
M
X
!
aj 2 sin(kl ) eIkl xi
j=1
!
M
N
X
df
1 X
(xi ) =
cl i
aj 2 sin(kl ) eikl xi
dx
l=N +1
j=1
(3.25)
By comparing Eq. (3.17) with Eq. (3.25), we can see that the coefficients of the
Fourier expansion of have now been modified. It should be clear that
66
&
-/.
*+,-'()!.
%"#
%
$"#
$
!"#
!
!"#
$"#
'()!
%"#
&
-2.
01*-' )!.!$
!!"#
!$
!$"#
!%
!"#
$"#
' )!
%"#
&
-0.
$!01*-'()!.
$"#
!"#
!"#
$"#
'()!
%"#
&
Figure 3.5: (a) Plot of l for point the central (2nd order), forward and backward
(1st order) differencing schemes. (b) Plot of l for the forward differencing scheme.
(b) Plot of l for the backward differencing scheme.
67
0.8
0.6
0.4
f(x)
0.2
0
0.2
0.4
0.6
0.8
0
Figure 3.6: Simulation carried out using the upwind differencing scheme at time
t = 2 with U = 1, N = 101 and the initial condition f (x, t = 0) = sin(2x).
exact solution
numerical solution.
68
1
kl =
M
X
!
aj 2 sin(kl )
j=1
(3.26)
For the 1st backward difference scheme, the imaginary part of the non-dimensional
effective wavenumber is
l = 1 cos(kl ).
(3.27)
Also show that for the 4th order central difference scheme, the nondimensional
effective wavenumber is
kl =
3.5
1
4
sin(kl ) sin(2kl )
3
6
(3.28)
In general, when the finite difference approximation is used to solve a linear partial
differential equation, the equation of motion can be written as a set of ordinary
differential equation
d
{f } = [A] {f } .
dt
(3.29)
Strictly speaking, the stability of the above system can be determined by just
finding the eigenvalues of [A]. The linear operator [A] usually contains information
regarding the boundary nodes. However, this is not usually convenient to find the
eigenvalues of [A]. From past experience, it is more more likely for the interior
nodes to go unstable thus there is really no need to worry about what is happening
close to the boundary points. Thus we will confine the stability analysis assuming a
periodic solution and thus confining ourselves to the operator at the interior points.
Consider the periodic function
N/2
X
fj =
69
cl eIkl xj
l=N/2+1
fj =
cl Ikl eIkl xj
l=N/2+1
0
fj =
fj
2
l
l
X 1
Ik xj
Ikl
Ikl
e
e l
=
cl
e
2
l
X
1
=
cl I
sin (kl ) eIkl xj
l
X
cl Ikl eIkl xj
=
where
1
sin (kl ) .
(3.30)
If we now use the 2nd order finite difference scheme to solve the convection equation,
kl =
f
f
= U
t
x
Discretising using finite difference scheme gives
X dcl
l
dt
eIkl xj =
U cl Ikl eIkl xj
l
Ikl xj
gives
dcl
= U cl Ikl
dt
dcl
1
= U cl I
sin (kl )
dt
70
Referring back to the model problem considered when analysing the stability of
the ordinary differential equations (dy/dt = y), the stability is determined by the
maximum magnitude of sin(kl ). Since 0 kl the maximum value of the
right hand side of the above equation is max = IU/. The stability is determined
by max h = IU h/. One can see that is purely imaginary so if we use the explicit
Euler time stepping scheme, the system will always be unstable.
Exercise 3.4: Using similar analysis as above, show that
1 4
1
dcl
= U cl I
sin(kl ) sin(2kl )
dt
3
6
(3.31)
if you use a higher order difference formula, say the 4th order central scheme
1
(fj2 8fj1 + 8fj+1 fj+2 )
12
to approximate the derivative (see section 3.1.4).
fj0 =
Equation (3.31) shows that the stability of your time stepping scheme is dependent on
1 4
1
IU
sin(kl ) sin(2kl ) .
3
6
Since the maximum value of
1
4
sin(kl ) sin(2kl )
3
6
occurs when kl = 1.7974 and the value of 34 sin(1.7974) 16 sin(2 1.7974) =
1.3722. Thus the value of max h for the 4th order central scheme is
h
1.3722.
In the discussion above, it was shown that the maximum value of max h for the 2nd
order central scheme is
IU
h
.
Note that for both cases the value of max h is purely imaginary. In Example 1.3 it
was shown that if max h is purely imaginary, then in order to keep the 4th order
Runge-Kutta scheme stable,
max h < 2.8284I
IU
Thus if we use the discretise the spatial domain using the 2nd order central scheme,
then
71
.
U
Whereas if you discretise the spatial domain using the 4th order central scheme,
then
h < 2.8284
2.8284
1.3722 U
< 2.0612 .
U
h <
So we have to take a slightly smaller time step when you use a higher spatial discretisation scheme. In general, this is usually the case.
3.6
Dispersion-Relation-Preserving Scheme
(3.32)
(3.33)
(3.34)
In order to ensure that the derivatives have better dispersion characteristics, a1 will
be chosen so as to minimize the integrated error, E defined as
Z
/2
E =
/2
/2
=
/2
4
9
i i2a1 sin() i a1 +
5
20
sin(2) i2
1
2
a1
5
15
2
sin(3) d
72
(3.35)
3584
992
2
84
a1
+
+
a1 = 0
375
1125 75
25
1 496 + 15
42 45 128
= 0.799266426974
a1 =
a1
Putting this value of a1 into Eqs. (3.32) and (3.33) will give you
a2 = 0.18941314
a3 = 0.02651995
3.7
j=1
j=1
00
2 fi
2!
2!
2!
2!
dx2
(M )3
(1)3
13
M3
d2 f
+
bM
+ ... + b1
+ b1 + ... + bM
3 2 (xi )
3!
3!
3!
3!
dx
4
+ O( )
00
bM + bM +1 + ... + b1 + b0 + b1 + ...bM 1 + bM
bM (M ) + ... + b1 (1) + b1 + ... + bM M
2
(1)2
12
M2
(M )
+ ... + b1
+ b1 + ... + aM
1
bM
2!
2!
2!
2!
(M )3
(1)3
13
M3
+ ... + b1
+ b1 + ... + aM
bM
3!
3!
3!
3!
..
.
= 0
= 0
= 0
= 0
.
= ..
The more equations we choose to satisfy, the higher the accuracy of the approximation. For example, if we only satisfy the first two equations, then the accuracy
of the approximation is O(). If we satisfy the first three equations, then the
approximation is O(2 ).
3.7.1
Lets now see what happens when we choose M = 2. The finite differencing schem
for the second derivative can be written as
00
fi =
1
(b2 f2 + b1 f1 + b0 f0 + b1 f1 + b2 f2 )
2
00
The equations that we need to satisfy in order for fi to approximate d2 f /dx2 (xi )
can be written as
b2 + b1 + b0 + b1 + b2
2b2 b1 + b1 + 2b2
1
1
2b2 + b1 + b1 + 2b2 1
2
2
4
1
1
4
b2 b1 + b1 + b2
3
6
6
3
2
1
1
2
b2 + b1 + b1 + b2
3
24
24
3
= 0
= 0
(3.37)
(3.38)
= 0
(3.39)
= 0
(3.40)
= 0
(3.41)
74
b2
b1
b0
b1
b2
=
=
=
=
=
0
0
1
2
1
fi =
b2
b1
b0
b1
b2
=
=
=
=
=
1
2
1
0
0
This is called the 1st order backward difference scheme for the second derivative.
If we now choose to satisfy Eqs. (3.37) to (3.41), then we will not have a free
parameter. The formula we get will be second order accurate and
00
fi =
b2 =
b1 =
4
3
b0 =
b1 =
1
12
5
2
4
3
b2 =
1
12
75
fi =
1
(fi1 + 16fi1 30fi + 16fi+1 fi1 ) .
122
(3.44)
This scheme is called the 4t h order centered difference scheme for the second derivative.
Exercise 3.5:
1. Show that the second order centered difference scheme for the second derivative is given by
00
fj =
1
(fj1 2fj + fj+1 )
2
(3.45)
(3.46)
3.8
3.8.1
(3.47)
Multidimensional problems
Steady problem
Lets now see how we can use differentiation schemes to solve the Poissons equation
k
2 2
+
= Q(x, y)
x2 y 2
(3.48)
2 Q(xi , yj )
k
76
2 Q(x1 , y1 )
k
2 Q(x2 , y1 )
k
2 Q(x3 , y1 )
k
2 Q(x1 , y2 )
k
2 Q(x2 , y2 )
k
2 Q(x3 , y2 )
k
2 Q(x1 , y3 )
k
2 Q(x2 , y3 )
k
2 Q(x3 , y3 )
k
4 1
0
1
0
0
0
0
0
1 4 1
0
1
0
0
0
0
0
1 4 1
0
1
0
0
0
1
0
0
4
1
0
1
0
0
0
1
0
1 4 1
0
1
0
0
0
1
0
1
4
1
0
1
0
0
0
1
0
0 4 1
0
0
0
0
0
1
0
1 4 1
0
0
0
0
0
1
0
1 4
1,1
2,1
3,1
1,2
2,2
3,2
1,3
2,3
3,3
2
Q(x1 , y1 )/k 1,0 0,1
2 Q(x2 , y2 )/k
2 Q(x3 , y2 )/k
77
(3.49)
One should then be able to solve Eq. (3.49) using solvers for linear systems. Because
there are so many zeros on the left hand side of Eq. (3.49), it is more common to
use iterative solvers to solve Eq. (3.49).
3.8.2
Unsteady problem
=
+
t
x2 y 2
(3.50)
(x, y = 1, t)
(x = 0, y, t)
(x, y = 0, t)
y
(x = 1, y, t)
x
= 300
= 20(1 y) + 300
= 0
= 0
If we use 2nd order finite difference scheme to represent the spatial derivatives, then
discretise form of Eq. (3.50) is
d
i+1,j 2i,j + i1,j i,j+1 2i,j + i,j1
i,j =
+
(3.51)
dt
2x
2y
where i,j is the approximated value of at point (xi , yj ).
At the points close to the boundaries, we have to implement the boundary conditions. Close to the east boundary,
N x1,j = N x,j
(3.52)
78
d
N 1,j =
dt x
(3.53)
The above is just a system of ODEs that one will have to solve in order to obtain
a time-dependent solution. One can use Euler, Runge-Kutta etc. to get a solution
to the set of equations. In order to carry out stability analysis, one has to put it
into matrix-vector form. Taking into account of the boundary conditions, the above
system of equations could be written into matrix-vector form (assuming N x = 5
and N y = 4 number of intervals in the x and y directions respectively) as follows.
d
dt
1,1
2,1
3,1
4,1
1,2
2,2
3,2
4,2
1,3
2,3
3,3
4,3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1,1
2,1
3,1
4,1
1,2
2,2
3,2
4,2
1,3
2,3
3,3
3,4
0,1 /2x
0,2 /2x
0
+
0
/
+
1,4 /2y
0,3
x
2,4 /2y
3,4 /2y
4,4 /2y
where
=
=
=
a =
b =
1
2x
2
2
2 2
x y
1
2y
1
2
2 2
x y
2
1
2 2
x y
One can then use can be the above system to perform stability analysis.
If you want to use implicit time-stepping scheme to compute the solution to Eq.
(3.51), things get a bit complicated. Lets use the Crank-Nicholson time-stepping
scheme
(n+1)
(n)
i,j
h
i,j
(n+1)
(n+1)
79
(n+1)
(n+1)
(n+1)
=
2
i+1,j 2i,j
2x
+
2
(n)
(n)
+ i1,j
i,j+1 2i,j
+
2y
(n)
(n)
(n)
(n+1)
+ i,j1
(n)
(n+1)
Rearranging and putting on the unknowns (i,j ) on the left hand side and all the
(n)
terms that are known (i,j ) on the right hand side gives
h
h
h h
h
h
(n+1)
(n+1)
(n+1)
(n+1)
(n+1)
i,j1
i1,j + 1 2 2 i,j
i+1,j
i,j+1
2
2
2
2
2y
2x
x y
2x
2y
h
h
h h
h
h
(n)
(n)
(n)
(n)
(n)
+
=
+
1
+
i,j+1
i,j1
i1,j
i,j
i+1,j
22y
22x
2x 2y
22x
22y
One can write the above system in matrix-vector form as
[A] (n+1) = {C}
where the matrix [A] is a sparse matrix consisting of lots of zeros. The vector {C}
(n+1)
(n)
which we know. For
is a linear function of i,j and the boundary values of i,j
Nx = 5 and Ny = 4, the matrix [A] would look similar to
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
a
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
80
where
h
22x
h h
1 2 2
x y
h
2
2y
3h
h
1
2
2x 22y
h
3h
1
22x 22y
=
=
=
a =
b =
In order to solve the system of equations, one can use point Jacobi iteration. We
can re-write
h h
1 2 2
x y
(n+1)
i,j
=
+
+
+
+
+
+
+
+
h
(n+1)
i,j1
22y
h
(n+1)
i1,j
2
2
x
h
(n+1)
i+1,j
2
2
x
h
(n+1)
i,j+1
22y
h
(n)
i,j1
2
2y
h
(n)
i1,j
22x
h h
(n)
1 + 2 + 2 i,j
x y
h
(n)
i+1,j
2
2
x
h
(n)
i,j+1
2
2y
(n+1)
(n+1)
The first four terms on the right hand side (i.e. terms containing i1,j , i+1,j ,
(n+1)
(n+1)
i,j1 , i,j+1 ) are unknowns, so we just have to guess their values in order to solve
(n+1)
(
for i,j+1 . The iterative formula one would use to iteratively solve for i,j n + 1) is
h h
1 2 2
x y
(n+1),r+1
i,j
81
=
+
+
+
+
+
+
+
+
h
(n+1),r
i,j1
22y
h
(n+1),r
i1,j
2
2
x
h
(n+1),r
i+1,j
2
2
x
h
(n+1),r
i,j+1
22y
h
(n)
i,j1
2
2y
h
(n)
i1,j
22x
h h
(n)
1 + 2 + 2 i,j
x y
h
(n)
i+1,j
2
2
x
h
(n)
i,j+1
2
2y
This scheme is called the point Jacobi scheme. One would use the formula above to
(n+1)
iteratively obtain the value for i,j
for the (r + 1)th iteration using the values of
(n+1)
(n+1)
(n+1)
(n+1)
i1,j , i+1,j , i,j1 , i,j+1 ) at the rth iteration.
If you want to use the Gauss-Seidel iteration scheme, one would use the latest
(n),r1
(n),r
value of i,j instead of the previous value i,j
.
3.8.3
XX
l
82
1
301
301
302
0.8
305
0.6
30
0.4
31
0
31
0.2
4
31
6
0
0.2
0.4
0.6
0.8
83
i1,j =
XX
i+1,j =
XX
i,j1 =
XX
i,j+1 =
XX
Substituting the above into Eq. (3.51) and compare coefficients of eIkl xi eIkm yj gives
1 cos (kl x ) 1 cos (kl y )
d
clm (t) = 2
+
clm (t)
dt
2x
2y
= clm (t)
In order for stability, we would like to fall within the stability diagram of the
particular numerical time-stepping scheme that you are using. If you are using
Euler time-stepping, then
2 > h > 0
The worst case scenario is when cos(kl x ) = cos(km y ) = 1. So if you time-step
using the explicit Euler scheme, then
2
2
1
+
.
h
2x 2y
Exercise 3.6: If one were to use 2nd order central difference (for spatial discretisation) and the Euler time stepping schemes to solve the diffusion equation, what
is the maximum value of h that one can use. Use, x = y = = 0.01 and
= 0.1. Compare your 1-d and 2-d results.
3.9
3.10
Test Problems
Euler equations
Using second order central difference with 4th order Runge-Kutta time stepping
with the following parameters dt = 0.1, x = y = 1, Mach number, M = 0.8 and
pressure forcing function
84
80
60
40
20
20
40
60
80
100
100
80
60
40
20
20
40
60
80
100
x
Figure 3.8: Contour plot of the solution to the Euler Equations at time t = 400.
The solution is computed using 2nd order central difference scheme and RungeKutta time stepping. Contour levels are at p = 0.1, 0.05, 0.01, 0.005, 0.001.
Note the dispersion error and also the reflected wave. The numerical solution was
obtained with dt = 0.1, x = y = 1.
sin(t)e((x+20)
2 +y 2 )/9
where = 0.03. The solution at time t = 400 is shown in Fig. 3.8. Note the high
wavenumber dispersion waves and the reflected waves.
The dispersion error could be eliminated by conducting the simulation with more
grid points. Figure 3.9 shows the results from the same simulation but carried out
with a better spatial resolution of x = y = 0.5 with dt = 0.02. The contour plot
shown is at t = 600. Note that no visible high wavenumber dispersive waves. Only
the reflected waves are the source of error. The effects the reflected waves from the
external boundaries could be lessen by computing the solution using a bigger domain.
Solution computed with x = y = 0.5, dt = 0.01 but with 300 x 300 and
300 y 300 is shown in Fig. 3.10.
85
100
80
60
40
20
20
40
60
80
100
100
80
60
40
20
20
40
60
80
100
x
Figure 3.9: Contour plot of the solution to the Euler Equations at time t = 600.
The solution is computed using 2nd order central difference scheme and RungeKutta time stepping. Contour levels equally spaced with 50 contour levels between
the maximum and minimum p values. The numerical solution was obtained with
dt = 0.02, x = y = 0.5.
86
300
200
100
!100
!200
!300
!300
!200
!100
100
200
300
Figure 3.10: Contour plot of p for the solution to the Euler Equations at time
t = 400. The solution is computed using 2nd order central difference scheme and
Runge-Kutta time stepping. The numerical solution was obtained with dt = 0.01,
x = y = 0.5. The solution was computed on a bigger domain, 300 x 300
and 300 y 300
Chapter 4
Differentiation: Unequally spaced
data
In order to find a formula to differentiate a function, consider the Taylor series
0
1 00
1 000
f (xi ) (xi+1 xi )2 + f () (xi+1 xi )3
2!
3!
(4.1)
where
xi+1 = xi + i
df
dx
and is somewhere in between xi and xi+1 . Equation (4.1) can be re-arranged to
obtain
0
f (x) =
f (xi ) =
f (xi+1 ) f (xi )
1 00
1 000
f (xi ) (xi+1 xi ) + f () (xi+1 xi )2
(xi+1 xi )
2!
3!
(4.2)
f (xi ) =
f (xi+1 ) f (xi )
(xi+1 xi )
(4.3)
Equation (4.3) is called the Forward Difference Scheme (FDS). The other terms
that have been neglected in Eq. (4.2) gives an approximation of the error in ap0
proximating f (xi ) with Eq. (4.3). The leading term in the truncation error for the
FDS approximation is
EF DS =
1 00
f (xi ) (xi+1 xi )
2!
87
88
Thus, in order to reduce the error, one would like to make (xi+1 xi ) as small as
possible.
The Taylor series expansion Eq. (4.1) could also be used to obtain an expression
for f (xi1 ).
0
(4.4)
where xi1 < 2 < xi . Equation (4.4) can be re-arranged to obtain an another
0
approximation for f (xi )
0
f (xi ) =
f (xi ) f (xi1 )
(xi xi1 )
(4.5)
1 00
f (xi ) (xi xi1 )
2!
Exercise 4.1: Subtract Eq. (4.4) from Eq. (4.1) and show that, after some
algebra, that the first derivative can be written as
0
00
f (xi+1 )f (xi1 )
(xi+1 xi )2 (xi xi1 )2
f
(x
)
i
xi+1 xi1
2!(xi+1 xi1 )
000
000
f () (xi+1 xi )3
f (2 ) (xi xi1 )3
3! (xi+1 xi1 )
3! (xi+1 xi1 )
f (xi ) =
f (xi ) =
f (xi+1 ) f (xi1 )
xi+1 xi1
(4.6)
This formula is called the Central Difference Scheme (CDS) and its leading order
error is given by
00
ECDS = f (xi )
89
Exercise 4.2:
For the case where all the xi s are equally spaced, i.e. xi+1 xi = = const, show
that Eqs. (4.3), (4.5), and (4.6) simplifies to the following expressions for the FDS,
BDS and CDS
f (xi+1 ) f (xi )
f (xi ) f (xi1 )
0
f (xi ) =
f (xi+1 ) f (xi1 )
0
f (xi ) =
2
Also show that the leading error term in the FDS and BDS is O(h) and the leading
error term in the CDS is O(h2 ).
0
f (xi ) =
f (x) = ex sin(10x)
This function is shown in Fig. 4.1. 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
(CDS, FDS or BDS) to calculate the derivatives, we would like to to confine most
of the grid points in the region 0 x 3 and not so many grid points in the region
x > 3. One method would be to describe place the grid points where
xi+1 = xi + i
i = ri1 and r > 1. Figures 4.2 and 4.3 shows the error in calculating the
derivatives using the CDS in the domain 0 x 5 using 30 grid points. The data
in fig. 4.2 was computed using r = 1.0 (i.e. equally spaced grid points) and the
data in fig. 4.3 was computed using r = 1.1. As you can clearly see, you can get
much better results by simply just making a small change in the location of the grid
points. Numerically, the mean error for r = 1.0 using the CDS scheme is 0.5411 and
for r = 1.1 is 0.2729.
4.1
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 xi1 (see Fig. 4.4).
df
df
(xi+1/2 ) dx
(xi1/2 )
d2 f
dx
(x
)
=
i
dx2
xi + 12 (xi+1 xi ) xi1 + 12 (xi xi1 )
(4.7)
Using Eq. (4.6) to approximate the first derivative, Eq. (4.7) can be written as
90
0.8
0.6
0.4
f(x)
0.2
0.2
0.4
0.6
0.8
0.5
1.5
2.5
x
3.5
4.5
f(x)
0.5
0
0.5
1
0.5
1.5
2.5
x
3.5
4.5
10
(b)
df(x)/dx
5
0
5
10
0.5
1.5
2.5
x
3.5
Error
10
4.5
(c)
10
0.5
1.5
2.5
x
3.5
4.5
Figure 4.2: Analysis of the error in the calculation for the derivatives for r = 1.0.
(a) shows the function (
) and the data points, xi ( ). The analytical derivative
of the function (
) and the calculated derivative using CDS is shown using the
symbols in (b). (c) shows the error, indicated by the symbols, calculated at every
grid point.
d2 f
(xi ) =
dx2
f (xi+1 )f (xi )
(xi1 )
f (xxi )f
xi+1 xi
i xi1
1
(xi+1 xi1 )
2
(4.8)
91
1
(a)
f(x)
0.5
0
0.5
1
0.5
1.5
2.5
x
3.5
4.5
10
(b)
df(x)/dx
5
0
5
10
0.5
1.5
2.5
x
3.5
Error
10
4.5
(c)
10
0.5
1.5
2.5
x
3.5
4.5
Figure 4.3: Analysis of the error in the calculation for the derivatives for r = 1.1.
(a) shows the function (
) and the data points, xi ( ). The analytical derivative
of the function (
) and the calculated derivative using CDS is shown using the
symbols in (b). (c) shows the error, indicated by the symbols, calculated at every
grid point.
xi-1
xi
xi-1/2
xi+1/2
xi+1
Figure 4.4: Figure showing the data points used to calculate the second derivative
92
ECDS
1
EBDS = f 00 (xi )i1
2
1 00
i1
i
= f (xi )
1
2
2
i1
where i = xi+1 xi
4.2
In many engineering problems, you will have to solve ordinary differential equations
that look something like
d2 y
dy
+ p(x) + q(x)y = r(x)
2
dx
dx
(4.9)
where
axb
We are also usually given the conditions at the boundaries
y(a) =
y(b) =
Exercise 4.4: Using CDS approximation for all the derivatives, show that Eq.
(4.9) can be approximated as
2
i
xi+1px
+
yi1
(xi+1 xi1 )(xi xi1 )
i1
2
+ qi y i
(4.10)
+ (xi+1 xi )(x
i xi1 )
i
+ xi+1px
+ (xi+1 xi12)(xi+1 xi ) yi+1 = ri
i1
where
yi = y(xi )
pi = p(xi )
qi = q(xi )
ri = r(xi )
93
(4.11)
Eq. (4.10) together with Eq. (4.11) represent a linear system that can be solved.
They can be represented by
[A]{X} = {C}
where
1 0
0
0
0
1 1 1
0
...
2
0
0 2 2
[A] =
..
.
.
.
0 0
..
..
..
.
0 0 N 2 N 2 N 1
0 0
0
0
y0
..
{X} =
.
..
..
yN
r
1
r2
{C} =
..
N
1
and
i
i
i
2
pi
+
=
+ qi
(xi+1 xi ) (xi xi1 )
pi
2
=
+
xi+1 xi1 (xi+1 xi1 ) (xi+1 xi )
94
Chapter 5
Galerkin Method
In this chapter, we will use Fourier series to solve a few classical partial differential
equations
5.1
Convection equation
=0
t
x
(5.1)
N/2
X
ck (t) eikx .
(5.2)
k=N/2+1
Z
0
dck i(kl)x
e
dx =
dt
95
N/2
X
k=N/1+1
Z
ck ik
|0
ei(kl)x
{z
}
=2kl
96
(5.3)
The equation above is on ordinary differential equation with t as the independent variable. The initial condition is obtained by Fourier transforming the initial
condition. The solution to Eq. (5.3) can be formally written as
ck (t) = ck (t = 0)eikt .
(5.4)
5.2
Burgers Equation
(5.6)
N/2
X
ak (t) eikx
(5.7)
k=N/2+1
All terms in Eq. (5.6), must be evaluated before one can solve it. Look at the
diffusion term on the right hand side first
u
1 X
=
ak (ik) eikx
x
N k
N
ikx
2u
1 X
2
a
k
e
=
k
x2
N k
(5.8)
Calculation of uu/x is a problem (will elaborate later), but for the moment,
it will be treated as just another periodic function and can be represented as
97
N
1 X
u
=
u
bk (t)eikx
x
N k
(5.9)
Substituting Eqs. (5.9), (5.8) and (5.7) into Eq. (5.6) gives
1 X dak ikx
1 X ikx
1 X
e +
k 2 ak eikx
bk e =
N k dt
N k
N k
Use Galerkin Formulation and multiply the above equation by (1/2)eimx gives
N
1
X
k=0
dak
dt
2
i(mk)x
dx +
N
1
X
2
i(mk)x
bk
dx =
k=0
N
1
X
k ak
ei(mk)x dx
k=0
Using the orthogonality condition will give the Burgers equation in Fourier space
dam
+ bm = k 2 am
(5.10)
dt
Equation (5.10) is just an ordinary differential equation and one can use RunggeKutta and other standard ode solvers to time step it forward with time. There is
another neater way of doing things. One can rearrange Eq. (5.10) as
dam
+ k 2 am = bm
dt
Multiplying both sides by an integrating factor
R
k2 dt
= ek
2t
gives
ek
2t
dam
2
2
+ ek t k 2 am = bm ek t
dt
d
2
2
am ek t = bm ek t
dt
(5.11)
am (t + h)em
2t
am (t)em
h
Re-arranging this equation will give us
2t
= bm (t)em
2 t
(5.12)
So to use Fourier Spectral Galerkin method to solve Burgers equation, do the following:
98
(1) Fourier transform u(x, t = 0) = sin(x) to get the Fourier coefficients, ak (t = 0).
(2) Calculate Fourier coefficients for u/x ck = ikak
(3) Inverse transform ck to get u/x(xi ) in physical space.
(4) Calculate u(xi )u/x(xi ) in physical space.
(5) Obtain bk by taking the Fourier transform of u(xi )u/x(xi ). This method is
commonly known as the pseudo spectral method.
(6) Use Eq. (5.12) to get am (t + h)
(7) Repeat steps (2) - (6) to step forward the solution in time.
5.3
Lets see what happen when we use pseudo spectral method to calculate the Fourier
coefficients of the product of two functions, u(x) v(x). Assume that both u(x)
and v(x) are periodic functions
N
1
u (xj ) =
N
2
X
u (p) eipxj
+1
p= N
2
N
v (xj ) =
1
N
2
X
(q) eiqxj
q= N
+1
2
(5.13)
The product of u and v can then be written as
u (xj ) v (xj ) =
1 XX
u (p) v (q) ei(p+q)xj
N2 p q
1 XXX
u (p) v (q) ei(p+qk)xj
N2 j p q
(5.14)
1 X
1
u (p) (q) +
N k=p+q
N
|
u (p) (q)
k=p+qN
{z
Errorterm
1 X
u (p) (q)
N k=p+q
Example
Let
u(x) = sin(3x)
v(x) = sin(5x)
Thus
u(x)v(x) = sin(3x) sin(5x) =
If N is big, this is what we will get
N
2 i
N2 i
u (k) =
0
N
2 i
N2 i
v (k) =
0
N
4 i
N4 i
u
cv (k) =
1
1
cos(2x) cos(8x)
2
2
if k = 3
if k = 3
otherwise
if k = 5
if k = 5
otherwise
if k = 2
if k = 8
otherwise
if N = 14
7i
7i
u (k) =
7i
7i
v (k) =
if k = 3
if k = 3
otherwise
if k = 5
if k = 5
otherwise
100
3.5 if k = 2
3.5 if k = 6
u
cv (k) =
0
otherwise
The mode at k = 6 is the aliasing error term of
1
N
u (p) v (q)
k=p+qN
The error mainly occurs because the FFT cannot differentiate between
uv(x) =
1
1
cos(2x) cos(8x)
2
2
uv(x) =
1
1
cos(2x) cos(6x)
2
2
and
Chapter 6
Collocation Method
In the collocation method, the approximated solution to the differential equation is
required to satisfy the differential equation at discrete collocation points. Thus it is
required to find the derivative at the collocation points.
6.1
ck =
N
1
X
f (xj )eikxj
(6.1)
j=0
1
f (xj ) =
N
N/2
X
ck eikxj
(6.2)
k=N/2+1
It was earlier shown that the derivative at the grid points, xj can be written as
(Du)l =
1
N
N/2
X
k=N/2+1
ck ikeikxl .
(6.3)
102
(Du)l
1
=
N
N/2
X
N
1
X
k=N/2+1 j=0
N 1
1 X
=
N j=0
N 1
1 X
=
N j=0
N/2
X
k=N/2+1
N/2
X
ikeik(xl xj ) f (xj )
k=N/2+1
(6.4)
remembering that
xj =
2j
N
gives
N 1
1 X
(Du)l =
N j=0
N/2
X
ike
2ik
(lj)
N
f (xj ).
(6.5)
k=N/2+1
If we define
dl,j
1
=
N
N/2
X
ike
2ik
(lj)
N
(6.6)
k=N/2+1
N
1
X
dlj f (xj ).
(6.7)
j=0
if l 6= j
if l = j
du
(x
dx
du
(x
dx
= x0 )
= x0 )
..
.
..
.
du
(x = xN 1 )
dx
d0,0
..
.
..
.
..
.
d0,1
..
.
...
dN 1,0
...
d0,N 1
..
.
..
.
..
.
dN 1,N 1
u(x = x0 )
u(x = x1 )
..
.
..
.
u(x = xN 1 )
Chapter 7
Some numerical examples
In this chapter, some examples of problems and solutions written by past and present
PhD. students are presented. It is intended to give the reader the vast majority of
problems that can be solved using the numerical methods presented in this set of
notes.
7.1
In this section, we will look at the problem where an oil droplet is suspended in
air. This same problem has been solved by Dombrovsky and Sazhin [2] so comparison could be made wit their data. The partial differential equation governing the
behaviour of temperature inside the liquid, l , is
1
l
2 l
= kl 2
r
.
(7.1)
l cp,l
t
r r
r
where kl = 0.14W/(mK) is the thermal conductivity of oil, cp,l = 2.83 kJ/(kg K)
is the specific heat capacity of the liquid and l = 600 kg/m3 is the density of the
oil. We will assume that the oil have constant properties. Lets also assume that
the droplet size is 20m so Eq. (7.1) is only valid for 0 < r < 2 106 . The initial
temperature of the droplet is set to be 300 K. Because of symmetry, the boundary
condition at r = 0 is
l (0)
=0
r
On the gas side, the governing equation for the temperature, g , is
g
1
2 g
g cp,l
= kg 2
r
.
t
r r
r
(7.2)
(7.3)
104
oil. Again, we will assume all properties are constant. Lets set the boundary of
the calculation to be 20 times the radius of the droplet so Eq. (7.1) is only valid for
2 106 < r < 4 105 . The initial temperature of the air is set to be 880K.
At the interface of the droplet, energy conservation requires that
g
l
= kg
r
r
Using one sided scheme to approximate the derivatives for Eq. (7.4) gives
kl
kl
(7.4)
s l,Nl 1
g,1 s
= kg
l
g
s =
kg l
kl g g,1
kg l
kl g
l,Nl 1 +
1+
(7.5)
where it is assumed that the temperature on the liquid side is represented by the
array [l,0 , l,1 , ....., l,Nl 1 , l,Nl ] and the temperature on the gas side is represented
as [g,0 , g,1 , ....., g,Ng 1 , l,Ng ]. Ng and Nl are the number of grid points (including
boundary points) in the gas and liquid side respectively. l and g are the grid size
in the gas and liquid side respectively. s = l,Nl = g,0 is the temperature on the
surface of the droplet.
The solution can be obtained by solving Eqs. (7.1) and (7.3) along with the
boundary conditions given by Eq. (7.5) and (7.2). g,i (i = 1..Ng ) is set at 880 K.
The solution can be computed using Runge-Kutta time stepping scheme using the
following numerical parameters g = 7.7551 106 , l = 2.2222 106 , h = 105 .
The solution at various times are shown in Fig. 7.1. The corresponding matlab code
is shown in Matlab Program 7.1 and 7.2
900
+!!
800
*!!
700
)!!
600
! (K)
,!!
7.1.
-.!/!!!!!!%0
(!!
400
'!!
300
!
"
1
!&
$%"!
900
900
800
800
700
700
600
600
! (K)
! (K)
t=0.001000 s
500
&!!
t=0.002000 s
500
1
r
2
!4
x 10
t=0.003000 s
500
400
400
300
300
0
1
r
!4
x 10
900
900
800
800
700
700
600
600
! (K)
! (K)
105
t=0.004000 s
500
1
r
t=0.005000 s
500
400
2
!4
x 10
400
300
300
0
1
r
2
!4
x 10
1
r
2
!4
x 10
Figure 7.1: Numerical solution for droplet heating problem. , temperature inside
the droplet, , air temperature and , temperature on the interface of the droplet.
106
rliquid=linspace(0,L,Nliquid);
rgas=linspace(L,20*L,Ngas);
rliquid=rliquid;
rgas=rgas;
Deltaliquid=rliquid(2)-rliquid(1);
Deltagas=rgas(2)-rgas(1);
Delta=Deltaliquid;
7.1.
107
ratio=(kg/kl)*(Deltaliquid/Deltagas);
fliquid=300*ones(size(rliquid));
fgas=880*ones(size(rgas));
fgas(1)=(fliquid(Nliquid-1)+ratio*fgas(2))/(1+ratio);
fliquid(Nliquid)=fgas(1);
%
%Plotting initial solution
%
t=0;
gg=plot(rliquid(1:Nliquid-1),fliquid(1:Nliquid-1),ko, rgas(2:Ngas),fgas(2:Ngas),ks,
set(gg(3),markerfacecolor,[0 0 0],markersize,10);
axis([0 20e-5 250 900]);
crap=sprintf(t=%f s,t);
text(1.5e-4,500,crap,fontsize,12);
xlabel(r);
ylabel(\phi)
%
% Starting time iteration
%
for j=1:nstep
foldgas=fgas;
k1=CalcRHS(nugas,Deltagas,rgas,foldgas);
temp=foldgas+k1*(h/2);
k2=CalcRHS(nugas,Deltagas,rgas,temp);
temp=foldgas+k2*(h/2);
k3=CalcRHS(nugas,Deltagas,rgas,temp);
temp=foldgas+k3*h;
k4=CalcRHS(nugas,Deltagas,rgas,temp);
fgas=foldgas+((1./6.)*k1+(1./3.)*(k2+k3)+(1./6.)*k4)*h;
foldliquid=fliquid;
k1=CalcRHS(nuliquid,Deltaliquid,rliquid,foldliquid);
108
%
% Plotting solution
%
figure(1)
gg=plot(rliquid(1:Nliquid-1),fliquid(1:Nliquid-1),ko, rgas(2:Ngas),fgas(2:Ngas),
set(gg(3),markerfacecolor,[0 0 0],markersize,10);
axis([0 20e-5 250 900]);
crap=sprintf(t=%f s,t);
text(1.5e-4,500,crap,fontsize,12);
xlabel(r);
ylabel(\phi (K))
figure(2)
plot(t,fgas(1),ko);
xlabel(t(s));
ylabel(\phi_s);
hold on
pause(0.001);
end
7.2.
dr=Delta;
N=length(f);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initializing matrix and allocating memory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d2f=zeros(N,1);
d2f(1)=0;
for i=2:N-1
d2f(i)=mu*((f(i+1)-2*f(i)+f(i-1))/dr^2+(2/r(i))*(f(i+1)-f(i-1))/(2*dr));
end
d2f(N)=0;
7.2
Blasius solution:
Contributed by Mr. M. Giacobello
To derived the Blasius boundary layer equation, the starting point is the steady
incompressible NavierStokes and continuity equations. For a freestream flow that
is aligned with the x-axis, these are given as
2
u
1 p
u 2u
u
+v
=
+
+
(7.6)
u
,
x
y
x
x2 y 2
2
v
v
1 p
v 2v
+v
=
+
+
,
u
(7.7)
x
y
y
x2 y 2
u v
+
= 0.
(7.8)
x y
For a boundary layer developing on a semi-infinite flat plate, these equations are
subject to the boundary conditions
u=v=0
for
x 0,
y = 0.
(7.9)
The Navier-Stokes equations can be simplified by noting that within the boundary layer u changes from zero at y = 0 to a value close to Uo at y = (x). Here (x)
is a measure of the boundary layer thickness and is defined as the distance from
110
the wall, to the point where u reaches 99% of the freestream velocity, Uo . In the
wall-parallel direction, the rate of change of u is much more gradual. By considering
the order of magnitude of each term in the Navier-Stokes and continuity equations
and retaining only the dominant terms, it can be shown that equation (7.6) and
(7.8) reduce to
u
u
2u
u
+v
= 2
x
y
y
u v
+
= 0.
x y
(7.10)
(7.11)
These equations must be solved within the boundary layer subject to the boundary
conditions
u=v=0
for
x 0,
y = 0.
(7.12)
y .
(7.13)
v0
as
To simplify the solution further, the stream function, , is introduced, such that
,
y
v= .
x
u=
(7.14)
(7.15)
(7.16)
(7.17)
(7.18)
f
(0) = 0,
(7.19)
(7.20)
7.2.
Equation (7.18) to (7.20) comprise a boundary value problem for f (). Due to its
non-linearity, equation (7.18) cannot be solved in closed form and Blasius resorted
to a series solution. Today it is much easier to solve this equation using a computer.
Here a fourth order RungeKutta method is used. RungeKutta integration methods are best suited to initial value problems. To solve this boundary value problem,
it is first recast as the initial value problem
1 2f
f
= 0,
2 2
f
(0) = 0,
f (0) =
3f
3
2f
() = tk .
2
(7.21)
(7.22)
(7.23)
(7.24)
Here the correct value of tk is converged-on iteratively using the secant method.
An initial guess of tk is made and equation (7.21) must be integrated toinfinity.
At this point the computed value of f / is compared to the specified boundary
value. If the agreement is not satisfactory another guess is made. This is known as
a shooting method.
In practice it is not necessary to integrate to infinity. Integrating from = 0
to 20 was found to be sufficient for an accurate solution. The numerical solution for
f , f / and 2 f / 2 versus is presented in figure 7.2. Once f () is obtained,
u() and v() may be calculated directly from their relationship to the dimensionless
stream function. That is
f
u() = Uo ,
1/2
1 Uo
f
v() =
f .
2
x
(7.25)
(7.26)
u()/Uo and v()/Uo are plotted versus y/(x) in figure 7.3 for Rex = Uo x/ = 1.
112
2.5
f
df/d!
d2f /d!2
2
1.5
0.5
3
!
7.2.
1
u/Uo
v/Uo
0.9
0.8
0.7
y/!(x)
0.6
0.5
0.4
0.3
0.2
0.1
0
0.2
0.4
0.6
0.8
U/Uo
Figure 7.3: Numerical solution for u()/Uo and v()/Uo versus y/(x) for Rex = 1.
114
function main
close all;
clear all;
clc;
error_tol = 1e-5;
nu = 1.0;
x = 1.0;
Uo = 1.0;
output_tag = 1;
% 1 == output, 0 == no output.
N = 800;
kk = 1;
M = 100;
a
b
= 0.0;
= 20.0;
% domain boundries.
ya
dya
dyb
=
=
=
% boundary conditions.
f(0) = 0.
%
f(0) = 0.
%
f(infty) = 1.0
= (b-a)/N;
0.0;
0.0;
1.0;
tkm2 = (dyb-dya)/(b-a);
% grid spacing.
% first guess of d2f/dx2.
7.2.
for ii = 1:N+1
eta(ii) = a + (ii-1)*h;
end
w1(1) = ya;
w2(1) = dya;
w3(1) = tkm2;
[w1,w2,w3] = RKutta(w1(1),w2(1),w3(1),a,h,N);
ykm2 = w2(N+1);
tkm1 = tkm2 +(dyb-ykm2)/(b-a); % second guess of d2f/dx2.
while (kk <= M)
w1(1) = ya;
w2(1) = dya;
w3(1) = tkm1;
[w1,w2,w3] = RKutta(w1(1),w2(1),w3(1),a,h,N);
ykm1 = w2(N+1);
tk = tkm1-(ykm1-dyb)*(tkm1-tkm2)/(ykm1-ykm2);
tkm2 = tkm1;
ykm2 = ykm1;
tkm1 = tk;
kk = kk + 1;
if(abs(w2(N+1)-dyb)<error_tol)
% if error is small plot to screen.
figure;
plot(eta,w1,-x,eta,w2,-o,eta,w3,-+)
axis([0.0 6.0 0.0 2.5])
xlabel(eta);
%title(Blasius ZPG Boundary Layer Profile);
grid on;
legend(f,df/deta,d^2f /deta^2);
string = sprintf(Convergence after %d iterations using secant method,kk);
disp(string);
break;
elseif(kk==M)
string = sprintf(Solution did not converge....);
116
end
%-------------------------------------------------------------%
%
return to physical space.
%
%
Note:
u = dPsi/dx is calculated using 2nd order FD.
%
%
The finite difference has been checked for
%
%
a known Psi.
%
%
v = 0.5*sqrt(nuUo/x)*(eta*f-f).
%
%-------------------------------------------------------------%
for ii = 1: N+1
y(ii)
= sqrt(nu*x/Uo)*eta(ii);
Psi(ii) = sqrt(nu*Uo*x)*w1(ii);
%Psi(ii) = y(ii)^3;
end
% The streamfunction.
%ii = 1;
% forward difference.
%u(ii) = (-3*Psi(ii) + 4.0*Psi(ii+1) - Psi(ii+2))/(2*(y(ii+1)-y(ii)));
%for ii = 2: N
% central difference.
%u(ii) = (Psi(ii+1)-Psi(ii-1))/(2*(y(ii+1)-y(ii)));
%end
%ii = N+1;
% backward difference.
%u(ii) = (Psi(ii-2)-4.0*Psi(ii-1)+3*Psi(ii))/(2*(y(ii)-y(ii-1)));
end
for ii = 1: N+1
u(ii) = Uo*w2(ii);
v(ii) = 0.5*sqrt(nu*Uo/x)*(eta(ii)*w2(ii)-w1(ii));
end
figure;
plot(u,y,-x, v,y,-o)
legend(u,v);
xlabel(U);
ylabel(y);
axis([0.0 1.1 0.0 10.0])
%title(Blasius ZPG Boundary Layer Profile in Physical Coordinates);
grid on;
%-------------------------------------------------------------%
7.2.
%
Finally nornalise by U_99.
%
%-------------------------------------------------------------%
for ii = 1: N+1
if(u(ii)>=0.99*Uo)
y99 = y(ii)
u99 = u(ii);
break;
end
end
figure;
plot(u/Uo,y/y99,-x, v/Uo, y/y99,-o)
legend(u/Uo,v/Uo);
%title(Normalised Blasius ZPG Boundary Layer Profile);
xlabel(U/Uo);
ylabel(y/delta(x));
axis([0.0 1.1 0.0 1.0])
grid on;
if (output_tag == 1)
out1 = sprintf(f_vs_eta.plt);
out2 = sprintf(u_v_vs_y.plt);
fid1 = fopen(out1,a);
fid2 = fopen(out2,a);
fprintf(fid1,variables="eta","f", "df", "d2f");
fprintf(fid2,variables="y","u", "v");
for ii = 1: N+1
fprintf(fid1,%10.8g %10.8g %10.8g %10.8gn, eta(ii), w1(ii), w2(ii), w3(ii));
fprintf(fid2,%10.8g %10.8g %10.8g n, y(ii)/y99, u(ii), v(ii));
end
fclose(fid1);
fclose(fid2);
end
118
for ii = 1:1:N
eta(ii) = a + (ii-1)*h;
K11 = h*rhs1(eta(ii), w1(ii), w2(ii), w3(ii));
K12 = h*rhs2(eta(ii), w1(ii), w2(ii), w3(ii));
K13 = h*rhs3(eta(ii), w1(ii), w2(ii), w3(ii));
K21 = h*rhs1(eta(ii)+h/2,w1(ii)+K11/2,w2(ii)+K12/2,w3(ii)+K13/2);
K22 = h*rhs2(eta(ii)+h/2,w1(ii)+K11/2,w2(ii)+K12/2,w3(ii)+K13/2);
K23 = h*rhs3(eta(ii)+h/2,w1(ii)+K11/2,w2(ii)+K12/2,w3(ii)+K13/2);
K31 = h*rhs1(eta(ii)+h/2,w1(ii)+K21/2,w2(ii)+K22/2,w3(ii)+K23/2);
K32 = h*rhs2(eta(ii)+h/2,w1(ii)+K21/2,w2(ii)+K22/2,w3(ii)+K23/2);
K33 = h*rhs3(eta(ii)+h/2,w1(ii)+K21/2,w2(ii)+K22/2,w3(ii)+K23/2);
K41 = h*rhs1(eta(ii)+h,w1(ii)+K31,w2(ii)+K32,w3(ii)+K33);
K42 = h*rhs2(eta(ii)+h,w1(ii)+K31,w2(ii)+K32,w3(ii)+K33);
K43 = h*rhs3(eta(ii)+h,w1(ii)+K31,w2(ii)+K32,w3(ii)+K33);
w1(ii+1)
w2(ii+1)
w3(ii+1)
end
= w1(ii) +(K11+2*(K21+K31)+K41)/6;
= w2(ii) +(K12+2*(K22+K32)+K42)/6;
= w3(ii) +(K13+2*(K23+K33)+K43)/6;
Bibliography
[1] K. Atkinson. Elementary Numerical Analysis. John Wiley and Sons, 1985.
[2] L. Dombrovsky and S. Sazhin. A parabolic temperature profile model for heating
of droplets. J. Heat Transfer, 125:535537, 2003.
[3] P. Moin. Fundamentals of Engineering Numerical Analysis. Cambridge University Press, 2001.
[4] C. Tam, J. C. Webb, and Z. Dong. A study of the short wave components in
computational acoustics. J. Comp. Acoustics., 1:130, 1993.
119