How does one test a numerical ODE solver implementation?

Question

I'm about to start working on a software library of numerical ODE solvers, and I'm struggling with how to formulate tests for the solver implementations. My ambition is that the library, eventually, will include solvers for both nonstiff and stiff problems, and at least one implicit solver (more or less on par with the capabilities of the ode routines in Matlab), so the test methodology needs to reflect the various types of problems and criteria for different solvers.

My problem now is that I don't know where to begin with this testing. I can think of a few different ways to test the output of an algorithm:

• Test for a problem that has an analytical solution, and check that the numerical solution is within tolerance levels for all the returned data points. This requires knowledge of a number of analytical problems which exhibit all the properties that I want the different solvers to work with (stiffness, implicit problems etc), which I don't have, at least not off the top of my head.

  This method tests the results of a solver method. Thus, there is no guarantee that the solver actually works, just that it works for the given test problem. Therefore, I suspect a large number of test problems is needed to confidently verify that the solver works.

• Manually calculate the solution for a few time steps using the algorithms I intend to implement, and then do the same with the solvers and check that the results are the same. This requires no knowledge of the true solution to the problem, but in turn requires quite a lot of hands-on work.

  This method, on the other hand, only tests the algorithm, which is fine by me - if someone else has proven that 4^th order Runge-Kutta works, I don't feel a desperate need to. However, I do worry that it will be very cumbersome to formulate test cases, as I don't know a good method to generate the test data (except maybe by hand, which will be a lot of work...).

Both the above methods have serious limitations for me with my current knowledge - I don't know a good set of test problems for the first one, and I don't know a good method of generating test data for the second.

Is there other ways to verify numerical ODE solvers? Are there other criteria on the implementations that should be verified? Are there any good (free) resources on testing ODE solvers out there¹?

EDIT:
Since this question is very broad, I want to clarify a little. The test suite I want to create will fill two main purposes:

1. Verifying that the solvers work as expected, for the problems they're intended to solve. In other words, a solver for non-stiff problems is allowed to go bananas on a stiff problem, but should perform well on non-stiff problems. Also, if there are other solvers in the library that offer higher accuracy, it might not be necessary to enforce very accurate results - just "accurate enough". Thus, part of my question is what tests should be used for what solvers; or, at least, how one should reason to decide that.

2. Sanity test upon installation of the library. These test need not (should not) be elaborate or time-consuming; just the very basics that can be run in under 5 seconds, but that will alert the user if something is off-the-charts weird. Thus, I also need a way to construct tests that are very simple, but that still tell me something about the state of the library.

---

^{1 Yes, I've been Googling my eyes out, but most of what I find is lecture notes with very trivial examples, with the notable exception of the CWI ODE test set from Bari which I don't know if, or how, I could use for my purposes, since it treats much more sophisticated solvers than the ones I want to test...}

Answer 1

This is a very broad question and I am going to give you some things to think about (some are already included in your post, but they are repeated here for completeness).

Scope of Problems

• You need to define the interface of how to specify problems.
• Are you going to allow parameters that can be fixed or can vary for solutions?
• Are you going to allow perturbation parameters to slightly perturb problems and see if they are still solvable (for example, an $\epsilon$ parameter to be defined anywhere) in a specific problem?
• Are you going to allow infinite precision?
• Are you going to test for speed and sensitivity to numerical precision?
• Have you chosen two (maybe more) libraries that already exist to compare results?
• How will you choose stopping criteria, will you use various methods and let the user select or define their own?
• Are you going to measure error using various measures and allow the user to turn those on and off?
• Have you looked at the professional packages like Computer-Algebra-Systems (CAS) and understand all of the options they allow?
• Are you going to allow displaying of results and/or comparisons and/or plots?

Problem Recommendations

• You need to write a test specification defining the source of problems, the scope of how problems were tested, capturing results and metrics of running the routines.
• I would certainly look to other libraries already out there for the problems they are using (maybe test files).
• I would go to college libraries and go through books on ODEs and pull out problems of all types, those with known closed form or numeric only solutions.
• Case 1: We want as many variations of closed form solution problems as we can get in order to compare exact versus numerical results.
• Case 2: I would go to every numerical analysis book I can find and capture the worked examples and duplicate them. I would additionally capture the problem sets, particularly the ones that have some pathology that exist in most books (sensitivity to this or that types).
• Case 3: I would go to different branches of applied math like Physics, Ecology, Biology, Economics, et. al and capture problems from each of those domains to validate that your specification language for problems allows for such examples.
• Case 4: I would research papers/journals that contain the most useful examples where the particular author had to modify a particular approach to account for some pathology or weirdness or hardness.
• Case 5: Search the web for additional examples. For stiff, see the references here and peruse them ALL to ferret out test problems. Here are some MATLAB examples to peruse.

This is not unique. If you look at the book "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by Dennis and Schnabel, Appendix B, "Test Problems", you can see how they did it. After developing one of the most beautiful set of algorithms write ups I have ever seen, they threw a collection of problems at it that made if go nuts. You had to tweak here and there! They included five very different and pathological problems that strained the capabilities of the solvers. This has taught me that we can continue to throw problems at algorithms that they are incapable of handling for a host of reasons. Note, they even borrowed this set of problems from More', Garbow and Hillstrom (you can also look up that reference and perhaps there are others you can use as a guide).

In other words, this is not a trivial task. You need Known-Answer-test cases that always allow you to test the validity of updates and don't break things. That is, a repeatable and extensive set of problems from low to high, from easy to hard, from possible to impossible, ... You also need a collection of problems that your solvers cannot handle in order to truly understand its limitations.

Answer 2

One sanity check I run against my ODE solvers is to simply check it on smaller linear systems via exactly computing the exponential of the system's matrix. i.e. given

$$\frac{d\mathbf{u}}{dt} = \mathbf{Au} $$

check the error in

$$\exp{(t\mathbf{A})}\mathbf{u}_0 - \hat{\mathbf{u}}(t) $$

where $\hat{\mathbf{u}}(t)$ is your calculated solution. This way you don't need to know any real analytic solution (assuming you have some equivalent of 'expm'), but you can still get an idea of whether your solver is doing what it should be.

Just don't calculate the exponential with one of your time-steppers (i.e. dubious method number 6 :) http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf )

Answer 3

You might try investigating the "Method of Manufactured Solutions", which is a general methodology used to test implementation of codes which solve PDEs (it can be used to find both mathematical and coding errors). I imagine it could be adapted to work for solving ODEs, if your solution methodology is general enough.

http://prod.sandia.gov/techlib/access-control.cgi/2000/001444.pdf