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

Question

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

stomachman · Accepted Answer

さてさて、第３弾です。

Lotka-Volterra(L-V)方程式の軌道を描くだけなら、もし
f(x) = x - ln(x) の逆関数、あるいは、
F(x) = exp(x) / x の逆関数（F(x) = exp(f(x))だから、Fの逆関数が分かればfの逆関数も分かる。）
を求めることができれば、わざわざ差分方程式を使わなくても、保存則Sだけで計算できます。（どうしてか分かりますよね？）また、差分方程式を安定に計算するのにも、Sを使うから、この逆関数の計算が必要です。

　既に、よく使われるニュートン法による計算法を示しましたが、xの値が1に近いときには計算が不安定になるという性質を持っていて、二通りの計算法を切り替える必要がありました。ここでは「もう少し遅いけれど、収束が安定している方法」を紹介します。これだと、切り替えが（多分）不要です。
　実は、面白いことに、最近同じ問題「F(x) = exp(x) / x の逆関数」が質問としてUPされているんです。(http://oshiete1.goo.ne.jp/kotaeru.php3?q=17648)
　以下はその質問の回答として移動中の時間を使って書いたんだけど、UPするまえに質問者がCLOSEしてしまって宙に浮いていたもの。（いやー、こんなに早く日のめを見るとは思わなかった。）
以下の(2)(3)の漸化式を使って計算すれば、F(x)の逆関数x(F)の値が出せます。L-V方程式ではx>0かつF(x)>e についてだけFの逆関数が求められれば十分です。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
(1) まず、F(x)のグラフを描いてみて下さい。表計算ソフト(excelとか)を使うと簡単です。
F(x)=exp(x)/xはx<0 では負, x>0 では正です。x->+0では∞、x->-0では-∞、x->∞では∞、x->-∞では0であることが分かります。F(x)の導関数（xによる微分）は
F'(x) = ((x-1)/(x^2))exp(x)
だから、F(x)は x=1 だけで極値（極小）を持つことが分かります。この極小値はF(1) = e　（e=exp(1)です）。F(x)はx>0では下に凸の関数です。これでF(x)のグラフが描けたことでしょう。
　明らかに、F(x)の逆関数 x(F)はF<0およびF>=eの範囲でだけ定義され、e>F>=0の範囲には値がありません。 x(F)はF<0では１価関数、F>=e では２価になります。一方の解はx>=1, もう一方は1>=x>0の範囲にある(x=1の時のみ二つの解が一致）という訳です。

(2) 逆関数は簡単な数式で表せません。まずx>1の方の解だけ考えることにしましょう。
（以下ではFの代わりにyと書いています。）y>eの時の、y = exp(x)/xの逆関数x=G(y) (ただしx>1となる枝）は、以下に示す漸化式x[0], x[1],...の極限です。（ln()は自然対数、exp()は指数関数です。）
x[0]= ln(y)
x[1]= ln(y) + ln(x[0])
:
x[n+1]= ln(y) + ln(x[n])
:
強いて一つの式にまとめると
x =G(y)= ln(y)+ln(ln(y)+ln(ln(y)+ln(ln(y)+ ......
と表せる。これが一つの答です。（他にも表し方はある。）

(3) では0<x<1におけるx(y)の漸化式は?
x[0]=1/y
x[n+1]=exp[n]/y
です。まとめるとどうなるか、はご自分で。

(4) すぐに信用しちゃいけません。(2),(3)をFに代入してF(G(y))=y になる事を確かめてみてください。またG(F(x))=xもやってみてください(面白いですよ）。数値的には、表計算ソフトでも(1)(2)のやり方でx[0], x[1],...が簡単に計算できますから確かめてみてください。（なお、もっと速い計算法としては「ニュートン法」があります。）

(5) では(2)の導出のあらすじを書いときましょう。（(3)はそのバリエーションでできます。）
x>=１という条件で考えます。y = exp(x)/x をln(y)=x-ln(x)と変形しておきます。
x>1ではxもln(x)も単調増加なのに、ln(y)=x-ln(x)が単調増加である。だから、ln(x)の項よりxの項の方が優位である(増え方が速い)ことが分かります。そこでとりえあえずln(x)の項を無視した第０次近似として
x[0]= ln(y)
を考えます。誤差（正解との差）は
x-x[0] = (ln(y)+ln(x))-x[0] = (ln(y)+ln(x))-ln(y)=ln(x)
だけありますから、x[0]にあとln(x)を加えれば正解になるわけです。しかしあいにくxが分からないので、代わりにx[0]を使って誤魔化します。これで第１次近似ができました。
x[1] = ln(y) + ln(x[0])
もちろん「x[0]よりx[1]の方が誤差の絶対値が小さくなる」という事を証明しなくちゃいけませんが、めんどくさいのでパス。さて、ここから数学的帰納法を使います。
n(>0)のときに
x[n]= ln(y) + ln(x[n-1])
であるとすると、誤差は
x-x[n] = (ln(y)+ln(x))-x[n] =ln(x)-ln(x[n-1])
従って、ln(x)-ln(x[n-1])をx[n]に加えたいのですが、ln(x)の代わりにln(x[n])で誤魔化して第n+1次近似を作ります。
x[n+1]= x[n]+ln(x[n])-ln(x[n-1]) = ln(y)+ ln(x[n])
（ここでもホントはx[n+1]の方がx[n]よりも正解のxに近いという事を証明しなくちゃいけません。）
　このように、誤差をどんどん修正していくように漸化式を構成する訳です。修正のやり方を変えればいろんな漸化式が得られ、ニュートン法もそのひとつです。

stomachman · Answer

さてと、もういっちょ行きますか。第２弾。
（やっぱり、高校の微分法の教科書か参考書は用意して、微分法だけでも読んで置いてくださいね。）
●記号がいちいちめんどくさいので、以後　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法の意味についてはまた別途質問してね。

stomachman · Answer

と、言うわけで、第一弾。（分からないところは多分TCMさんが教えてくれますよね。）

●Lotka-Volterraのモデル
アドリア海の魚の発生数の変動を説明するために考案されたモデル。現象を非常に単純化しているが、個体数変化曲線の複雑な挙動をある程度説明する。

餌になる魚1と、それを捕食する魚2の２種類を考える。個体数をn1, n2とする。
d(n1)/dt  = (単位時間当たりの魚1の増殖）- (単位時間当たりの魚1の死滅数）
である。
(単位時間当たりの魚1の増殖）= (α1) n1  
とモデル化する。これだけだと、魚1は指数関数的に増えることになる。（n1 = A exp(α1 n1)が解だから。）(α1)は増殖率。
(単位時間当たりの魚1の死滅数）= αn1 n2
とする。つまり、魚１と魚２の数の積が、それらが出会って捕食が行われる頻度を近似すると考える。魚１が寿命で死ぬのは考慮しない。
実際には、豊富な餌があって、ある程度以下のn1であれば飢え死にする心配はなく、また喰われずに寿命を全うできるほど生き延びることもない、ということだと思えばよい。
次に、魚２
d(n2)/dt = (単位時間当たりの魚2の増殖）- (単位時間当たりの魚2の死滅数）
魚２は魚１がいないと暮らしていけない。喰った内の一部の物質なりエネルギーなりが増殖に使われるのだから、どれだけ魚１を喰ったか、によって増殖数が決まる。従って、魚2の増殖は(単位時間当たりの魚1の死滅数）に比例していると考えられる。
 (単位時間当たりの魚2の増殖）=βn1 n2
魚２は一定の比率で死ぬ。個体数(population)当たりの死亡率を2(κ2)とすると
 (単位時間当たりの魚2の死滅数）= 2(κ2)n2

以上から、
d(n1)/dt = (α１)(n1) - α(n1)(n2)
d(n2)/dt =β(n1)(n2) - 2(κ2)(n2)
これがLotka-Volterra方程式。(α１),α,β,(κ2)という４個の定数、t,(n1),(n2)という３個の変数を含んでいる連立常微分方程式である。

●この方程式を無次元化（計測する単位に依らない形に）する。
新しい変数(n1*)、(n2*)、(t*)と、新しい定数aを使って
(n1) = 2(κ2)(n1*) /β
(n2) = (α1)(n2*)/α
t = (t*)/(α1)
とおくと、無次元化したLotka-Volterra方程式は
d(n1*)/d(t*) = (n1*) - (n1*)(n2*)
d(n2*)/d(t*) =-a {(n2*) - (n1*)(n2*)}
である。

●この方程式において、保存則
((n1*)-ln(n1*))+((n2*)-ln(n2*))/a = 一定
が成り立つ。
●この方程式の定常状態 （d(n1*)/d(t*) もd(n2*)/d(t*)もゼロ）は
(n1*) = 1
(n2*) = 1
である。個体数がこの状態にあれば、以後、個体数の変化はない。
　しかしこの定常状態が安定かどうかは別問題。つまり、偶然に個体数に僅かな変動が生じたとするとき、自然に定常状態に戻るのかどうか、が最も興味ある問題である。

●n1*とn2*をそれぞれ横軸、縦軸にとったグラフを描く。これを「n1*-n2*相平面上の軌道」という。
d(n1*)/d(t*) = (n1*) - (n1*)(n2*)
d(n2*)/d(t*) =-a {(n2*) - (n1*)(n2*)}
を差分方程式にする、つまり時刻k(Δt*)におけるn1*を(n1*)[k]、n2*も同様に表し、
d(n1*)を(Δn1*)[k], d(n2*)を(Δn2*)[k], d(t*)を(Δt*)で置き換えてやると、

(Δn1*)[k] = (Δt*){(n1*)[k] - (n1*)[k] (n2*)[k]}
(Δn2*)[k] = -a(Δt*) {(n2*)[k] - (n1*)[k] (n2*)[k]}
(n1*)[k+1]  = (n1*)[k]+(Δn1*)[k]
(n2*)[k+1]  = (n2*)[k]+(Δn2*)[k]

となる。これをk=0,1,2,....と計算すればよい。
パラメータaはシステム（系）の条件を決めるものである。これは計算の初めに決めて置かねばならない。
(n1*)[0]と(n2*)[0]は初期条件あるいは初期値と言い、適当に与えなくてはならない。
Δt*は十分小さくとるが、いろいろ試して適当な大きさを決めなくてはならない。
また、数値の精度を変えてみる（多倍長計算。doubleを使うかsingleをつかうか。FORTRANならREAL*4かREAL*8か）ことによって、計算の結果があまり変わらないようにしないと、何をやっているのか分からんことになる。
この計算自体はEXCEL位で十分である。せいぜいVBで良い。
計算を続けて、どんどん相平面上の軌道を描いていくと、パラメータと初期条件が適当なら、閉軌道になるであろう。これが周期解である。
●もうひとつ経時変化のグラフを描く。こんどは横軸がt*で縦軸が(n1*)、(n2*)である。個体数が時間と共に変わる様子が描かれる。
●周期解が得られたら、わざと(n1*)か(n2*)を少し変化させてみる。（パラメータは同じにして、初期値を少し変えればよい。）それでもほとんど同じ周期軌道を描くなら、この周期解は「安定」である。さもなくば不安定である。

●この方程式の軌道は安定であるが、漸近安定(asymptotic orbital stability)ではない。つまり、二つの（相平面上で）近接した軌道上の、互いに近接した2点を初期値としてそれぞれ計算をしてみると、経時変化のグラフではどんどん波の山の位置がずれていってしまう。

●系のパラメータが変化すると、軌道が大きく変化することがある。たとえば、魚１の餌が不足して、増殖率(α1)が低下した、とするとどうなるか。

まあ、こういった所がシミュレーションである。
さて、次の数学的ステップは、ポアンカレ-ベンディクソンの定理だとか、リヤプノフ指数だとかを使って、系の安定性を数学的に解析することなんだけど.....どうしようか。というのも、これを追っていっても、遙か昔に確立された理論をなぞっていくだけのことになってしまう。大変な努力がいる割に、分かっている人には答は自明である。

むしろ、卒論としては、一応研究なんだから、 生物学が得意なら、系をもうちょっとリアルにしてみるというようないたずらの方が面白いだろう。

●大昔は漠然と、微分方程式を数値計算のために近似するのが差分方程式（いわばまがいもの）だ、と思われていた。ところが実は、差分方程式は、微分方程式では説明できないようなモデルを表現するということが認識されたのも、生物学が最初だったように思う。
　繁殖が年１回しか起こらないような動物の捕食関係においては、差分方程式
(Δn1) = (α１)(n1)[k] - α(n1)[k] (n2)[k]
(Δn2) =β(n1)[k](n2)[k] - 2(κ2)(n2)[k]
でないと、本当の姿が見えてこない。Δt=１年である。微分方程式で表した場合には、n1,n2の個体数の変化が他方の個体数に直ちに反映するが、１年単位の繁殖を考えると、大繁殖や絶滅が簡単に起こることが分かる。

●(単位時間当たりの増殖）、 (単位時間当たり死滅）を単純に個体数の積で表していていいのか？ここをもっと精密化したらどんな現象が見えてくるだろうか？確率の知識がないと旨いモデルが作れないかもしれないが、面白いバリエーションである。

●拡散方程式がどう繋がるのか、というのは多分、個体数が空間的にむらになったとき、そのむらが解消していくプロセスを表現するモデル、ということだと思われる。
この発想で行くと、実空間の中に個体がいっぱいいて、その場所場所で捕食が起こり、個体数にむらが発生する、というような、シミュレーションの話（つまらないゲームのようなものだ。）になる。こっちをやれと言われているのだろうか？

●やっぱり、もう少し、卒論のテーマの詳細をUPして貰わないと....

stomachman · Answer

補足を見ました。
さぼってた、なんて言ってごめんね。

　近頃、大学は学生の学力低下を問題視していますが、教える方のレベルもpoorなのか...でも勉強は自分でするもんです。

こうなったら、「教えてgoo」をとことん利用しましょう。
「何でもいいので教えて！」と言われたら、その道の大家でもないかぎりビビリますよ。もそっと、卒論の中身を具体的にUPして下されば、わらわらと回答が集まるやも知れません。

....講義になんか全然出なかったstomachmanでした。

TCM · Answer

いやあ皆さん、手厳しい(>.<)　とりあえず、botanさんの努力を応援してあげましょうよ。ちょっと説明しますが、目をつぶっておいてくださいね。まずは拡散方程式から。

２．拡散方程式
　拡散方程式は熱伝導方程式ともいい、定常状態の場合、
　　∂φ/∂t＝p∇2φ
という形をした偏微分方程式です。針金の熱伝導とかの解析で説明されることもありますね。∇2の２は、ここでは表示できませんが、本当は上付きで空間の２階偏微分演算子です。φが求めたい変数で、pは拡散の程度を表す係数です。空間をx軸の１次元とすると、
　　∂φ/∂t＝p(∂2φ/∂x2)
となります。やはり、ここでも２は本当は上付きです。他にも説明したいことはありますが、とりあえずここまで。

１．差分法
　差分法というのは、微分方程式の微分を差分に置換することによって元の微分方程式を解く手法です。上式ならば、まず図１のように間隔d（デルタ）ｘの格子（１次元なので格子とは言えませんが）を考え、節点iの時刻ndtでの値φiをφ(i,n)とすると、
　{φ(i,n+1)-φ(i,n)}/dt＝p{φ(i+1,n)－2φ(i,n)＋φ(i-1,n)}/dx^2
　　　　　　　　　　　　　　　　　　　　d:デルタ
と離散化します。そして、全ての節点についてこの方程式を作り、系全体の連立１次方程式を組み立て、節点１と節点ｍの値や時刻ゼロでの各節点の値など境界条件を与えて解きます。そして時刻をdtだけ進めながら、所望の解を得るまで同じ手順を繰り返します。

　１－－２－－３－－４･･･ｉ－－i+1･･･m-1－－ｍ
　｜dｘ｜
　　　　　　図１．１次元の差分格子

　しかし実際には、差分スキームの安定性、適合性、収束性などの検討をしたり、離散化を工夫したり、解の妥当性の検証など気をつけなければいけないことが多々あります。偏微分や行列がわからないとかなり苦しい戦いになるでしょう。ひょっとすると先輩の残したプログラムがあって、それを改良したり検証を行ったりという研究テーマなのかもしれませんね。で、それはLotka-Volterraモデルを使った生態系シミュレーションの研究だったりするのでしょう。botanさんにはがんばってもらうしかないですが、確かに数学がネックになりそうですね。
　差分法の参考書は多々ありますが、「数値計算プログラミング」森正武（岩波書店）が比較的簡単かもしれません。ただ言語がFORTRANなので、まずは大きな書店で自分に合った本を探すことをお勧めします。

追記：Lotka-Volterraモデルについてはよく知りませんが、方程式は拡散方程式に付加項がついた形をしていますね。

stomachman · Answer

deeppurple14さんてば、人聞きの悪い(^o^)
参考書だけ教えてパニックさせるんじゃ可哀想だってんで、凄くサービスしたつもりなんですけど。
　実際、Lotka-Volterra方程式のキモは何か、なんて、stomachmanだってまだ勉強中てな所ですから。

deeppurple14 · Answer

stomachmanさんがかなり冷たく突き放しておられますが、私も同じような考えです。
今のテーマをなんとかしようと考えるより、テーマを変えてもらった方が良いと思います。

数値解析のお薦めの本はいろいろ知っています。でも、それらを読んで理解するには、基礎となる数学（特に解析学）をかなり高いレベルでマスターしていないと厳しいです。しかも、差分の方は比較的初歩的ですが、拡散方程式はいろいろありますが簡単に理解できるものではないです。
しかも大学の解析学を勉強するには高校の微積は完璧に理解していないと無理です。

結論として、これは不可能です。もし抜け出せないのなら、stomachmanさんの言うようにごまかしで切り抜けるというのはできるかもしれない。

ぜったいテーマ変えた方が良いと思いますが。

stomachman · Answer

困りましたねえ。もうちょっと質問者の置かれた状況に関する補足をお願いします。
　本来なら文系の頭ならではのテーマを自分で選ぶべきでした。
　いつまでに仕上げるのか、が一番気になります。来春卒業予定というのなら、今までさぼってたんだからしょうがない。大急ぎで、もっとずっと易しい問題に切り替えるべきです。もし本当に「数値解析プログラム」というだけなら、何も積分方程式なんか相手にしなくても、数学自体は中学レベルでも理解できてしかも実用性もある、というテーマだってあるんじゃないでしょうか。
　やっぱり微分方程式・積分方程式の数値解析でなきゃだめだ、という事だとすると...
「数学が苦手」と諦めているようでは、いきなり専門書を読んでも５ページ目までに挫折するのは目に見えています。まずは、微分法・積分法の一番基礎の概念だけでも把握しないとダメでしょう。「文章題」：つまり現象のモデルを微分方程式や積分方程式で表す。逆に数式を文章に戻す、という技術は最低限必要です。偏微分の概念が必要ですね。高校数学でやる紙の上での微積分の計算の練習はこの際省くしかないでしょうが、結局、高校数学の教科書が一番早道かも知れません。
　どうしても今すぐに、というのなら、最後の手段として：どこかで出来合いのプログラムを探してきて、適当な境界条件を入れていろいろ計算してみました。幾つかのプログラムの性能を比較してみました。数値計算の誤差を検討するために、同じ計算を多倍長計算の精度を変えて行って、どこまで一致するかを確認しました。誤差のグラフを描いてみました。で勘弁して貰うしかなさそう。これだけでも徹夜の連続になるでしょう。（でもこの程度で「情報系で数値解析プログラムを専攻しました」などとは絶対に言って欲しくないです。計算の中身が理解できていない以上、内容的には小学生の夏休みの自由研究と幾らも変わらないレベルですから。）

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

さてさて、第３弾です。

さてと、もういっちょ行きますか。

と、言うわけで、第一弾。

補足を見ました。

いやあ皆さん、手厳しい(>.<) とりあえず、botanさんの努力を応援してあげましょうよ。

deeppurple14さんてば、人聞きの悪い(^o^)

この回答への補足

stomachmanさんがかなり冷たく突き放しておられますが、私も同じような考えです。

困りましたねえ。

この回答への補足

似たような質問が見つかりました

このQ&Aを見た人はこんなQ&Aも見ています

関連するカテゴリからQ&Aを探す

このQ&Aを見た人がよく見るQ&A

デイリーランキングこのカテゴリの人気デイリーQ&Aランキング

マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング

　いやあ皆さん、手厳しい(>.<)　とりあえず、botanさんの努力を応援してあげましょうよ。