- 追加された行はこの色です。
- 削除された行はこの色です。
* 常微分方程式の求解 [#ua9e43bb]
ここでは、n 階の常微分方程式は一階の常微分方程式に変形し、u(t) = f(u,t) という形式として取り扱うことにしておく．
さて、常微分方程式の求解には、大雑把に言って
- 一段法: ある時点 t のデータ u(t) が有れば、Δt だけ経った次の時間 t+Δt の近似値 u(t+Δt) が計算できる方法
- 一段法: u(t) の近似値を計算するために、u(t-Δt) のデータだけあれば済む方法．
- 多段法: u(t) の近似値を計算するために、u(t-Δt) だけではなく、u(t-2Δt), u(t-3Δt) など、複数の過去データを必要とする方法
の二種類がある．
** 一段法 [#re808e45]
簡単だが安定性と精度に難ありの Euler 法と、頑健な Runge-Kutta 法の二つを授業では説明した．
以下、プログラムを載せておこう．
*** サンプルプログラム: Euler 法 1 [#n03818b6]
人口問題に対するかなりシンプルなモデル方程式である
du/dt = u(1-u) という常微分方程式に対して Euler 法でのプログラム．
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# Euler method program for the initial value problem (1).
include Math
##### constants #####
# calculation t from 0 to Tmax.
Tmax = 10.0
# Delta t
Dt = 1.0/10
# Initial value of u
IniU = 0.4
##### functions #####
# Right hand side of ODE
def f(u)
return u*(1.0-u)
end
# update of the approximate solution
def new_by_euler(u)
return u + Dt * f(u)
end
# result
def print_result(u,t)
printf("%.6f %.6f\n", t, u)
end
##### main program #####
# initialization
u = IniU
print_result(u, 0.0)
# main loop
for i in (1..(Tmax/Dt).ceil )
u = new_by_euler(u)
print_result(u, i*Dt)
end
}}
*** サンプルプログラム: Euler 法 2 [#vaac1276]
サメと魚のような、被食者-捕食者の関係があるときのその種の数の変動に関するシンプルなモデル常微分方程式に対して Euler 法でのプログラム．
方程式は 魚相当の種の匹数を u(t), サメ相当の種の匹数を v(t) とすると
&br;
du/dt = (2-v)u,
&br;
dv/dt = (2u-3)v, &br;
という連立一階常微分方程式である．
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# Euler method program for the initial value problem (2).
include Math
##### constants #####
# calculation t from 0 to Tmax.
Tmax = 10.0
# Delta t
Dt = 1.0/100
# Initial value of u
IniU = [4.0, 1.0]
##### functions #####
# Right hand side of ODE
def f(u)
a = u[0]
b = u[1]
return [ (2.0-b)*a, (2.0*a-3.0)*b ]
end
# update of the approximate solution
def new_by_euler(u)
v = f(u)
return [ u[0] + Dt*v[0], u[1] + Dt*v[1] ]
end
# result
def print_result(u,t)
printf("%.4f %.4f %.4f\n", t, u[0], u[1])
end
##### main program #####
# initialization
u = IniU
print_result(u, 0.0)
# main loop
for i in (1..(Tmax/Dt).ceil )
u = new_by_euler(u)
print_result(u,i*Dt)
end
}}
*** サンプルプログラム: Runge-Kutta 法 1 [#r95c9463]
同様に
du/dt = u(1-u) という常微分方程式に対して Runge-Kutta 法でのプログラム．
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# Runge--Kutta method program for the initial value problem (1).
include Math
##### constants #####
# calculation t from 0 to Tmax.
Tmax = 10.0
# Delta t
Dt = 1.0/10
# Initial value of u
IniU = 0.4
##### functions #####
# Right hand side of ODE
def f(u)
return u*(1.0-u)
end
# update of the approximate solution
def new_by_rk(u)
f1 = f(u)
f2 = f(u + 0.5*Dt*f1)
f3 = f(u + 0.5*Dt*f2)
f4 = f(u + Dt*f3)
return u + Dt*(f1/6.0 + f2/3.0 + f3/3.0 + f4/6.0)
end
# result
def print_result(u,t)
printf("%.6f %.6f\n", t, u)
end
##### main program #####
# initialization
u = IniU
print_result(u, 0.0)
# main loop
for i in (1..(Tmax/Dt).ceil )
u = new_by_rk(u)
print_result(u, i*Dt)
end
}}
*** サンプルプログラム: Runge-Kutta 法 2 [#ye83643e]
同様に
&br;
du/dt = (2-v)u,
&br;
dv/dt = (2u-3)v, &br;
という連立一階常微分方程式に対する Runge-Kutta 法でのプログラム．
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# Runge--Kutta method program for the initial value problem (2).
include Math
##### constants #####
# calculation t from 0 to Tmax.
Tmax = 10.0
# Delta t
Dt = 1.0/100
# Initial value of u
IniU = [4.0, 1.0]
##### functions #####
# Right hand side of ODE
def f(u)
a = u[0]
b = u[1]
return [ (2.0-b)*a, (2.0*a-3.0)*b ]
end
# update of the approximate solution
def new_by_rk(u)
f1 = f(u)
f2 = f([ u[0]+0.5*Dt*f1[0], u[1]+0.5*Dt*f1[1] ])
f3 = f([ u[0]+0.5*Dt*f2[0], u[1]+0.5*Dt*f2[1] ])
f4 = f([ u[0]+Dt*f3[0], u[1]+Dt*f3[1] ])
a = u[0] + Dt * (f1[0]/6.0 + f2[0]/3.0 + f3[0]/3.0 + f4[0]/6.0)
b = u[1] + Dt * (f1[1]/6.0 + f2[1]/3.0 + f3[1]/3.0 + f4[1]/6.0)
return [ a,b ]
end
# result
def print_result(u,t)
printf("%.4f %.4f %.4f\n", t, u[0], u[1])
end
##### main program #####
# initialization
u = IniU
print_result(u, 0.0)
# main loop
for i in (1..(Tmax/Dt).ceil )
u = new_by_rk(u)
print_result(u,i*Dt)
end
}}
** 多段法 [#r0fc4e7e]
一段法の Runge-Kutta 法と同様に精度を(理論的には)自由にあげられる多段法として、線形多段階法を授業で説明した．
線形多段階法は、計算精度をあげても函数計算としての計算量が増えないという大変すぐれた性質がある．
ただし、解の函数がある程度滑らかであることが前提とされることや、過去の一定のデータの保存が必要なこと、時間離散幅 Δt をその場その場で変更しようとすると大変なこと、というあたりが代わりに要求される．
これについてもプログラム例をあげておこう．
*** サンプルプログラム 1 [#d604cf70]
上と同様に
du/dt = u(1-u) という常微分方程式に対する線形多段階法でのプログラム．
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# Linear multistep method with 4 stages to the initial value problem 1).
include Math
##### constants #####
Tmax = 10.0
Dt = 1.0/10
IniU = 0.4
##### functions #####
# Right hand side of ODE
def f(u)
return u*(1.0-u)
end
# update of the approximate solution (Runge-Kutta)
def new_by_rk(u)
f1 = f(u)
f2 = f(u + 0.5*Dt*f1)
f3 = f(u + 0.5*Dt*f2)
f4 = f(u + Dt*f3)
return u + Dt*(f1/6.0 + f2/3.0 + f3/3.0 + f4/6.0)
end
# predictor
def predictor(old_u,old_f)
return old_u + Dt * ((55.0/24.0)*old_f[0] - (59.0/24.0)*old_f[1] + (37.0/24.0)*old_f[2] - (3.0/8.0)*old_f[3])
end
# corrector
def corrector(old_u,old_f,new_f)
return old_u + Dt * ((3.0/8.0)*new_f + (19.0/24.0)*old_f[0] - (5.0/24.0)*old_f[1] + (1.0/24.0)*old_f[2])
end
# result
def print_result(u,t)
printf("%.6f %.6f\n", t, u)
end
##### main program #####
# initialization
u = IniU
old_f = [f(u)]
print_result(u,0.0)
# we make more three initial values by Runge-Kutta method
for i in (1..3)
u = new_by_rk(u)
old_f.unshift(f(u))
print_result(u,i*Dt)
end
# main loop by a LM with PECE cycle
for i in (4..N = (Tmax/Dt).ceil)
new_u = predictor(u,old_f) # P: prediction
new_f = f(new_u) # E: Evaluation
new_u = corrector(u,old_f,new_f) # C: Correction
new_f = f(new_u) # E: Evaluation
u = new_u
old_f.unshift(new_f)
old_f.delete_at(4)
print_result(u,i*Dt)
end
}}
*** サンプルプログラム 2 [#jfc36af8]
同様に
&br;
du/dt = (2-v)u,
&br;
dv/dt = (2u-3)v, &br;
という連立一階常微分方程式に対する線形多段階法でのプログラム．
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# Linear multistep method with 4 stages to the initial value problem 2).
include Math
##### constants #####
Tmax = 10.0
Dt = 1.0/100
IniU = [4.0, 1.0]
##### functions #####
# Right hand side of ODE
def f(u)
a = u[0]
b = u[1]
return [ (2.0-b)*a, (2.0*a-3.0)*b ]
end
# update of the approximate solution (Runge-Kutta)
def new_by_rk(u)
f1 = f(u)
f2 = f([ u[0]+0.5*Dt*f1[0], u[1]+0.5*Dt*f1[1] ])
f3 = f([ u[0]+0.5*Dt*f2[0], u[1]+0.5*Dt*f2[1] ])
f4 = f([ u[0]+Dt*f3[0], u[1]+Dt*f3[1] ])
a = u[0] + Dt * (f1[0]/6.0 + f2[0]/3.0 + f3[0]/3.0 + f4[0]/6.0)
b = u[1] + Dt * (f1[1]/6.0 + f2[1]/3.0 + f3[1]/3.0 + f4[1]/6.0)
return [ a,b ]
end
# predictor
def predictor(old_u,old_f)
a = old_u[0] + Dt * ((55.0/24.0)*old_f[0][0] - (59.0/24.0)*old_f[1][0] + (37.0/24.0)*old_f[2][0] - (3.0/8.0)*old_f[3][0])
b = old_u[1] + Dt * ((55.0/24.0)*old_f[0][1] - (59.0/24.0)*old_f[1][1] + (37.0/24.0)*old_f[2][1] - (3.0/8.0)*old_f[3][1])
return [a,b]
end
# corrector
def corrector(old_u,old_f,new_f)
a = old_u[0] + Dt * ((3.0/8.0)*new_f[0] + (19.0/24.0)*old_f[0][0] - (5.0/24.0)*old_f[1][0] + (1.0/24.0)*old_f[2][0])
b = old_u[1] + Dt * ((3.0/8.0)*new_f[1] + (19.0/24.0)*old_f[0][1] - (5.0/24.0)*old_f[1][1] + (1.0/24.0)*old_f[2][1])
return [a,b]
end
# result
def print_result(u,t)
printf("%.4f %.4f %.4f\n", t, u[0], u[1])
end
##### main program #####
# initialization
u = IniU
old_f = [f(u)]
print_result(u,0.0)
# we make more three initial values by Runge-Kutta method
for i in (1..3)
u = new_by_rk(u)
old_f.unshift(f(u))
print_result(u,i*Dt)
end
# main loop by a LM with PECE cycle
for i in (4..(Tmax/Dt).ceil)
new_u = predictor(u,old_f) # P: prediction
new_f = f(new_u) # E: Evaluation
new_u = corrector(u,old_f,new_f) # C: Correction
new_f = f(new_u) # E: Evaluation
u = new_u
old_f.unshift(new_f)
old_f.delete_at(4)
print_result(u,i*Dt)
end
}}
// ━┃┏┓┛┗┣┳┫┻╋
// コマンドライン入力は「行頭を > で始める」.
// コマンドライン出力は「行頭をブランクで始める」.
// 実習アイコン
// &ref(/materials/notes.png);
// 注意アイコン
// &ref(/materials/warning.png);
// Link アイコン
// &ref(/materials/JNorth_arrow-right-sm.png);
// OK アイコン
// &ref(/materials/OK.png);
// NG アイコン
// &ref(/materials/NG.png);
// サンプルアイコン
// &ref(/materials/Gnome-Preferences.png);
// 大文字での強調
// CENTER:&size(24){''ほげほげ''};
// programu source 表記
// #highlighter(language=ruby,number=on,cache=on){{}}
// パイプ
// ¦#contents
