What algorithm/option in DifferentialEquations should I choose to solve efficiently a large and stiff system of ODE (a system of 100 to 10_000 say) when I am only interested in the value of the stationary solution (i.e I know that system of ODE converges to a constant vector and I just want to find this constant).
CVODE_BDF(linear_solver = :GMRES) with the solver options save_everystep=false,save_start=false is likely what you’re looking for.
This blog post explains a workflow for doing exactly this: http://www.stochasticlifestyle.com/solving-systems-stochastic-pdes-using-gpus-julia/
If you know that there is a unique stationary solution and you are only interested in it, maybe you try a Newton method to find the equilibrium directly.
1 Like
True. But that doesn’t capture initial condition dependence if that’s a thing in this model. FWIW if you use the steady state stuff then this is simplified, and then switching between rootfinding and a dynamical form is trivialized:
http://docs.juliadiffeq.org/latest/solvers/steady_state_solve.html#SteadyStateDiffEq.jl-1
# Newton method via nlsolve choice, defaults to NLsolve.jl
sol = solve(prob,SSRootfind())
# Dynamical solve to steady-state, automatically exits when `norm(du)< ...`
sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)
2 Likes
DynamicSS is exactly what I was looking for. Thanks a lot!