Let $$ \begin{matrix} y'(t) = -150y(t)+49-150t, t\in[0,1]\\ y(0) = 1/3+0.1 \end{matrix} $$ 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 (ie. Euler's forward/backward, trapezoidal, Runge-Katta, Heun). Can anyone help me determine the stability for the Heun method?

Heun's method: $\phi(t,u,h) := \frac{1}{2}[f(t,u)+f(t+h,u+hf(t,u))]$


To get a practical feeling for this, just solve your problem for various step sizes from $0$ to $1/3$, 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 "h=%.4e 150*h=%.4e, y(1/3)=%.12g"%(h,150*h, Heun(f,t0,y0,1./3,h))

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

While the stability region would suggest $0<h<\frac2{150}$ are valid choices, the table above gives only sensible results for about $h<\frac{1.5}{150}$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.