# How to represent an Ornstein-like stochastic ODE in Julia?

**URL:** <https://discourse.julialang.org/t/how-to-represent-an-ornstein-like-stochastic-ode-in-julia/95987>\
**Category:** General Usage\
**Tags:** ode, sde\
**Created:** [March 13, 2023, 7:50am UTC](https://discourse.julialang.org/t/how-to-represent-an-ornstein-like-stochastic-ode-in-julia/95987 "2023-03-13T07:50:27Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![tomkimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomkimpson/32/39168_2.png) [@tomkimpson](https://discourse.julialang.org/u/tomkimpson)\
**Post date:** [March 13, 2023, 7:50am UTC](https://discourse.julialang.org/t/how-to-represent-an-ornstein-like-stochastic-ode-in-julia/95987/1 "2023-03-13T07:50:27Z")

</div>

## 1. The equation I want to solve

I have some vector stochastic ODE that has an almost [Langevin form](https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process), but with an additional time dependence:

\frac{dy}{dt} = -\gamma (y - a - bt) + b + \xi(t)

where a,b,\gamma are constant parameters and \xi(t) is a white noise process characterised by variance \sigma (i.e. \langle \xi(t),\xi(t') \rangle = \sigma^2 \delta(t-t') )

This can also be written as  
dy = \gamma (\mu(t) - y) dt + \xi dt = \gamma (\mu - y) dt + dW\_t

for \mu(t) = a + b t + b and W\_t a Wiener process

## 2. My attempt to implement in Julia

I want to integrate this equation in Julia.

A MWE might be something like

```julia
using DifferentialEquations, Noise, Plots 

#Vector of parameters
γ = [1e-13, 1e-14, 1e-15]
a = [300,400,500]
b = [-1e-15,2e-15,3e-15]

#Initial conditions on y 
y0 = [300,400,500]

#Define the ODEs 

f(du,u,p,t) = (du .= -γ.*u .+ γ.*(a .+ b*t) .+ b)
g(du,u,p,t) = (du .= 1.0)

#Specify the noise process to use 
noise = WienerProcess(0., 0.) # WienerProcess(t0,W0) where t0 is the initial value of time and W0 the initial value of the process

#Time over which to integrate 
t = range(0.0, 3e8, length=500)
tspan = (first(t),last(t))

#Setup SDE problem and solve
prob = SDEProblem(f,g,y0,tspan,tstops=t,noise=noise)
solution = solve(prob,EM())

#Plot it for interest
i = 1
solution_i = solution[i,:]
plt = plot(t,solution_i)
display(plt)

```

### 3. My question

How do I include \sigma as a parameter which can be set?

[From the docs](https://docs.sciml.ai/DiffEqNoiseProcess/stable/noise_processes/) it suggests that the variance of the Wiener process is given by `dt`. I understand I could do something like `noise.dt = σ` or set the stepsize in the integrator accordingly (dropping `tstops`) but this then becomes a very long integration depending on choice of \sigma, and I only care about the value of the solution at `tstops`.

Another option is to include it in `g(du,u,p,t)`, but then I am unsure how this “interacts” with the `WienerProcess`.

Thanks in advance for any help

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 13, 2023, 10:52am UTC](https://discourse.julialang.org/t/how-to-represent-an-ornstein-like-stochastic-ode-in-julia/95987/2 "2023-03-13T10:52:37Z")

</div>

> [@tomkimpson](#):
>
> How do I include \sigmaσ\sigma as a parameter which can be set?

> [@tomkimpson](#):
>
> Another option is to include it in `g(du,u,p,t)`, but then I am unsure how this “interacts” with the `WienerProcess`.

When X ~ N(0,dt), sigma\*X ~ N(0,sigma^2 \* dt) (using N(mean, variance)). So yes, `g(du,u,p,t) = (du .= p.σ)` or however you want to write it would set the variance.