Understanding Heun's method

Question

I want to have an intuitive understanding of Heun's method when applied to a 2nd order ODE rather than rely on plug and chug techniques.

The iteration is calculated by the formula $y_{n+1}=y_n+\frac{h}4(k_1+3k_2)$

where $k_1=f(x_n,y_n)$ ,$k_2=f(x_n+\frac{2h}{3},y_n+\frac{2hk_1}{3})$

Approximate $y(0.1)$ with step-size $h=0.1$ given $y''=x^2+y-xy'$ and $y(0)=1,y'(0)=0$

Now I can convert a second order ODE to a first order by substituting $u=y', u'=y''$, so $x_0=0,y_0=1,u_0=0$ to give $u'=x^2+y-xu$

The principle behind Heun's method is to use the average of the two slopes, with $k_1$ and $k_2$ denoting the slopes at the initial value and at the first iteration.

Now the slope of the graph is $y'$ and in this example $y'=u$ so finding the value of u at the respective points gives us the graph. We know that $u(0)=0$ therefore the slope at $(0,1)=0$ so $k_1=0$, $k_2$=$f(\frac{1}{16}$,$\frac{15}{16})$

My only issue is that it seems to be a circular problem. Before I can work out $u$ I have to know $u'$, but $u'$ is obtained based on the value of $u$.

Ralston Graph

Answer 1

You have everything together, you need just to get used to the idea to apply the method to a system of ODE, or to an ODE for a vector-valued function. $$ \begin{align} y'&=u\\ u'&=f(x,y,u)=x^2+y−xu \end{align} $$ which can be implemented as

dx = x[j+1]-x[j]

k1y = dx*u[j]
k1u = dx*f(x[j],y[j],u[j])

k2y = dx*(u[j]+2/3*k1u)
k2u = dx*f(x[j]+2/3*dx, y[j]+2/3*k1y, u[j]+2/3*k1u)

y[j+1] = y[j] + (k1y+3*k2y)/4
u[j+1] = u[j] + (k1u+3*k2u)/4