[Math] Use Heun’s method without iteration to solve a second order differential equation, Chapra Numerical methods 7th edition chapter 25 p.3

numerical methods

I've spent a day trying to solve this problem and I can't figure out how:

Use the (b) Heun (without iteration) method to
solve
$$\dfrac{d^2y}{dt^2}-0.5t+y=0$$
where $y(0)=2$ and $y'(0)=0$. Solve from x=0 to 4 using h=0.1.

Someone who really knows how to think please just show the first two or three steps. If you don't see the circular dilemma, or think I'm referring to the order of the error please don't respond. The answer I've posted below. You are starting from the second derivative and have to use $y''$ slope to find $y'$ which you then use to find $y^j_{i+1}$ which would have been used to plug into $y'$ to average the slopes to find $y_{i+1}$ but we started from $y''$. For a reminder a link to heun's rule I linked below too. I assume anyone who can answer this and get the right numbers wouldn't need a reminder for the easy first derivative version. Thankfully there's NO iteration. Just average the slopes once.

https://drive.google.com/file/d/0B0TrO0YEHpxxWDFaRkJoZWFyclk/view?usp=sharing

https://en.wikipedia.org/wiki/Heun%27s_method

Best Answer

To break the circle, you do exactly as indicated in the partial solution, you transform into a system of order 1. Set $z=y'$, then $$ \frac d{dt}\binom yz =\binom{z}{0.5·t-y} $$ This now has the required form $u'=f(t,u)$ and can be used in a generic Heun solver

import numpy as np

def ODEfunc(t,u):
    y,z = u
    return np.array([ z, 0.5*t-y ])



def HeunStep(f,t,u,h):
    k1 = h*f(t, u)
    k2 = h*f(t+h, u+k1)
    return u+0.5*(k1+k2)

def HeunTest(f,a,b,x0,h):
    t = a
    x = x0.copy()
    while t < b+0.5*h:
        print t, x
        x = HeunStep(f,t,x,h)
        t=t+h


HeunTest(ODEfunc,0,4,np.array([2.,0.]),0.1)

Also note that you can easily obtain the exact solution $$ y(t)=2\cos(t)-0.5\sin(t)+0.5t $$ by standard methods to compare against.

Related Question