Numeric Analysis 3 - ODE¶

In [1]:

from math import *
import matplotlib.pyplot as plt
%matplotlib inline 

First-Order ODEs¶

Initial value problem

$$y'=f(x,y), \qquad y(x_0) = y_0$$

at the equidistant points on the x-axis

$$x_1 = x_0 + h, \quad x_2 = x_0 + 2h, \quad x_3 = x_0 + 3h, \quad \cdots$$

where the step size $h$ is a fixed number, for instance, 0.2 or 0.1 or 0.01

We compute approximate numeric values of $y(x)$

$$y_1, \quad y_2, \quad y_3, \quad \cdots$$

By the Taylor series

$$y(x+h)=y(x)+hy'(x)+ \frac{h^2} 2 y''(x) + \cdots$$

For small $h$ the higher powers $h^2, h^3, \cdots$ are very small.

$$y(x+h) \approx y(x)+hy'(x) = y(x)+hf(x,y)$$

The Euler method is

$$y_{n+1} = y_n + hf(x_n,y_n) \qquad (n=0,1, \cdots )$$

[ Example ] Obtain numerical solution with $h = 0.2$

$$y' = y + x, \qquad y(0)=0$$

The exact solution is $y = e^x -x -1$

---

In [2]:

from math import *

def f(x,y):   # y'
    return y + x

h = 0.2
N = 6
x = [0]*N
y = [0]*N   # initial value y[0] = 0

# print(x,y)

for n in range(N-1):     # n = 0, 1, 2, ... 4
    x[n] = n*h
    x[n+1] = (n+1)*h
    y[n+1] = y[n] + h * f( x[n], y[n] )


print(x,y)

yEuler = y     # 나중에 비교를 위해서 보관

[0.0, 0.2, 0.4, 0.6000000000000001, 0.8, 1.0] [0, 0.0, 0.04000000000000001, 0.12800000000000003, 0.27360000000000007, 0.4883200000000001]

Improved Euler Method. Predictor, Corrector¶

The predictor is calculated by

$$y_{n+1} ^* = y_n + hf(x_n,y_n)$$

The corrector is then calculated by

$$y_{n+1} = y_n + \frac 1 2 h \left[ { f(x_n,y_n) + f(x_{n+1},y_{n+1}^* ) } \right]$$

[ Example ] Obtain numerical solution with $h = 0.2$

$$y' = y + x, \qquad y(0)=0$$

The exact solution is $y = e^x -x -1$

In [3]:

from math import *

def f(x,y):   # y'
    return y + x

h = 0.2
N = 6
x = [0]*N
y = [0]*N   # initial value y[0] = 0
ypre = [0]*N

for n in range(N-1):
    x[n] = n*h
    ypre[n+1] = y[n] + h * f( x[n], y[n] )  # predictor
    x[n+1] = (n+1)*h
    y[n+1] = y[n] + h/2 *( f(x[n],y[n]) + f(x[n+1],ypre[n+1]) ) # corrector

    
# 비교
print('수치해를 정답과 비교 \n')

print('     x       Euler    Improved Euler   Exact')
for n in range(N):
    yexact = exp(x[n]) - x[n] - 1
    print('%9.6f   %9.6f   %9.6f   %9.6f' %( x[n], yEuler[n], y[n], yexact ) )

    
yImEuler = y

수치해를 정답과 비교 

     x       Euler    Improved Euler   Exact
 0.000000    0.000000    0.000000    0.000000
 0.200000    0.000000    0.020000    0.021403
 0.400000    0.040000    0.088400    0.091825
 0.600000    0.128000    0.215848    0.222119
 0.800000    0.273600    0.415335    0.425541
 1.000000    0.488320    0.702708    0.718282

Runge–Kutta Method¶

In [4]:

# Runge-Kutta method

from math import *

def f(x,y):   # y'
    return y + x

h = 0.2
N = 6
x = [0]*N
y = [0]*N   # initial value y[0] = 0

for n in range(N-1):
    x[n] = n*h
    k1 = h* f( x[n], y[n] )
    k2 = h* f( x[n]+h/2, y[n]+k1/2 )
    k3 = h* f( x[n]+h/2, y[n]+k2/2 )
    k4 = h* f( x[n]+h, y[n]+k3 )
    x[n+1] = (n+1)*h
    y[n+1] = y[n] + 1/6 *( k1 + 2*k2 + 2*k3 + k4 )

    
# 비교
print('수치해를 정답과 비교 \n')

print('     x   Improved Euler  Runge-Kutta   Exact')
for n in range(N):
    yexact = exp(x[n]) - x[n] - 1
    print('%9.6f   %9.6f   %9.6f   %9.6f' %( x[n], yImEuler[n], y[n], yexact ) )

수치해를 정답과 비교 

     x   Improved Euler  Runge-Kutta   Exact
 0.000000    0.000000    0.000000    0.000000
 0.200000    0.020000    0.021400    0.021403
 0.400000    0.088400    0.091818    0.091825
 0.600000    0.215848    0.222106    0.222119
 0.800000    0.415335    0.425521    0.425541
 1.000000    0.702708    0.718251    0.718282

Methods for Systems and Higher Order ODEs¶

Initial value problems for first-order systems of ODEs

$$y'_1 = f_1 (x,y_1,y_2,\cdots,y_m), \qquad y_1(x_0) = y_{10}$$ $$y'_2 = f_2 (x,y_1,y_2,\cdots,y_m), \qquad y_2(x_0) = y_{20}$$ $$\vdots \qquad \qquad \qquad \qquad \vdots$$ $$y'_m = f_m (x,y_1,y_2,\cdots,y_m), \qquad y_m(x_0) = y_{m0}$$

Initial value problem for an m th-order ODE

$$y^{(m)} = f \left( {x,y,y',y'', \dots ,y^{(m-1)}} \right)$$

with initial conditions $y(x_0)=K_1, y'(x_0)=K_2, \cdots , y^{(m-1)}(x_0)=K_m$

The higher-order ODE can be rewritten as a system of the first-order ODEs by setting

$$y_1=y, \quad y_2=y', \quad y_3=y'', \,\, \cdots , \,\, y_m = y^{(m-1)}$$

Then we obtain the system of first-order ODEs

$$y'_1 = y_2$$ $$y'_2 = y_3$$ $$\vdots$$ $$y'_{m-1} = y_m$$ $$y'_m = f(x,y_1,y_2,y_3,\cdots,y_m)$$

with initial conditions $y_1(x_0)=K_1, y_2(x_0)=K_2, \cdots , y_m(x_0)=K_m$

Euler Method for Systems¶

This method will not be accurate enough for practical purposes.

For first-order systems of ODEs, the Euler method prescribes

$$\mathbf y_{n+1} = \mathbf y_n + h \mathbf f( x_n, \mathbf y_n)$$

in components $$y_{1,n+1} = y_{1,n} + h f_1 (x,y_{1,n},y_{2,n},\cdots,y_{m,n}),$$ $$y_{2,n+1} = y_{2,n} + h f_2 (x,y_{1,n},y_{2,n},\cdots,y_{m,n}),$$ $$\vdots \qquad \qquad \qquad \qquad \vdots$$ $$y_{m,n+1} = y_{m,n} + h f_m (x,y_{1,n},y_{2,n},\cdots,y_{m,n}),$$

[Example] Second-Order ODE. Mass–Spring System

$$y'' +2y' + 0.75y =0, \quad y(0)=3, \quad y'(0)=-2.5$$

Obtain the numerical solution with step $h=0.2$

(exact solution) $y=2 e^{-0.5x} + e^{-1.5x}$

---

We let $$y_1=y, \quad y_2=y'$$

The the second-order ODE becomes

$$y'_1 = y_2, \qquad \qquad \,\,\,\,\,\, y_1(0) = 3$$ $$y'_2 = -2y_2 - 0.75y_1, \qquad y_2(0) = -2.5$$

In [5]:

# Euler method

from math import *

def f1(x,y1,y2):            # y1' = f1(x,y1,y2)
    return y2

def f2(x,y1,y2):            # y2' = f2(x,y1,y2)
    return -2*y2 - 0.75*y1

h = 0.2
N = 6
x =  [0]*N
y1 = [0]*N   
y2 = [0]*N

# initial value 
y1[0] = 3
y2[0] = -2.5

for n in range(N-1):     # n = 0,1,2, ... 4
    x[n] = n*h
    x[n+1] = (n+1)*h
    y1[n+1] = y1[n] + h* f1( x[n], y1[n], y2[n] )
    y2[n+1] = y2[n] + h* f2( x[n], y1[n], y2[n] )

   
# 비교
print('수치해를 정답과 비교 \n')
print('     x       Euler       Exact')

for n in range(N):
    yexact = 2*exp(-0.5*x[n]) + exp(-1.5*x[n])
    print('%9.6f   %9.6f   %9.6f' %( x[n], y1[n], yexact ) )


yEuler = y1  # 나중을 위해 저장

수치해를 정답과 비교 

     x       Euler       Exact
 0.000000    3.000000    3.000000
 0.200000    2.500000    2.550493
 0.400000    2.110000    2.186273
 0.600000    1.801000    1.888206
 0.800000    1.552300    1.641834
 1.000000    1.349050    1.436191

Runge–Kutta Methods for Systems¶

At each step n, four auxiliary quantities are calculated

$$\begin{align} \mathbf k_1 &= h \mathbf f(x_n, \mathbf y_n) \\ \mathbf k_2 &= h \mathbf f(x_n + \frac 1 2 h, \mathbf y_n + \frac 1 2 \mathbf k_1 ) \\ \mathbf k_3 &= h \mathbf f(x_n + \frac 1 2 h, \mathbf y_n + \frac 1 2 \mathbf k_2 ) \\ \mathbf k_4 &= h \mathbf f(x_n + h, \mathbf y_n + \mathbf k_3 ) \end{align}$$

New values at the next step n+1 are calculated by

$$\mathbf y_{n+1} = \mathbf y_n + \frac 1 6 (\mathbf k_1 + 2 \mathbf k_2 + 2 \mathbf k_3 + \mathbf k_4 )$$

[Example] Second-Order ODE. Mass–Spring System

$$y'' +2y' + 0.75y =0, \quad y(0)=3, \quad y'(0)=-2.5$$

Obtain the numerical solution with step $h=0.2$

(exact solution) $y=2 e^{-0.5x} + e^{-1.5x}$

---

In [6]:

# Runge-Kutta method

from math import *

def f1(x,y1,y2):            # y1' = f1(x,y1,y2)
    return y2

def f2(x,y1,y2):            # y2' = f2(x,y1,y2)
    return -2*y2 - 0.75*y1

h = 0.2
N = 6
x =  [0]*N
y1 = [0]*N   
y2 = [0]*N

# initial value 
y1[0] = 3
y2[0] = -2.5

for n in range(N-1):     # n = 0,1,2, ... 4
    x[n] = n*h
    x[n+1] = (n+1)*h

    k11 = h* f1( x[n], y1[n], y2[n] )
    k12 = h* f2( x[n], y1[n], y2[n] )
    
    k21 = h* f1( x[n]+h/2, y1[n]+k11/2, y2[n]+k12/2 )
    k22 = h* f2( x[n]+h/2, y1[n]+k11/2, y2[n]+k12/2 )
    
    k31 = h* f1( x[n]+h/2, y1[n]+k21/2, y2[n]+k22/2 )
    k32 = h* f2( x[n]+h/2, y1[n]+k21/2, y2[n]+k22/2 )
    
    k41 = h* f1( x[n]+h, y1[n]+k31, y2[n]+k32 )
    k42 = h* f2( x[n]+h, y1[n]+k31, y2[n]+k32 )
    
    y1[n+1] = y1[n] + 1/6 *( k11 + 2*k21 + 2*k31 + k41 )
    y2[n+1] = y2[n] + 1/6 *( k12 + 2*k22 + 2*k32 + k42 )
    
    
# 비교
print('수치해를 정답과 비교 \n')
print('     x        Euler     Runge-Kutta   Exact')

for n in range(N):
    yexact = 2*exp(-0.5*x[n]) + exp(-1.5*x[n])
    print('%9.6f   %9.6f   %9.6f   %9.6f' %( x[n], yEuler[n], y1[n], yexact ) )

수치해를 정답과 비교 

     x        Euler     Runge-Kutta   Exact
 0.000000    3.000000    3.000000    3.000000
 0.200000    2.500000    2.550512    2.550493
 0.400000    2.110000    2.186302    2.186273
 0.600000    1.801000    1.888238    1.888206
 0.800000    1.552300    1.641866    1.641834
 1.000000    1.349050    1.436221    1.436191