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

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性OPアンプPICCPA常微分方程式モンテカルロ解析ecaljodeトランジスタ状態密度インターフェースDOS定電流スイッチング回路PDS5022半導体シェルスクリプトレベルシフト乱数HP6632AR6452AI2C可変抵抗分散関係トランジスタ技術ブレッドボード温度解析反強磁性確率論バンドギャップセミナー数値積分熱設計非線形方程式ソルババンド構造絶縁偏微分方程式ISO-I2CLM358フォトカプラ三端子レギュレータカオスLEDシュミットトリガGW近似A/Dコンバータ発振回路PC817C直流動作点解析USBマフィンティン半径数値微分アナログスイッチTL43174HC4053カレントミラーサーボ量子力学単振り子チョッパアンプ補間2ちゃんねる開発環境bzqltyFFT電子負荷LDAイジング模型BSch基本並進ベクトルブラべ格子パラメトリック解析標準ロジックアセンブラ繰り返し六方最密充填構造SMPコバルトewidthFET仮想結晶近似QSGW不規則合金VCAMaximaGGA熱伝導cygwinスレーターポーリング曲線キュリー温度スイッチト・キャパシタ失敗談ランダムウォークgfortran抵抗相対論位相図スピン軌道相互作用VESTA状態方程式TLP621ラプラス方程式TLP552条件分岐NE555LM555TLP521マントル詰め回路MCUテスタFXA-7020ZR三角波過渡解析ガイガー管自動計測QNAPUPSWriter509ダイヤモンドデータロガー格子比熱熱力学awkブラウン運動起電力スーパーセル差し込みグラフ第一原理計算フェルミ面fsolveCIFxcrysden最大値最小値ubuntu最適化平均場近似OpenMPシュレディンガー方程式固有値問題井戸型ポテンシャル2SC1815TeX結晶磁気異方性OPA2277非線型方程式ソルバフラクタルFSM固定スピンモーメントc/agnuplotPGA全エネルギーfccマンデルブロ集合縮退正規分布キーボード初期値interp1multiplotフィルタ面心立方構造ウィグナーザイツ胞L10構造半金属二相共存ZnOウルツ鉱構造BaOSIC重積分磁気モーメント電荷密度化学反応クーロン散乱岩塩構造CapSenseノコギリ波デバイ模型ハーフメタルRealforceフォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt合金等高線凡例軸ラベル線種シンボルトラックボールグラフの分割MAS830LPIC16F785トランス入出力CK1026PC直流解析パラメータ・モデル等価回路モデル不規則局所モーメント関数フィッティング日本語ヒストグラムTS-112ExcelGimp円周率TS-110LMC662片対数グラフ三次元specx.fifortUbuntu文字列疎行列不純物問題ジバニャン方程式ヒストグラム確率論マテリアルデザインP-10境界条件連立一次方程式AACircuit熱拡散方程式HiLAPW両対数グラフ陰解法MBEナイキスト線図負帰還安定性Crank-Nicolson法EAGLE最小二乗法

最新コメント
リンク

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