メタンが燃える反応では以下の変化が起こります。
反応式の係数を求めるスクリプトを作りましょう。
/* 7.1.4m, 7.1.8m */ reactants : [[C, 1, H, 4], [O, 2]]$ products : [[C, 1, O, 2], [H, 2, O, 1]]$ atoms : []$ lexp : 0$ numvars : 1$ for r : 1 thru length(reactants) do ( m : reactants[r], for j : 1 step 2 thru length(m) do ( if not member(m[j], atoms) then atoms : endcons(m[j], atoms), lexp : lexp + a[numvars]*m[j + 1]*m[j] ), numvars : numvars + 1 )$ rexp : 0$ for r : 1 thru length(products) do ( m : products[r], for j : 1 step 2 thru length(m) do ( rexp : rexp + a[numvars]*m[j + 1]*m[j] ), numvars : numvars + 1 )$ numvars : numvars - 1$ varlist : []$ for j : 2 thru numvars do varlist : endcons(a[j], varlist)$ e : []$ for j : 1 thru length(atoms) do ( e : endcons(diff(lexp, atoms[j]) = diff(rexp, atoms[j]), e) )$ a[1] : 1$ sol : solve(e, varlist)$ print(atoms)$ print(lexp)$ print(rexp)$ print(e)$ print(sol)$
であることが計算できました。
上のスクリプトは、a[1]=1を仮定したため、他の係数が分数になる場合があります。これを整数化するスクリプトも用意しておきましょう。
/* 7.1.6m */ den : 1$ for j : 2 thru numvars do ( k : at(a[j], sol[1]), if not integerp(k) then den : lcm(den, denom(k)) )$ for j : 1 thru numvars do ( a[j] : at(a[j] * den, sol[1]), print(j, " ", a[j]) )$
参考文献
- Crandall, Mathematica―理工系ツールとしての (アジソン ウェスレイ・トッパン情報科学シリーズ) p.187-p.190