頭は文系なんですが・・・

質問

高校で数学をとっていなかったんですが（商業科なので）、大学の卒業論文でなぜだかテーマが「数値解析プログラム」になってしましました。情報系の学科に所属してはいますが、数学は大の苦手で、中学以来といってもいいくらいです。どなたか、
１．差分法の概説
２．拡散方程式
３．Ｌｏｔｋａ－Ｖｏｌｔｅｒｒａ方程式（生物）
について何でもいいので教えてください。お願いします。良い本など、何でもいいです！！

回答No.7

stomachman

さてと、もういっちょ行きますか。第２弾。
（やっぱり、高校の微分法の教科書か参考書は用意して、微分法だけでも読んで置いてくださいね。）
●記号がいちいちめんどくさいので、以後　x= n1*、y= n2*、z = t*　と書くことにします。またdx/dzというのは
dx
--
dz
という微分の記号のことですよ。これで書くと、無次元化したLotka-Volterra方程式(L-V方程式)は
dx/dz = x(1-y)
dy/dz = -ay(1-x)
ですね。
●この式を見て（中学生でも）分かるように、
＃　dx/dz =0となるのはx=0かy=1のときである。
＃　dy/dz =0となるのはx=1かy=0のときである。
しかし、絶滅しない限りx=0やy=0は起こらないので、
＃　「dx/dz =0となるのはy=1のときで、dy/dz =0となるのはx=1のとき。」と言える。
さらに、dy/dx = (dy/dz)/(dx/dz) （高校の教科書より）であるから、
＃　dy/dx = -ay(1-x)/(x(1-y))
同様にして、
＃　dx/dy = -x(1-y)/(ay(1-x))
である。従って、
＃　dy/dx = 0となるのはx=1かy=0のときである。
＃　dx/dy = 0となるのはx=0かy=1のときである。
ということが分かる。dy/dx = 0とは相平面の軌道がx軸と平行になる位置、dx/dy = 0とはy軸と平行になる位置である。

●保存則は
＃　S=(x-ln(x))+(y-ln(y))/a = 一定
証明：
Sをzで微分すると、
dS/dz = dx/dz - d(ln(x))/dz + (1/a)dy/dz - (1/a)d(ln(y))/dz
ここで　d(ln(x))/dz = (1/x) dx/dz　であるから、
dS/dz = dx/dz - (1/x) dx/dz + (1/a)dy/dz - (1/(ay))dy/dz
= (1- 1/x) dx/dz + (1/a)(1-1/y) dy/dz
ここにL-V方程式を代入すると、
dS/dz = (1- 1/x)x(1-y) + (1/a)(1-1/y)(-ay(1-x))
= (x-1)(1-y) - (y-1)(1-x) = 0
これで dS/dz = 0が証明できた。
　この保存則を満たすような曲線が、相平面上の軌道、すなわちL-V方程式の解である。（微分方程式の解というのは、数値ではなく、曲線あるいは関数です。）しかし、この保存則は当然の事ながらzを含んでいない。したがって、これだけでは時系列でx,yの変化を描くことはできない。

●x-ln(x) は(x>0)において、x=1のとき最小値1を取る。（グラフを描いてみると分かります。）従って
＃　S≧(1+1/a)である。
＃　x=1のときdy/dx = 0となる。このとき、y-ln(y)=a(S-1) である。 この方程式は２つの解 ymax,yminを持つ。グラフを描いてみると分かります。）
＃　y=1のときdx/dy = 0となる。このとき、x-ln(x)=S-1/a である。この方程式は２つの解xmax,xminを持つ。つまり、
＃　相平面上の軌道は、y=ymin, y=ymax, x=xmin, x=xmaxの４本の直線で作られる長方形に内接している。長方形との接点は(x,y) = (1,ymin), (1,ymax), (xmin,1),(xmax,1)の４点である。

●保存則を使った計算
dx/dz = x(1-y)
dy/dz = -ay(1-x)
を単純に差分化して計算すると、(x*, y*を次の時刻における値として）
Δx = x(1-y)Δz、x* = x+Δx
Δy = -ay(1-x)Δz、y* = y+Δy
であるが、通常、誤差が累積して渦巻きになってしまう。（多分、やってみると分かります。）そこで、初期値を決めたときに、あらかじめS=(x-ln(x))+(y-ln(y))/aを計算しておく。以後、このSをずっと使って、
(1) Δx = x(1-y)Δz、x* = x+Δx はそのまま使い、y* は、保存則を使って計算する。
(2) Δy = -ay(1-x)Δz、y* = y+Δy はそのまま使い、x* は、保存則を使って計算する。
という二つの方法を使うと良い。

まず、(2) の場合の計算法を示す。
step1) Δy = -ay(1-x)Δz を使って、次のyの値を計算し、これをy*とする。
step2) S=(x*-ln(x*))+(y*-ln(y*))/a
を満たすようなx*を求めたい。すなわち、C = (y*-ln(y*))/a - S　とおき、関数f(x)を
f(x) = (x-ln(x)) +C
と定義すると、f(x*) = 0 になるx*を求めたい。そこで、この方程式をNewton法で解く。まずf(x)をxで微分した物をf'(x)と書くと、
f'(x) = d(f(x))/dx = 1 - (1/x)
である。
Newton法step-1) Δx = x(1-y)Δzを使って、次のxの値を計算し、これをx0とする。
Newton法step-2) x1 = x0 - f(x0)/f'(x0), x2=x1-f(x1)/f'(x1),....を計算すると、数回～10回程度で収束する。これをx*とする。

(1)の場合も、同様のやり方でy*を計算できる。（導出できますか？）

＃　xが1に近いときには、(1)の方法を使い、yが1に近い時には(2)の方法を使うと良い。
なぜなら、xが1に近いとf'(x)がほとんど0になってしまうので、分母がほとんど0のわり算を行うことになり、f(x)/f'(x)がとんでもなく大きくなる。

この方法を使うと、渦巻きにならないで周期軌道が描ける。

●Newton法を利用して、ymax,ymin,xmax,xminを求めることができる。

●ここまで来れば、初期値をいろいろ変えたり（と言っても、どちらか一方を1にして、他方だけ動かせばよい。）aを変えたりして、色々な軌道を描けるようになりました。また軌道を描く前に、軌道に外接する長方形が分かるようになった訳です。

●Newton法の意味についてはまた別途質問してね。