微分方程式 1階常微分方程式の数値解法

常微分方程式の数値解法

微分方程式を

とする。初期値をとして、におけるyの値yiまでは計算済みとする。

【オイラー法】

この式により初期値(x0, y0)から出発して、逐次yi の値が求まる。これがオイラー法である。

【実習】

1次反応の速度式: につきオイラー法をおこなうと、

となる。これを計算して解析解と比較する。ただしyの初期値は1、・k=1とする。<rk.xls>

 

 

 

【ルンゲ−クッタ(Runge-kutta)法】

におけるyi+1を計算するために、

1.中点におけるyの値 y

オイラー法により計算する。

また、中点の勾配y’の計算をおこなう。

2.中点の値の再計算。

勾配の再計算。

3.におけるyの値の概略値をつくる。

それをもとにの概略値を求める。

4.以上4回計算された勾配を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