3
$\begingroup$

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)

enter image description here

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

enter image description here

Do you have any idea?

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ @davidhigh : If you follow the links you will find that this method is based on two 1995 papers where methods of claimed orders up to 5 were constructed for SDE. The question is if the methods are valid in some way. Apparently there is no strong convergence, but perhaps there is higher order convergence in some statistical measures. $\endgroup$ Commented May 30, 2020 at 10:00
  • $\begingroup$ @LutzLehmann So the comparison is valid? $\endgroup$ Commented May 30, 2020 at 10:22
  • 1
    $\begingroup$ For those who read it later journals.aps.org/pre/abstract/10.1103/PhysRevE.74.068701 $\endgroup$ Commented May 30, 2020 at 16:35

1 Answer 1

Reset to default
2
$\begingroup$

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.

plots of Euler approximations

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.

plot of numerical solutions with different step sizes

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}

enter image description here

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.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.