Numerical Partial Differential Equations 1

Homework 1 (Due: October 13, 2005)

• Problem 1
  Chapter 2, page 27, problem 2 of the "Flaherty notes."

• Problem 2
  Chapter 2, page 27, problem 3 of the "Flaherty notes."

• Problem 3
  Consider the initial value problem from Section 2.2 of the "Flaherty notes,"
  y' = -µ (y - t²) + 2t,
  
  0 < t < T, y (0) = y₀ .

  Recall that the exact solution is y (t) = y₀ e^{-µ t} + t².

  Set µ = 520, T = 2 and y₀ = 4.

  A. For this IVP, write (and present) four short MATLAB scripts implementing the following fixed step-size methods:

  • the forward Euler method (2.1.2);
  • the backward Euler method (2.2.3);
  • the trapezoidal rule (2.2.7);
  • "the predictor-corrector" trapezoidal rule (2.2.8).

  Note that because the right-hand side of the ODE is linear in y, you do not need to solve a nonlinear equation to obtain y_n+1 from y_n in the two implicit schemes. Use this fact in your codes.

  B. Stability

  The "test equation" for the above ODE is y' = -µ y.

  The results about the regions of absolute stability (that were derived in the class and in Problem 2, above) should give you an indication what the needed step-size to achieve absolute stability is for each of the above methods. Use your MATLAB codes with various time steps (e.g., of the form 2^-k ) to see whether the results of your numerical experiments correspond to the theory. Present your results in a graphical form. Plot the exact and the numerical solution in the same picture. Also plot the global error e_n or its absolute value. Discuss the results.

  As part of this problem, run your methods for the same equation on the interval [a, b] = [1,2] with the exact initial condition y₀ = y (1) = 4 e^-520 + 1. You should see that even though on this interval [a, b] the solution is essentially t², the forward Euler method still requires a small step-size to get reasonable results. What about the other methods?

  (As an illustrative example, you may want to change your codes so that they solve the equation y' = 2t on the interval [1,2] with y₀ = 1. Observe that in this case the numerical methods do not exhibit any of the instability features.)

  C. Order of Convergence

  For each of the four methods, establish numerically the order of convergence by computing the maximum of the global errors |e_n | and observing its behavior as h converges to 0. Do the results correspond to the theory? Present the results in a graphical form or in a table.

  For the trapezoidal rule, you should notice that in most of the interval [0,2] the global error is almost negligible. What is the order of the error? Can you explain why?