# Order of accuracy of Numerov's method

I used Numerov's method to solve the following ODE $$\ddot{y}+\omega^2(t)y=0$$ where $$\omega (t)$$ is a $$\tanh$$ function. However, when I looked at the order of accuracy of the method (I let $$\omega(t)$$ be a constant) by halving the step size. The error only halved instead of reduced by a factor of $$2^4$$ as it supposed to be since Numerov's method is a fourth order method. I also tried letting $$\omega(t)$$ be the $$\tanh$$ function I used then calculate the order of accuracy using $$p=\log_2\left(\frac{y(h)-y(0.5h)}{y(0.5h)-y(0.25h)}\right)$$ where $$h$$ is the step size, it gives $$p=1$$. What is the reason for this disagreement?

• Numerov's method is a multistep method. How do you compute $x_1$ from $x_0$? Commented Mar 28, 2019 at 21:51
• See math.stackexchange.com/q/3064879/115115, where the 4th order of the method is confirmed via error plots. The last somewhat demonstrates the sensitivity towards the $x_1$ computation. Commented Mar 28, 2019 at 22:10
• I calculated the $y_1$ by $y(\Delta t) = y(0) + y'(0) - \frac{(\Delta t)^2}{2} \omega^2(t)y(0)$. Commented Mar 29, 2019 at 6:36
• This is $O(Δt^3)$, due to the double summation from the second order difference quotient the error initially behaves like a quadratic function, thus it can be expected that it contains a term $t^2Δt$ aside some term $tΔt^2$, giving a first order global error. Commented Mar 29, 2019 at 7:14

For $$\ddot y(t)+\tanh^2(t)y(t)=0$$ one finds that the solution close to zero has a power series expansion $$y(t)=\sum c_kt^k$$ that starts with $$[2c_2+6c_3t+12c_4t^2+20c_5t^3+...]+[t^2+\tfrac23t^4+...][c_0+a_1t+c_2t^2+...]=0\\ c_2=0\\c_3=0\\12c_4=-c_0\\20c_5=-c_1\\...\\ \implies y(t)=\bigl[1-\tfrac12t^4\pm...\bigr]y_0+\bigl[1-\tfrac1{20}t^4\pm...\bigr]t\dot y_0.$$ This means that the error in setting $$y(Δt)\approx y_1=y(0)+\dot y(0)Δt+\frac12\ddot y(0)Δt^2$$ is not the generally expected $$O(Δt^3)$$ but for this case $$O(Δt^4)$$.
The global error $$e(t_k)=y_k-y(t_k)$$ of the Numerov method has two main contributions, the general one from the discretization error and secondly the error in the first step. Their leading terms are $$e(t)=a(t)Δt^p+b(t)Δt^4$$ where in the case $$p<4$$ the functions $$a,b$$ satisfy the differential equations \begin{align} \ddot a(t)+\tanh^2(t)a(t)&=0,&~~a(0)&=0,~~a(Δt)Δt^{p+1}=y_1-y(Δt)+O(Δt^{p+2})\\ \ddot b(t)+\tanh^2(t)b(t)&=-r_6(t),&~~b(0)&=0,~~\dot b(0)=0 \end{align} The function $$r_6(t)Δt^6+...=y(t+Δt)-2y(t)+y(t-Δt)-\frac1{12}[f(t+Δt,..)+10f(t,..)+f(t-Δt,..)]Δt^2$$ is the local truncation error of the exact solution. As $$a$$ has a power series expansion like $$y$$, one finds that $$a(Δt)Δt^3=[1-\frac1{20}Δt^4\pm\dots]Δt^4\dot a(0)=-\frac1{12}y(0)Δt^4$$ so that in leading order and for small $$t$$ one has $$a(t)=-\frac1{12}y(0)\cdot t$$ and the dominant error term is $$a(t)Δt^3$$.
In conclusion, with a correct implementation of the method you should find an empirical error order $$3$$ for the given problem with the quadratic first step. In general, if the third derivative at zero is non-zero, that method assembly should only give global error order $$2$$.
To remove any visible influence of the error in the first term one needs to compute $$y_1$$ with an order 4 method or better, as then also $$p=4$$, so that both error components have the same power.