Solver RK4 still surprising behavior

[DaluS]

I had trouble with the homemade rk4, and decided to verify if our solvers were working well. In consequence, I coded a very simple couple ODE (dampened harmonic oscillator) and compared results with its analytical solution.

The Homemade RK4 still have a weird behavior. Here is a difference between the numeric and theoric value :

The slope is 1/2, it should be 4, this system surprisingly behave as a system of order 1/2.

It seems to me that the code is right, so I do not know why it is not working https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

def _rk4(dydt_func=None, dt=None, y=None, t=None):
    """
    a traditional RK4 scheme, with:
        - y = array of all variables (all ode)
        - dt = fixed time step
    """
    dy1_on_dt = dydt_func(t, y)
    dy2_on_dt = dydt_func(t + dt/2., y + dy1_on_dt * dt/2.)
    dy3_on_dt = dydt_func(t + dt/2., y + dy2_on_dt * dt/2.)
    dy4_on_dt = dydt_func(t + dt, y + dy3_on_dt*dt)
    return (dy1_on_dt + 2*dy2_on_dt + 2*dy3_on_dt + dy4_on_dt) * dt/6.

[Didou09]

Will try to investigate this in a few weeks (mid-December probably) if the problem is still there by then :-)

[DaluS]

Yes ! I'll focus on using the rk8 from scipy for the moment :)

[DaluS]

I coded a Lorenz system to check the qualitative system, for the moment both the rk2-scipy and the homemade rk4 are wrong :) (converging toward one of the point instead of being a strange attractor)

For both cases, a fork of the cockpit to test them will be joined to the next pull request

[DaluS]

I went a bit further on the subject :

• an Rk1 has been added for the very-simple case
• Coded a convergence test (see testConvergence_DampOsc and both Rk1/Rk4 are not power-law which imply that either the convergence test is wrong or that the problem is somewhere else
• testConvergence_DampOsc is showing that the first slope estimation is wrong, it might give insight on the problem

The way I see it : The initial condition on the first derivative is not respected, even if well noted in get_summary.

Maybe there is an error in the order of resolution ?

[Didou09]

So I had a look using the DampOscillator with Perfect preset (because it is easy to see any problem) and it appears to me that the homemade rk4 behaves correctly if the time step is fine enough.

Below I compare:

• erk4-scipy (variable time step optimized by scipy)
• erk4-homemade for dt = 0.01 s, 0.001 s and 0.0001 s

We can see that the homemade output converges towards the scipy result as dt decreases. I put the code so you can reproduce it ( I reduce Tmax to avoid very large computations).

In [79]: hub = pgm.Hub(model='DampOscillator', preset='Perfect', verb=False)

In [80]: hub.set_dparam(key='Tmax', value=4, verb=False)

In [81]: hub.run(solver='eRK4-scipy')

In [82]: dax = hub.plot(key='Energy', color='k')

In [83]: hub.run(solver='eRK4-homemade')

In [84]: dax = hub.plot(key='Energy', color='r', dax=dax)

In [85]: hub.set_dparam(key='dt', value=0.001, verb=False)

In [86]: hub.run(solver='eRK4-homemade')

In [87]: dax = hub.plot(key='Energy', color='b', dax=dax)

In [88]: hub.set_dparam(key='dt', value=0.0001, verb=False)

In [89]: hub.run(solver='eRK4-homemade')

In [90]: dax = hub.plot(key='Energy', color='g', dax=dax)

This is an aspect that we had already discussed here, in PR #95 On my opinion, it shows that we should:

• Recommend eRK4-scipy by default
• Learn how to optimize the time steps size to improve eRK4-homemade
• Never use eRK4-homemade without doing a convergence study with respect to the time step size

[Didou09]

I am still not sure why eRK8-scipy gives wrong results though

[Didou09]

Solved, see PR #130

[Didou09]

I m closing this issue as it was solved