自己触媒反応に現れる多根を持つ２変数の非線形方程式の簡便な数値解法

Amih SAGAN, 長嶋 雲兵, 森 義仁

Return

1 はじめに

　自己触媒反応の反応物とその滞留時間との関係[1, 2]など、化学における多くの問題は、多くの変数の高次の複雑な非線形方程式系で記述されることが知られている。非線形方程式の解法は、或区間に関する解の存在を統計的に推定する方法[3]や古典的なNewton-Raphson法を基礎とする方法[4]など数多く提唱されているが、それらの多くは一般的に多根をもつ非線形方程式つまり一つのxに対していくつかのyの解があるような非線形方程式に適用することが難しい。そこで、我々は２変数の非線形方程式F(x,y)=0を満足するxとyの関係が連続で単調な多価関数である場合、つまり一つのxに対していくつかのyの解があり、それらが一筆書きで連続的にたどれる場合の解を数値的に得る簡便な算法を開発したので報告する。
　次節でその算法を簡単に説明し、３節では、２つの簡単な例を用いて本方法の特徴と限界及び本法で導入されたパラメータの選び方を説明する。さらに本方法は微分方程式の初期値問題の簡便な数値解法として知られているオイラー法における時間幅の調節の指針をあたえるので、４節で光照射によって振幅と周期が変動する振動反応[5]を例にとり、その有効性を説明する。そして最後に本方法の特徴をまとめる。

2 算法

　２変数の任意の非線形関数F(x,y)を考える。

F(x, y)=0(1)

　式(1)はxとyのある種の依存関係を示していると考えることができる。この場合xとyは、例えば自己触媒反応の反応物とその滞留時間などである。
　式(1)をxで微分すると、

∂F/∂x + ∂F/∂y ・dy/dx　=0(2)

を得る。これを変形して、

dy/dx = -(∂F/∂x)/(∂F/∂y) = g(x, y)(3)

を得る。式(1)の解であるxとyの関係は、もし特別な解（x0, y0）が既知であるならば、それを初期値として式(2)または(3)で与えられる微分方程式の初期値問題を解けばよい。
　非常に簡便な微分方程式の初期値問題の解法としてオイラー法がある。オイラー法は簡単であるが、近似精度が悪く非常に細かな時間刻み幅を用いないと十分な精度が得られないことが知られているが、それを用いるとすると、式(1)の数値解は、

xi+1=xi+Δx(4)
yi+1=yi + g(xi, yi)Δx(5)

で与えられる。ここでΔxは時間刻み幅であり、i = 0,1,2.....である。式(4)と(5)はyがxの１価単調連続な関数の場合に適用可能である。これを簡単にyがxの多価の単調で連続な関数の場合に拡張するには、g(x, y)の絶対値があるしきい値Gより大きくなったときに積分の軸をｘからｙに変更すればよい。すなわち|g(x, y)| > G　（ここで軸の回転のしきい値Gは、適当な値である。この値の取り方については次節で議論する。）となった場合に軸を変更する。
　軸を変更した時の式は、

yl+1 = yl +Δy(6)
xl+1 = xl +ξlΔy/g(xl, yl)(7)

となる。ここでξlは多価の単調で連続な関数に拡張するために新たに導入した、xの進行方向を変更させるパラメータであり、

ξl= 1 　 　 　 (xl > xl-1)
　 = -1　 　 　 (xl < xl-1)

である。このξlをもちいた、式(4)と(5)に対応する|g(x, y)| < G の場合の式は以下のようになる。

xi+1=xi+1 +ξi+1Δx(8)
yi+1=yi+1 +ξi+1g(xi+1, yi+1)Δx(9)

　式(6)と(7)の組と式(8)と(9)の組を交互に用いることで、yがxの多価の単調で連続な関数の場合の数値解を得ることができる。
　ここで、通常の微分方程式の解法という観点から式(7)をながめて見ると、この場合、常にξl=1である。ここで式(7)のΔy をTyと書くとすると、Tyをパラメータとする刻み幅自動調節機能付きオイラー法の式を得ることができる。

Δxi= Ty/|g(xi, yi)| (10)
xi+1=xi+Δxi
yi+1=yi + g(xi, yi)Δxi

　式(10)はｘの刻み幅を直前の微分の値によって調整する事を示している。

3 簡単な例

　ここでは、x2+y2-1=0とx4+y4-2y(y2-x2)=0というよく知られた解析的な関数について本方法の適用例を示し、本方法の性質およびパラメータＧの取り方を説明する。

3. 1 x2+y2-1=0 (Circle) の場合

　x2+y2-1=0は、よく知られた非線形方程式であり、この場合yは -1<x<1でxの２価の関数となる。
　解くべき微分方程式は

dy/dx= g(xi, yi) = -x/y

である。
　Figure 1に、初期値(0.0,-1.0)　Δx=0.02でG=1.0の場合を示した。この場合、オイラー法の数値計算上の誤差のために閉じた円にはならない。もちろん刻み幅を十分小さくすれば徐々に正しい円になっていく。


 Figure 1. The case of x2+y2-1=0. Initial value: x=0.0, y=-1.0, Dx=0.02, G=1.0


Figure 2. Calculated dy/dx for the case of x2+y2-1=0

　Gの依存性を見るためにG=10.0と100.0で計算を実行してみたが、どちらでもほぼ同様な結果が得られた。ここには具体的には示さないが、Gの違いによる誤差は、オイラー法の数値計算上の誤差に比べると平均して１桁以上小さい。これは、Figure 2に示すように計算された微係数dy/dxの変化がx=0付近で比較的緩やかで、x=1.0と x=-1.0付近で急激に変化して、x=1.0と x=-1.0で発散するためである。ただし数値計算上は、発散する領域ではすでに軸の変換が行われていおり、ｙ軸方向に積分が行われるため、x=1.0または-1.0付近でのオーバーフローは起こらない。
　この結果を微分方程式の数値解法という観点からみると、-1.0< x <1.0でG>1.0の時すなわちxが-1.0から-0.5付近に近づくにつれて刻み幅Δxが自動的に大きくなっていき、また0.5付近から1.0になるにつれて小さくなっていく様に見える。

3. 2 x4+y4-2y(y2-x2)=0の場合

　x4+y4-2y(y2-x2)=0は、Figure 3に示すような形を示すことが知られている。この場合yは、１価から４価の関数となる。


 Figure 3. The case of x4+y4-2y(y2-x2)=0. Initial value: x=0.0, y=2.0, Dx=0.02, G=1.0

　この場合、解くべき微分方程式は

dy/dx= g(xi, yi) =-(2y3+2xy)/(2x3-3x2+y2)

である。Figure 3は初期値 x=0.0,y=2.0, Δx=0.02でG=1.0の場合を示した。
　この場合も先の円の場合と同様、オイラー法の数値計算上の誤差のためにサークルは閉じない。Figure 3とFigure 4からも判るようにこの場合dy/dxは６カ所で無限大に発散する。この関数は、先の円に比べて微係数の変化の度合いが大きい。　
　この場合も、Gの依存性を見るためにG=10.0と100.0で計算を実行してみた。Figure 5とFigure 6にそれぞれの場合の計算結果を示す。Figure 5、Figure 6に示すようにG=10.0および100.0では期待した解が得られていない。またここには例示しないがG<1.0の場合も正しい解が得られない。それは、式(5)と式(7)から判るようにG=g(x, y)=1.0の場合は軸変換前と後でg(x, y)と1/g(x, y)の値に大きな違いがないにもかかわらず、例えばG=10.0の場合にでも両者にはオーダー２の違いが生じる。そのためオイラー法がもつ数値的誤差が悪影響をおよぼし、正しく解をたどることができなくなっているためである。つまり軸変換のパラメータGの値として1.0が最適である。先の円の場合はGの依存性が無いようにみられたが、これは偶然であったことが判る。


 Figure 4. Calculated dy/dx for the case of x4+y4-2y(y2-x2)=0


Figure 5. The case of x4+y4-2y(y2-x2)=0. Initial value: x=0.0, y=2.0, Dx=0.02, G=10.0


Figure 6. The case of x4+y4-2y(y2-x2)=0. Initial value:x=0.0, y=2.0, Dx=0.02, G=100.0

4 光照射によって振幅と周期が変動する振動反応の例

　森ら[5]は、硫酸酸性、イオン触媒共存下で、マロン酸の臭素酸による酸化反応（BZ反応）について、混合溶液の酸化還元電位の振動の振幅と周期の光照射による変動を観測し、それを詳細に解析した。その振動の振幅と周期は、光照射のタイミングと強度によって変動する。
　我々は、よく知られている振動反応のモデルであるオレゴネータ[1]を若干修正して、実験事実を矛盾無く説明するモデルを構築した。以下の反応モデルがそれらを説明する事がわかった。

dx/dt = s(y - xy + qx2 + A exp(-a(ti-t)2)(c-z))
dy/dt= (fz - y - xy)/s(11)
dz/dt= w(x - z + A exp(-a(ti-t)2)(c-z))

　ここで、x, y, zはそれぞれ無次元化されたHBrO2、Br-、Ce4+の量であり、cは光励起するCe4+の総量、s, f, wは通常のオレゴネータで用いられるパラメータである。　dx/dtとdz/dtの式に現れるガウス型関数A exp(-a(ti-t)2)(c-z)の項が光照射の影響を表現しており、強度　A、と半値幅aのガウス型の光が 反応開始後ti秒に入射する。詳細は[5]にゆずるが、ガウス型の光の照射によって瞬間的にHBrO2(x)とCe4+(z)が増加するというモデルである。このモデルでは、光照射はBr-の変動に間接的にしか影響しない。この問題は変数が、x, y, z, tという４変数の場合であるが、われわれの方法はこのような場合でも通常の連立一次常微分方程式の解法と同様に解くことができる。今回はzの増減に注目するためにｚとtの振る舞いに注目する。
　このようなオレゴネータ型の微分方程式の解は、微分の値がマイナス無限大から０を経て、無限大近くまで急激に変化する櫛形または鋸歯形の解を与えることが知られている。このような場合に通常のオイラー法を適用すると数値解が発散してしまいオーバーフローが起こったり誤差の蓄積が大きくなったりして、数値的には安定に解くことができないことが知られている。有効な解を数値的に得るためには、時間刻み幅を10-8程度に固定したオイラー法でも十分でなく、また、たとえ固定刻み幅のルンゲクッタ法を用いても、時間刻み幅を非常に小さくして積分を行わなければならない。また、固定刻み幅では、本例のように振動周期が変動する場合、正しいピーク位置および強度を得ることが難しい。さらに、積分の刻み幅を小さくすることは、計算量の増大を招き、また数値計算上の誤差の蓄積も多くなる。
　そのため通常は、ルンゲクッタの４次または５次の公式と、ブルリッシュとストアの方法[6]やフェールベルグの方法[7]またはメルルッツィらの方法[8]などいくつかある時間幅自動調節機能を組み合わせて積分を実行することが必要である。しかし、当然であるが、それらは計算時間がかかることが知られている。
　ある初期条件での性能評価の例としてTable 1に、式(11)のうちｚの量の振動周期の変化に注目した場合の計算された第一ピークの位置、強度、その際の微分計算の回数およびその時点での時間刻み幅の例を示した。この計算例は、t=0から微係数が負から0となり、t=290あたりから急激に正の値をとり、t=297あたりのピーク位置でまた急激に微係数の符号が変わるプロセスを追っている。
　ここで比較に用いた方法は、ルンゲ・クッタ４次法にフェールベルグ（RKF）[7]とメルルッツィ（RKM）[8]による自動時間調節刻み法を加えたもの、および本法（SNM）である。初期データとしての時間幅は1.0×10-5である。またSNMで、式(10)の時間調節パラメータTyは1.0とした。Table 1の数値は、初期条件その他によって変化するため大まかな性能評価を与えるにすぎないが、性能を比較する目安にはなる。
　Table 1から判るように、ピーク位置及び強度はRKFとRKMおよびSNMで大きな差は見られない。オイラー法の離散化誤差がO(Δx2) で累積誤差がO(Δx)であり、ルンゲ・クッタ４次法のそれらが、それぞれO(Δx5) とO(Δx4)であることを考慮すると、ルンゲ・クッタ法とオイラー法の結果がこの程度の精度で一致するのは驚きである。この一致は偶然かもしれない。

Table 1 Position and Value of the first Peak, No. of dx/dt calculation, and Dt at the Peak.
RKFRKMSNMdifference/retio*
1st Peak(RKF-SNM)(RKM-SNM)
Position (sec.)297.34297.34297.341.0×10-71.0×10-7
Value (×104)4.0091824.0091834.0092062.4×10-52.3×10-5
No. of dx/dt calc.120102541835004240300102.98*4.55*
Dt at the Peak1.3633×10-31.8723×10-37.34×10-7
RKF: Runge-Kutta- Fehlberg: 4th and 5th order[7].
RKM: Runge-Kutta-Merluzzi-Brosilow: 4th order[8].
SNM: This work.

　この例では、最適な時間刻み幅が初期値1.0×10-5から徐々に大きくなり、また減少するプロセスとなっている。最大ピーク位置での時間刻み幅は、RKFとRKMは初期入力より大きくなっているのにくらべ、SNMは逆に小さくなっていることが判る。ちなみに結果の詳細は割愛するが、刻み幅1.0×10-5の固定刻み幅のルンゲ・クッタ４次法または５次法では、4.00×10+4以上の高さを持つピークは見いだせなかった。また、ピーク位置もTable 1に掲げた方法による結果とは異なっていた。
　計算性能を評価する場合、計算の実行時間を比較することが多いが、実行時間はプログラムの実装により大きく違うため、微分計算のためのサブルーチンの呼び出し回数で性能を評価することにする。[9]　今回用いたプログラムは非常に素朴に書かれたものであり、微分計算の部分はどの方法も共通である。微分計算サブルーチンの呼び出しの回数は、時間刻み幅がルンゲクッタ法とオイラー法で全く同じ場合、オイラー法の場合はルンゲクッタ４次法の1/4となる。本例の場合RKF、RKMに比較するとSNMのそれは、それぞれ1/3または1/5程度になっている。これは、SNMがdx/dt=0の近傍で時間刻み幅をRKFやRKMなどにくらべ大きくとっており、逆にピーク近傍のdx/dtの絶対値が大きくなる領域では、小さくとっていることを示している。またRKFはSNMの1/3程度の微分計算回数であるので、RKFがSNMに比べ効率のよい時間刻み幅の選択を行っていることを示している。
　先にのべたように最大ピーク位置での時間刻み幅は、RKFとRKMは初期データより大きくなっているのにくらべ、SNMのそれは逆に小さくなっていることと、ピーク位置とその値の計算精度がほぼルンゲクッタ法と同等となっていることを考慮すると、SNMの時間刻み幅調節機能はRKFやRKMに比べ、柔軟になっていることがうかがえる。SNMの時間刻み幅調節機能の妥当性に関しては、より多くの例で確かめる必要がある。

5 まとめ

　開放系の自己触媒反応などに現れる多根を持つ２変数の非線形方程式F(x, y)=0の簡便な数値解法とその応用例を報告した。本方法は、そのような非線形方程式の解法を、既知の解の一つを初期値とする微分方程式の初期値問題の解法に帰着するため、xとyの関係が連続で単調な多価関数である場合にのみ有効である。本稿では微分方程式の初期値問題の解法としてのオイラー法を積分軸変換パラメータGと積分方向パラメータξlを導入して拡張した。オイラー法を用いていくつかの例を解いた経験から、積分軸変換パラメータGは1.0と取ることが最適であることがわかった。
　また本方法は微分方程式の初期値問題において、簡便な数値解法として知られているオイラー法における時間刻み幅の調節の指針を与える。本方法をルンゲクッタ４次法にフェールベルグとメルルッツィによる自動時間調節刻み法を加えたものと比較したところ、偶然にも、本方法とそれらには計算精度上の大きな差は見られなかった。また微分の計算回数で性能の評価をすると、本方法はそれらに比較してほぼ３倍程度高速であることがわかった。
　今後の課題として、SNMの時間刻み幅調節機能の妥当性に関して、より多くの例に適用して確かめる必要がある。また、xとyの関係がなだらかにつながっておらず、むしろ離散的な場合の簡便な解法の開発が必要である。また本方法を多変数の場合に拡張することも今後の課題としてあげられる。

　有益な議論と示唆をいただいた、北海道教育大学の小原　繁助教授、図書館情報大学の長谷川秀彦講師、電子技術総合研究所の関口智嗣主任研究官および建部修見研究官、物質工学工業技術研究所の三浦信明博士に感謝する。

参考文献

[ 1] 三池　秀敏、森　義仁、山口　智彦, 非平衡系の科学III, 講談社サイエンティフィック, 東京 (1994).
[ 2] P. Gray, K. Scott and K. Stephen, Chemical Oscillations and Instabilities, Clarendon Press, Oxford (1990).
[ 3] B. Jones, M. Banerjee, J. Jones, Computer Journal, 27, 184 (1984).
[ 4] R. H. Perry and C. H. Chilton, Chemical Engineers Handbook, Chapter 2, McGraw-Hill, New York (1973).
[ 5] Y. Mori, N. Kaminaga, N. Okazaki, U. Nagashima, and I. Hanazaki, J. Phys. Chem., (1998), submitted for publication.
[ 6] R. Bulisch and J. Stoer, Numer. Math., 8, 1 (1966).
[ 7] E. Fehlberg, Computing, 6, 61 (1970).
[ 8] P. Merluzzi and C. Brosilow, Computing, 20, 1 (1978).
[ 9] 今回性能評価に用いたプログラムを用いてSUN Ultra II( 177MHz)上での計算の実時間を測定すると、SNMは両者に比べて７倍から８倍速いが、これは、RKFおよびRKMはSNMに比較してプログラムが複雑になっているためであろうと考えられる。 プログラムのチューニングを十分すれば、この差は微分計算の回数の比である３倍程度に落ち着くものと期待できる。
Return