Possible issue in RungeKutta::step() stage-time calculation

Issues
1

Hi all,

I am sitting trying to implement a fully coupled RK with Coulomb forces and noticed in the step function of runge_kutta.h something I can’t quite understand, due it being a RungeKutta utility with a derivative that depends explicitly on time.

In runge_kutta.h, the stage state and stage time are calculated as:


for(int j = 0; j < i; ++j) {

yt += h_ * tableau_(i, j) * k.row(j);

tt += tableau_(i, j);

}

k.row(i) = function_(tt, yt);

Should the second line instead be:


tt += h_ * tableau_(i, j);

For classical RK4 with h_ = 0.05 ns, the current code passes callback times


t, t + 0.5, t + 0.5, t + 1

while I would expect


t, t + 0.025 ns, t + 0.025 ns, t + 0.05 ns

The cited Fehlberg formulation evaluates the derivative at x_0 + \alpha_k h, which appears to support multiplying the tableau coefficient by the step size.

Is the first argument of StepFunction intended to be the absolute independent-variable value, or is it meant to represent a normalized within-step stage coordinate?

All the best,
Rickard

2

It seems as this is not something that affects the propagation, since the first argument in

[&](double, const Eigen::Vector3d& cur_pos) -> Eigen::Vector3d

seems to be deliberately unnamed and ignored. But I would appreciate confirmation on that.
Though it might seems small, I am also interested in why the mixing of unit-less tableu values and time units was implemented.

All the best,
Rickard

3

Hi @rickard ,

I completely agree with you:

a) Adding unit-less numbers to the time stamp does not make sense. I’ll have a more thorough look at this, but as you say, it seems like here the tableau value should be multiplied by the time step h_.

b) As of now, this likely has not been spotted as it does not have any effect due to the first argument of the listed function options do not use this. Nevertheless it should be fixed.

Thank you very much for reporting this and digging into it! Would you prefer to file a merge request to our repository yourself, else we can implement it.

Cheers
Paul

4

Hi @pschutze ,

No problem, I found it by accident during debugging of an implementation of first time-step refinement for the initial Coulomb expansion.

Since my current branch contains many unrelated research changes and I’m not so comfortable with Git to safely isolate this into a clean merge request, perhaps the best way is for you maintainers to implement it.

All the best,
Rickard

5

Hi @rickard ,

no problem, we can take care of this! :slightly_smiling_face:

Cheers
Paul