# Convergence-test for ODE approximates wrong limit

I am trying to numerically solve a differential equation but I am having trouble getting the convergence test to run properly. The problem is as follows:

Consider an ODE

$$y'(t) \enspace = \enspace f(t,y) \quad , \qquad y(0) \enspace = \enspace y_0 \quad . \tag{1}$$

Let now $$y^{h}$$ be a numerical solution to $$(1)$$ with step size $$h$$. My numerical solver is of order $$\mathcal{O}(h^n)$$. As far as I understand, I can test convergence of the solution by verifying the limit

$$\lim_{h \longrightarrow 0} \frac{ || y^{h/2} - y^{h} || }{ || y^{h/4} - y^{h/2} || } \enspace = \enspace 2^n \quad .$$

In order to keep the code for this question nice and simple and clean, I will use a standard Euler-Forward approach of order $$\mathcal{O}(h^2)$$ and do not care about efficiency etc. of the implementation itself. My working example will be the differential equation

$$y'(t) \enspace = \enspace - \sin(t) \quad , \qquad y(0) = 1$$

with known exact solution

$$y(t) \enspace = \enspace \cos(t)$$

Here is the broken-down and simplified code in PYTHON (Note that I am using grid point number $$N$$ instead of step size $$h$$):

import numpy as np
import matplotlib.pyplot as plt

def Solve( N ) :

y = np.zeros(N)
y[0] = 1

t = np.linspace( 0, 2*np.pi, num=N )
h = t[1] - t[0]

def func( x ) :
return -np.sin(x)

for i in range(1,len(t)) :
y[i] = y[i-1] - h*np.sin(t[i-1])

return y

def main() :
N = 10
for i in range( 10 ):
a = Solve( N*2**i )
b = Solve( N*2**(i+1) )
c = Solve( N*2**(i+2) )

res = np.linalg.norm( b[::2] - a ) / np.linalg.norm( c[::2] - b )

plt.scatter( N*2**i, res )

main()


This should by what I stated above converge to $$2^2 = 4$$. However, it rather converges to some different value (probably $$\sqrt{2}$$, see picture) and I have no clue why. Any ideas?

• To get step size $$h=1/N$$, you need $$N$$ steps, which gives $$N+1$$ nodes in the time subdivision. Currently you compare sequences with step sizes $$\frac1{N-1}$$, $$\frac1{2N-1}$$ and $$\frac1{4N-1}$$. These are not exact multiples, and so in the difference of the sequences you compare values at slightly different points.
• You have indeed a built-in $$\sqrt2$$, as in the denominator the vector for the norm has twice as many as the one in the numerator. So the limit is indeed $$\frac{2^1}{\sqrt2}$$ The easiest way to repair this is to use the same points below as above
               res = np.linalg.norm( b[::2] - a ) / np.linalg.norm( c[::4] - b[::2] )