Numerical Differential Equations Homework 3
(Due: Jul. 26, 2007) Consider the initial value problem
y0 =−µ(y − t2) + 2t, 0 < t < T, y(0) = y0. Recall that the exact solution is
y(t) = y0e−µt+ t2. Set µ = 520, T = 2 and y0 = 4.
A For this IVP, write (and present) three short MATLAB scripts implementing the following fixed step-size methods:
– the forward Euler method;
– the backward Euler method;
– the trapezoidal rule;
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 yn+1 from yn in the two implicit schemes. Use this fact in your codes.
B Stability The ”test equation” for the above ODE is
y0 =−µy.
The results about the regions of absolute stability 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 en 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 y0 = y(1) = 4e−520+ 1. You should see that even though on this interval [a, b] the solution is essentially t2, 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 y0 = 2t on the interval [1, 2] with y0 = 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 |en| 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?
