Solving a system of second order PDEs using Runge Kutta in C

Question

I have a problem solving a system of differential equations using the Runge Kutta algorithm. So far I have rewritten the second order PDE into a set of two coupled equations where

    f(L1,L2) = L2
    g(L1,L2) = A*(B*L1-C*L2-D)

are the two equations and A, B, C and D are constants. In order to get the value for the next step, I proceeded as follows for each time step dt:

    k1 = f(L1,L2)
    l1 = g(L1,L2)

    k2 = f(L1 + 0.5 * dt * k1,L2 + 0.5 * dt * l1 )
    l2 = g(L1 + 0.5 * dt * k1,L2 + 0.5 * dt * l1 )

    k3 = f(L1 + 0.5 * dt * k2,L2 + 0.5 * dt * l2 )
    l3 = g(L1 + 0.5 * dt * k2, L2 + 0.5 * dt * l2 )

    k4 = f(L1 + dt * k1,L2 +  dt * l1 )
    l4 = g(L1 + dt * k1,L2 + dt * l1 )

Where I use the values for L1 and L2 of the current time step and calculate the coefficients iteratively.

As a result I get L1 and L2 by summing up and weighting the coefficients at the end. My problem is, that the whole algorithm becomes unstable after 4 time steps.

Does anybody know if the realization is technically correct? Thanks!

L1 and L2 are values. for the first step, I used the starting values here. — MichaelScott, Nov 10 '13 at 13:19
And by *values* you mean constants, right? What starting values do you use? Can you show the output of your integration for a couple of steps? — Jonas Bötel, Nov 10 '13 at 16:49

score 0 · Answer 1 · answered Nov 10 '13 at 09:40

0

just a guess, since you don't say what dt value you use: keep it small as possible, because

local truncation error is on the order of O(h^5), while the total accumulated error is order O(h^4).

(cited from this wikipedia article, dt plays the h role).

answered Nov 10 '13 at 09:40

CapelliC

59,646
5
47
90

Thanks for your answer. The thing is that it does not depend on dt (time inkrement) in my case. There is no big difference in the values if I put dt=0.001 or 0.01. That is why I think it is wrong...My question is, if that is the correct way to implement Runge Kutta since I have no idea how to test it differently... – MichaelScott Nov 10 '13 at 09:46

score 0 · Answer 2 · answered Nov 10 '13 at 09:46

Two things:

Runge-Kutta in general is not stable. It is just "more stable" than Euler's. Depending on the condition of your differential equation and the dt it might just not be sufficient. Does a smaller dt help?

I'm missing the notion of t in your definition of f and g. Assuming L1 and L2 are not constant in t you better pass it through f and g. Like f(t,L1,L2). This forces you to think about those coefficient calculations, where you now need to pass another t' accordingly. This will lead to the evaluation of L1 and L2 at midpoints.

score 0 · Answer 3 · answered Jun 15 '18 at 07:17

In the computation of k4, l4 you used k1, l1 to compute the state where the ODE function is evaluated, while the method demands to use k3, l3. This error will reduce the order of the method, probably to 2, and the dependence of the error on the step size will reflect that.

Solving a system of second order PDEs using Runge Kutta in C

3 Answers3