Numerical solution of differential equationsEdexcel A-Level Further Maths: Revision notes

Many differential equations cannot be solved exactly. A numerical method starts from a known point $(x_0,y_0)$ and moves forward in equal steps of length $h$, so $x_n=x_0+nh$. At each step the derivative at $x_n$ is replaced by a difference approximation that uses nearby values $y_{n-1}$, $y_n$, $y_{n+1}$. The notation $y_n$ means the approximation to $y(x_n)$. Smaller steps usually give more accurate results but need more calculation. The three approximations on the specification are $\left(\frac{dy}{dx}\right)_n=\frac{y_{n+1}-y_n}{h},\quad\left(\frac{dy}{dx}\right)_n=\frac{y_{n+1}-y_{n-1}}{2h},\quad\left(\frac{d^2y}{dx^2}\right)_n=\frac{y_{n+1}-2y_n+y_{n-1}}{h^2}.$

Substitute $\left(\frac{dy}{dx}\right)_n=\frac{y_{n+1}-y_n}{h}$ into $\frac{dy}{dx}=f(x,y)$ and rearrange: $y_{n+1}=y_n+h\,f(x_n,y_n).$ Example: $\frac{dy}{dx}=x+y$, $y(0)=1$, $h=0.2$. Then $y_1=1+0.2(0+1)=1.2$, $y_2=1.2+0.2(0.2+1.2)=1.48$ and $y_3=1.48+0.2(0.4+1.48)=1.856$. The exact solution is $y=2e^x-x-1$, so $y(0.6)=2.044$ and the percentage error in $y_3$ is $\frac{2.044-1.856}{2.044}\times100=9.2\%$. Each step uses the gradient at the start of the interval, so errors build up.

The approximation $\left(\frac{dy}{dx}\right)_n=\frac{y_{n+1}-y_{n-1}}{2h}$ is more accurate, because it is centred on $x_n$. Substituting into $\frac{dy}{dx}=f(x,y)$ gives $y_{n+1}=y_{n-1}+2h\,f(x_n,y_n).$ It needs two starting values, $y_0$ and $y_1$, because every step looks back two values. Find $y_1$ by the forward difference method. Example: $\frac{dy}{dx}=x-y$, $y(0)=1$, $h=0.1$. First $y_1=1+0.1(0-1)=0.9$. Then $y_2=y_0+0.2(x_1-y_1)=1+0.2(0.1-0.9)=0.84$ and $y_3=y_1+0.2(x_2-y_2)=0.9+0.2(0.2-0.84)=0.772$.

Replace $\frac{d^2y}{dx^2}$ by $\frac{y_{n+1}-2y_n+y_{n-1}}{h^2}$. For $\frac{d^2y}{dx^2}=f(x,y)$ this gives $y_{n+1}=2y_n-y_{n-1}+h^2f(x_n,y_n).$ If the equation also contains $\frac{dy}{dx}$, replace it by the central difference $\frac{y_{n+1}-y_{n-1}}{2h}$ as well, then collect the $y_{n+1}$, $y_n$ and $y_{n-1}$ terms. Example: $\frac{d^2y}{dx^2}=-y$ with $h=0.1$ gives $y_{n+1}-2y_n+y_{n-1}=-0.01y_n$, so $y_{n+1}=1.99y_n-y_{n-1}$.

The second-order formula also needs two values, $y_0$ and $y_1$, but the initial conditions give $y_0$ and $\frac{dy}{dx}$ at $x_0$. Use the central difference at $n=0$: $\frac{y_1-y_{-1}}{2h}=\left(\frac{dy}{dx}\right)_0,$ which brings in the extra value $y_{-1}$. Write the recurrence formula with $n=0$, and solve the two equations for $y_1$ and $y_{-1}$. Example: $y_0=0$, $\frac{dy}{dx}=1$, $h=0.1$, and $y_{n+1}=1.99y_n-y_{n-1}$. Then $y_1-y_{-1}=0.2$ and $y_1=-y_{-1}$, so $y_1=0.1$, and $y_2=1.99(0.1)-0=0.199$. (The exact value is $\sin0.2=0.1987$.) A simpler start is $y_1=y_0+h\left(\frac{dy}{dx}\right)_0$.

Compare a numerical estimate with an exact value using the percentage error: $\frac{\text{exact}-\text{estimate}}{\text{exact}}\times100.$ Errors grow as more steps are taken. To improve an estimate, use a smaller step length $h$, so that each difference approximation is closer to the true derivative, or use the central difference rather than the forward difference. For $\frac{d^2y}{dx^2}+2\frac{dy}{dx}+y=0$, $y(0)=1$, $y'(0)=0$ and $h=0.5$, the estimate of $y(1)$ is $0.6875$ against an exact $0.7358$: an error of $6.6\%$.