微分方程式 １階常微分方程式の数値解法
常微分方程式の数値解法
微分方程式を
とする。初期値を
として、
におけるyの値yiまでは計算済みとする。
【オイラー法】
この式により初期値(x0, y0)から出発して、逐次yi
の値が求まる。これがオイラー法である。
【実習】
１次反応の速度式：
につきオイラー法をおこなうと、
となる。これを計算して解析解と比較する。ただしyの初期値は1、・k=1とする。<rk.xls>
【ルンゲ−クッタ(Runge-kutta)法】
におけるyi+1を計算するために、
１．中点
におけるyの値 y中を
オイラー法により計算する。
また、中点の勾配y’中の計算をおこなう。
２．中点の値の再計算。
勾配の再計算。
３．
におけるyの値の概略値
をつくる。
それをもとに
の概略値
を求める。
４．以上４回計算された勾配
を1:2:2:1の重みで平均して、最終の
とする。
【実習】上と同様、表計算でルンゲクッタ法をおこなう <rk.xls>
【実習】Basicプログラムによりルンゲクッタ法の積分を組み込んだ表により、微分方程式を積分することでレーリー式の解を求める。dis6.xls
雛形 dis6_temp.xls
留出率θを
とすると、先の微分方程式 
は
となる。これをθ=0でx=x0を初期値として積分する。プログラム中ではθ：X0, x: Y0(0)の変数となる。微分方程式の部分は
' [
DYDX=Y(0)' ]
DYDX =-(A *
Y0(0) / (1 + (A - 1) * Y0(0)) - Y0(0)) / (1 - X0)
である。
ルンゲクッタ法BASICプログラム
Option
Explicit
Dim Y(),
Y0(), YW(), YS(), YQ(), YD()
Dim A, NM,
YA, YB, Lin, YH0, YE0, YH1, X0, NK, NZ, YH, X, YER, E1, YE1, YDD, DYDX, TP1,
TP2
Sub
A_onClick()
'----
NM = 1
NM = NM - 1
ReDim Y(NM),
Y0(NM), YW(NM), YS(NM), YQ(NM), YD(NM)
'---- [a,b]
YA = 0
YB =
Cells(6, 2)
'----
Y0(0) =
Cells(2, 2)
YH0 =
Cells(7, 2)
YH1 = YH0
X0 = YA
NZ = 0
Lin = NZ +
11
Cells(Lin,
1) = X0
Cells(Lin,
2) = Y0(0)
Do While (X0
/ YB) <= 1#
NZ = NZ + 1
For NK = 0
To NM
YS(NK) =
Y0(NK)
Next
YH = YH1
X = X0
Call Sub1640
For NK = 0
To NM
YQ(NK) =
YW(NK)
Next
'1240------
YH = YH1 /
2#
X = X0
Call Sub1640
For NK = 0
To NM
YS(NK) =
YW(NK)
Next
X = X0 + YH
Call Sub1640
For NK = 0
To NM
YD(NK) = (YW(NK) - YQ(NK)) / 15
E1 = Abs(YD(NK)) / YH1 / YW(NK)
Cells(Lin,
4) = E1
Next
'1500-----
X0 = X0 +
YH1
For NK = 0
To NM
Y0(NK) =
YW(NK) + YD(NK)
Next
Lin = NZ +
11
Cells(Lin,
1) = X0
Cells(Lin,
2) = Y0(0)
Loop
End Sub
Sub
Sub1640()
'------Runge-Kutta
For NK = 0
To NM
Y(NK) =
YS(NK)
YW(NK) =
YS(NK)
Next
For NK = 0
To NM
Call SUB2000
YDD = DYDX *
YH
YW(NK) =
YW(NK) + YDD / 6#
Y(NK) =
YS(NK) + YDD / 2#
Next
X = X + YH /
2#
For NK = 0
To NM
Call SUB2000
YDD = DYDX *
YH
YW(NK) =
YW(NK) + YDD / 3#
Y(NK) =
YS(NK) + YDD / 2#
Next
For NK = 0
To NM
Call SUB2000
YDD = DYDX *
YH
YW(NK) =
YW(NK) + YDD / 3#
Y(NK) =
YS(NK) + YDD
Next
X = X + YH /
2#
For NK = 0
To NM
Call SUB2000
YDD = DYDX *
YH
YW(NK) =
YW(NK) + YDD / 6#
Next
End Sub
Sub
SUB2000()
'---------ﾊﾟﾗﾒｰﾀ
A = Cells(1,
2)
'---------ﾋﾞﾌﾞﾝ
ﾎｳﾃｲｼｷ
Select Case
NK
Case 0
' [ DYDX=Y(0)'
]
DYDX = -(A *
Y0(0) / (1 + (A - 1) * Y0(0)) - Y0(0)) / (1 - X0)
End Select
End Sub