MATH281-Project 3 Solved

Project Outline
The term project is divided into three section, each with a corresponding Jupyter workbook. Each workbook will ask you to investigate various aspects of numerical methods or PDEs

Take care with your graphs. Some things to consider are:

labeling (what is plotted on the axes? What are their units [if any]?)
legibility (labels, axis numbers, etc should be large enough to be easily read when included in the report
captions are a useful way to give a (brief) explanation of what is shown in a figure.
Workbook Questions [45 marks in questions + 5 general presentation = 50 points total]
In the third and final component of the term project, you will investigate the dynamical behaviour of three physical mechanisms represented by PDEs. You will apply the tools from the previous workbooks in order to numerically solve the PDEs, and then analyze the behaviour of the system in both physical- and spectral-space. The idea is to qualitatively understand the dynamical impact of each mechanism / mathematical term. Understanding how each component of a PDE impacts the governing dynamics is critical to gaining insight into the over-all behaviour of the system reflected by the PDE.

Assume periodic boundary conditions for each question in the workbook. The initial conditions and domain for each question are given by

u = u(x,t) u(x,0) = sin(4πx/Lx) + 0.25sin(8πx/Lx) x ∈ [−Lx/2,Lx/2] t ∈ [0,50]

Lx = 10

Note about Simulations keep in mind that for at least the first two problems, you know qualitatively what kind of solution to expect. And so you are expected to recognize if your solution explodes or gives results wildly outside of the expected behaviour. Analysing simulation results is only meaningful if the simulation results themselves are sound.

Note about PDES notice that between the three cases, we are only changing the order of the spatial derivative. The physical mechanisms described by the three PDEs are often called, in order, advection, diffusion, and dispersion.

Note about Grading in addition to 45 marks spread across the three questions, there are five additional marks for over-all presentation.

[15 marks]: Question 1: Linear Advection
                                                                 ∂tu + U0∂xu = 0;                                                 U0 = 0.1

Substitute the ansatz u = exp(i(kx − ωt)) into the PDE and solve for ω as a function of k[1].
Substitute back into the ansatz and write it as a product of a term that oscillates in time and a term that decays/grows in time. Note that one or neither of these terms may simply be the constant 1 (one).

Using your understanding of representing functions as Fourier sums, how / what does this analysis reveal about the behaviour of the system? What temporal behaviour does it describe? Does it depend on wavenumber k?

Using a second order finite difference scheme, compute a numerical solution to the PDE. Explain and justify your choice ofnumber of spatial gridpoints (Nx), and (ii) size of timestep (∆t).
Your choices should be such that the spatial features are cleanly resolved, and the time-step is sufficiently small to provide stability[2]. Do not assume that the values given in the sample workbook are sufficient.

Create appropriate plots to show the evolution of the system. Qualitatively describe how the system evolves. How do these results compare with your Fourier analysis from part (a)?

Compute (and plot) the root-mean-squared[3] (RMS) error as a function of time. What does this tell you about the validity of your simulation? For reference / comparison, compute the RMS of the initial conditions.
Compute and describe (and plot) the time evolution of the spectral power. How do these results compare to what you expect given your knowledge of the PDE? If there are any differences, can you explain / justify them using what you know about numerical methods?
(i) Using the PDE, obtain an ODE for E(t) := Rx u2(x,t)dx, the domain integral of u2. Is this result
consistent with your Fourier analysis in part (a)?

Compute (and plot) the spatial integral of E(t) from your numerical solutions. Is this consistent with your predictions? Why or why not?
[15 marks]: Question 2: Diffusion
                                                                     ∂tu = ν∂xxu;                                                    ν = 0.1

The CFL restriction that we use in advective problems doesn’t apply to diffusive processes, since information isn’t propagating in the same sense.
Consider again the linear advection equation, but this time with the derivatives replaced by finite differences:

.

Simple re-arranging then gives

,

which was the upper bound of our CFL restriction.[4]

Using a similar method, what result do you find for the diffusion equation?

Remembering that we are assuming a periodic grid and uniform spacing, derive the centred, secondorder finite difference scheme for second derivatives. Recall that you can use the error-curve plots from the previous two workbooks to verify the order of your scheme. Include the error curve (error vs. ∆x) as evidence of the validity of your scheme.
To do this, you will need to write the Taylor expansions for ui−1,ui,ui+1 out to the fourth-order term.

An appropriate code-snippet to compute this derivative would then be:

d2udx2 = ( np.roll(u, 1) - 2 * u + np.roll(u, -1) ) / ( delta x**2 )

Repeat part (a) from question 1.
Repeat part (b) from question 1.
Note that these simulations may take several minutes depending on your choices of ∆x and ∆t. You may find it beneficial to save your simulation results after computing them to avoid unnecessary recomputing.

Repeat part (e) from question 1.
[15 marks]: Question 3: Dispersion
                                                                    ∂tu = α∂xxxu;                                                    α = 0.1

Repeat parts (a) - (e) from question 2.

For part (b), you will need to write the Taylor expansions for ui−2,ui−1,ui,ui+1,ui+2 out to the fifth-order term. An appropriate code-snippet to compute this derivative would then be:

d3udx3 = (

5 * np.roll(u, 2) + np.roll(u, 1)
roll(u, -1) + 0.5 * np.roll(u, -2)
) / ( delta x**3 )

When analyzing your numerical simulations, you may find it beneficial to run additional simulations that include the sin terms from the initial conditions separately.

[1] this equation is known as the dispersion relation

[2] A good way to test your choice of ∆t is to run a second simulation with 0.5∆t (but the same Nx). If the results are nearly identical [and no obviously bad things are happening], then you can have some faith in your results.

r

[3] 3RMS =              meanh(utrue − unumerical)2i

[4] Recall that we needed to use a substantially smaller time step for the purpose of stability, but this gives us a reference point to work from, and tells us how to scale the time-step when the grid changes.