Time independent Runge Kutta integration of SDE

Question

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?

Answer 1

The first method is Euler-Maruyama with its strong convergence of order $0.5$. This can be seen in the difference plot on the right. The step sizes vary with a factor of $4$, while the error does not decrease that rapidly.

In the second method, if one inserts the constants actually used, then the steps are \begin{align} k_1&=f(x_t)Δt+g(x_t)\sqrt2ΔW_t\\ k_2&=f(x_t+k_1) Δt+g(x_t+k_1)\sqrt2ΔW_t\\ x_{t+Δt}&=x_t+0.5(k_1+k_2) \end{align} Now insert the functions for the geometric Brownian motion, $f(x)=\alpha x$, $g(x)=\beta x$, to get \begin{align} k_1&=αx_tΔt+βx_t\sqrt2ΔW_t\\&=(αΔt+β\sqrt2ΔW_t)x_t\\ k_2&=α(x_t+k_1) Δt+β(x_t+k_1)\sqrt2ΔW_t\\ &=(αΔt+β\sqrt2ΔW_t)(1+αΔt+β\sqrt2ΔW_t)x_t\\ x_{t+Δt}&=x_t+(αΔt+β\sqrt2ΔW_t)x_t+0.5(αΔt+β\sqrt2ΔW_t)^2x_t \end{align} In the limit $Δt\to 0$ and with $(ΔW_t)^2\simΔt$ this step equation converges to $$ dX_t=(α+β^2)X_t\,dt+β\sqrt2X_t\,dW_t $$ This is still a geometric Brownian motion, but with rather changed parameters. So this method is not even consistent.

And indeed, changing the true solution to these parameters,

Xtrue = Xzero*np.exp(alpha*t+2**0.5*beta*W)

gives a picture of converging plots, albeit with the same convergence speed as the first method.

For a method that uses the second-order Heun method on the regular side to implement a derivative-free version of the Milstein method of strong order 1, see for instance https://stackoverflow.com/questions/44820604/have-i-implemented-milsteins-method and the source links there.

The improvement over Euler-Maruyama of the Milstein method is an additional term in a Taylor-like expansion, $$ x_{t+Δt}=x_t+f(x_t)Δt+g(x_t)ΔW_t+0.5\,g'(x_t)g(x_t)\,((ΔW_t)^2-Δt). $$ Now factorize $((ΔW_t)^2-Δt)=(ΔW_t-\sqrt{Δt})(ΔW_t+\sqrt{Δt})$ and in the second term $ΔW_t=0.5((ΔW_t-\sqrt{Δt})+(ΔW_t+\sqrt{Δt}))$. Note that terms $Δt^{3/2}$ or $ΔW_tΔt$ or $ΔW_t^3$ are effectively zero in this situation, so this formula can be approximated by absorbing the derivative as \begin{align} x_{t+Δt}&=x_t+0.5\Bigl[f(x_t)Δt+g(x_t)(ΔW_t-\sqrt{Δt})\Bigr]\\ &\qquad +0.5\Bigl[f(x_t)Δt+g\bigl(x_t+g(x_t)(ΔW_t-\sqrt{Δt})\bigr)(ΔW_t+\sqrt{Δt})\Bigr] \\ &=x_t+0.5(k_1+k_2)\\ k_1&=f(x_t)Δt+g(x_t)(ΔW_t-\sqrt{Δt})\\ k_2&=f(x_t+k_1)Δt+g(x_t+k_1)(ΔW_t+\sqrt{Δt}) \end{align}

Now with an appeal to belief one could almost say that the error here reduces as the step sizes by factors of 4, while in the other plots it only reduces mainly with a factor of 2. For more precise empirical observations one would need a more precise measure of the error and compute this over a larger sample of paths.