I am trying to compare the result of numerical integration of time independent Runge_Kutta, github page for stochastic differential equations with the analytical solution.
True answer match the rk1_ti_step
. But I get different results for rk2_ti_step
.
I used the same random normal array and pass it to the integration function to compare the analytical and numerical results.
from numpy.random import normal
np.random.seed(2)
def rk1_ti_step(x, t, h, q, fi, gi, dW):
a21 = 1.0
q1 = 1.0
x1 = x
w1 = dW * sqrt(q1 * q / h)
k1 = h * fi(x1) + h * gi(x1) * w1
xstar = x1 + a21 * k1
return xstar
def fi(x):
return alpha * x
def gi(x):
return beta * x
T = 1
N = 2**10
dt = T / N
Xzero = 1.0
alpha = 2.0
beta = 1.0
t = np.arange(0, T, dt)
dW = np.random.normal(loc=0.0, scale=1.0, size=N)
W = np.cumsum(np.sqrt(dt) * dW)
Xtrue = np.zeros(t.size)
Xtrue[0] = Xzero
Xtrue[1:] = Xzero * exp((alpha - 0.5 * beta ** 2) * t[1:] + beta * W[1:])
Xem = np.zeros(N)
Xem[0] = Xzero
for i in range(1, t.size):
Xem[i] = rk1_ti_step(Xem[i-1], t[i], dt, 1.0, fi, gi, dW[i])
pl.plot(t, Xtrue, color="k", lw=0.8, label="True")
plt.plot(t, Xem, lw=0.8, label="rk1_ti")
emerr = np.abs(Xem[-1] - Xtrue[-1])
print("error is : ", emerr)
def rk2_ti_step(x, t, h, q, fi, gi, dW):
a21 = 1.0
a31 = 0.5
a32 = 0.5
q1 = 2.0
q2 = 2.0
x1 = x
w1 = dW * sqrt(q1 * q / h)
k1 = h * fi(x1) + h * gi(x1) * w1
x2 = x1 + a21 * k1
w2 = dW * sqrt(q2 * q / h)
k2 = h * fi(x2) + h * gi(x2) * w2
xstar = x1 + a31 * k1 + a32 * k2
return xstar
Do you have any idea?