data留学生编程讲解、辅导Python程序、Java，c++程序语言辅导
解析Haskell程序|解析C/C++编程
Lab assignment #8: Partial Differential Equations, Pt. I
Due Friday, November 6th 2020, 5 pm
First, try to sign up for room 1. If it is full, sign up for room 2 (and ask your partner to join you there
if they are in room 1).
General Advice
• Work with a partner!
• This lab’s topics revolve around computing solutions to PDEs using basic methods.
• Read this document and do its suggested readings to help with the pre-lectures and labs.
• Ask questions if you don’t understand something in this background material: maybe we can
explain things better, or maybe there are typos.
• Carefully check what you are supposed to hand in for grading in the section “Lab Instructions”.
• Whether or not you are asked to hand in pseudocode, you need to strategize and pseudocode
before you start coding. Writing code should be your last step, not your first step.
• Test your code as you go, not when it is finished. The easiest way to test code is with print('')
statements. Print out values that you set or calculate to make sure they are what you think
they are.
• Practice modularity. It is the concept of breaking up your code into pieces that as independent
as possible from each other. That way, if anything goes wrong, you can test each piece
independently.
• One way to practice modularity is to define external functions for repetitive tasks. An external
function is a piece of code that looks like this:
1
PHY407F Lab08 Due 6 November 2020, 5pm
def MyFunc(argument):
"""A header that explains the function
INPUT:
argument [float] is the angle in rad
OUTPUT:
res [float] is twice the argument"""
res = 2.*argument
return res
Place them in a separate file called e.g. MyFunctions.py, and call and use them in your answer
files with:
import MyFunctions as fl2 # make sure file is in same folder
ZehValyou = 4.
ZehDubble = fl2.MyFunc(ZehValyou)
Computational background
Boundary value problem methods (for Q1): The first class of PDEs that Newman discusses in
Chapter 9 are solutions of elliptic partial differential equations of which the Poisson equation (9.10)
is a classic example. He discusses the Jacob method, which is a relaxational method for solving
Laplace’s or Poisson’s equation, and then speedups to this method using overrelaxation and
Gauss-Seidel replacement. For this lab we would like you to focus on Gauss-Seidel, with or without
overrelaxation. To obtain the solution of a field φ(x, y) that satisfies Laplace’s equation, (1)
subject to boundary conditions (see Physics Background), overrelaxation is written [see (9.17)], (2)
where the left arrow indicates replacement of the left hand side with the right, and ω is a relaxation
parameter. For these methods, you can either pick a fixed number of iterations or can choose a
level of convergence for the solution. In Q1a, we will be using the latter approach.
Stream plots show streamlines of a vector field F~. That is, at each location, F~(~x) is tangential to
the local streamline. For example, for an electric field E~, streamlines are electric field lines. Matplotlib’s
streamplot function (https://matplotlib.org/3.2.1/api/_as_gen/matplotlib.pyplot.
streamplot.html) can do such a plot, and colour-code the lines with the scalar of your choice. See
example below for electric field lines due to the presence of two point charges ±4π²0, separated by
a distance d = 2. It is adapted from the example at the bottom of the web page above.
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-2, 2, 100) # does not include zero
y = np.linspace(-1, 1, 50)
X, Y = np.meshgrid(x, y)
R1 = ((X-1)**2 + Y**2)**.5 # 1st charge located at x=+1, y=0
R2 = ((X+1)**2 + Y**2)**.5 # 2nd charge located at x=-1, y=0
2 of 9
PHY407F Lab08 Due 6 November 2020, 5pm
V = 1./R1 - 1./R2 # two equal-and-opposite charges
Ey, Ex = np.gradient(-V, y, x) # careful about order
fig = plt.figure(figsize=(6, 3))
strm = plt.streamplot(X, Y, Ex, Ey, color=V, linewidth=2, cmap='autumn')
cbar = fig.colorbar(strm.lines)
cbar.set_label('Potential $V$')
plt.title('Electric field lines')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.axis('equal')
plt.tight_layout()
plt.show()
Time/space PDE problems For Q2-Q3, we will introduce and implement some techniques for
solving time-dependent partial differential equations. For simplicity, we will only consider examples
with one space dimension, thus variables will depend on x and t. Concepts include finite
differences, von Neumann stability analysis, the forward-time centred-space (FTCS) numerical
scheme, and the Lax-Wendroff numerical scheme.
Note that the Lax-Wendroff numerical scheme introduced in this lab is not discussed in the Newman
textbook, but online resources exist, simply starting with https://en.wikipedia.org/wiki/
Lax%E2%80%93Wendroff_method.
A large class of time-dependent PDEs can be written in flux-conservative form. In one space dimension, (3)
where, in general, ~u and F~ are vectors consisting of a set of multiple fields. The variable we are
solvoing for is ~u, while F~ is a given function that can depend on ~u and on spatial derivatives of ~u.
The simplest scheme to discretize this system is the forward-time centred-space (FTCS) scheme, in
which one uses a forward difference for the time derivative, and a centred difference for the spatial
derivative. If, for example, the problem is one dimensional and ~u and F~ have just one component
field,(4)
where superscripts like n refer to the time step index and subscripts like j refer to the spatial index
(see figure ).
Figure 1: Notations for spatial grid.
3 of 9
PHY407F Lab08 Due 6 November 2020, 5pm
Substituting these approximations into eqn.(5)
However, as you will show below, this very simple discretization is unconditionally unstable for typical
wave equations. You will implement an additional scheme, the Lax-Wendroff method, which
is second-order accurate in time and has improved stability. This method will be introduced in the
Q3.
Animations in Matplotlib Animations are possible with Matplotlib. Here is a sample code that
will animate the curve sin(x − t).
from numpy import arange, pi, sin
from pylab import clf, plot, xlim, ylim, show, pause, draw
t = arange(0, 4*pi, pi/100) # t coordinate
x = arange(0, 4*pi, pi/100) # x coordinate
for tval in t:
clf() # clear the plot
plot(x, sin(x-tval)) # plot the current sin curve
xlim([0, 4*pi]) # set the x boundaries constant
ylim([-1, 1]) # and the y boundaries
draw()
pause(0.01) #pause to allow a smooth animation
The script simply clears each frame and then replots it, with a little pause to make it run smoothly.
This is not an ideal method but gets the job done. A “cleaner” method would be to follow the
guidelines of
https://matplotlib.org/api/animation_api.html,
but it does represent a bit more work.
Physics background
Electrostatics. In empty space, the electrostatic potentialV (x, y) satisfies Laplace’s equation ∇
2V =
0. In Q1, we will use relaxation methods as in Newman’s Chapter 9 to solve forV given some boundary
conditions.
Shallow Water Equations (for Q2-3). The shallow water equations are a simplification of the full
Navier-Stokes equations, which describe the motion of fluid flow (see § 3.1 of Atmospheric and
Oceanic Fluid Dynamics by Geoffrey K. Vallis, available on the UofT library website at http://
go.utlib.ca/cat/11621853, for a basic introduction; we will not consider the Earth’s rotation
and you can therefore ignore f , the Coriolis parameter). They are appropriate for describing the
motion in shallow layers of fluid under the influence of gravity. By “shallow” we mean that the
depth of the fluid is small compared to the relevant horizontal length scales of the fluid. Thus, in
some cases, the ocean can be considered a shallow fluid, as long as we are interested in phenomena
whose horizontal length scales are much larger than the depth of the ocean (on average, around
4 km). In particular, it is common to use the shallow water equations for modelling the tides, and
4 of 9
PHY407F Lab08 Due 6 November 2020, 5pm
Fig. .
Figure 2: Variable definitions for the shallow water system. Figure shows a cross-section on the
(x, z) plane. The upper dashed line represents the equilibrium surface, the lower dashed line the
z = 0 origin, H is the average (in time and space) water column height, the top solid line represents
the free surface of the fluid at z = η(x,t), and the bottom solid line represents the fixed bottom
topography at z = ηb(x). From Vallis (2017, fig. 3.1).
for simulating tsunamis in the ocean. In this lab, we will attempt to do the latter, albeit in a highly
simplified setting.
For simplicity, we will restrict ourselves to the 1D version of the shallow-water equations, which
are (see also fig. 2 for notations)(6)
where u is the fluid velocity in the x directions, h = η−ηb is the water column height, η the altitude
of the free surface, ηb is the altitude of the bottom topography, which is typically a known function,
and g is the acceleration due to gravity. Note that all the variables are functions of x and t, but not
z (i.e. every column of fluid moves together). This is a consequence of the assumptions that lead
to the shallow water equations.
In particular, these equations support wave solutions. These waves are non-dispersive, their phase
and group speeds over flat bottom being equal to p
g H, where H is the water column height at
rest.
Questions
1. [25% of the lab] Electrostatics and Laplace’s equation (adapted from Newman, ex. 9.3, p.
417)
(a) Consider the simple model of an electronic capacitor sketched in fig. 3, consisting of
two flat metal plates enclosed in a square metal box.
For simplicity let us model the system in two dimensions. Using the Gauss-Seidel method
without overrelaxation, write a program to calculate the electrostatic potential in the
5 of 9
PHY407F Lab08 Due 6 November 2020, 5pm
10 cm 6 cm
+1 V −1 V
0 V
2 cm 6 cm 2 cm
Figure 3: Capacitor for Q1 (from Newman, p. 417).
box on a grid of 100×100 points, where the walls of the box are at voltage zero and the
two plates (which are of negligible thickness) are at voltages ±1V as shown. Have your
program calculate the value of the potential at each grid point to a precision of 10−6
volts and then make a contour plot of the potential (filled or not), and a stream plot
of the electric field lines, colour-coded by the value of the electric potential. You may
overlay them on the same figure (careful about colour schemes and readability!) or plot
two separate figures.
Note: This exercise is similar to Newman’s Ex. 9.1, which you don’t have to do. In effect,
the capacitor plates are part of the boundary condition in this case: they behave the same
way as the walls of the box, with potentials that are fixed at a certain value and cannot
change.
HAND IN PLOTS AND WHICHEVER COMMENTS ARE RELEVANT.
(b) Now, implement the same method with overrelaxation. Try with ω = 0.1 and ω = 0.5.
What do you notice?
SUBMIT YOUR WRITTEN ANSWERS, AND CODE.
2. [50% of the lab] Simulating the shallow water system, Part I
(a) Problem set up: Show that the 1D shallow water Equations (6) can be rewritten in the
flux-conservative form of eqn. 3, with ~u = (u,η). (7)
Now use the FTCS scheme (eqn. 5) to discretize the 1D shallow water equations. You
should have two equations of the form u
n+1
j
= ... and η
n+1
j
= ... where the right-hand
sides depend on u and η at the timestep n. Assume ηb(x) is a known function (i.e.,
ηb(xj) = ηb,j
is given).
HAND IN YOUR DERIVATION.
6 of 9
PHY407F Lab08 Due 6 November 2020, 5pm
(b) FTCS Implementation: Implement the 1D shallow water system with the FTCS scheme
in Python. Use the following set up (perhaps simulating waves in a very shallow one
dimensional flat bathtub!):
• Take your domain as x = [0,L], with L = 1 m, and use a grid spacing as shown in
fig., with J = 50, so ∆x = 0.02m.
• Take g = 9.81 m s−2
• Assume ηb = 0, H = 0.01 m (flat bottom topography).
• As boundary conditions, you should impose u(0,t) = u(L,t) = 0, i.e., rigid walls.
Note that you will have to treat the grid points at the boundaries separately, because
you cannot use a centred difference approximation here, since you do not
have values for u or η at j = −1 or j = J + 1. To get around this issue, use forward
and backward differences at each of the boundary points to approximate the spatial
derivative of the fluxes:
with A = 0.002m, µ = 0.5 m, σ = 0.05 m, and where 〈〉 is the x-average operator,
ensuring that H retains its definition of the free surface altitude at rest. This represents
a small gaussian peak at the centre of the domain
See the background material for an example of how to create an animation out of a
series of lineplots, which will be helpful for visualizing your results. As part of your
write up, you should include plots of η(x,t) at t = 0 s, t = 1 s and t = 4 s.
Note 1: you may (actually, should) find that your solution eventually diverges. Do not be
alarmed, this doesn’t necessarily mean you have an error in your code!
Note 2: In the finite difference scheme, you need to use only the previous time step values
to update the variable. The easiest way to do this is by defining second arrays for u and η
(e.g., called u_new and eta_new) and then your update lines in the loop can look like:
u_new[j] = u[j] + ... # (all in terms of u and eta)
eta_new[j] = eta[j] + ... # (all in terms of u and eta)
eta = np.copy(eta_new)
u = np.copy(u_new)
It is important to use the np.copy command in the last 2 lines, otherwise python just
updates eta at the same time as eta_new in the next iteration. This is because Python
treats the equality sign as assignment to a pointer when the object on the right-hand side
is a list or Numpy array.
HAND IN YOUR CODE, PLOTS, EXPLANATORY NOTES.
7 of 9
PHY407F Lab08 Due 6 November 2020, 5pm
(c) Do a von Neumann stability analysis of the FTCS scheme for the 1D shallow water equations
(as described in § 9.3.2, pp. 425-430 of the Newman’s textbook). Note that this type
of stability analysis can only be applied to linear equations, so first you need to show
that the linearized eqns (6)a,b about (u,η) = (0,H) are, with constant topography ηb,
Then you can closely follow the development in the textbook for the stability analysis
of the wave equation (pp. 429–430) to show that the magnitude of the eigenvalues λ is,
for both eigenvalues,
Do you expect the FTCS scheme to be stable? Why or why not?
HAND IN YOUR DERIVATION.
3. [25% of the lab] Simulating the shallow water system, Part II
We can do much better than the FTCS scheme. As one example (of many possible computational
schemes), in this question you will implement the Two-Step Lax-Wendroff scheme.
There is a slightly higher computational cost for this scheme, but it leads to improved stability
and accuracy.
The first step is to calculate the values at the “half-steps” tn+1/2 and xj+1/2 using the Lax
method, i.e.,
Once you have the values of u and η at the half steps, you can calculate the fluxes F
n+1/2
j±1/2 at
the half steps. Then, calculate the values at the next full timestep (and at the “normal” grid
points) using the simple FTCS scheme. (13)
(a) Implement this scheme, and run it using the same set up described above. Again, save
figures at t = 0, 1, 4 s.
HAND IN YOUR PLOT(S) AND BRIEF COMMENTS.
(b) As a final step, you should allow for variable bottom topography. This should be a relatively
straightforward change to your code. Instead of having ηb be a constant, it will be
an array defined on your grid. Thus, when you calculate the fluxes F
n
j
you need to make
sure you are using the value of the topography at the grid point j (in fact, for the LaxWendroff
scheme you will need ηb at the half-steps j +1/2, so you should just estimate
it at these points as the mean between the two nearest integer points).
You should now be able to simulate a 1D tsunami. Try the following setup.
• Domain: x = [0, 1], J = 150, so ∆x = 1/150.
• Average altitude of free surface H = 0.01 m,
8 of 9
PHY407F Lab08 Due 6 November 2020, 5pm
0.0 0.2 0.4 0.6 0.8 1.0
x [m]
0.000
0.005
0.010
b(x) [m]
Figure 4: Bottom topography for tsunami simulation.
• Use a bottom topography that goes from ηb = 0 on the side of the abyssal ocean,
has a sharp change in depth in the middle of the domain, simulating the presence
of a continental shelf break, to end up at ηb = H −4×10−4 m (see Figure 4), namely,
ηb(x) =
ηbs
2
{1+tanh[(x − x0)α]}, (14)
with ηbs = H −4×10−4 m, α = 8π m−1
and x0 = 0.5 m.
• Boundary conditions: u(0,t) = u(L,t) = 0, i.e., rigid walls.
• Timestep: ∆t = 0.001; number of iterations: 5,000 at least.
• Initial conditions (wide gaussian bump at far left boundary, i.e., in deep ocean),
with A = 2×10−4 m, σ = 0.1m, and where 〈〉 is the same as in eqn. (9).
You should see a shallow, wide bump starting at the left of the domain move to the
right, towards the continental shelf. What happens to the waveform when it reaches
the shallower waters (in terms of its shape, height, speed, etc.)?
Make plots of η(x,t) (you might find it helpful to also plot the bottom topography on
the same figure) at t = 0, t = 1, t = 2, and t = 4 s.
HAND IN YOUR CODE, PLOTS, EXPLANATORY NOTES. EVEN THOUGH THE CODE IS A
MODIFICATION OF THE PREVIOUS CODE, PLEASE DO NOT TRY TO COMBINE THEM INTO ONE,
AND HAND IN SEPARATE SCRIPTS FOR Q2 AND Q3, FOR MARKING PURPOSES.
9 of 9