Scilabで化学反応

私の練習のためにScilabで微分方程式による物理現象のモデル化(PDF)の問題を順番に解いてきました。最早、基本的な常微分方程式ならその物理的な背景を理解していなくてもルーチンワークでプログラムを書くだけで計算できてしまいます。

今回は、可逆反応と非可逆反応の反応が進む様子をシミュレーションします。

001_20130802112714a35.jpg

Fig.1: 化学反応は反応物質の数に依存して反応速度が異なる。また、可逆反応と非可逆反応では反応速度を表す微分方程式の形が異なる。画像は(c)GLOBALHAWK90



非可逆反応


非可逆な化学反応とは、反応する化学物質
非可逆反応の微分方程式とその解析解をプロットするScilabプログラムを作成します。

常微分方程式:
\frac{\mathrm{d}x}{\mathrm{d}t} = k(a-x)(b-x)

解析解:
x(t) = \begin{cases}ab\frac{\exp(gt) -\exp(-gt)}{a\exp(gt)-b\exp(-gt)} & (a > b) \\a - \frac{1}kt-\frac{1}{a}} & (a = b)\end{cases}

プログラムはunreversible_sce.txtです。

002_201308021127143ec.png

Fig.2: 非可逆反応の計算結果。


clear;

// ************************************
//
// 共通部分
//
// ************************************
// *** 入力パラメータ ***
a = 1;
k1 = input("k1 = ");
k2 = input("k2 = ");
// 時間ベクトル
T = linspace(0,10);

// ************************************
//
// 数値解
//
// ************************************
// 解くべき微分方程式
function dx = unreversible(t,x)
dx = k1 * (a - x) ^ 2 - k2 * x;
endfunction
// 初期値
x0 = 0;
// 微分方程式の数値解
X = ode(x0,0,T,unreversible);

// ************************************
//
// 解析解
//
// ************************************
b1 = (2 * a + k2 / k1 + sqrt(4 * a * k2 / k1 + (k2 / k1) ^ 2)) / 2;
b2 = (2 * a + k2 / k1 - sqrt(4 * a * k2 / k1 + (k2 / k1) ^ 2)) / 2;
d = sqrt(4 * a * k2 / k1 + (k2 / k1) ^ 2);
Xa = b2 * (exp(k1 * d * T) - 1) ./ (exp(k1 * d * T) - b2 / b1);

// ************************************
//
// プロット
//
// ************************************
plot(T,X,'or');
plot(T,Xa,'-g');
xtitle(strcat(["Chemical Reaction: ","a = ",string(a),", k1 = ",string(k1),", k2 = ",string(k2)]));
xlabel("$t [s]$");
ylabel("$N_x$");
legend(["Numerical","Analytical"],4);


可逆反応


可逆反応の微分方程式と解析解をプロットするScilabプログラムを作成します。

常微分方程式:
\frac{\mathrm{d}x}{\mathrm{d}t} = k_1 (a-x)^2 - k_2 x

解析解:
x(t) = b_2 \frac{\exp(k_1 \Delta b t)-1}{\exp(k_1 \Delta b t) - \frac{b_2}{b_1}}

ただし
b_1 = \frac{2a+\frac{k_2}{k_1} + \sqrt{4a\frac{k_2}{k_1}+\left(\frac{k_2}{k_1} \right)^2}}{2}

b_2 = \frac{2a+\frac{k_2}{k_1} - \sqrt{4a\frac{k_2}{k_1}+\left(\frac{k_2}{k_1} \right)^2}}{2}

\Delta b = b1 - b2 = \sqrt{4a\frac{k_2}{k_1}+\left(\frac{k_2}{k_1} \right)^2}


プログラムはbalanced_sce.txtです。

003_20130802112714957.png

Fig.3: 可逆反応の計算結果。


clear;

// ************************************
//
// 共通部分
//
// ************************************
// *** 入力パラメータ ***
a = 1;
b = input("b = ");
k = input("k = ");
// 時間ベクトル
T = linspace(0,10);

// ************************************
//
// 数値解
//
// ************************************
// 解くべき微分方程式
function dx = balanced(t,x)
dx = k * (a - x) * (b - x);
endfunction
// 初期値
x0 = 0;
// 微分方程式の数値解
X = ode(x0,0,T,balanced);

// ************************************
//
// 解析解
//
// ************************************
g = (a - b) / 2 * k;
if a == b then
Xa = a - 1 ./ (k .* T + 1 / a);
else
Xa = a * b * (exp(g * T) - exp(- g * T)) ./ (a * exp(g * T) - b * exp(- g * T));
end

// ************************************
//
// プロット
//
// ************************************
plot(T,X,'or');
plot(T,Xa,'-g');
xtitle(strcat(["Chemical Reaction: ","a = ",string(a),", b = ",string(b),", k = ",string(k)]));
xlabel("$t [s]$");
ylabel("$N_c$");
legend(["Numerical","Analytical"],4);


常微分方程式の数値解法


以下の通り、これまで微分方程式による物理現象のモデル化(PDF)に倣って、Scilabで常微分方程式の数値解法を扱ってきました。



2013年の更新はこれにて終了です。特にひねりは無いですが、来年は量子力学の問題から行こうと思います。
ただし、量子力学以降は単純な常微分方程式ではありません。常微分方程式は(ある意味何も考えずに)常微分方程式ソルバodeを使うだけで計算できてしまいますが、偏微分方程式は全自動で解いてくれるソルバが存在しないので少し考える必要があります。

元PDFには、(遍微分方程式ではない)普通の常微分方程式の問題として電子回路関係が残っていますが、LTspiceを使えば極めて簡単に計算できる回路をわざわざプログラム書きたくないという思いがあり、やる気がおきないです。

それでは、良いお年を。(と、予約投稿なので、書いている今はまだ7月なのですが。)

参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 常微分方程式 ode 化学反応 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoCCPAOPアンプPIC強磁性常微分方程式モンテカルロ解析トランジスタode状態密度DOSインターフェースecaljスイッチング回路定電流PDS5022半導体シェルスクリプト乱数レベルシフトHP6632A温度解析可変抵抗I2Cブレッドボード分散関係トランジスタ技術R6452A数値積分反強磁性バンドギャップ確率論セミナー絶縁偏微分方程式非線形方程式ソルババンド構造熱設計カオスA/DコンバータISO-I2Cフォトカプラ三端子レギュレータシュミットトリガLEDGW近似LM358アナログスイッチ数値微分TL43174HC4053マフィンティン半径発振回路サーボ直流動作点解析カレントミラーPC817CUSB単振り子bzqlty開発環境BSch2ちゃんねる電子負荷イジング模型LDAチョッパアンプ量子力学補間アセンブラFFTブラべ格子標準ロジックパラメトリック解析基本並進ベクトルewidthキュリー温度QSGWGGA失敗談MaximaSMPTLP621スイッチト・キャパシタ熱伝導コバルト相対論スピン軌道相互作用六方最密充填構造繰り返しFETランダムウォークcygwingfortran不規則合金状態方程式ラプラス方程式抵抗スレーターポーリング曲線位相図格子比熱マントルデータロガー自動計測ダイヤモンドガイガー管QNAPUPS固有値問題条件分岐井戸型ポテンシャルシュレディンガー方程式詰め回路MCU第一原理計算起電力熱力学スーパーセルVCALM555仮想結晶近似awkTLP521NE555ubuntufsolveブラウン運動OpenMPVESTA最大値テスタ差し込みグラフFXA-7020ZRWriter509三角波TLP552平均場近似最適化最小値過渡解析LMC662トランスPIC16F785CapSenseMBEナイキスト線図CK1026フィルタP-10負帰還安定性EAGLEAACircuit2SC1815OPA2277PGAノコギリ波縮退非線型方程式ソルバL10構造fcc面心立方構造結晶磁気異方性TeX全エネルギー固定スピンモーメントFSMウィグナーザイツ胞interp1ヒストグラム確率論マテリアルデザインspecx.fジバニャン方程式等高線初期値フェルミ面正規分布c/agnuplotBaO岩塩構造ルチル構造ウルツ鉱構造ZnO重積分SIC二相共存スワップ領域リジッドバンド模型半金属合金multiplotハーフメタルデバイ模型edeltquantumESPRESSOフォノンifortUbuntuマンデルブロ集合キーボードRealforce関数フィッティングフラクタルクーロン散乱CIF化学反応三次元最小二乗法日本語直流解析PCトラックボールExcelTS-110パラメータ・モデル等価回路モデルTS-112疎行列文字列HiLAPW両対数グラフ片対数グラフ熱拡散方程式陰解法境界条件連立一次方程式Crank-Nicolson法グラフの分割軸ラベルヒストグラム不規則局所モーメント入出力円周率Gimp凡例線種シンボルMAS830L

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ