SciML / StochasticDelayDiffEq.jl

Stochastic delay differential equations (SDDE) solvers for the SciML scientific machine learning ecosystem
MIT License
25 stars 10 forks source link

Solving stiff SDDE problems #1

Closed anirudhpammi closed 4 years ago

anirudhpammi commented 4 years ago

I tried using this package to solve a stiff problem, I am only able to use the EM or LambaEM solvers which result in some instability. I want to use stiff solvers like ImplicitEM or ImplicitRKMil. Can you help me out? Thanks!

ChrisRackauckas commented 4 years ago

Minimal example?

anirudhpammi commented 4 years ago

using BenchmarkTools, PyPlot, DifferentialEquations, StochasticDelayDiffEq

function yamada_norm!(du,u,p,t)

no noise ODE

b1,b2,mu1,mu2,s=p;
du[1] = (b1*mu1) - ( (b1*u[1])*(1+u[3]) ) ;
du[2] = (b2*mu2) - (  (b2*u[2])*( 1 + (s*u[3]) )  );
du[3] = ( u[1] - u[2] - 1 )*u[3] ;

end

function yamada_dde(du,u,h,p,t)

b1, b2, mu1, mu2, s, kappa, tau = p;
histI = h(p,t-tau)[3];
du[1] = (b1*mu1) - ( (b1*u[1])*(1+u[3]) ) ;
du[2] = (b2*mu2) - (  (b2*u[2])*( 1 + (s*u[3]) )  );
du[3] = ( u[1] - u[2] - 1 )*u[3] + kappa*(histI);

end

function yamada_sigDDE(du,u,h,p,t) b1, b2, mu1, mu2, s, kappa, tau = p; du[1] = b1*0.01; du[2] = 0; du[3] = 0; end

function yamada_sigODE(du,u,p,t)

noise component of ODE

b1, b2, mu1, mu2, s = p;
du[1] = b1*0.01;
du[2] = 0;
du[3] = 0;

end

b1 = 0.01; b2 = 0.02; mu1 = 2.4; mu2 = 2.2; s = 5; kappa = 0.05; tau = 1100.0; RT_max = 100; tspan_hist = (0.0,tau); tspan_dde = (tau,tau*RT_max) lags = [tau]; u0_hist = [mu1; mu2;3]; u0_dde = [mu1; mu2;0]; p_hist = [b1,b2,mu1,mu2,s]; p_dde = [b1,b2,mu1,mu2,s,kappa,tau];

u0_hist = [mu1; mu2;3];

p_hist = [b1,b2,mu1,mu2,s]; p_dde = [b1,b2,mu1,mu2,s,kappa,tau];

prob_SDE = SDEProblem(yamada_norm!,yamada_sigODE,u0_hist,tspan_hist,p_hist); sol_sde = solve(prob_SDE,ImplicitEM());

u0_dde = [mu1; mu2;0]; h(p,t) = sol_sde(t);

prob = SDDEProblem(yamada_dde, yamada_sigDDE, u0_dde, h, tspan_dde, p_dde; constant_lags=tau);

sol = solve(prob,ImplicitEM());

HTSykora commented 4 years ago

Hi,

SROCK methods are working, and they are also stabilized methods. The necessary pull request that make the implicit methods work are no yet included, but they are already done. So until this package is not integrated into the package, I recommend, that you use SROCK methods.

If you desperately need implicit methods, hit me up on the Julia slack channel, and I'll help you make the implicit method work, but unfortunately I will be available only from Monday.

anirudhpammi commented 4 years ago

Hi, Thanks, the S-ROCK methods work and I don't seem to have the instabilities now. While I don't need the implicit method desperately, could you please tell me if I will have any significant difference/speedup if I use ImplicitEM() v/s SROCK1()? Thanks!

ChrisRackauckas commented 4 years ago

The SROCK are better if you only have large real eigenvalues because they are stable (for real eigenvalues) while being explicit. You only need to go implicit when you have large complex eigenvalues.

ChrisRackauckas commented 4 years ago

Looks like this is all good now.