Skip to article frontmatterSkip to article content

ODE with periodic solution

from pylab import *

Consider the ODE system

x′=−y,y′=xx(0)=1,y(0)=0\begin{gather} x' = -y, \qquad y' = x \\ x(0) = 1, \qquad y(0) = 0 \end{gather}
(1)

The exact solution is

x(t)=cos⁡(t),y(t)=sin⁡(t)x(t) = \cos(t), \qquad y(t) = \sin(t)
(2)

This solution is periodic. It also has a quadratic invariant

x2(t)+y2(t)=1,∀tx^2(t) + y^2(t) = 1, \qquad \forall t
(3)
theta = linspace(0.0, 2*pi, 500)
xe = cos(theta)
ye = sin(theta)

Forward Euler scheme

def ForwardEuler(h,T):
    N = int(T/h)
    x,y = zeros(N),zeros(N)
    x[0] = 1.0
    y[0] = 0.0
    for n in range(1,N):
        x[n] = x[n-1] - h*y[n-1]
        y[n] = y[n-1] + h*x[n-1]
    return x,y
h = 0.02
T = 4.0*pi
x,y = ForwardEuler(h,T)

plot(xe,ye,'r--',label='Exact')
plot(x,y,label='FE')
legend()
grid(True)
gca().set_aspect('equal')
<Figure size 640x480 with 1 Axes>

The phase space trajectory is spiralling outward.

Backward Euler scheme

xn=xn−1−hyn,yn=yn−1+hxnx_n = x_{n-1} - h y_n, \qquad y_n = y_{n-1} + h x_n
(4)

Eliminate yny_n from first equation and solve for xnx_n to get

xn=xn−1−hyn−11+h2x_n = \frac{x_{n-1} - h y_{n-1}}{1 + h^2}
(5)

and now yny_n can also be computed.

def BackwardEuler(h,T):
    N = int(T/h)
    x,y = zeros(N),zeros(N)
    x[0] = 1.0
    y[0] = 0.0
    for n in range(1,N):
        x[n] = (x[n-1] - h*y[n-1])/(1.0 + h**2)
        y[n] = y[n-1] + h*x[n]
    return x,y
h = 0.02
T = 4.0*pi
x,y = BackwardEuler(h,T)

plot(xe,ye,'r--',label='Exact')
plot(x,y,label='BE')
legend()
grid(True)
gca().set_aspect('equal')
<Figure size 640x480 with 1 Axes>

The phase space trajectory is spiralling inward.

Trapezoidal scheme

xn=xn−1−h2(yn−1+yn),yn=yn−1+h2(xn−1+xn)x_n = x_{n-1} - \frac{h}{2}(y_{n-1} + y_n), \qquad y_n = y_{n-1} + \frac{h}{2}(x_{n-1} + x_n)
(6)

Eliminate yny_n from first equation

xn=(1−14h2)xn−1−hyn−11+14h2x_n = \frac{ (1-\frac{1}{4}h^2) x_{n-1} - h y_{n-1} }{1 + \frac{1}{4}h^2}
(7)

and now yny_n can also be computed.

def Trapezoid(h,T):
    N = int(T/h)
    x,y = zeros(N),zeros(N)
    x[0] = 1.0
    y[0] = 0.0
    for n in range(1,N):
        x[n] = ((1.0-0.25*h**2)*x[n-1] - h*y[n-1])/(1.0 + 0.25*h**2)
        y[n] = y[n-1] + 0.5*h*(x[n-1] + x[n])
    return x,y
h = 0.02
T = 4.0*pi
x,y = Trapezoid(h,T)

plot(xe,ye,'r--',label='Exact')
plot(x,y,label='Trap')
legend()
grid(True)
gca().set_aspect('equal')
<Figure size 640x480 with 1 Axes>

The phase space trajectory is exactly the unit circle.

Multiply first equation by xn+xn−1x_n + x_{n-1} and second equation by yn+yn−1y_n + y_{n-1}

(xn+xn−1)(xn−xn−1)=−h2(xn+xn−1)(yn+yn−1)(x_n + x_{n-1})(x_n - x_{n-1}) = - \frac{h}{2}(x_n + x_{n-1})(y_n + y_{n-1})
(8)
(yn+yn−1)(yn−yn−1)=+h2(xn+xn−1)(yn+yn−1)(y_n + y_{n-1})(y_n - y_{n-1}) = + \frac{h}{2}(x_n + x_{n-1})(y_n + y_{n-1})
(9)

Adding the two equations we get

xn2+yn2=xn−12+yn−12x_n^2 + y_n^2 = x_{n-1}^2 + y_{n-1}^2
(10)

Thus the Trapezoidal method is able to preserve the invariant.

Partitioned Euler

This is a symplectic Runge-Kutta method. xx is updated explicitly and yy is updated implicitly.

xn=xn−1−hyn−1,yn=yn−1+hxnx_n = x_{n-1} - h y_{n-1}, \qquad y_n = y_{n-1} + h x_{n}
(11)

The update of yy uses the latest updated value of xx.

def PartEuler(h,T):
    N = int(T/h)
    x,y = zeros(N),zeros(N)
    x[0] = 1.0
    y[0] = 0.0
    for n in range(1,N):
        x[n] = x[n-1] - h * y[n-1]
        y[n] = y[n-1] + h * x[n]
    return x,y
h = 0.02
T = 4.0*pi
x,y = PartEuler(h,T)

plot(xe,ye,'r--',label='Exact')
plot(x,y,label='PE')
legend()
grid(True)
gca().set_aspect('equal')
<Figure size 640x480 with 1 Axes>

We get very good results, even though the method is only first order. We can also switch the implicit/explicit parts

yn=yn−1+hxn−1,xn=xn−1−hyny_n = y_{n-1} + h x_{n-1}, \qquad x_n = x_{n-1} - h y_n
(12)

Classical RK4 scheme

def rhs(u):
    return array([u[1], -u[0]])

def RK4(h,T):
    N = int(T/h)
    u = zeros((N,2))
    u[0,0] = 1.0 # x
    u[0,1] = 0.0 # y
    for n in range(0,N-1):
        w = u[n,:]
        k0 = rhs(w)
        k1 = rhs(w + 0.5*h*k0)
        k2 = rhs(w + 0.5*h*k1)
        k3 = rhs(w + h*k2)
        u[n+1,:] = w + (h/6)*(k0 + k1 + k2 + k3)
    return u[:,0], u[:,1]

We test this method for a longer time.

h = 0.02
T = 10.0*pi
x,y = RK4(h,T)

plot(xe,ye,'r--',label='Exact')
plot(x,y,label='RK4')
legend()
grid(True)
gca().set_aspect('equal')
<Figure size 640x480 with 1 Axes>

RK4 is a more accurate method compared to forward/backward Euler schemes, but it still loses total energy. The trajectory spirals inwards towards the origin.

Numerical Analysis
ODE System
Numerical Analysis
Shooting method for BVP