In [2]:
# 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
# Brownian paramator (p) = 1
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
p = 1
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 + p * 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 [4]:
# choose the paramator of Brownian
# Brownian paramator (p) = 0.1
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
p = 0.1
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 + p * 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 [5]:
# choose the paramator of Brownian
# Brownian paramator (p) = 0.5
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
p = 0.5
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 + p * 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 [6]:
# choose the paramator of Brownian
# Brownian paramator (p) = 1.5
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
p = 1.5
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 + p * 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 [8]:
# choose the paramator of Brownian
# Brownian paramator (p) = 1.75
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
p = 2
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 + p * 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 [10]:
# choose the paramator of Brownian
# Brownian paramator (p) = 2
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
p = 2
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 + p * 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 [ ]: