Order of accuracy of Numerov's method

Question

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?

Answer 1

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.