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.
Here is the code I am using:
γ = 34.0
D = 69.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?
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)
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:
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:
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.
You might want to join our Slack if you have more questions.