Numeric Analysis 3 - ODE¶
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
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
For small $h$ the higher powers $h^2, h^3, \cdots$ are very small.
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 $
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 # 나중에 비교를 위해서 보관
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 $
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
Runge–Kutta Method¶
# 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 ) )
Methods for Systems and Higher Order ODEs¶
Initial value problems for first-order systems of ODEs
Initial value problem for an m th-order ODE
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
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
# 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 # 나중을 위해 저장
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}$
# 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 ) )
댓글 없음:
댓글 쓰기