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カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoCOPアンプPICCPA強磁性常微分方程式モンテカルロ解析odeトランジスタ状態密度インターフェーススイッチング回路ecaljPDS5022DOS定電流半導体シェルスクリプト乱数レベルシフトHP6632Aブレッドボード分散関係温度解析R6452Aトランジスタ技術I2C可変抵抗反強磁性セミナー数値積分確率論偏微分方程式バンド構造非線形方程式ソルババンドギャップ絶縁熱設計シュミットトリガLEDA/Dコンバータ三端子レギュレータLM358ISO-I2CGW近似カオスフォトカプラマフィンティン半径TL431数値微分PC817Cアナログスイッチ直流動作点解析発振回路USBサーボカレントミラー74HC4053パラメトリック解析LDAbzqltyチョッパアンプ量子力学FFT2ちゃんねるアセンブラBSch開発環境電子負荷ブラべ格子イジング模型補間基本並進ベクトル標準ロジック単振り子キュリー温度繰り返しMaxima状態方程式失敗談相対論スピン軌道相互作用FETランダムウォーク熱伝導六方最密充填構造コバルトewidthTLP621GGAQSGW不規則合金位相図抵抗SMPcygwinラプラス方程式スレーターポーリング曲線gfortranスイッチト・キャパシタ詰め回路TLP552三角波格子比熱TLP521条件分岐LM555MCUNE555QNAPマントルテスタ過渡解析FXA-7020ZRダイヤモンドデータロガーガイガー管自動計測Writer509UPSシュレディンガー方程式ブラウン運動awk差し込みグラフ熱力学平均場近似仮想結晶近似VCAfsolve井戸型ポテンシャルVESTA起電力スーパーセルOpenMP第一原理計算ubuntu固有値問題L10構造OPA2277interp12SC1815fccウィグナーザイツ胞面心立方構造フィルタジバニャン方程式ヒストグラム確率論マテリアルデザインspecx.f等高線正規分布PGAフェルミ面非線型方程式ソルバ初期値固定スピンモーメントスワップ領域ルチル構造リジッドバンド模型edeltquantumESPRESSO岩塩構造BaOSIC二相共存ZnOウルツ鉱構造フォノンデバイ模型c/aノコギリ波全エネルギーFSMTeXgnuplotmultiplotハーフメタルCapSense半金属合金結晶磁気異方性Ubuntu文字列入出力TS-110TS-112疎行列Excel直流解析ヒストグラム円周率不規則局所モーメントトラックボールPC等価回路モデルパラメータ・モデルキーボードRealforce三次元マンデルブロ集合フラクタル化学反応重積分縮退日本語最小二乗法関数フィッティングGimpMAS830LHiLAPW熱拡散方程式両対数グラフナイキスト線図負帰還安定性陰解法Crank-Nicolson法P-10クーロン散乱境界条件連立一次方程式片対数グラフEAGLEPIC16F785LMC662トランスシンボルCK1026線種凡例MBEAACircuitグラフの分割軸ラベルifort

最新コメント
リンク

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