# Calculating the analytical solution of an Ornstein-Uhlenbeck SDE
**URL:** https://discourse.julialang.org/t/calculating-the-analytical-solution-of-an-ornstein-uhlenbeck-sde/37359
**Category:** Numerics
**Tags:** diffeq, sde
**Created:** [April 10, 2020, 5:50pm UTC](https://discourse.julialang.org/t/calculating-the-analytical-solution-of-an-ornstein-uhlenbeck-sde/37359 "2020-04-10T17:50:49Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@williamjsdavis](https://discourse.julialang.org/u/williamjsdavis)
#### Post date: [April 10, 2020, 5:50pm UTC](https://discourse.julialang.org/t/calculating-the-analytical-solution-of-an-ornstein-uhlenbeck-sde/37359/1 "2020-04-10T17:50:49Z")
I want to compare different SDE solvers against the analytical solution for the [Ornstein-Uhlenbeck process](https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_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:
```plaintext
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](https://docs.sciml.ai/latest/tutorials/sde_example/#Using-Higher-Order-Methods-1), 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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 10, 2020, 6:04pm UTC](https://discourse.julialang.org/t/calculating-the-analytical-solution-of-an-ornstein-uhlenbeck-sde/37359/2 "2020-04-10T18:04:00Z")
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:
```julia
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](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](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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 10, 2020, 6:04pm UTC](https://discourse.julialang.org/t/calculating-the-analytical-solution-of-an-ornstein-uhlenbeck-sde/37359/3 "2020-04-10T18:04:25Z")
You might want to join our Slack if you have more questions.