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]


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-29-71a84621a139> in <module>()
      3 # 1回小さくなりだすとどんどん小さくなる
      4 for i in range(2050, 2080):
----> 5     print y_check[i]

NameError: name 'y_check' is not defined

In [ ]: