反応ネットワークの化学平衡を求める 

平衡性を持つ化学反応は，通常整方程式系として直接モデル化することができる．これらの問題は保存と反応の方程式により定義される．以下の例題は燃焼化学のものであり，車のエンジンなどの燃焼室で燃料を燃やすことをモデル化している．このモデルで，さまざまな段階で燃焼室にどのような化学成分や化学物質が存在するかを調べる．調べる反応系は，，，，，，で与えられる．

各成分の量に対する変数を導入する．は単独では存在しないので，ここではを化学成分として扱う．

各化合物の量に対する変数も導入する．

In[1]:=

X vars = {Subscript[x, 1], Subscript[x, 2], Subscript[x, 3], Subscript[ x, 4], Subscript[x, 5], Subscript[x, 6], Subscript[x, 7], Subscript[x, 8], Subscript[x, 9], Subscript[x, 10]};

閉化学系を想定する場合，量は保存される．これにより保存方程式がその化学および代数表記で表示される．

In[2]:=

X conservation = {Subscript[x, 2] + 2 Subscript[x, 6] + Subscript[x, 9] + 2 Subscript[x, 10] == Subscript[T, H], Subscript[x, 3] + Subscript[x, 8] == Subscript[T, CO], Subscript[x, 1] + Subscript[x, 3] + 2 Subscript[x, 5] + 2 Subscript[x, 8] + Subscript[x, 9] + Subscript[x, 10] == Subscript[T, O], Subscript[x, 4] + 2 Subscript[x, 7] == Subscript[T, N]};

平衡方程式を与える．

In[3]:=

X equilibrium = {Subscript[x, 5] == Subscript[R, 1] Subscript[x, 1]^2, Subscript[x, 6] == Subscript[R, 2] Subscript[x, 2]^2, Subscript[x, 7] == Subscript[R, 3] Subscript[x, 4]^2, Subscript[x, 8] == Subscript[R, 4] Subscript[x, 1] Subscript[x, 3], Subscript[x, 9] == Subscript[R, 5] Subscript[x, 1] Subscript[x, 2], Subscript[x, 10] == Subscript[R, 6] Subscript[x, 1] Subscript[x, 2]^2};

パラメータ値を与える．

In[4]:=

X pars = {Subscript[T, O] -> 5 10^-5, Subscript[T, CO] -> 3 10^-5, Subscript[T, H] -> 1 10^-5, Subscript[T, N] -> 1 10^-5, Subscript[R, 1] -> 24.528, Subscript[R, 2] -> 22.206, Subscript[R, 3] -> 47.970, Subscript[R, 4] -> 24.942, Subscript[R, 5] -> 22.120, Subscript[R, 6] -> 46.989};

実数上で系を解き，非負の解だけを選ぶ．

In[5]:=

X rsols = NSolve[Join[conservation, equilibrium] /. pars, vars, Reals];

In[6]:=

X Thread[vars -> SelectFirst[vars /. rsols, VectorQ[#, NonNegative] &]]

Out[6]=

高精度解を計算する．

In[7]:=

X rsols = NSolve[ SetPrecision[Join[conservation, equilibrium] /. pars, 20], vars, Reals, WorkingPrecision -> 20];

In[8]:=

X Thread[vars -> SelectFirst[vars /. rsols, VectorQ[#, NonNegative] &]]

Out[8]=