0

I have a problem with approximating ROC for numerical algorithms for solving ODE's. It's known, that the global error for Euler's method is of order $O(h)$, and for 4th order Runge Kutta it's $O(h^4)$. Using linear regression for log-log correspondence between error and the step size however returns a variety of results for different equations (e.g. pendulum eq or trivial $y'=y$).

Is it a mistake to expect exact (sort of) analytical results for various problems? Should I pick large points for error evaluation instead of choosing e.g. $t = 10$? Or maybe I'm using the wrong formula? I've chosen the one from wikipedia - the $p$ in discretization methods https://en.wikipedia.org/wiki/Rate_of_convergence.

Lutz Lehmann
  • 126,666
jerry
  • 1
  • Please give some more of what you tried and what you got. The numerical methods for a given floating point data type have a working interval or sweet point where they behave as you expect. Above that point, for larger step sizes, the higher order error terms become noticeable, below that point floating point noise accumulates to a size that dominates the method error. Thus your log-log plot should look like a V shape if you have a high enough sampling density in your step sizes. Euler should be linear for h=1e-7 .. 1e-3, RK4 for 1e-3 .. 1e-1. – Lutz Lehmann Oct 23 '18 at 10:01
  • Longer answer with table demonstrating the leading error term in https://math.stackexchange.com/q/1238995/115115, https://math.stackexchange.com/q/1786896/115115 – Lutz Lehmann Oct 23 '18 at 10:32
  • Your advice on seeking the so called "sweet spot" worked perfectly on Euler's method. The problem still lies with RK4. The error does not converge to 0 linearly in the pendulum case - it oscillates. The method conserves the oscillation property, thus despite a large stepsize the error does not explode. I'm able to find the stopping time and stepsize vector such that the '4' is achieved, but the loglog plot is just unsatisfactory mess, as it includes error oscillations. If I exclude them, I achieve ROC above 10. I'd like to achieve a neat linearity as in the Euler case if it's possible. – jerry Oct 23 '18 at 18:11
  • For a detailed error analysis I would need the code, especially how the time loop is constructed. It is easy to get a "naive" implementation that produces a random error of size O(h) by not meeting exactly the final time. – Lutz Lehmann Oct 23 '18 at 20:02
  • If you use a loop like t,y = t0, y0; while t+h<tf: t,y = t+h, methodstep(sys,t,y,h); you need to finish this by t,y = tf, methodstep(sys,t,y,tf-t). To avoid almost-zero steps, use something like while t+1.01*h<tf:... – Lutz Lehmann Oct 23 '18 at 21:33

0 Answers0