Numerical Solution of Partial Differential Equations (446-1)

Project 1

Due November 6, 2000

I. Various Ways for the Two-Way Wave Equation

Consider the wave equation

¶t u = 1 r(x) ¶x v,      ¶t v = r(x) c(x)2 ¶x u

(1)

with periodic boundary conditions on a domain [0,1]. Note that the density r and the velocity c(x) will be allowed to depend on the position x. This could, for instance, describe the propagation of sound waves in a medium with space-dependent density and elastic constant.

Solve (1) numerically using different schemes as described below. In each case use as an initial condition a pulse with width D

u(x,0) = A e-[((x-0.5)2)/(2D2)],        v(x,0) = A c r(x) e-[((x-0.5)2)/(2D2)]

(2)

with A = 0.3. Let the solution evolve from t = 0 to t = 1. In each time step measure the total energy of the wave which is given by

E = ó õ 1 0  r(x) c(x)2 u(x)2+ 1 r(x) v(x)2  dx.

(3)

In order to get a good idea of what each scheme is doing it is very useful to plot the temporal evolution of the solution during the computation. In addition to the numerical solution plot also

ue(x,t) = A e-[((x+ct-0.5)2)/(2D2)],        ve(x,t) = A c re-[((x+ct-0.5)2)/(2D2)],

(4)

which would be the exact solution in an infinite domain, in the graphics window. Pay attention to the proper implementation of the periodic boundary conditions for this (essentially) exact solution. Note that it is useful to define functions uex and vex to calculate the exact solution. Keep also track of the temporal evolution of the total energy E.

1. Implement the (2,2)-MacCormack scheme to solve (1) for r = 2 and c = 1 with width D = 0.03. For the simulations choose Dt close to the stability limit, i.e. l = 0.95 fixed.

   Make a graph of three solutions at t = 1 with different Dx which demonstrates that the dispersion disappears with Dx ® 0. Judged from your simulations what is the sign of the phase error?

   Decrease Dx in factors of 2 until you obtain the final state at t = 1 to `graphical accuracy', i.e. when the whole solution is plotted (unzoomed) decreasing Dx by a factor of 2 does not show any change. Note that this value of Dx need not be determined precisely, just proceed in factors of 2. Demonstrate your result in a plot overlaying the relevant solutions.

   Plot the temporal evolution of the energy E for a reasonable choice of l and Dx. Judged from your simulations what is the sign of the phase error?

   Similarly, determine the Dx₁ that gives an accuracy of at least 1%. For this define the error E₂ by comparing the numerical solution with the exact solution u_e(x,t = 1) where the error (squared) is defined as

   
   

   E22 = N å j = 1  (ujn-ue(xj,t = 1))2 N å j = 1  ue(xj,t = 1)2 .

   (5)

   Here u_jⁿ stands for the numerical solution at t = n·Dt and position x = j·Dx.

   Plot E₂(Dx) double-logarithmically for l = 0.95 over a range in Dx which demonstrates the correct convergences of your implementation of the scheme.

2. Implement the (2,4)-MacCormack scheme to solve (1) for r = 2 and c = 1. Use again as initial condition a pulse with width D = 0.03.

   Measure the error E₂(Dx) over a range of Dx for various values of l. Plot the error double-logarithmically plotting the results for E₂(Dx) for the different values of l into the same graph. Choose the range in Dx and the values of l such that the plot reflects the convergence properties expected from the truncation error of the scheme.

   Based on the theoretical truncation error, for which l as a function of Dx do you expect the scheme to run most efficiently, i.e. for given accuracy and maximal time t_max, for which l do you expect the computational effort to be smallest? Check whether this expectation is borne out in your simulation, i.e. compare the computational effort of runs that were obtained with different values of l but delivered essentially the same accuracy.

   For a reasonable choice of l and Dx, plot the temporal evolution of the energy E. Judged from your simulations what is the sign of the phase error?

3. Implement the (4,4)-Runge-Kutta scheme to solve (1) for r = 2 and c = 1. Use again as initial condition a pulse with width D = 0.03. Measure the error E₂(Dx) over a range of Dx for various values of l. Plot the error double-logarithmically to demonstrate that your code exhibits the correct order of convergence, again plotting the results for E₂(Dx) for the different values of l into the same graph. How does the error depend on l at fixed Dx?

   For a reasonable choice of l and Dx, plot the temporal evolution of the energy E. Judged from your simulations what is the sign of the phase error?

II. Reflections on a Rope

Consider now the wave equation (1) with variable density r(x) and/or phase velocity c(x).

In each of the cases below, investigate (1) numerically using (2,2)- and (2,4)-MacCormack, and (4,4)-Runge-Kutta. As initial condition use again a pulse with width D.

1. Study the reflection off a small inhomogeneity in the density.

   
   

   r(x) = 2        0 < x £ 0.3        and        0.32 < x £ 1, r(x) = 2.01        0.3 < x £ 0.32.

   Launch a narrow pulse with D = 0.005 at x = 0.5 and measure the reflected signal u(x = 0.5,t). From this signal one can determine the structure of the inhomogeneity. In particular, its location X and width L are obtained by comparing the delay of its arrival time at x = 0.5 with the known velocity of pulses in this non-dispersive wave equation. Vary Dx over a range of values, Dx = 0.01, 0.005, 0.0025, to determine X and L from the reflected signal. To obtain these values, determine the timing of the relevant extrema of u(x = 0.5,t). Use a polynomial interpolation to determine the maxima accurately (e.g. with polyfit). Do this for each of (2,2)- and (2,4)-MacCormack, and (4,4)-Runge-Kutta.
   Use the scheme that gives the best results in these runs to resolve the reflected pulse graphically (when zoomed in to show only the reflected pulse at the time when it returns to x = 0.5).
   Note: Similar methods using sound waves can be used for exploration of structures in the earth (earth quakes, location of oil fields ...) or in the sonography of tissue.

2. Consider now a large jump in the density

   
   

   r(x) = 4        0 < x £ 0.36,        r(x) = 2        0.36 < x £ 1.

   (6) Start with a pulse of width D = 0.03 at x = 0.5 and evolve it to t = 0.4. How do you have to choose l?
   Decrease Dx to achieve graphical accuracy of the solution. Show this in one plot for each scheme. Do you notice a difference in the performance of the MacCormack and the Runge-Kutta scheme? If so, what is its origin?
   To see the effect of multiple scattering (reflection) follow the evolution of the pulse up to t = 4 using (2,2)-MacCormack and Dx = 0.005 and plot u(x) at two representative intermediate times and at t = 4. Plot the temporal evolution of the energy E for this run.

3. Now investigate waves in a stratified medium in which the wave speed depends on space. Use as a model the above rope with wave speed

   
   

   c(x) = 1 + 0.9 sin(2 px).

   (7)

   The density is now assumed constant, r = 2.
   Run the above schemes with a pulse of width D = 0.03 up to t = 1. How should you choose l? Can you achieve graphical resolution with each of the schemes? Determine u(x = 0.75,t = 1) with 1% accuracy using the scheme that is best suited for this problem. Comment on the difference in the performance of the three schemes. What is its origin?