Movatterモバイル変換

[0]ホーム
URL:

Solving Stochastic Differential Equations(SDE) in R with diffeqr

Chris Rackauckas

2024-12-04

1D SDEs

Solving stochastic differential equations (SDEs) is the similar toODEs. To solve an SDE, you usediffeqr::sde.solve and givetwo functions:f andg, wheredu = f(u,t)dt + g(u,t)dW_t

de<- diffeqr::diffeq_setup()f<-function(u,p,t) {return(1.01*u)}g<-function(u,p,t) {return(0.87*u)}u0<-1/2tspan<-list(0.0,1.0)prob<- de$SDEProblem(f,g,u0,tspan)sol<- de$solve(prob)udf<-as.data.frame(t(sapply(sol$u,identity)))plotly::plot_ly(udf,x = sol$t,y = sol$u,type ='scatter',mode ='lines')
geometric_sdes
geometric_sdes

Systems of Diagonal Noise SDEs

Let’s add diagonal multiplicative noise to the Lorenz attractor.diffeqr defaults to diagonal noise when a system of equations is given.This is a unique noise term per system variable. Thus we generalize ourprevious functions as follows:

f<-function(u,p,t) {  du1= p[1]*(u[2]-u[1])  du2= u[1]*(p[2]-u[3])- u[2]  du3= u[1]*u[2]- p[3]*u[3]return(c(du1,du2,du3))}g<-function(u,p,t) {return(c(0.3*u[1],0.3*u[2],0.3*u[3]))}u0<-c(1.0,0.0,0.0)tspan<-c(0.0,1.0)p<-c(10.0,28.0,8/3)prob<- de$SDEProblem(f,g,u0,tspan,p)sol<- de$solve(prob,saveat=0.005)udf<-as.data.frame(t(sapply(sol$u,identity)))plotly::plot_ly(udf,x =~V1,y =~V2,z =~V3,type ='scatter3d',mode ='lines')

Using a JIT compiled function for the drift and diffusion functionscan greatly enhance the speed here. With the speed increase we cancomfortably solve over long time spans:

tspan<-c(0.0,100.0)prob<- de$SDEProblem(f,g,u0,tspan,p)fastprob<- diffeqr::jitoptimize_sde(de,prob)sol<- de$solve(fastprob,saveat=0.005)udf<-as.data.frame(t(sapply(sol$u,identity)))plotly::plot_ly(udf,x =~V1,y =~V2,z =~V3,type ='scatter3d',mode ='lines')

stochastic_lorenz Let’s see how much faster the JIT-compiledversion was:

>system.time({for (iin1:5){ de$solve(prob    ) }})   user  system elapsed146.400.75147.22>system.time({for (iin1:5){ de$solve(fastprob) }})   user  system elapsed1.070.101.17

Holy Monster’s Inc. that’s about 145x faster.

Systems of SDEs with Non-Diagonal Noise

In many cases you may want to share noise terms across the system.This is known as non-diagonal noise. TheDifferentialEquations.jlSDE Tutorial explains how the matrix form of the diffusion termcorresponds to the summation style of multiple Wiener processes.Essentially, the row corresponds to which system the term is applied to,and the column is which noise term. Sodu[i,j] is theamount of noise due to thejth Wiener process that’sapplied tou[i]. We solve the Lorenz system with correlatednoise as follows:

f<- JuliaCall::julia_eval("function f(du,u,p,t)  du[1] = 10.0*(u[2]-u[1])  du[2] = u[1]*(28.0-u[3]) - u[2]  du[3] = u[1]*u[2] - (8/3)*u[3]end")g<- JuliaCall::julia_eval("function g(du,u,p,t)  du[1,1] = 0.3u[1]  du[2,1] = 0.6u[1]  du[3,1] = 0.2u[1]  du[1,2] = 1.2u[2]  du[2,2] = 0.2u[2]  du[3,2] = 0.3u[2]end")u0<-c(1.0,0.0,0.0)tspan<-c(0.0,100.0)noise_rate_prototype<-matrix(c(0.0,0.0,0.0,0.0,0.0,0.0),nrow =3,ncol =2)JuliaCall::julia_assign("u0", u0)JuliaCall::julia_assign("tspan", tspan)JuliaCall::julia_assign("noise_rate_prototype", noise_rate_prototype)prob<- JuliaCall::julia_eval("SDEProblem(f, g, u0, tspan, p, noise_rate_prototype=noise_rate_prototype)")sol<- de$solve(prob)udf<-as.data.frame(t(sapply(sol$u,identity)))plotly::plot_ly(udf,x =~V1,y =~V2,z =~V3,type ='scatter3d',mode ='lines')
noise_corr
noise_corr

Here you can see that the warping effect of the noise correlations isquite visible! Note that we applied JIT compilation since it’s quitenecessary for any difficult stochastic example.

[8]ページ先頭
©2009-2025 Movatter.jp