Using SDE for the first time

New to Julia
1

When I solve problems that stochastic effects are presents, I usually discretize the domain and then I create an interpolation vector with random numbers as my stochastic source for the dynamics. But now that I am learning Julia, I notice that there are tons of good packages for these systems, and I think its good practice to start using them, so, lets consider this pictorial example:

Suppose we want to solve an electrical resistor with thermal noise, i.e.,

\dot{U} + \frac{U}{RC} = F_0

where F_0 obeys the following fluctuation-dissipation relation

<F_0(t)F_0(t’)> = \frac{2k_bT}{RC^2}*\delta(t - t’)

Regardless of what each variable is, what’s the most efficient way of solving this using DifferentialEquations.jl ? The reference that I am trying to follow is Stochastic Differential Equations · DifferentialEquations.jl (sciml.ai) but I’m having problems

2

What did you try?

3

I think as a start it would be good to translate from physics to the more mathematical language of SDEs. Can you identify drift function, diffusion coefficient?

4

You are completly right, I will put the code here, it will be easy to discuss.
The non-stochastic problem I’m doing this way

#----------------------------------------------------
R = 1.0e3 # resistance [Ω]
C = 1.0e-9 # capacitance [F]

function circuit!(du,u,p,t)
du[1] = -u[1]/(R*C)
end

uₒ = [10.0]
Δt = (0.0, 10RC)
prob_ode = ODEProblem(circuit!, uₒ, Δt)
@time sol = solve(prob_ode)
plot(sol.t, sol[1,:])
#----------------------------------------------------

If there is anything already here that I could do to optimize things, I’d appreciate that. Then, going to the stochastic case (and following Stochastic Differential Equations · DifferentialEquations.jl (sciml.ai) we would have something like:

#----------------------------------------------------
kb = 1.380649e-23 # bolztmann constant [J/K]
R = 1.0e3 # resistance [Ω]
C = 1.0e-9 # capacitance [F]
T = 300.0 # temperatura [K]

function circuit!(du,u,p,t)
du[1] = -u[1]/(R*C)
end

function fluctuation(du,u,p,t)
du[1] = 2kbT/(RCC)
end

uₒ = [10.0]
Δt = (0.0, 10RC)
prob_sde = SDEProblem(circuit!, fluctuation, uₒ, Δt)
@time sol = solve(prob_sde)
plot(sol.t, sol[1,:])
#----------------------------------------------------

We do have the exponential profile from the ODE case, which is an agreement, however, my problems are:
1- what exactly the function ‘fluctuation’ is doing? At the tutorial, I understand that it is a Wieger process, if so, how can I chance it (if one desires to)
2- I’m almost sure that putting “2kbT/(RCC)” in the ‘fluctiation’ function defined stills do not guarantee the fluctuation-dissipation theorem <F_0(t)F_0(t’)> = \frac{2 kb T}{RC^2}*\delta(t - t’). In other words, I’m confused about the implementation of the noise variable.

One thing that I sure might do is to compute the mean and variance of my system, and check if it matches the theory, but this is really expensive because, in these examples shown, the ODE case was done in 3seconds, but the SDE case in 77seconds, and I would have to simulate lots of points so…

6

You want a square root in your fluctuation function to be in the right units and adding noise to it shouldn’t help things, because by Donsker’s theorem it just acts as further scaling of the noise