Mathematics 172 – Numerical Analysis 2 Discussion 1 – First Numerical Methods for First Order IVP’s February 2, 1999

Background

A numerical approach to solving the standard first order initial value problem:

\{ \begin{matrix}y′ = f(t, y) \\ y(t_{0}) = α \end{matrix}

produces as output a table of (approximate) function values y_{j} for the solution function of the IVP, for instance at a collection of equally-spaced t-values t_{1} = t_{0} + h, t_{2} = t_{0} + 2h:

y_{j} ≐ y(t_{j}) = y(t_{0} + jh),

where y_{j} are the approximate values, and y(t_{j}) are the exact values of the solution function at t = t_{j} = t_{0} + jh, j = 1, …, n.

In class last Thursday, we recalled a very basic method of this type called Euler’s Method: Starting from y_{0} = α, t_{0}, we construct a sequence y_{k} (one step at a time) by the rule

(1)  y_{j} = y_{j−1} + f(t_{j−1}, y_{j−1})h

Today we will consider some of the properties of this method, and related methods.

Discussion Questions

A) Consider the Exponential Growth/Decay equation y′ = ky, and the initial condition y(0) = α.

1) With step size fixed at h, what do you get after n steps of Euler’s Method? (what is y_{n}?)

2) Next, consider a fixed t̃ and let h = t̃/n (so that for all n ≥ 1, n steps take you to the same final value t_{n} = t̃). What happens as you let n → ∞ (so h → 0). Does

lim_{n→∞} y_{n}

agree with the value of the exact solution of the IVP at t = t̃?

B) What is the error after one step of Euler’s method on the general equation y′ = f(t, y)? That is, what is |y(t_{1}) − y_{1}|? Can you give a O(h^{ℓ}) estimate for the one-step error? (Hint: think Taylor!)

C) Euler’s method is hardly the “last word” in numerical methods for ODE’s. In fact it is a rather poor (i.e. inaccurate) method. Thinking in terms of the slope field of the equation y′ = f(t, y), we can see one reason immediately. Note that f(t, y) will usually be changing continuously as we move from t = t_{j−1} to t = t_{j}. Hence using the constant slope f(t_{j−1}, y_{j−1}) everywhere between t = t_{j−1} and t = t_{j} may not do a very good job of “going with the flow” of the slope field. Here is an idea for a first way to try to improve the situation:

• First we can find

y_{j}^{(Euler)} = y_{j−1} + f(t_{j−1}, y_{j−1})h

as in Euler’s Method.

• Next, we can “sample” the slope field again at the point (t_{j}, y_{j}^{(Euler)})

• If h is relatively small, then there is a good chance f is either strictly increasing or strictly decreasing along the exact solution curve between t = t_{j−1} and t = t_{j}. So it is at least plausible to suppose that the average rate of change of y between t = t_{j−1} and t = t_{j} is somewhere between f(t_{j−1}, y_{j−1}) and f(t_{j}, y_{j}^{(Euler)}).

Use the idea above to devise a modification to Euler’s method, using the two slope field values

f(t_{j−1}, y_{j−1}), f(t_{j}, y_{j}^{(Euler)})

that should be more accurate than Euler’s method itself. Give an iteration formula for your method like (1) above.

D) One test for whether your method is really better than Euler’s method is to find a O(h^{ℓ}) estimate for its one-step error. Do this for your method. Is your method really superior? Why? In working on this problem, the following facts may be useful:

• Taylor expansions for functions of two variables:

\begin{aligned}f(t_{0} + h, y_{0} + k) & = f(t_{0}, y_{0}) + \frac{∂f}{∂t} (t_{0}, y_{0}) · h + \frac{∂f}{∂y} (t_{0}, y_{0}) · k+ \\ & \frac{1}{2} (\frac{∂^{2}f}{∂t^{2}} (t_{0}, y_{0}) · h^{2} + 2\frac{∂^{2}f}{∂t∂y} (t_{0}, y_{0}) · hk + \frac{∂^{2}f}{∂y^{2}} (t_{0}, y_{0}) · k^{2}) + ⋯ \end{aligned}

(these are the first and second degree terms – the expansion keeps going, but the higher terms are not relevant here).

• If y(t) is a solution of the differential equation y′ = f(t, y), then y′′(t) can be computed by differentiating both sides:

y′′(t) = \frac{d}{dt}y′(t) = \frac{d}{dt}f(t, y(t)).

• Chain rule: If f(x, y) is a function of two variables, and a vector function of one variable (u(t), v(t)) is substituted into f, then the derivative of the resulting function of t is:

\frac{d}{dt}f(u(t), v(t)) = \frac{∂f}{∂x} (u(t), v(t)) · u′(t) + \frac{∂f}{∂y} (u(t), v(t)) · v′(t)

Assignment

Group writeups due Tuesday, February 9.