オイラー法を使ったpython

Question

Q。

という問題で、オイラー法を使って解く問題なのですが、分からないのでご教授お願い致します。

import math
import matplotlib.pyplot as plt
import pandas as pd

r = 0.01
k = 30
n = 1
t = 2

while r > 0:
    r1 = r - n*t
    k1 = k - n*t
    print(r1,k1)

こんな感じで書いてみたんですが、同じ値しか出てきませんでした。

どうかよろしくお願いいたします。

Comments

[srsnsts]

解くべき微分方程式の記述がないようですが。

[ngm_programmer]

dn/dt = r(1-n/k)n
です。大変申し訳ございません。よろしくお願いいたします。

[srsnsts]

ご返答ありがとうございます。
微分方程式了解です。
回答を書くのでしばらくお待ち下さい。

[ngm_programmer]

すみません。
ありがとうございます。
自分では、何が何だかわからないので、差し支えなければ、詳細に御教授頂けますと幸いです。
よろしくお願いいたします。

[srsnsts]

ご返答ありがとうございます。
> 自分では、何が何だかわからないので、差し支えなければ、詳細に御教授頂けますと幸いです。
まぁ最初は誰でもそうですよ。
ちょっと待ってて下さい。

[ngm_programmer]

御親切にありがとうございます。
心折れかけていたので、、、。
よろしくお願いいたします。

[srsnsts]

> Q:ケージでネズミを飼っているとする。1匹のネズミは複数匹出産するので個体数は増えるが、個体数が増えるほ
> ど餌や生活空間の取り合いで成長率が個体数に比例して減少すると考えられる。このとき、ネズミの個体数の時間
> 変化を計算し、グラフに表せ（横軸を時間、縦軸を個体数とする）。簡単のため、ネズミは死なないし、雌雄の区別な
> く単体で継続して出産できるものとする。
> 解くべき微分方程式は次のようになる。r = 0.01 [/日], k = 30 [匹] とし、最初 1匹だったネズミの個体数を1日後
> までプロットせよ。

この部分を質問から削除したのはなぜですか？
あとから見た人がどんな質問なのかわからなくなるじゃないですか。
もとに戻していただけませんか？
このようなことをしたら、以後あなたの質問に対して回答しませんよ？

Answer 1

解くべき微分方程式は

dn/dt = r(1 - n/k)n

ですね。この式の左辺を離散化します。

微少時間 dt の間に、n が next_n に変化するとします。

すると微分方程式の左辺はオイラー法に従うと

(next_n - n)/dt と離散化できます。

つまり、

(next_n - n)/dt = r(1 - n/k)n

となります。この式を next_nについて解くと、

next_n = n + dt * r * (1 - n/k) * n

となります。

初期値として n = 1.0, r = 0.01, k = 30.0 を代入すると、next_n を計算することができます。

next_n を求めた後、n に 求めた next_n　の値を代入し、繰り返し next_n の値を求めればよいです。

これをpythonでコード化します。

依頼を受けてグラフを表示できるように修正しました。

import matplotlib.pyplot as plt

r = 0.01
k = 30.0
n = 1.0
dt = 1.0

t = 0
x = [t]
y = [n]
while t <= 1000:
	next_n = n + dt * r * ( 1.0 - n / k ) * n
	t = t + 1
	x.append(t)
	y.append(next_n)
	n = next_n

plt.plot(x,y)
plt.show()

このプログラムで計算すると、900秒ぐらいでnが30に収束することがわかります。
参考にしてみて下さい。

Comments

[ngm_programmer]

ご回答いただきまして、誠にありがとうございます。
とても勉強になりました。
オイラー法の当てはめ方と今回の計算方法を理解することができました。
ありがとうございます。

これをグラフにするにはどのようにしたらよろしいでしょうか。。。
申し訳ございません。よろしくお願いいたします。

[srsnsts]

> オイラー法の当てはめ方と今回の計算方法を理解することができました。
> ありがとうございます。
それはなにより。

> これをグラフにするにはどのようにしたらよろしいでしょうか。。。
> 申し訳ございません。よろしくお願いいたします。
あ〜僕はgnuplot使ったんだけど、
gnuplotって知ってる？
それともパイソニスタっぽいやり方のほうがいい？

[ngm_programmer]

ご連絡遅くなってしまい申し訳ございません。

現在使用しているものがGoogle Colaboになっております。。。

[srsnsts]

ブラウザで動くPython環境みたいな感じかな？

[srsnsts]

僕が提示したプログラムをPythonだけでグラフ表示できるようにしてみるね。