Method of Line with DifferentialEquations.jl


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:


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.


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:

# 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)


DynamicSS is exactly what I was looking for. Thanks a lot!