In Matlab, we might be tempted to implement this method with periodic boundary conditions on the grid xi = ih for i = 1,2,··· ,m + 1 with h = (b − a)/(m + 1) within a time-stepping loop as
U(1) = U(1) − 0.5 ∗ nu ∗ (U(2) − U(m + 1))...
+0.5 ∗ nu^2 ∗ (U(m + 1) − 2.0 ∗ U(1) + U(2));
U(2 : m) = U(2 : m) − 0.5 ∗ nu ∗ (U(3 : m + 1) − U(1 : m − 1))...
+0.5 ∗ nu^2 ∗ (U(1 : m − 1) − 2.0 ∗ U(2 : m) + U(3 : m + 1));
U(m + 1) = U(m + 1) − 0.5 ∗ nu ∗ (U(1) − U(m))...
+0.5 ∗ nu^2 ∗ (U(m) − 2.0 ∗ U(m + 1) + U(1));
Note that periodic boundary conditions mean that U0 = Um+1, so we do not solve for U0 in the equations above. However, this implementation is incorrect. Can you find the bug?
2. Consider the numerical scheme:
(a) With which system of PDEs is this scheme consistent?
1
Math 232: Numerical Analysis II Homework 9
(b) Find the modified equations for this scheme.
(c) Using von Neumann stability analysis find conditions on k and h sufficient for this scheme to be stable.
3. The m-file advection LW pbc.m implements the Lax-Wendroff method for the advection equation on 0 ≤ x ≤ 1 with periodic boundary conditions.
(a) Observe how this method behaves with m + 1 = 50,100,200 grid points. Change the final time to tfinal = 0.1 and use the m-files error table.m and error loglog.m to verify second order accuracy.
(b) Modify the m-file to create a version advection up pbc.m implementing the upwind method and verify that this is first order accurate.
(c) Keep m fixed and observe what happens with advection up pbc.m if the time step k is reduced, e.g. try k = 0.4h, k = 0.2h and k = 0.1h. When a convergent method is applied to an ODE we expect better accuracy as the time step is reduced and we can view the upwind method as an ODE solver applied to a Method of Lines system. However, you should observe decreased accuracy as k → 0 with h fixed. Explain this apparent paradox. (Hint: What ODE system are we solving with more accuracy? You might also consider the modified equation given by (10.44) in LeVeque).
