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:
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?