ExerciseODEsadvanced

Exercise: stability of explicit and implicit Euler

Full analysis of y′=λyy'=\lambda y: the amplification factor of each method, the explicit method's stability condition h<−2/λh<-2/\lambda, the unconditional stability of the implicit one and a numerical check with λ=−10\lambda=-10.

Amplification factor

Study the model problem y′=λyy'=\lambda y with λ<0\lambda<0 (decaying Malthus model): the exact solution y(t)=yaeλty(t)=y_a e^{\lambda t} tends to zero, and a reasonable method should reproduce that behaviour.

  1. Explicit Euler applied to f(t,y)=λyf(t,y)=\lambda y multiplies the solution by a constant factor at each step:

    yk+1=yk+hλyk=(1+hλ) yk  ⇒  yk=(1+hλ)k yay_{k+1}=y_k+h\lambda y_k=(1+h\lambda)\,y_k\;\Rightarrow\; y_k=(1+h\lambda)^k\,y_a

  2. The numerical solution decays if and only if ∣1+hλ∣<1|1+h\lambda|<1, that is, if h<−2λh<-\frac{2}{\lambda}. With larger hh, (1+hλ)k(1+h\lambda)^k oscillates with growing amplitude: the method is unstable even though it is consistent.

  3. Implicit Euler has a closed-form solve in this linear problem:

    yk+1=yk+hλyk+1  ⇒  yk+1=yk1−hλy_{k+1}=y_k+h\lambda y_{k+1}\;\Rightarrow\; y_{k+1}=\frac{y_k}{1-h\lambda}

  4. Since λ<0\lambda<0, the denominator 1−hλ=1+h∣λ∣1-h\lambda=1+h|\lambda| exceeds 1 for every h>0h>0: the amplification factor stays below 1 and the method is unconditionally stable.

Numerical check

With λ=−10\lambda=-10, y(0)=1y(0)=1 on [0,2][0,2] (exact solution y=e−10ty=e^{-10t}), the explicit method's stability condition is h<0.2h<0.2, i.e. N>10N>10. The table shows the maximum error of both methods:

NENE_N explicitENE_N implicit
281.00000.0909
4256.00000.1599
825.62890.2036
160.53650.1579
320.16030.0922
Maximum error against the exact solution y=e−10ty=e^{-10t} on [0,2][0,2].

Solving the implicit step in a nonlinear case

ExampleImplicit Euler on the logistic equation

Set up the implicit Euler step for y′=y(1−y)y'=y(1-y) and solve for yk+1y_{k+1}.

  1. The implicit step is yk+1=yk+h yk+1(1−yk+1)y_{k+1}=y_k+h\,y_{k+1}(1-y_{k+1}). Rearranging yields a quadratic equation in yk+1y_{k+1}:

    h yk+12+(1−h) yk+1−yk=0h\,y_{k+1}^2+(1-h)\,y_{k+1}-y_k=0

  2. For positive solutions take the plus-sign root:

    yk+1=−(1−h)+(1−h)2+4h yk2hy_{k+1}=\frac{-(1-h)+\sqrt{(1-h)^2+4h\,y_k}}{2h}

  3. Here each step's equation has a closed-form solution; in general there is none and g(yk+1)=0g(y_{k+1})=0 is solved with the Newton-Raphson method.

The extra cost per step (solving an equation) is the price of the implicit method's unconditional stability.