読者です 読者をやめる 読者になる 読者になる

読み書きプログラミング

日常のプログラミングで気づいたことを綴っています

(39) 化学反応式

メタンが燃える反応では以下の変化が起こります。

反応式の係数を求めるスクリプトを作りましょう。

/* 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])
  )$