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 )


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?

enter image description here


1 Answer 1

  • The Euler method has global error order 1, not 2

  • 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] )
  • $\begingroup$ Thank you for your answer! I accepted it gladly, because point 1 and point 3 were very enlightening. However, I dont quite understand point 2. I am taking the limit h -> 0, so this discrepancy shouldnt have any influence on the outcome of the convergence test. Was this just intended as a general comment? $\endgroup$
    – Octavius
    Nov 20 at 15:59
  • $\begingroup$ That is correct, the influence of this error goes to zero. It is just a question of style to avoid such systematic errors. $\endgroup$ Nov 20 at 16:16
  • 1
    $\begingroup$ Also, the time scales do not match. This is without consequence for order 1 methods, but can degrade the computation to order 1 also for higher order methods. $\endgroup$ Nov 20 at 16:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.