ロジスティック方程式の数値解法

ロジスティック方程式の数値解法

 

ネットで次のような問題を見つけた。

 

問題　次の微分方程式（ロジスティック方程式）

　　

を、t∈[0,1]の範囲で初期条件

　　

および

　　

のもとにオイラー法を適用して、次のことを調べよ。

（１）　t=20の値を解析解と比較せよ。

（２）　独立変数の差分間隔をΔt=1.0とした場合、どのような解が得られるか調べよ。

（３）　この方程式はどのような物理現象を表しているか説明せよ。

（４）　この方程式の差分方程式とカオスの関係について説明せよ。

 

問題の題意が多少曖昧であるが、面白そうなので、解いてみることにした。

 

オイラー法とは、微分方程式

　　

の解を

　　

と近似し、x₀=x(t₀)を出発点とし、逐次的に計算することで微分方程式の数値解を求める方法。

ここで、はにおける微分方程式の数値解である。

差分間隔と一定の場合、差分方程式は

　　

となる。

本問の場合、

　　

とおけば、

　　

という差分方程式が得られ、初期条件を計算の出発点に、上の漸化式を用いて、

と次々と計算すればよい。

 

微分方程式

　　

は、変数分離法を用いると、

　　

となる。

x₀=x(0)とすると、初期条件から

　　

となり、解は次のようになる。

 　　

したがって、t→∞のとき、初期値x₀に無関係に

　　

に収束する。

t=20のとき、e²⁰は非常に大きい数になるので、

　　

となることが容易に予想できるだろう。

 

Δt=1、x₀=0.1として、表計算ソフトを用いて計算した結果は次の通り。

なお、Δt=1のとき、

　　

両辺を１で引くと、

　　

となるので、Δt=1のとき、これを用いて計算してもよい。

計算結果をグラフにすると、次のようになる。

 

 

参考までに、４次のルンゲクッタを用いて解いた結果を次に示す。

 

Δt=1、x₀=2.0のときの計算結果は、次の通り。

この条件で、オイラー法を用いて、微分方程式の近似解を求めると、t=1以降のxの値がすべて０となり、t→∞のときの極限値が0になってしまい、解析解の挙動を捉えていないことがわかる。

このことは、x₀=2のとき、⑨より、x₀−1=2−1=1となり、

　　

となることからわかる。

また、x₀>2のとき、x₀−1>1となり、xは単調に減少し、

　　

となることも容易に理解できるだろう。

 

参考までに、同一の計算条件で、ルンゲ・クッタ法を用いて、解いた結果は次の通り。

 

 

（１）と（２）はこれでよしということにしよう。

 

（３）だが、ネムネコは物理屋さんじゃないので、この方程式がどのような物理現象を表しているかについては答えられない。

この件については、ddt³さんに丸投げする(^^ゞ

 

ロジスティック方程式は、マルサスの人口論にある

「人口は幾何級数的に増加する」

つまり、

　　

という人口増加のモデルに対する反論、または、その修正として作られたモデル。

マルサスのモデルだと、

　　

となり、t→∞のとき、x→∞になってしまう。

こんな解は非現実的だということで、人口、生物の個体数の増加を抑制する、

　　

という項を加え、

　　

としたものをロジスティック方程式という。

ここで、rは内的自然増加率、Kは環境収容力などと呼ばれる。

本問の場合は、r=1、K=1である。

この微分方程式を解くと、

　　

となり、

　　

なぜ、Kを環境収容力と呼ぶのか、理解してもらえるのではないだろうか。

 

（４）　いったい、何をもってカオスと呼ぶのか不明なので、何とも答えようがないのだけれど、たとえば、x₀=0.1、Δt=2とすると、次のような振動解が得られるってことをいいたいのか。Δt=2のときは振動しながらも1に収束するが、2<Δt≦3のときは収束せずにΔtが増加するとカオス的な状況を呈するようになる。

それとも、Δx=1のとき、オイラー法では、t→∞の時に得られる定常解が1、0、−∞と初期値x₀の値によって変わるということを言いたいのだろうか。

はたまた、時間進展とともに、微小な擾乱（この場合、オイラー法の打切誤差）が時間進展とともに増大し（打切誤差の蓄積のため）、積もり積もって、破局・カタストロフィーをもたらすというバタフライ効果を言いたいのか(・・?

 

 

続き

 

差分間隔Δt=0.1、初期条件x₀=0.1、x₀=2.0で、Euler法を用いて解いた結果を下に示す。

この計算に使用したプログラムは、次の通り。

 

parameter(n=200)
real x(0:n), y(0:n), y1(0:n)

x=0.; y=0.
y(0)=0.1; y1(0)=y(0)
dx=0.1

c = y(0)/(y(0)-1)
do i=1,n
    x(i)=x(0)+dx*i
end do

call Euler(x,y,dx,n)
call Runge_Kutta(x,y1,dx,n)
do i=0,n
    write(*,100) x(i),y(i),y1(i),exact(x(i),c)
end do

100 format(f8.4,1x,2(f8.5,1x), f8.5)
end

function f(x,y)
f = y-y*y
end

subroutine Euler(x,y,dx,n)
real x(0:n), y(0:n)

do i=1, n
    y(i) = y(i-1)+dx*f(x(i-1),y(i-1))
end do
end

subroutine Runge_Kutta(x,y,dx,n)
real x(0:n), y(0:n)

do i=1,n
    dk1=dx*f(x(i-1),y(i-1))
    dk2=dx*f(x(i-1)+dx/2.,y(i-1)+dk1/2.0)
    dk3=dx*f(x(i-1)+dx/2.,y(i-1)+dk2/2.0)
    dk4=dx*f(x(i),y(i-1)+dk3)
    y(i)=y(i-1)+(dk1+2*dk2+2*dk3+dk4)/6.0
end do
end

function exact(x,c)
exact=c*exp(x)/(c*exp(x)-1)
end

コメント 2

　ddt^3です。

　ロジスティクス方程式の物理的意味ですか。ちょっと脚色しますね(^^)。


　まずはマルサスモデル。

　ある世代tの個体数をx(t)とした時、その半分が女の子だとして、女の子は全員残りの雄どもとつがいになって、子供をみんなα人産むと仮定します。そうするとある世代の人口増加は、x/2×α＝αx/2です。
　一方、その世代tに老衰して死ぬジジババは（だんだん自分もその範疇に・・・(^^;)）、死亡率をβとしてβxだとします。

　そうするとdx/dt＝(α/2－β)x。もし富国強兵で「産めよ殖やせよ」が国策なら、α/2－β＞0となり、人口は幾何級数的に増えるはず(^^)。

　でも上記は、一夫一婦制が前提ですよね？。人間って、やりたいと思えば何でもやる生物だから乱交かも知れないじゃない。乱交の場合、雄x/2がランダムウォークして、女の子x/2がランダムウォークして出会う確率は？。
　じつは化学反応論の標準則（経験則）として、その頻度を与える法則がある。それは、2つの化学物質の反応速度は、それぞれの濃度の積だ。x/2×x/2＝x^2/4．

　という訳で乱交の場合、dx/dt＝α/4･x^2－βx。幾何級数以上に人口は増える！。やってらんねぇ～ぜ。だから道徳や婚姻制度があるのかしら？(^^;)。
　一夫多妻も多夫一妻という地域もあるけれど、とにかく乱交を認める社会は地球上にはない。そうするとx/2×x/2のどちらかが固定人数になるので（そういう羨ましい人達は小数なのだ！(^^;)）、けっきょくx(t)の挙動は、一夫一婦とはそんなに違わないものになりそうだ。

　その前提でマルサスモデルを考えると、やっぱりdx/dt＝γx。


　・・・しかし幾何級数で人口が増えても問題は山積みなのだ。地球には有限の資源しかない。少なくとも利用可能な資源はとても有限だ。そうすると幾何級数の爆発的速度で増える人間どもの間には、資源を求める争いが起こる。

　いや、その前に殺人が起こる。人間は、邪魔だと思えば何でも破壊するし殺す生物である。だから資源云々以前に、人口が増えれば殺人は激増する。化学反応論の標準則に従おう。人類xがランダムウォークして、人類xがランダムウォークして出会う確率は、x×x＝x^2だ。それに比例して殺人が起こる。

　しかし人類は協力し合える生物でもあるはずだ。それはその通りだ。でもそれはあくまで身内に限られる。他宗教を信じてる奴らなんて人間じゃねぇ～！、って考えは確かにある。そこに資源を求める競合が重なる。

　資源の競合度合いもx^2に比例するはずだ。そして資源の競合に負けた方は、資源を簒奪された方は、餓死するとかそんな選択肢しか残らない。それは間接的な殺人だ。
　他宗派，他民族は人間じゃないから、勝った方は絶対に敗者を助けない。従って人口は、x^2に比例する形で減少もする。これが環境収容力なのだと思う。

　よってロジスティック方程式は、平均的には恐らく正しい気がする。平均とは、エルゴード定理の意味において。

by ddtddtddt (2018-05-10 18:44) 

こんばんは。

ランダムウォークですか(^^)

何でも、このモデルは、流体力学の「流体からうける抵抗は、速度の２乗に比例する」という経験則と関係があるらしいですね。
重力gを受けて落下する場合、運動方程式は
　m(dv/dt)=mg - kv²
となる。
――速度Vが小さい時は、抵抗は−vに比例するけれど、速度が大きくなると、−v²に比例する――
どうやら、これがモデルのヒントになったらしい。

アナロジー的な発想というやつですよね。
なぜ、そのアナロジーが成立するかはわからない、メカニズムはわからないけれど、現象をよく表すモデルならば・・・。


by nemurineko (2018-05-10 19:27)