 # Calculating the analytical solution of an Ornstein-Uhlenbeck SDE

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.

(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()")
``````

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?

1 Like

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: