2

I am trying to solve this system of SDEs in R, dX = (-X*a + (Y-X)b + h)dt + g dW and dY = (-Ya + (X-Y)*b)dt for time [0,200], a=0.30, b=0.2, g=1 and h is 1 for time [50,70] and 0 otherwise. I have the following code for the system for constant h=1:

library(diffeqr) ##for solving ODEs SDEs DDEs..
library(deSolve)
library(Julia)  ##used  for the sde 
library(JuliaCall) 
library(data.table)
library(minpack.lm)


tryCatch(system2(file.path(JULIA_HOME, "julia"),
                 "-E \"try println(JULIA_HOME) catch e println(Sys.BINDIR) end;\"", stdout = TRUE)[1],
         warning = function(war) {},
         error = function(err) NULL)

julia_setup(JULIA_HOME ="/home/tim/.local/share/R/JuliaCall/julia/v1.5.0/bin")

##calling Julia
JuliaCall:::julia_locate(JULIA_HOME = "/home/tim/.local/share/R/JuliaCall/julia/v1.5.0/bin")


## setup of the packages
de <- diffeqr::diffeq_setup() 
julia <- julia_setup()


#u0 is the starting value for the vector
#func1 is the ODE TO SOLVE
#f is the drift part in the sde
#g is the diffusion part in the sde
#p are the parameters associated to the sde
# noise_rate_prototype is a sparse matrix by making a dense matrix and setting some values as not zero

#starting value for the three variables
u0<-c(0,0)


## parameters
a=0.33;
b=0.2;
g=1;
h=1;

## parameters of the sde 
p<-c(a, b, g, h)

## drift part in the SDE
f <- JuliaCall::julia_eval("
function f(dy,y,p,t)
  dy[1] = -y[1]*p[1] + (y[2]-y[1])*p[2] + p[4]
  dy[2] = -y[2]*p[1] + (y[1]-y[2])*p[2]
end")



## diffusion part in the SDE
g <- JuliaCall::julia_eval("
function g(dy,y,p,t)
  dy[1,1] = p[3]
  dy[2,1] = 0
end")



## matrix associated to the noise
noise_rate_prototype <- matrix(rep(1,2), nrow = 2, ncol = 1)

tspan <- c(0.,200)

JuliaCall::julia_assign("u0", u0); ## initial value
JuliaCall::julia_assign("tspan", tspan); ## grid
JuliaCall::julia_assign("p", p); ## parameters
JuliaCall::julia_assign("noise_rate_prototype", noise_rate_prototype); ## matrix associated to the noise

## call the SDE problem
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)))
    
#plot the 2 figures, 
plot(sol$t,udf$V1,type="l",col="purple", ylab="Y1") #L
lines(sol$t,udf$V2,col="red",ylab="Y2") #R

I am not sure how to include the time dependence of h, I tried to do it by including a forcing variable but I got a lengthy Julia error (stating that h is not defined- code below) . So perhaps thats not the way to go. I'd appreciate any help regarding this.

The chunk of the code I changed to include the time dependent parameters is this (the new or changes lines have an arrow to their right):

## drift part in the SDE
f <- JuliaCall::julia_eval("
function f(dy,y,p,t)
  MU <- signal(t)                                    # <-
  dy[1] = -y[1]*p[1] + (y[2]-y[1])*p[2] + MU         # <- 
  dy[2] = -y[2]*p[1] + (y[1]-y[2])*p[2]
end")

signal <- approxfun(x = c(0, 35, 125, 200),          # <- 
                    y = c(0, 1,  0,  0), 
                    method = "constant", rule = 2)

## diffusion part in the SDE
g <- JuliaCall::julia_eval("
function g(dy,y,p,t)
  dy[1,1] = p[3]
  dy[2,1] = 0
end")



## matrix associated to the noise
noise_rate_prototype <- matrix(rep(1,2), nrow = 2, ncol = 1)

tspan <- c(0,200)

JuliaCall::julia_assign("u0", u0); ## initial value
JuliaCall::julia_assign("tspan", tspan); ## grid
JuliaCall::julia_assign("p", p); ## parameters
JuliaCall::julia_assign("signal", signal);                   # <-
JuliaCall::julia_assign("noise_rate_prototype", noise_rate_prototype); ## matrix associated to the noise

## call the SDE problem
prob <- JuliaCall::julia_eval("SDEProblem(f, g, u0,tspan,p, noise_rate_prototype=noise_rate_prototype, signal=signal)");      # <-

The rest of the code remained the same.

14thTimeLord
  • 363
  • 1
  • 14
  • What's the lengthy Julia error? A stacktrace would help narrow down what went wrong. – BatWannaBe Nov 06 '21 at 05:59
  • I am editing the post to include that information. – 14thTimeLord Nov 06 '21 at 07:00
  • 1
    You didn't add the stack trace (the actual message that is printed out when the error happens), but something stands out already: `MU <- signal(t)` was inserted into the Julia code, but it's R syntax. Amend it to Julia syntax `MU = signal(t)`. Never mixed Julia and R so I'm not sure if you can even make a Julia function use an R function like that from an R script, but if it errors, you can just edit the post and include the updated stacktrace. – BatWannaBe Nov 06 '21 at 08:24
  • 1
    I apologise, still getting myself familiarised with the terminology. But you were right, the syntax was causing the error that h (or MU in the second chunk of code) isn't defined. It's running now :) – 14thTimeLord Nov 06 '21 at 12:05
  • 1
    Pleased to hear it works now. Very cool how the languages can use each other's user-defined functions both ways like that – BatWannaBe Nov 06 '21 at 14:11
  • It would be also cool if the OP would formulate the solution explicitly as an answer. Then it is more visible and contributes to the knowledge base. SO guidelines allow (and in my opinion) even recommend this. – tpetzoldt Nov 10 '21 at 06:57
  • Yes I agree, the code is a complete way to solve SDEs with time dependent parameters (after the correection) so if @batwannabe you could put your comment as an answer, i'll approve it :) – 14thTimeLord Nov 11 '21 at 10:53

0 Answers0