Calculating the analytical solution of an Ornstein-Uhlenbeck SDE

Specific Domains Numerics
,
1

I want to compare different SDE solvers against the analytical solution for the Ornstein-Uhlenbeck process. Currently I am doing this by solving numerically first (e.g. Euler-Maruyama), saving the noise, and calculating the integral solution, e.g.

Screen Shot 2020-04-10 at 10.39.52
Screen Shot 2020-04-10 at 10.39.52806×120 7.09 KB
(1)

Here is the code I am using:

using DifferentialEquations
using Plots

γ = 34.0
D = 69.0
μ= 5.3
u₀=7.0
f(u,p,t) = -γ*(u - μ)
g(u,p,t) = sqrt(D)
dt = 0.01
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u₀,span)

sol = solve(prob,EM(),dt=dt,save_noise=true)

# Analytical solution
ex = exp.(-γ*sol.t)
inNoise = [0; diff(sol.W.u[1:end-1])]
anaSol = u₀*ex + μ*(1 .- ex) + sqrt(D)*ex.*cumsum(inNoise.*exp.(γ*sol.t))

plot(sol.t,anaSol,label = "Ana")
plot!(sol,label = "EM()")

Screen Shot 2020-04-10 at 10.29.56
Screen Shot 2020-04-10 at 10.29.561204×790 109 KB

This is the result I want, but it involves saving the noise and calculating Weiner increments for the integration. I know that SDEProblem can be used to add a defined analytical solution to an SDE, but I’m not sure the solution (1) can be cast in this form. Is there a more elegant way to calculate the analytical solution?

2

Cool to see someone’s interested in this. You might want to dig a little deeper into the dev tools for this. You can do:

analytical = TestSolution(t,u)
errsol = appxtrue(sol,analytical)

and now errsol has the errors calculated, plots analytical solutions in the normal way, etc. You can see that in a lot of benchmarks we use TestSolution when we need alternative ways to calculate analytical results, so all of the benchmarking functions allow it, like in here:

https://benchmarks.sciml.ai/html/NonStiffODE/Pleiades_wpd.html

I am not sure we added in a hook for W-dependent TestSolution into work-precision diagrams, so that probably needs a PR to DiffEqDevTools. The closest to what you’re trying to do are the strong stochastic Lotka-Volterra benchmarks:

https://benchmarks.sciml.ai/html/NonStiffSDE/LotkaVolterraSDE.html

where there it’s using a low tolerance setup for calculating an “analytical solution” but of course that’s highly costly and makes the benchmark take awhile.

3

You might want to join our Slack if you have more questions.