逐次反応 −微分方程式−

逐次反応のシミュレーション

最も簡単な逐次反応は,
A→R,  r1=k1CA
R→S,  r2=k2CR
である.各成分の濃度を表わす微分方程式は次のように書ける.
dCA/dt = -k1CA
dCR/dt =k1CA - k2CR
dCS/dt = k2CR

 初期濃度をCA=1,CR=0,CS=0としてk1=k2,  k1>k2,   k1<k2,の3つの場合について各成分の濃度の径時変化を横軸にk1×時間,縦軸に濃度をとって,右図のような様式で示しなさい. (A4方眼紙にて提出)

常微分方程式の数値解法

 常微分方程式は以下のようにRunge-Kutta法で数値計算がおこなわれる.
 微分方程式を y '=f (x ,y ), 初期値を y (x0 )=y とする.
x=xiにおけるyの値yiまでは計算ずみとする)
次の値 xi+1 = xi+h における yi+1 を計算するには,

常微分方程式解法プログラム

プログラム上では t→X,CA→Y(0),CR→Y(1),CS→Y(2),と変数を割り当てる.

(BEEP が鳴ってx=?での結果がでたら RET をおす)

1000 '====レンリツ ジョウビブン ホウテイシキ ========
1010 ' RUNGE-KUTTA METHOD FOR SYSTEMS OF DIFFERENTIAL EQUATIONS
1020 READ NM
1030 NM = NM - 1
1040 DIM Y(NM), Y0(NM), YW(NM), YS(NM), YQ(NM), YD(NM)
1050 READ YA, YB
1070 FOR NK = 0 TO NM
1080 READ Y0(NK): NEXT NK
1110 READ YH0, YE0
1130 YH1 = YH0: X0 = YA
1150 FOR NK = 0 TO NM: YS(NK) = Y0(NK): NEXT NK
1180 YH = YH1: X = X0
1200 GOSUB 1640
1210 FOR NK = 0 TO NM: YQ(NK) = YW(NK): NEXT NK
1240 YH = YH1 / 2: X = X0
1260 GOSUB 1640
1270 FOR NK = 0 TO NM: YS(NK) = YW(NK): NEXT NK
1300 X = X0 + YH
1310 GOSUB 1640
1320 YER = 0
1330 FOR NK = 0 TO NM
1340 YD(NK) = (YW(NK) - YQ(NK)) / 15
1350 E1 = ABS(YD(NK)) / YH1 / YW(NK)
1360 IF E1 > YE0 THEN 1400
1370 IF E1 > YE0 / 100 THEN YER = 1
1380 NEXT NK
1390 GOTO 1500
1400 '------シュウソク シナイ
1410 YH1 = YH1 / 2
1420 PRINT " x ノ キザミハバヲ"; YH1; "ニ カエマス"
1450 FOR NK = 0 TO NM: YQ(NK) = YS(NK): YS(NK) = Y0(NK): NEXT NK
1490 GOTO 1240
1500 '------シュウソク
1510 X0 = X0 + YH1
1520 WAIT: BEEP 3
1530 FOR NK = 0 TO NM
1540 Y0(NK) = YW(NK) + YD(NK)
1550 PRINT " y"; NK; "="; Y0(NK);
1560 NEXT NK
1565 PRINT USING "x=##.# y0=#.### y1=#.### y2=#.### RETヲオセ"; X0;Y0(0);Y0(1);Y0(2)
1570 IF YER = 1 THEN 1620
1580 IF YH1 = YH0 THEN 1620
1590 YH1 = YH1 * 2
1600 PRINT " x ノ キザミハバヲ"; YH1; "ニ モドス"
1620 IF X0 < YB THEN 1150
1630 END
1640 '------Runge-Kutta
1650 FOR NK = 0 TO NM: Y(NK) = YS(NK): YW(NK) = YS(NK): NEXT NK
1690 FOR NK = 0 TO NM
1700 GOSUB 2000
1710 YD = DYDX * YH
1720 YW(NK) = YW(NK) + YD / 6
1730 Y(NK) = YS(NK) + YD / 2
1740 NEXT NK
1750 X = X + YH / 2
1760 FOR NK = 0 TO NM
1770 GOSUB 2000
1780 YD = DYDX * YH
1790 YW(NK) = YW(NK) + YD / 3
1800 Y(NK) = YS(NK) + YD / 2
1810 NEXT NK
1820 FOR NK = 0 TO NM
1830 GOSUB 2000
1840 YD = DYDX * YH
1850 YW(NK) = YW(NK) + YD / 3
1860 Y(NK) = YS(NK) + YD
1870 NEXT NK
1880 X = X + YH / 2
1890 FOR NK = 0 TO NM
1900 GOSUB 2000
1910 YD = DYDX * YH
1920 YW(NK) = YW(NK) + YD / 6
1930 NEXT NK
1940 RETURN


2000 '---------ビブン ホウテイシキ
2010 '----ホウテイシキ ノ カズ
2020 DATA 3
2030 '----セキブン クカン [a,b]
2040 DATA 0,30
2050 '----ショキチ x=a ニオケル Y(0),Y(1),Y(2)...ノ アタイ
2060 DATA 1.0,0,0
2070 '----xノ キザミハバh, ゴサ ノ キョヨウ ゲンカイ E
2080 DATA 2, 0.01
2090 '----EQUATIONS( Y,X,N デハジマル ヘンスウ ハ ツカワナイ)
2095 K1 = .1: K2 = .1
2100 ON NK + 1 GOTO 2110, 2130, 2150
2110 ' [ DYDX=Y(0)' ]
2120 DYDX = -K1 * Y(0)
2125 RETURN
2130 ' [ DYDX=Y(1)' ]
2140 DYDX = K1 * Y(0) - K2 * Y(1)
2145 RETURN
2150 ' [ DYDX=Y(2)' ]
2160 DYDX = K1 * Y(0) - (K1 * Y(0) - K2 * Y(1))
2165 RETURN