Excelで簡単!微分方程式

「化学工学資料のページ」に戻る
伊東研


 

1.微分方程式とは方程式を微分したもの

方程式と微分方程式(常微分方程式)は同じ関数の別の表示法である。しかし微分すると絶対値の情報が欠けるので,微分方程式には初期値が必要となる。

【例題1】微分方程式で方程式を再現する<rk01.xlsm>

 

2.Excel VBAによる「常微分方程式解法シート」

自然現象や物理現象の数式モデルの多くは微分方程式となる。微分方程式が1階の が与えられれば,初期値から出発して数値積分できる。 このような数値積分の手法としてはルンゲ・クッタRunge-Kutta法が標準的である。

また,多階の常微分方程式も1階の正規形常微分方程式の連立式に書き直せるので,初期値がすべて与えられていればRunge-Kutta法により解くことができる。

ここでは正規形の連立常微分方程式を解くための,「微分方程式解法シート」を提供する。シート内にはRunge-Kutta法 およびRunge-Kutta-Fehlberg法のVBAプログラムを組み込んである。しかし,VBAプログラム中の変数とシート上のセルの値を対応させておくことで,式,定数をVBAプログラム中に書く必要がないようにした。微分方程式はシート上のセル内にセル座標で記述する形式なので,いろいろな常微分方程式の解法に使える。

> 参考 (同じようにExcelのセルに記述する形式で常微分方程式を解く方法を考案している例)
舞鶴高専 伊藤教授 講義資料 微分方程式熊本大学 趙教授 Excel VBAによる数値計算法についての研究POLYMATH ODE_Solver| 北陸電力 
ミックスオイラー法Excelで使える数値計算ツールSUITEXL

3.生物学の微分方程式−生態系のモデル−

【例題2】生態系のモデル<rk02.xlsm>

生態系における被食者と捕食者の生存闘争モデルとして有名なLotka-Volterraモデルは次式の連立常微分方程式である。

ここではx をウサギ(被食者)の個体数,y をキツネ(捕食者)の個体数とする。t は時間で,両式の第2項(xy )が個体群間の相互作用を表す「交差項」である。ウサギが1000匹,キツネが100匹の状態から,このモデルにより個体数の変化を予測せよ。a =0.01, b =0.05, c =0.0001とする。

図が「微分方程式解法シート」である。セルB1に微分方程式の数2を書く。定数をG2:G4に記しておく。B5に式(1)を,C5に式(2)を記述する。その際,式中の変数x, y はセルB3, C3を用いる。時間t の積分区間と刻み幅をB7:B9に設定する。x, y の初期値1000, 100をB12,C12に記入して,ボタン”Runge-Kutta”をクリックすると連立常微分方程式の積分が実行される。結果がシートの13行以下に出力された。 計算結果を図中のグラフで示す。初期はキツネが増え,ウサギが減少するが,ウサギが減少しすぎるとキツネも減少し,その後ウサギは増加する。このような生態系の周期変化が,簡単なモデル式により表現されていることがわかる。

 

4.化学の微分方程式−逐次反応−

【例題3】逐次反応 <rk03.xlsm>

反応物AがRを経てSになる簡単な逐次反応:

では,各成分の濃度を表わす連立常微分方程式が次のように書ける。

ここでCi は各成分の濃度,r は反応速度,k1, k2 は反応速度定数である。k1 = 0.2,k2 = 0.2,初期濃度をCA= 1, CR = 0, CS = 0 として濃度の時間変化を示せ。

図が「微分方程式解法シート」である。定数をG2:G3に,上式をB5:D5に記述する。積分区間・刻み幅,初期値(B12:D12)を設定し,ボタンクリックで積分を実行する。図中のグラフに各成分の濃度変化の様子を示す。

 

5.量子化学の微分方程式−Schrödinger 方程式−

【例題】Schrödinger 方程式 <schroed.xlsm>

量子化学の基礎式,シュレーディンガ−方程式Schrödinger equationは波動関数wave function φ に関する2階の常微分方程式である。(定常状態。左辺の定数 2m/h2は1とする。)

 

E が固有値であり,波動関数の解は固有値E のとびとびの値で存在する。ここではV に井戸型ポテンシャルwell-type potential:

を考えて解いてみる。 (参考資料:物理のかぎしっぽ:井戸型ポテンシャルの数値解)基礎式を次式の正規形の連立常微分方程式 :

として解く。固定刻み幅のRunge-Kutta法による。

 

 

6.物理の微分方程式−粒子の飛跡−

【例題4】粒子の飛跡<rk04.xlsm>

重力gのみ考慮して,射出された物体の軌跡を求める。水平方向をx 座標,垂直方向をy 座標とし,各座標の粒子の位置[m]をx, y,粒子の速度[m/s]をux, uyとし,時間をt [s]とすると粒子の運動方程式は次式である。

【参考】粒子の飛跡シミュレーション<fluid3.xlsm>
床面上0.2 mの高さの点から,密度ρ p = 2500 kg/m3, 粒径D p = 0.0003 mの球形粒子が初速度u 0 = 15 m/s, 迎角 45度(=α0)で静止空気中に放出された。このときの粒子の飛跡を求め,床面に落ちる位置を決定せよ。ただし,空気密度ρ f = 1.2 kg/m3, 空気粘度μ f = 0.000018 Pa・sである。

 水平方向をx 座標,垂直方向をy 座標とし,各座標の粒子の位置[m]をx, y,粒子の速度[m/s]をux, uyとし,時間をt [s]とすると粒子の運動方程式は次式である。

また,ux, uyの初期値は共にu0sinα0=u0cosα0=10.61 m/sである。よって,x, y, ux, uy に関する連立常微分方程式(11.5)〜(11.8)を,初期条件t =0でx = 0, y =0.2, ux=uy= 10.61により解く問題となる。

図が「微分方程式解法シート」である。諸定数およびu, Re, CDの計算式をセルG2:G9に記入・記述する。このときux, uyはD3, E3とする。方程式数(B1)を4とし,B5:E5に微分方程式(11.5)〜(11.8)を記述する。積分区間と刻み幅を設定し,初期値をB12:E12に書く。ボタンクリックでRunge-Kutta法による積分が実行される。その結果から得られた粒子の飛跡(x-y)を図中に示す。また床面に落ちる位置(y =0 )は計算結果よりx =0.1735 mである。

 

7.物理の微分方程式−主人と犬問題−

【例題】主人と犬問題 the master-and-dog problem<dog.xls>

主人はy 軸上を(0,0)から出発して速度v で歩く。犬は主人の方向に向かって速度w で走る。犬の位置を(x ,y )  として,犬の軌跡は次式の2階の微分方程式となる。

 ここでy ’ =z とすると,2階の微分方程式が以下の連立常微分方程式となる。 刻み幅可変のFehlberg法により解く。

 

 

.機械の微分方程式−バネの振動−

単振動の微分方程式は重力とフックの法則より、

 である。x は位置,おもりの質量をM, kをバネ定数として,ω2 =k/M である。

この2階の常微分方程式をルンゲ・クッタ法により解く。正規形とするため速度v (=dx/dt)を用いると,上式は

のような連立常微分方程式となる。初期条件は、
x =1, v =0 (距離1引っ張って離す場合) または
x =0, v =1 (静止状態からおもりに速度を与える場合)
などが与えられる。

 【例題5】単振動 <rk05.xlsm>

表計算+BASICで振動の微分方程式の解を求める。これはバネ振動のシミュレーションと言える。初期条件を変えて振動の様子をみる。

 

さらに、図のようにダンパーを加えた系を考える。

力のバランスは、
(質量Mによる慣性力)+(ダンパーの制動力)+(ばねの力)=(外力) より、

である。これは速度v を位置 で書き換えると、v =dx/dt なので,x に関する2階の常微分方程式:

 である。 (1)式は以下の正規形連立常微分方程式となる。

【例題6】バネマスダンパ系<rk06.xlsm>

9.電気回路の微分方程式−回路の過渡現象−

図のようなR(抵抗), L(コイル), C(コンデンサ)からなる回路で,スイッチを入れたときの過渡現象を調べる。直列回路の関係式より,

である。これは電荷 q と電流 i に関する次式の連立常微分方程式となる。

【例題7】R-L-C直列回路の過渡現象<rk07.xlsm>

 

10.制御

15.1  液面制御系のモデルと動特性

化学工学ではタンクの液面レベル系が制御のモデルとしてよく使われる。底面積A [cm2]のタンクに液面レベルh[cm]で液が入っており,これに流入流量q [cm3/s]で液が流入し,タンク底からバルブを通して流出している。タンク底からの流出流量は液面レベルh とバルブの抵抗R に依存し,(流出流量)=h/R (注)である。この1槽モデルを直列に接続したのが,右図の2槽モデルである。

(注:=√h/R が現実的である(「トリチェリの定理」)が,ここでは簡単のためh に比例するとした。このようにすると,機械工学におけるバネ−ダンパ系の問題および電気工学におけるRC回路問題とモデル式が同じになる。)

 制御では値を制御したい変数すなわち被制御変数と、これを操作するため変数すなわち操作変数がある。液面レベル系では被制御変数が液面レベルh ,操作変数が流入流量 q である。こ

1槽モデルにおいて,物質収支((タンク内液量変化)=(流入流量)―(流出流量))からh q の関係は次式である。

 この形式のモデルを1次遅れ系という。

制御のためには先ず操作変数の特定の変化に対する被制御変数の変化(応答)を知る必要がある。これをプロセスの動特性という。ここでは系が定常状態(qss, h (0))にあるところから,時間t ≧0 で操作変数q qd にステップ変化した場合のh の変化,すなわちステップ応答を考える。

 【例題15.1】1次遅れ系のステップ応答<contr1.xlsm>

 1槽モデルでqss =20 cm3/s,h (0)=20 cmの定常状態であった。qqd=30 cm3/sにステップ変化させた場合の応答を求めよ。

 式(15.1)中のq を定数qdとして,h に関する以下の微分方程式を初期条件h (0)=20 で解く問題となる。

図が「微分方程式解法シート」である。本章では以下全てこのシート上で計算をおこなう。定数をセルG2:G4に書き,B5に式(15.4)を,B7:B9に積分区間・刻み幅を記述する。初期値をB12に設定し,ボタンクリックで積分を実行する。結果を図中のグラフで示す。1次遅れ系のステップ応答は指数関数となる。

15.2 比例制御

このような動特性を持つプロセスに対して,フィードバック制御系を構成し,液面高さの目標値をhsetに設定変更した場合の応答を求める。制御では被制御変数h と設定値hsetとの差すなわち偏差e

 をもとに操作変数q を操作する。比例制御では操作変数を偏差e に比例させて操作する。

 Kc を比例ゲインという。初期状態(t <0)では操作変数qss,被制御変数h (0)の定常状態にあるとして,この式を微分すると次式である。

【例題15.3】1次遅れ・比例制御系の応答<contr3.xlsm>

 1槽モデルでqss =20 cm3/s,h (0)=20 cmの定常状態であった。t ≧0で新たなhの設定値hset=30を与えて,比例制御によるh の応答を示せ。Kc =2とする。

1次遅れ系のモデル(式(15.1))と比例制御の基礎式(式(15.9)) による以下の連立常微分方程式を解く問題となる。

図がこれを解いたシートである。セルG4:G6に定常値qss, 設定値hsetKcの値を設定する。B5,C5に式(15.10), (15.11)を記述し,B7:B9に積分区間,刻み幅,B12:C12に変数の初期値を入れ,ボタンクリックで積分を実行する。得られたq, hの時間応答を図中のグラフに示す。比例制御では設定値hset=30に至らずに定常値になる。これがオフセットである。

 

15.3 比例・積分制御

比例・積分制御では偏差eの時間t =0からの積分値を計算し,これが0になるまで操作変数を変える。

TI が積分時間である。上式を微分した形式が次式である。

【例題15.5】1次遅れ・比例・積分制御系の応答<contr5.xlsm>

1槽モデルでqss=20 cm3/s,h (0)=20 cmの定常状態であった。t ≧0で新たなhの設定値hset =30を与えて,比例・積分制御によるh の応答を示せ。

モデル(式(15.1))と制御式(式(15.16))から次式の連立常微分方程式を解く問題となる。

図のシートでG2:G7に定数を入れ,B5:C5に式(15.17),(15.18)を記述する。ここで式(15.18)中の(dh/dt) 項はB5を用いる。積分区間,初期値を設定してボタンクリックで積分を実行する。得られたq, h の時間応答を図中のグラフに示す。積分を加えることでオフセットを無くすことができる。