Stability and convergence for Heun's method

Question

Let \begin{align} y'(t) &= -150y(t)+49-150t, ~~~t\in[0,1]\\ y(0) &= 1/3+0.1 \end{align} I've know the solution: $y(t) = 0.1·e^{-150t}-t+1/3$.

I'm testing a couple different numerical methods to determine which is the best (i.e., Euler's forward/backward, trapezoidal, Runge-Kutta, Heun). Can anyone help me determine the stability for the Heun method?

Heun's method: $\begin{align}y_{n+1}&=y_n+h\phi(t_n,y_n,h)\\[0.3em]\phi(t,u,h) &= \tfrac{1}{2}[f(t,u)+f(t+h,u+hf(t,u))]\end{align}$

Answer 1

Stability of explicit Runge-Kutta methods around a fixed point or a slow moving equilibrium is obtained when $Lh\sim 1$ with $L$ the $y$-Lipschitz constant and $h$ the step size. The exact constant on the right side depends on the order of the method, usually it is between 1 and 5. For first and second order methods that is only a guideline as the stability region does not contain a semi-circle around the origin in the negative complex half-plane.

---

For the Heun method and any other second-order two-stage explicit Runge-Kutta method, the stability region is given by $$ |0.5z^2+z+1|\le1, $$ with $z=\lambda h$ and $\lambda$ ranging over all the eigenvalues of the Jacobian near the equilibrium.

In the given equation the one eigenvalue is real, and the inequality $-1\le 0.5x^2+x+1\le 1$ is satisfied over $[-2,0]$. Thus qualitatively correct results, bounded with a trend to falling, can be expected for $150h<2$. For solutions that fall nearly as fast as the exact solution one needs a substantially smaller step size $h$. (Which is not really visible in the experiment as the true value is practically zero, 2e-23, while the floating-point accuracy gives zero at level 1e-16. One would have to try at $t=0.1$ or $t=0.2$.)

---

To get a practical feeling for this, just solve the given task for various step sizes over the integration interval from $0$ to $1/3$, and compare how the values fall to zero, as $y(1/3)=0.1e^{-50}$ is a really small number, relative to the initial value and floating point precision.

def Heun(f,t0,y0,tf,h):
    def HeunStep(t,y,h): 
        k1 = h*f(t,y); 
        k2 = h*f(t+h,y+k1); 
        return y+0.5*(k1+k2);
    t, y = t0, y0
    while t+1.1*h<tf: t, y = t+h, HeunStep(t,y,h)
    return HeunStep(t,y,tf-t)

def f(t,y): return 49-150*(t+y)

t0,y0 = 0.0,1./3+0.1
for h in np.linspace(3.5,0.5,16)/150:
    print(f"h={h:10.4e} 150*h={150*h:10.4e}, y(1/3)={Heun(f,t0,y0,1./3,h):18.12g}")

with the results

    h=2.3333e-02 150*h=3.5000e+00, y(1/3)=3382845.45109
    h=2.2000e-02 150*h=3.3000e+00, y(1/3)=1820496.14682
    h=2.0667e-02 150*h=3.1000e+00, y(1/3)=558708.429462
    h=1.9333e-02 150*h=2.9000e+00, y(1/3)=79763.1519923
    h=1.8000e-02 150*h=2.7000e+00, y(1/3)=9204.09082254
    h=1.6667e-02 150*h=2.5000e+00, y(1/3)=1648.41784102
    h=1.5333e-02 150*h=2.3000e+00, y(1/3)=37.6118108918
    h=1.4000e-02 150*h=2.1000e+00, y(1/3)=0.740437945724
    h=1.2667e-02 150*h=1.9000e+00, y(1/3)=0.00432803787737
    h=1.1333e-02 150*h=1.7000e+00, y(1/3)=1.0688786734e-05
    h=1.0000e-02 150*h=1.5000e+00, y(1/3)=1.14794370277e-08
    h=8.6667e-03 150*h=1.3000e+00, y(1/3)=5.57823318786e-12
    h=7.3333e-03 150*h=1.1000e+00, y(1/3)=2.8015784137e-15
    h=6.0000e-03 150*h=9.0000e-01, y(1/3)=2.34187669257e-17
    h=4.6667e-03 150*h=7.0000e-01, y(1/3)=-1.38777878078e-17
    h=3.3333e-03 150*h=5.0000e-01, y(1/3)=2.60208521397e-17

As the stability region suggests, $0<h<\frac2{150}$ are valid choices in the table above, giving results close to zero. Sensible results that are practially zero as floating point numbers (and relative to the initial value) are only obtained for about $h<\frac{1.5}{150}$.