1. No category

Lab # 3
Computational Physics III (Phys581)
Rachid Ouyed
Due April 14, 2015 (at the end of class)
[Total: 100 marks including 20 marks for the report]
• The latexed report is worth 20% of the total mark. Only well documented and neatly
presented reports will be worth that much! It is not sufficient to give only numerical
results and show plots, you should also discuss your results. Include complete figure
captions, introduction and conclusion sections.
• Your report must be in a two-column format.
• Your report should include your Fortran code and Gnuplot scripts in an Appendix
using the verbatim command.
• If applicable, animations should be shown to the teacher and/or TA before handing
in the lab report. The animation should be well documented and contains necessary
information (student name, assignment number, run time etc ...). Basically, the
information should be included in each frame before they are put together.
• You must name your report using the names (last names only): student1-student2phys581-lab#.pdf.
• Procedure for Handing in your lab report (see instructions on phys581
website):
• You must check with your TA that your report was received and is readable BEFORE you leave the lab.
1
1
Warm Up [16 Marks]
First, go over
→ all of the material covered in Chapter 12 in Ouyed&Dobler on pjl.
→ all of the material covered in Chapter 19 in Numerical Recipes (f77) or equivalently in Chapter 20 in Numerical Recipes (f95).
• [1.0 Marks] Describe in a few paragraphs the direct method for solving time-independent
PDEs.
• [1.0 Marks] Describe in a few paragraphs the Iterative Solution Methods for solving
time-independent PDEs.
• [1.0 Marks] Describe in a few paragraphs the Jacobi iteration method and the GaussSeidel iteration method for solving time-independent PDEs. Comment briefly on key
differences.
• [3.0 Marks] The voltage u(x, y) inside a condensator with parallel plates separated
by a distance L obeys the Laplace’s equation
uxx + uyy = 0,
(1)
0 ≤ x ≤ L 0 ≤ y ≤ 1.
(2)
in the region
The boundary conditions are given by:
2πx
)
L
2πx
u(x, 1) = −100 sin(
)
L
u(0, y) = u(L, y) = 0
u(x, 0) = 100 sin(
(3)
(4)
(5)
(i) Use the Jacobi and Gauss-Seidel (GS) methods to solve for the problem on a
64 × 64 grid with L = 1 and convergence criteria = 1x10−8 ; is defined by the
maximum difference between any two grid points in consecutive iterations.
(ii) Use gnuplot to plot the resulting solution side-by-side on a 2-panel figure.
(iii) Which one converged faster? It means how many steps were necessary for the
Jacobi and Gauss-Seidel methods to converge?
Hint: Use a convergence criteria defined by the maximum difference between any
two grid points in consecutive iterations.
2
• [1.0 Marks] Give the mathematical expressions and the order of the following PDEs:
a) Poisson’s equation
b) Laplace’s equation
c) 3-Dimensional Wave equation
d) Heat equation in one spatial dimension
e) Wave equation in one spatial dimension
f) Ginzburg-Landau equation
g) Korteweg-de Vries equation
• [1.0 Marks] Classify (elliptic, parabolic or hyperbolic) the following PDEs:
a) Utt = 2Uxx
b) Uxx + Uyy = 0
c) Ut − U xx = 0.
• [1.0 Marks] Complete table 1.
PDE
Ux
Ux
Ux
Uxx
Ut
Ut
Ut
Utt
FD approximation
Type
Spatial discretization
—–
forward
—–
backward
—–
central
—–
symmetric
Time discretization
—–
forward
—–
backward
—–
central
—–
symmetric
Order
—–
—–
—–
—–
—–
—–
—–
—–
Table 1: PDEs: Type and Order
• [2.0 Marks] Take the function U (x) = x2 and find the first order forward FD
approximation to Ux (3) using step size h = 0.1 and for h = 0.05. Compare to the
exact formula and conclude.
• [1.0 Marks] Repeat the question above (i.e. U (x) = x2 ) using the central FD approximation and show that it is second order accurate.
• [1.0 Marks] There are several finite difference schemes which include (see Table 2)
the: (i) explicit schemes: Forward Euler, Upwind, Leap-Frog, and Lax-Wendroff; (ii)
the implicit schemes: Backward Euler, and Crank-Nicolson.
3
Complete Table 2 for the 1D advection equation; ut = cux with c > 0. In the
last column describe the stability (in terms of the COURANT NUMBER r =
c∆t/∆x).
Name
FD Scheme
Forward Euler
Upwind
Leap-Frog
Lax-Wendroff
—
—
—
—
Backward Euler
Cranck-Nicolson
—
—
Order of accuracy
Explicit
—
—
—
—
Implicit
—
—
CFL stability restriction
—
—
—
—
—
—
Table 2: PDEs: Explicit and Implicit methods
• [1.0 Marks] What are the key differences between implicit FD methods and explicit
FD methods for PDEs?
• [1.0 Marks] Explain the purpose of the Von Neumann stability analysis for FD
schemes of PDEs (see §12.5.1 in Ouyed&Dobler).
• [1.0 Marks] Use von Neumann stability analysis to show that the Crank-Nicolson
scheme is unconditionally stable.
2
FDs in the Numerical Recipes [2 Marks]
You will find the material covered in Chapter 19 in Numerical Recipes (f77) or
equivalently in Chapter 20 in Numerical Recipes (f95) helpful.
Apply the FD methods the for solution of flux conservative PDEs of the type
∂u(x, t) ∂f (u)
+
=0,
∂t
∂x
(6)
∂f
using (consider the two cases: ∂f
∂x > 0 and ∂x < 0):
(i) [0.5 Marks] The FTCS (Forward Time Centered Space) method (Numerical Recipes
Eq. 20.1.11). Why is this method rarely (in fact, never) used?
(ii) [0.5 Marks] The Lax Method (Numerical Recipes; Chapter 20.1.2). Note, and
describe, the small difference compared to FTCS.
(iii) [0.5 Marks] The Upwind Method ((Numerical Recipes Eq. 20.1.27).
4
(iv) [0.5 Marks] The Lax-Wendroff scheme (Numerical Recipes chapter 20.1.4). Note:
there is a misprint in 20.1.37; there is a factor 1/2 missing in the second term on the
right-hand side.
3
Popular cases [28 Marks]
3.1
The advection equation [10 Marks]
One can think of the terms in the advection equation as representing a concentration u(x, t)
(e.g. a river pollutant) and c the speed of the flow (e.g. a river) so that in 1-D: ut = cux
with c > 0 or c < 0.
• [1.0 Marks] Give the exact analytical solution of the 1D advection equation (ut = cux
with c > 0) using the method of characteristics.
• [1.0 Marks] If you set a pulse as the initial function u0 what happens to it when c > 0
and when c < 0?
• [8.0 Marks] Apply the FD schemes listed in Table 2 to the 1D advection equation for
the following 4 cases (each representing different initial conditions):
u(x, 0) = sin (2x) .
(7)
(
0 if x < 0
u(x, 0) =
1 if x ≥ 0
(8)


0 if x < 0
1
u(x, 0) =
∆x if x = 0


0 if x > 0
(9)
u(x, 0) = exp (−x2 ) .
(10)
You are asked to find the solution for the 7 cases (each with different set of parameters) given in Table 3.
(i) Using gnuplot, and for each case, show the solution for 0 ≤ t ≤ 2.0 and −2.0 ≤
x ≤ 2.0. Use multiplot to better display your result.
(ii) Comment on your results. Is the solution behaving as expected for different values
of r in each case? Explain.
5
Trial
1
2
3
4
5
6
7
c
0.5
0.5
0.5
0.5
0.5
0.5
0.5
∆x
0.04
0.02
0.0137
0.0101
0.0099
0.02
0.02
∆t
0.02
0.02
0.02
0.02
0.02
0.01
0.04
r
0.25
0.5
0.728
0.99
1.11
0.25
1
Table 3: FD parameters for the advection equation
3.2
The diffusion equation [6 Marks]
The pure diffusion equation, Ut −βUxx = 0, has an analytical solution for special conditions
that is useful for checking numerical schemes; β is called the diffusion coefficient1 . In a
pure diffusion case (i.e. no advection), the initial pulse (e.g. a pollutant) diffuses into the
surrounding fluid (e.g. a river).
For x in the interval [0, 1] with initial condition, U (0, x) = sin(πx) and boundary
conditions, U (t, 0) = U (t, 1) = 0 for all t the solution is:
2
U (t, x) = e−βπ t sin (πx) .
(11)
(i) [1.0 Marks] Replace the time derivative by a first order forward difference and the
spatial derivative by a symmetric difference and display the resulting explicit FD equation
β∆t
for Uin+1 in terms of r = (∆x)
2.
(ii) [1.0 Marks] After some manipulation show that a von Neumann stability analysis
yields:
G = (1 − 2r) + 2r cos(k∆x) .
(12)
where k is the spatial wave number and G is the amplification factor (see e.g. equation
12.137 in Ouyed&Dobler on pjl).
(iii) [1.0 Marks] Show that for stability one must impose
∆t < ∆tmax =
(∆x)2
2β
(13)
(iv) [1.0 Marks] Before you proceed, check that the initial condition satisfies the FD
equation.
(v) [1.0 Marks] Run the FD equation for 300 time steps with 51 grid points with β = 1
and ∆t = 0.9∆tmax .
(vi) [1.0 Marks] Plot the initial condition, the solution from the FD (i.e. after 300 time
steps) and compare it to the analytical solution. Conclude.
1
The presence of a second order spatial derivative often indicates a diffusive process.
6
3.3
The advection-diffusion equation [6 Marks]
Now, we combined advection and diffusion together and get an advection-diffusion equation
ut + cux = βuxx ,
(14)
where in this case any initial concentration is advected at speed c while also diffusing into
the surrounding fluid at a rate given by β.
• [2.0 Marks] Discretize the advection-diffusion equation using a first order forward
difference for the time derivative, and a first oder backward difference for the first
space derivate (assuming c > 0) and a (second order) symmetric difference for the
second space derivative. Write down the resulting FD as un+1
= .... in terms of
i
r = c∆t/∆x and s = β∆t/(∆x)2 .
• [1.0 Marks] Use the Von Neumann stability analysis to show that ∆t ≤ ∆tmax =
(∆x)2 /(c∆x + 2β).
• [2.0 Marks] Solve the FD equation for c = 0.5, β = 0.1 for N = 101 grid points and
for ∆t = 0.9∆tmax and run it up to t = 57 seconds. Using gnuplot, and in the same
graph, plot the initial pulse and the pulse at t = 57 seconds.
• [1.0 Marks] How many time-steps were required to reach t = 57 seconds and how far
has the peak of the initial pulse has moved (in meters).
3.4
Conservative versus non-conservative forms [6 Marks]
Conservative forms of PDEs admit shock capturing. Conserved equations do not
have “source terms”. The source term usually is a driving force, a sinking term, or a
dissipative/friction term.
Imagine a force Fx pushing on a mass dm (with cross section area A), making it move in
the x-direction with velocity u(x, t). Then we can write
dux
dux
∂ux ∂x ∂ux
∂ux
∂ux
Fx = dm ax = dm
= (ρAdx)
= (ρAdx)
+
= (ρAdx)
+ ux
.
dt
dt
∂t
∂t ∂x
∂t
∂x
(15)
2
When Fx = 0 one gets the inviscid (i.e. without viscosity) Burgers’ equation:
∂u
∂u
+u
=0.
∂t
∂x
2
(17)
With viscosity Burgers’ equation is
∂u(x, t)
∂u(x, t)
∂ 2 u(x, t)
=0
+ u(x, t)
−ν
∂t
∂x
∂x2
where ν is the kinematic viscosity.
or equivalently
7
ut + uux − νuxx = 0 ,
(16)
where we dropped the subscript “x” since we are concerned with 1-dimensional (1D) Burgers’ equation.
• [1.0 Marks] Discuss the difference between the advection equation you saw earlier
and the Burgers’ equation above.
• [1.0 Marks] Rewrite equation above in its conservative form ∂u/∂t + ∂f (u)/∂x = 0
and give the expression for f (u).
• [2.0 Marks] Discretize the conservative form and the non-conservative form with the
Upwind method and the Lax-Wendroff method and display the FD forms.
• [2.0 Marks] How do the two forms of the equation (for the two FD schemes) behave
near x = 0 for a grid in −1.0 ≤ x ≤ 1.0 (chose a grid resolution ∆x that guarantees
stability when you integrate in time). Comment and conclude. Hint: Look at the
velocity profile for rather small time t in order to capture the onset of shock waves.
4
General Applications [18 Marks]
Many PDEs of interest contain second order (and higher) partial derivatives.
Hyperbolic equations describe many time-dependent (transient) phenomena (e.g.
fluid flow) and are characterised by a speed of propagation of information (often called
“wave speed”).
4.1
Fluid mechanics: time-dependent cases [8 Marks]
An equation frequently used as a model for fluid mechanics is the general Burgers’ equation
∂u
∂u
∂2u
+u
−ν 2 =0 .
∂t
∂x
∂x
(18)
Here, you will be solving it using the two-step Lax-Wendroff method (an explicit method)
for ν = 10−2 /π and with the following initial condition
u(x, 0) = − sin πx,
−1 < x < 1 .
(19)
Set the boundary conditions as u(−1, t) = u(1, t) = 0.
• [2.0 Marks] Rewrite the Burgers’ equation above into its conservative form ∂u/∂t +
∂f (u, ν)/∂x = 0 and give the function f (u, ν).
8
• [2.0 Marks] What order is the accuracy of the two-step Lax-Wendroff? Give the
stability condition of this method. Also, give the FD form of the conservative and
non-conservative Burgers’ equation above.
• [2.0 Marks] Here you are asked to display your results in a 2-panel figure using
multi-plot: Use the left panel for the conservative case and the right panel for the
non-conservative case. Plot u(x, t) versus the spatial position x (with −1.0 ≤ x ≤ 1.0)
at t = 0.0, 0.01, 0.5, 1.0 and 2.0 s. Discuss your figure. In particular, what happens
to the solution at x = 0 as time progresses using both schemes?
• [2.0 Marks] Discuss the accuracy and stability of the explicit two-step Lax-Wendroff
scheme in the context of the Burgers’ equation for the conservative and non-conservative
2
schemes. What does the viscous term (ν ∂∂xu2 ) adds when you compare your results to
the non-viscous case (i.e. ν = 0) you explored earlier?
4.2
Quantum mechanics: time-independent cases [6 Marks]
An equation frequently used as a model for Quantum mechanics is Schr¨odinger’s equation
∂2u ∂2u
+ 2 − αu = 0 .
∂x2
∂y
(20)
Here you will be comparing two methods to solve it: (i) The direct method; (ii) The
weighted Jacobi method. Use the following boundary conditions:
u(0, y) = u(1, y) = u(x, 0) = u(x, 1) = 1 .
(21)
You will be solving equation above for 0 ≤ x ≤ 1.0.
• [2.0 Marks] For α = 1.0, apply the direct method in a 50 × 50 grid and find the
solution. Using Gnuplot, plot u(x, y) versus (x, y). In the same 3D plot, project the
solution into the (x, y) plane and show it as contours.
• [2.0 Marks] Use the weighted Jacobi method in a 50 × 50 grid to explore the effect of
the weighting factor ω on the solution (do runs for ω = 0.1, 0.25, 0.5, 0.75, and 0.9).
Using Gnuplot, plot u(x, y) versus (x, y). In the same 3D plot, project the solution
into the(x, y) plane and show it as contours.
• [2.0 Marks] Redo the above for α = 1000. What do you conclude?
4.3
Econo-physics: time-dependent cases [4 Marks]
An equation frequently used as a model for stock option pricing is the Black-Scholes equation
1
∂2V
∂V
∂V
+ σ 2 S 2 2 − rV + rS
=0,
(22)
∂t
2
∂S
∂S
9
where V is the price of the option as a function of stock price S and time t, r is the risk-free
interest rate, and σ is the volatility of the stock.
Below you are asked to solve for V (S, t) in the region defined by 0 ≤ S ≤ 20 and for
0 ≤ t ≤ 100. The boundary conditions (i.e. the final stock value) are
(
0,
0 ≤ S < 10
V (S, t = 100) =
(23)
S − 10,
10 ≤ S ≤ 20
with volatility σ = 0.1 and return r = 0.05.
• [2.0 Marks] Chose an explicit FD method (justify your choice) to solve for the BlackScholes equation above. Give the resulting FD equation.
• [2.0 Marks] Using Gnuplot, plot V (S, t) versus S for t = 0.0, 50, 75, 85 and 100 s.
Discuss your figure. In particular, what happens to the solution at S = 10 as time
progresses? Describe any other behavior of the function you may notice.
5
The Riemann Problem/Solvers [16 Marks]
A Riemann problem consists of a step at x = 0 in parameter q with constant states
ql (left) and qr (right) on the left and right side.
(i) I encourage you to visit the VH-1 site:
http://wonka.physics.ncsu.edu/pub/VH-1/index.php
In particular read over their discussion of the Riemann problem at:
http://wonka.physics.ncsu.edu/pub/VH-1/bproblems.php
(ii) For this part of the lab, you will find a very useful paper at:
http://arxiv.org/pdf/1303.3999v1.pdf
• [2.0 Marks] Among the most interesting and difficult problems in computational fluid
dynamics is the simulation of discontinuities like shock fronts. Simple finite difference
schemes cannot handle this type of singular behavior. Explain why?
• [7.0 Marks] A special case of a Riemann problem in Eulerian hydrodynamics is when
the initial state has zero velocity (v = 0) on both sides of the dividing point, but the
pressure has a jump. If one solves such problems using numerical hydrodynamics one
finds that the self-similar solution that follows from such a problem typically has 5
regions which we shall call region 1,2,3,4,5.
(i) Write a program to find the exact solution to this Riemann problem (i.e. the Sod
shock tube problem) at t = 5000 for γ = 7/5, ρl = 105 , Pl = 1, ρr = 1.25 × 104 and
Pr = 0.1. Here the subscript “l” stands for “left side of the Sod tube” while “r”
10
stands for“right side of the Sod tube” with ρ and P being the density and pressure,
respectively; γ is the adiabatic index.
(ii) Using Gnuplot, plot the solution and identify the 5 regions described above.
Your plot should show the density (upper panel), the pressure (middle panel) and
the energy (lower panel) at t = 5000. Describe these regions and how the physical
quantities evolve in space and time for the parameters given.
• [2.0 Marks] Among the most popular Riemann solvers (i.e. FD methods making
use of the Riemann idea) are the Godunov schemes and the van Leer schemes.
Describe these in details and in particular why and how the second order van Leer
method is superior. You will find a nice introduction to the MUSCL3 scheme at :
http://en.wikipedia.org/wiki/MUSCL_scheme
The second-order van Leer MUSCL schemes are actually one of the most
popular high order scheme for fluid dynamic computations.
• [5.0 Marks] For the same parameters given in the Sod problem above (v = 0) with
γ = 7/5, ρl = 105 , Pl = 1, ρr = 1.25 × 104 and Pr = 0.1, write a code to compute the
numerical solution a t = 5000 using the :
(a) Explicit Euler time-stepping scheme with a CFL number of 0.5 on a grid of 100
uniform cells
(b) Godunovs method.
(c) van Leer Flux Vector-Splitting (FVS)4 method.
(d) Plot the density, pressure and energy profiles and compare these three numerical
solutions with the exact solution you solved for above. The analytical solution should
be shown as solid lines while the different FD schemes using symbols.
(e) Make an observation about the quality of the numerical solutions obtained by
these three methods.
3
4
MUSCL stands for Monotonic Upstream-Centered Scheme for Conservation Laws (van Leer, 1979)
Split flux such that
F (U ) = F + (U ) + F − (U ) ,
(24)
with F + due to waves moving to right, (positive speed) and F − due to waves moving to left (negative
speed.)
11