What are some tips on developing a problem-specific ODE solver?

Question

I have a small system of stiff ODEs describing a chemical reaction. The right-hand side is quite complicated, as well as the Jacobian. This equation will be solved many times with different initial conditions. Given that I've already selected a stiffness-oriented method, is it realystic to create more effective code than some existing solver? For example, is it worth trying to invest into optimizing LU-decomposition for my small system by loop unrolling, or whatever?

In this thread some general ideas are given, but maybe there are some tips about stiff ODE solution?

Update: there are only three equations. I'm using Fortran and compared the efficiency of RADAU (implicit RK), RODAS (Rosenbrock) and DLSODE (BDF method) codes. RODAS seems to work a bit faster.

Answer 1

I assume you have already verified that your system is stiff? Otherwise you might be paying a significant performance penalty for using an implicit solver.

The next thing I suggest is to gather more data from the ODE solvers you are using. Typically this kind of data is available in some additional output variables or by turning on some kind of verbose output flag.

1. Total number of time steps in the solution.
2. Minimum and maximum time step sizes.
3. Number of times in the solution where the time step is changed.
4. Number of times the Jacobian is updated and factored.

You say your RHS is quite complicated so I'm guessing your problem may be very nonlinear. The algorithm will change time step size either because the problem is very nonlinear and the algorithm is having difficulty obtaining a converged solution or because it didn't meet the specified accuracy requirements (Have you experimented with different accuracy tolerances?) If there are a large number of Jacobian formations and factorizations, those routines might be candidates for optimization.

You might want to implement a simple fixed time-step backward Euler solver and see how the solution time compares with the off-the-shelf solvers.

Finally, if you able to get some insights about the characteristics of your ODE system from the above experiments, you could try to profile your code. If you are running on Linux, you can use gprof or valgrind/callgrind (http://valgrind.org/docs/manual/cl-manual.html) to obtain data about where the time is being spent.