Runge–Kutta of 4th order

Question

Lets consider I have a system of 4 ODE: dX/dt = F(X), where X is a vector(4-demension) and F: R^4 -> R^4. F is called vectorDE_total_function, and I'm trying to calculate the solution using RK-4.

def solvingDES():
    previous_vector = np.array ([theta_1, omega_1, theta_2, omega_2]);
    for current_time in time:
        temp_vector = previous_vector;
        RK_vector = np.array([0.0,0.0,0.0,0.0]);
        for c in [6,3,3,6]:
            RK_vector = vectorDE_total_function(previous_vector + c * RK_vector/6) * time_step;
            temp_vector += RK_vector / c;
        previous_vector = temp_vector;
        current_time += 1;

And it looks like I'm wrong somewhere, but I'm not sure where. Is it seems legit?

score 1 · Accepted Answer · answered Dec 11 '16 at 08:47

1

It is a strange, seldom seen way to implement it and it only works for classical RK4, other Runge-Kutta methods would not work like that. But the general idea seems correct.

You have a common error in an usually unexpected place. Setting

temp_vector = previous_vector;

and later

previous_vector = temp_vector;

you do not make a copy of the vectors, but make both object references share the same vectors. Use

temp_vector = previous_vector.copy();

or

previous_vector = temp_vector[:];

to force copying of the vector data.

answered Dec 11 '16 at 08:47

Lutz Lehmann

25,219
2
22
51

I'm not sure then what '=' does. Is it binding, rather than an assignment? – Mikhail Tikhonov Dec 11 '16 at 11:07
The vector variables are references, or pointers, to the actual memory objects. Thus what `=` is doing is copying the pointers. Probably, but not surely, with some reference counting. -- `+=` also is acting on the object pointed to on the left side, thus unlike `a=a+b` where the new `a` is a new object, in `a+=b` `a` is the same object as before. Thus `temp_vector += RK_vector / c;` does not break the pointer coincidence and thus also "changes" `previous_vector` ("…" because no further action takes place). – Lutz Lehmann Dec 11 '16 at 11:29
So no, it is not a binding, as that keeps two objects synchronized via some event handling mechanism, i.e., a change in object A causing a method of object B to be called. Nothing of that sort here. – Lutz Lehmann Dec 11 '16 at 11:42
I mostly use C-code, so.. You claim that something like b=a, a+=b cause b to change also? I thought b is passed by reference(i.e. not the exact memory address, but the value) – Mikhail Tikhonov Dec 11 '16 at 12:16
Yes, I had exactly that happen in some RK4 code. `ylast=ycurr; ycurr+=h*k;` where `k` is the averaged slope of a RK step, and then wondered why the resulting graph had these strange oscillations in the interpolation from the last to the current value. It needed an explicit addition `ycurr=ycurr+h*k;` to get the correct behavior. – Lutz Lehmann Dec 11 '16 at 13:43
The code, RK4 with adaptive time step and cubic Bezier interpolation, is here: http://stackoverflow.com/a/40866829/3088138 – Lutz Lehmann Dec 11 '16 at 13:48
It helped a little, now I got the oscillations, but they're still strange anyway. I guess that's the code for F-function problem. Thanks, i'll try to fix it. – Mikhail Tikhonov Dec 11 '16 at 14:19

Runge–Kutta of 4th order

1 Answers1