Numerical Integration of ODEs - Consistency

Question

EDIT: The answers below were helpful, but didn't get at the core of the problem, which I think has to do with the impulse aspect of things. I've since re-branded this question at Runge Kutta with Impulse . Cheers.

So I'm developing a method for studying ODEs that is similar in style to 'impulse response' but which is more general, though I have yet to find something as analytically beautiful as a Greens Function for it. Right now, everything is numerical/computational.

Having never done this before, I basically copy and pasted Runge-Kutta equations from Wikipedia into my simple python script. The results look beautiful, but I'm concerned that there are inconsistencies. Basically I'm trying to measure something I call 'energy residence time' using simulations of the ODEs, but I get different results for different step sizes and I'm not sure things are converging at smaller step sizes. Are there standard ways I can test my code/algorithm to make sure its legit, and hone in on where my problems lie?

If you want more info on the method itself (and the notion of energy residence time), I'd love to discuss!

Answer 1

I would recommend taking known examples (for example, those which have known closed form solutions so you can compare), run those through your code and then compare error bounds.

Additionally, I would look for examples that have peculiarities and difficulties like pathological examples and use those to help test the robustness of your code. This can include everything from singularities, dependence on initial conditions or really hard examples.

Lastly, you should have some idea about the computational complexity of what you've done. Are there more efficient ways, are there things that can be done to improve results, ...

Regards

Answer 2

For a Runge-Kutta method, if your code converges at some step size $h$, it is guaranteed to converge for some step size $\tilde{h} < h$.

Explicit Runge-Kutta methods have a region of absolute convergence in the complex plane. Effectively, depending on the order of the method, convergence is guaranteed if $h\lambda$ is inside the region of absolute stability. The region of absolute stability is simply connected (i.e. any closed curve inside the region is homotopic to the point -1), so shrinking $h$ does not remove you from this region.

If you are getting different results when you shrink your step size, then I would expect you are getting round-off error, which would lead me to question your results at larger step sizes.

Alternatively, I might suspect that your smaller step size solutions are more accurate, and that you may indeed have a stiff system, and hence your larger step size solutions are converging to the solution of a non-stiff ODE instead.

You could apply an adaptive RK method that would avoid this step size issue entirely, but it will not avoid the problem caused by stiffness.

I would suggest posting your ODE system so that the issues can be diagnosed.