Has anyone used Kasdin's Runge-Kutta SDE algorithm?

General Usage
1

Dear all, I am doing SDE now. I know Julia has fancy differentialequations.jl. It has several SDE solver like SRIW1 or something like that, which uses adaptive Runge-Kutta (RK).

I also found there is a fortran version of fix step RK method for SDE, by John Burkardt,

https://people.math.sc.edu/Burkardt/f_src/stochastic_rk/stochastic_rk.html

It is based on Kasdin’s method,

  1. Jeremy Kasdin,
    Runge-Kutta algorithm for the numerical integration of stochastic differential equations,
    Journal of Guidance, Control, and Dynamics,
    Volume 18, Number 1, January-February 1995, pages 114-120.
  2. Jeremy Kasdin,
    Discrete Simulation of Colored Noise and Stochastic Processes and 1/f^a Power Law Noise Generation,
    Proceedings of the IEEE,
    Volume 83, Number 5, 1995, pages 802-827.

Eqns is simple,

image
image1061×870 17 KB

Although John’s version is for just for one SDE which means it is a scaler version, it is trivial to extend it to make the vector version to solve SDEs. For example, a 2D version (so 2 SDEs) can be like this in Fortran,

subroutine rk4_ti_step_vec_mod ( x, t, h, q, fi, gi, n, xstar )
use random
implicit none
interface				   !! Declare function        
    function fi ( x ) result (f)
    use constants
    real ( kind = r8 ) :: f(2)
    real ( kind = r8 ) :: x(2)
    end function
    function gi ( x ) result (g)
    use constants
    real ( kind = r8 ) :: g(2)
    real ( kind = r8 ) :: x(2)
    end function      
end interface
integer ( kind = i8 ) n
real ( kind = r8 ) h
real ( kind = r8 ) k1(n)
real ( kind = r8 ) k2(n)
real ( kind = r8 ) k3(n)
real ( kind = r8 ) k4(n)
real ( kind = r8 ) q

real ( kind = r8 ) t
real ( kind = r8 ) t1
real ( kind = r8 ) t2
real ( kind = r8 ) t3
real ( kind = r8 ) t4
real ( kind = r8 ) w1
real ( kind = r8 ) w2
real ( kind = r8 ) w3
real ( kind = r8 ) w4
real ( kind = r8 ) x(n)
real ( kind = r8 ) x1(n)
real ( kind = r8 ) x2(n)
real ( kind = r8 ) x3(n)
real ( kind = r8 ) x4(n)
real ( kind = r8 ) xstar(n)

real ( kind = r8 ), parameter :: a21 = 2.71644396264860_r8 &
,a31 = - 6.95653259006152_r8 &
,a32 =   0.78313689457981_r8 &
,a41 =   0.0_r8 &
,a42 =   0.48257353309214_r8 &
,a43 =   0.26171080165848_r8 &
,a51 =   0.47012396888046_r8 &
,a52 =   0.36597075368373_r8 &
,a53 =   0.08906615686702_r8 &
,a54 =   0.07483912056879_r8

real ( kind = r8 ), parameter, dimension(4) :: qs = [ 2.12709852335625_r8 &
    ,2.73245878238737_r8 &
    ,11.22760917474960_r8 &
    ,13.36199560336697_r8 ]

real ( kind = r8 ) :: qoh
real ( kind = r8 ) :: warray(n,4)
integer (kind = i8) :: i

qoh = q / h

do i =1,4
    warray(:,i) = gaussian(n)*sqrt(qs(i)*qoh)
enddo

!t1 = t
x1 = x
k1 = h * ( fi ( x1 ) + gi ( x1 ) * warray(:,1) ) 

!t2 = t1 + a21 * h
x2 = x1 + a21 * k1
k2 = h * ( fi ( x2 ) + gi ( x2 ) * warray(:,2) )

!t3 = t1 + ( a31  + a32 )* h
x3 = x1 + a31 * k1 + a32 * k2
k3 = h * ( fi ( x3 ) + gi ( x3 ) * warray(:,3) )

!t4 = t1 + ( a41  + a42 + a43 ) * h
x4 = x1 + a41 * k1 + a42 * k2 + a43 * k3 ! the a43 * k3 seems need to be added.
k4 = h * ( fi ( x4 ) + gi ( x4 ) * warray(:,4) )

xstar = x1 + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4

return
end subroutine rk4_ti_step_vec_mod

An improved version can be like below ( mostly suggested by veryreverie from stackoverflow,
algorithm - Is there room to further optimize the stochastic_rk Fortran 90 code? - Stack Overflow),

subroutine rk4_ti_step_vec_mod1 ( x, t, h, q, fi_vec, gi_vec, n, xstar )
    use random
    implicit none
    interface				   !! Declare function        
        function fi_vec ( x ) result (f)
        use constants
        real ( kind = r8 ) :: f(2)
        real ( kind = r8 ) :: x(2)
        end function
        function gi_vec ( x ) result (g)
        use constants
        real ( kind = r8 ) :: g(2)
        real ( kind = r8 ) :: x(2)
        end function      
    end interface       
    integer(kind = i8), intent(in) :: n
    real(kind=r8), intent(in) :: x(n)
    real(kind = r8), intent(in) :: t  
    real(kind = r8), intent(in) :: h
    real(kind = r8), intent(in) :: q
    !real(kind = r8), external :: fi_vec
    !real(kind = r8), external :: gi_vec
    real(kind = r8), intent(out) :: xstar(n)
    real ( kind = r8 ), parameter :: alphas(4) = [ 0.47012396888046_r8 &
    ,0.36597075368373_r8 &
    ,0.08906615686702_r8 &
    ,0.07483912056879_r8 ]
    real(kind = r8), parameter :: as(4,4) = reshape([ &
        &  0.0_r8,              0.0_r8,              0.0_r8,              0.0_r8, &
        &  2.71644396264860_r8, 0.0_r8,              0.0_r8,              0.0_r8, &
        & -6.95653259006152_r8, 0.78313689457981_r8, 0.0_r8,              0.0_r8, &
        &  0.0_r8,              0.48257353309214_r8, 0.26171080165848_r8, 0.0_r8 ], [4,4])
    real(kind = r8), parameter :: qs(4) = [  2.12709852335625_r8, &
        &  2.73245878238737_r8, &
        & 11.22760917474960_r8, &
        & 13.36199560336697_r8 ]
    real(kind = r8) :: ks(n,4)
    real(kind = r8) :: ts(n,4)
    real(kind = r8) :: ws(n,4)
    real(kind = r8) :: xs(n,4)
    real(kind = r8) :: warray(n,4)
    real(kind = r8) :: ak(n)
    integer(kind = i8) :: j,i
    real(kind = r8) :: qoh
    qoh = q / h
    xstar = x    
    do j=1,4
        warray(:,j) = gaussian(n)*sqrt(qs(j)*qoh)
        !ts(j) = t + sum(as(:j-1,j)) * h      
        do concurrent (i=1:n)
            ak(i)=dot_product(as(:j-1,j), ks(i,:j-1))  
        enddo        
        !ak=zero
        !do i=1,j-1
        !    ak = ak + as(i,j)*ks(:,i)
        !enddo
        xs(:,j) = x + ak  
        ks(:,j) = h * (  fi_vec(xs(:,j)) + gi_vec(xs(:,j))*warray(:,j)  )
        xstar = xstar + alphas(j)*ks(:,j)     
    enddo
return
end subroutine rk4_ti_step_vec_mod1

Now, just a simple question, has anyone used Kasdin’s method?
How is it compared with those SDEs solvers in the differentialequations.jl?
Thank you very much!

PS
I used Kasdin’s method in Fortran with very a simple model for just one SDE, f=-x, g=1. Time is in [0,1], x(0)=20. As a bench mark test, I set step size = 1/1000, and run it for 10000 times. So totally call John’s RK subroutine for 10^7 times, and it took my laptop 1.1 seconds using Fortran.

subroutine rk4_ti_step ( x, t, h, q, fi, gi, seed, xstar )

!*****************************************************************************80
!
!! RK4_TI_STEP takes one step of a stochastic Runge Kutta scheme.
!
! Discussion:
!
! The Runge-Kutta scheme is fourth-order, and suitable for time-invariant
! systems in which F and G do not depend explicitly on time.
!
! d/dx X(t,xsi) = F ( X(t,xsi) ) + G ( X(t,xsi) ) * w(t,xsi)
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 20 June 2010
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Jeremy Kasdin,
! Runge-Kutta algorithm for the numerical integration of
! stochastic differential equations,
! Journal of Guidance, Control, and Dynamics,
! Volume 18, Number 1, January-February 1995, pages 114-120.
!
! Jeremy Kasdin,
! Discrete Simulation of Colored Noise and Stochastic Processes
! and 1/f^a Power Law Noise Generation,
! Proceedings of the IEEE,
! Volume 83, Number 5, 1995, pages 802-827.
!
! Parameters:
!
! Input, real ( kind = 8 ) X, the value at the current time.
!
! Input, real ( kind = 8 ) T, the current time.
!
! Input, real ( kind = 8 ) H, the time step.
!
! Input, real ( kind = 8 ) Q, the spectral density of the input white noise.
!
! Input, external real ( kind = 8 ) FI, the name of the deterministic
! right hand side function.
!
! Input, external real ( kind = 8 ) GI, the name of the stochastic
! right hand side function.
!
! Input/output, integer ( kind = 4 ) SEED, a seed for the random
! number generator.
!
! Output, real ( kind = 8 ) XSTAR, the value at time T+H.
!
implicit none

real ( kind = 8 ) a21
real ( kind = 8 ) a31
real ( kind = 8 ) a32
real ( kind = 8 ) a41
real ( kind = 8 ) a42
real ( kind = 8 ) a43
real ( kind = 8 ) a51
real ( kind = 8 ) a52
real ( kind = 8 ) a53
real ( kind = 8 ) a54
real ( kind = 8 ), external :: fi
real ( kind = 8 ), external :: gi
real ( kind = 8 ) h
real ( kind = 8 ) k1
real ( kind = 8 ) k2
real ( kind = 8 ) k3
real ( kind = 8 ) k4
real ( kind = 8 ) q
real ( kind = 8 ) q1
real ( kind = 8 ) q2
real ( kind = 8 ) q3
real ( kind = 8 ) q4
real ( kind = 8 ) r8_normal_01
integer ( kind = 4 ) seed
real ( kind = 8 ) t
real ( kind = 8 ) t1
real ( kind = 8 ) t2
real ( kind = 8 ) t3
real ( kind = 8 ) t4
real ( kind = 8 ) w1
real ( kind = 8 ) w2
real ( kind = 8 ) w3
real ( kind = 8 ) w4
real ( kind = 8 ) x
real ( kind = 8 ) x1
real ( kind = 8 ) x2
real ( kind = 8 ) x3
real ( kind = 8 ) x4
real ( kind = 8 ) xstar

a21 = 2.71644396264860D+00
a31 = - 6.95653259006152D+00
a32 = 0.78313689457981D+00
a41 = 0.0D+00
a42 = 0.48257353309214D+00
a43 = 0.26171080165848D+00
a51 = 0.47012396888046D+00
a52 = 0.36597075368373D+00
a53 = 0.08906615686702D+00
a54 = 0.07483912056879D+00

q1 = 2.12709852335625D+00
q2 = 2.73245878238737D+00
q3 = 11.22760917474960D+00
q4 = 13.36199560336697D+00

t1 = t
x1 = x
w1 = r8_normal_01 ( seed ) * sqrt ( q1 * q / h )
k1 = h * fi ( x1 ) + h * gi ( x1 ) * w1

t2 = t1 + a21 * h
x2 = x1 + a21 * k1
w2 = r8_normal_01 ( seed ) * sqrt ( q2 * q / h )
k2 = h * fi ( x2 ) + h * gi ( x2 ) * w2

t3 = t1 + a31 * h + a32 * h
x3 = x1 + a31 * k1 + a32 * k2
w3 = r8_normal_01 ( seed ) * sqrt ( q3 * q / h )
k3 = h * fi ( x3 ) + h * gi ( x3 ) * w3

t4 = t1 + a41 * h + a42 * h + a43 * h
x4 = x1 + a41 * k1 + a42 * k2
w4 = r8_normal_01 ( seed ) * sqrt ( q4 * q / h )
k4 = h * fi ( x4 ) + h * gi ( x4 ) * w4

xstar = x1 + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4

return
end

2

No, it’s not trivial to do that. The reason is because the iterated integrals do not have an analytical solution, and so you need to do Wiktorsson approximation which is not an easy calculation. Wiktorsson approximations are implemented in StochasticDiffEq.jl so getting to order 1 is possible. If you want to know what 3rd order looks like, you can dig into https://arxiv.org/pdf/2003.14184.pdf on page 60. Note that it expands to pages and pages of calculations in many cases.

So what’s going on here? Kasdin’s paper isn’t the only one to make this mistake, Wilke had a paper like this as well:

“4th order”. So why are we always talking about 1st order and 1.5 order? The fact that they are using integer orders should scare you because stochastic ODEs have 1/2-order expansion terms in the stochastic Taylor series. Indeed, the main authors of this field have a really nice paper which comments on this last one describing what is going wrong in the calculation:

So no, you cannot just take RK4 and add a few random numbers to it and get a 4th order method: it’s missing terms from the Taylor expansion that fundamentally do not exist in the deterministic Taylor expansion and dropping those terms means you still have a lower order scheme. In particular, you cannot have a first order scheme without having a dW^2 - dt term in the calculation (i.e. Milstein’s correction) which is a clear indicator that this will be at most 1/2-order strong 1st weak. If you cross reference the classical RK coefficients to Rossler’s schemes:

https://epubs.siam.org/doi/abs/10.1137/09076636X?journalCode=sjnaam

then you can see that all classical RK methods do satisfy the 1/2-order conditions, which then means that taking an ODE method and sticking random numbers on it is a (relatively bad) way to create a 1/2-order method, not a 4th order method, for SDEs.

The SDE literature is filled with landmines. Proceed with caution.

10 Likes
3

Thank you very much for your detailed response!
Alright, I kind of understand it a little better.
So, it looks like some guys claim interger order precision like 4th order so, they missed an 1/2 order term from the stochastic noise. Therefore their method is really like a 1/2 order precision instead of 4th order as they claimed.

That explains why I got a little confused when I checked your papers and you are discussing 1st and 1.5th order instead of 4th order.

Now, perhaps a really stupid and layman question, so, for the precision, let us say,
For some 1.5 order precision method, there is a factor, say S, in front of the the 1.5th order, so, S * O(1.5).
Now, Kasdin has a 1/2 order precision, also there is a factor K, so, K * O(1/2).

How big is the ratio S/K ?
Uhm, perhaps an O(1.5) method with big factor S is somewhat similar with an O(1/2) method with a small factor K?

But anyway, I believe your guys O(1) , O(1.5) methods could be better.

4

Remember how Taylor series coefficients have a 1/k! in front? There’s something similar in the stochastic case, and so you generally get much better coefficients just by default from a higher order scheme. This isn’t always the case because you can have RK methods optimize the coefficient, but it’s “generally the case”, and these “accidental 1/2-order methods” aren’t calculating the right coefficient anyways so they aren’t in the class of things heavily optimizing this choice.

Though there are some notable cases. The trapezoidal method for SDEs is 1/2-order, but known to have a much lower coefficient than the implicit Euler scheme which is also 1/2-order, and so it’s still suggested in some cases.

2 Likes
5

Thank you very much!

By the way, have you guys perhaps tested some factors in front of the orders ( like, you know, S * O(1.5), K * O(0.2) , the factors S and K) for different SDE methods and equations?

It might be very interesting to do some comparison perhaps.

In quantum Monte Carlo, we always have similar problem, the time step delta t dependence (due to Trotter Formula), final true result needs to be extrapolated to delta t = 0 case.
Usually O(delta t)^2 works better than O(delta t), but it might be a little bit more expensive.
Sometimes an O(delta t) is not bad, it depends on the model (Hamiltonian) we are using.
Sometimes we might need to go to some very very tiny time step to see the advantage of O(delta t)^2 or even higher order methods.

I think the high order methods could beat lower order methods when the integrand is vertically changed rapidly and a lot. Perhaps for many pharmacology situations when the integrand is relatively smooth, perhaps the difference will not be very significant.

Anyway, great job in DIfferentialEquations.jl !