In [26]:
# ordinary differential equation
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
time = 5000
# analytical solution
x = np.arange(0,time, 1)
y_anal = [2] *time
#numerical solution
steps = [0.1, 0.01, 0.001, 0.0001]
for num, j in enumerate(steps):
value = 1
step = j
y_num = [0] * time
y_num[0] = value
for i in range(time - 1):
update_value = step * value * (2 - value) + value
y_num[i + 1] = update_value
value = update_value
fig = plt.figure(figsize=(10, 20)) # horizontal, vertical
plt.subplot(4, 1, num+1)
plt.plot(x, y_num)
plt.plot(x, y_anal)
plt.ylim(0,3)
In [31]:
# euler-maruyama method
# plot by variable steps
# average _ 100 times
import pylab
from random import normalvariate as norm
from random import seed
from math import sqrt
x = np.arange(0,time, 1)
steps = [0.1, 0.01, 0.001, 0.0001]
time = 10000
average_num = 100
for num, j in enumerate(steps):
step = j
y_all = []
for l in range(average_num):
y_num = [0] * (time)
value = 1
y_num[0] = value
for i in range(time - 1):
update_value = value + step * (2 - value) * value + sqrt(step) * value * norm(0, 1)
y_num[i + 1] = update_value
value = update_value
y_all.append(y_num)
y = np.array(y_all)
np.reshape(y, (average_num, time))
y_ave = [np.average(y[:, k]) for k in range(time) ]
fig = plt.figure(figsize=(10, 20)) # horizontal, vertical
plt.subplot(4, 1, num+1)
plt.plot(x, y_ave)
plt.ylim(0,3)
In [29]:
# step = 0.1の時に途中で消える問題
# 2000期目ぐらいで小さくなりだす
# 1回小さくなりだすとどんどん小さくなる
for i in range(2050, 2080):
print y_check[i]
In [ ]: