Understanding coupled Runga Kutta derivation

[ognik]

I understand Runga Kutta 4th order for use with a single ODE, but am not confident I have correctly derived how to apply RK4th to a coupled or 2nd order system. As example I have dy/dt=p; dp/dt=-4pi²y
I recognize this as SHM with w=2pi so yes, I don't need RK to solve it, but I want to use this simple situation to understand how to apply RK. The '1st part' seems logical enough taking f(t_n,y_n,p_n)=-4pi²y and dt=h.
Then k₁=h*f(t,y_n,p_n) - effectively f(y_n); kj₂=h*f(y_n+k₁/2); k₃=hf(y_n+k_2/2); k₄=h*f(y_n+k3) and then p_n+1=p_n+(k₁+2k₂+2k₃+k₄)/6.
I choose the origin at t=0 and y(0) = 1 for simplicity. I integrate dp/dt=-4pi²y to give p=-4pi²yt. So at t=0, p₀=0.
Thus far I am feeling comfortable - but please correct/elucidate the above as appropriate.
-----------------------
Now I have dy/dt=p and clearly p_n depends on y_n and t_n. This is where I am not very confident. So dy=pdt and I'm thinking this must be the first 'slope' point, so dy=p_n*h = j₁, and I have p₀=0. I want j(1,2,3,4) to find y_n+1. I am still a little confident so far :-)
But what of j₂=h(p_n+?/2). It seems to me that I must use k₁/2 and go on as j₃=h(p_n+k₂/2) etc. with y_n+1=(j₁+2j₂+2j₃+j₄)/6 - because intuitively that 'couples' y and p.
But I can also see that j_i already includes p - so maybe it should be j₃=h(p_n+j₂/2) etc. ?
I would really like to understand the above, not just get the answer please?

 

[DrClaude]

I The '1st part' seems logical enough taking f(t_n,y_n,p_n)=-4pi²y and dt=h.

But you have 2 f's:
$$
\begin{align}
f_1 ( t_n, y_n, p_n) &= -4 \pi^2 y_n \\
f_2 ( t_n, y_n, p_n) &= p_n
\end{align}
$$

Then k₁=h*f(t,y_n,p_n) - effectively f(y_n); kj₂=h*f(y_n+k₁/2);

You have to modify both ##y_n## and ##p_n## at each step.

A simple way to see this may be to take it as a vector equation,
$$
\mathbf{y}_n = \begin{pmatrix} y_n \\ p_n \end{pmatrix}
$$
and
$$
\mathbf{f}(t_n, \mathbf{y}_n) = \begin{pmatrix} -4 \pi^2 y_n \\ p_n \end{pmatrix}
$$
and apply the single-variable RK to it.

 

[ognik]

I was in fact trying to apply RK to both y_n and p_n at each step, so when I said '1st part' above, that would refer to the bottom component of your f_n vector - and I use k_i to find p_n+1 = p_n + (k1+2k2+2k3+k4)/6. so k₁=h*(-4pi²)*y_n etc.
------------------------------
So unless I have that wrong, its just simultaneously finding the value of y_n at each step that I am unclear on. Using j_i to find y_n+1 I am fairly sure that j1=h*p_n - is that right? Then y_n+1=y_n + (j1+2j2+2j3+j4)/6
What I am not seeing clearly is if j₂=h*(p_n + j₁) etc.?

 

[ognik]

Hi again, having satisfied myself that I have the right approach to apply Runge Kutta 4th order to a coupled system, I wrote a Fortran program to test how I should apply it in practice.

But the program is wrong - the values rapidly diverge. Use a debugger, I can see that the k_i and j_i values I am calculating are the problem, but I can't see why - unless either the textbook has the wrong formula (unlikely) or I am somehow misunderstanding the text. I have uploaded the fortran program - please could someone look at how I am coding the k calculations in the first loop - and tell my what the blindingly obvious is that I am missing? Thanks.

 

[the_wolfman]

You're not applying RK4 correctly to the system of equations. You have to advance both equations at the same time. You're advancing one equation then the other.

Treat the system of ODEs as a vector equation [itex] \frac{d \vec x}{d t}= \vec f\left( \vec x \right)[/itex]
The first step of RK4 is [itex] \vec k_1= \Delta t \vec f \left( \vec x^n \right)[/itex]
The second step is [itex] \vec k_2= \Delta t \vec f \left( \vec x^n+\vec k_1/2 \right)[/itex]
The third step is [itex] \vec k_3= \Delta t \vec f \left( \vec x^n+\vec k_2/2 \right)[/itex]
The forth step is [itex] \vec k_4= \Delta t \vec f \left( \vec x^n+\vec k_3 \right)[/itex]
Then the final step is to sum everythin [itex] \vec x^{n+1} = \vec x^n + \vec k_1/6 + \vec k_2/3 + \vec k_3/3 + \vec k_4/6[/itex]

Every step is a vector equation. You have a second order system, so you have to update 2 equations at each step.

 

[ognik]

Hi - I'm not quite following, I thought I was advancing both eqtns for each step=h? I have to admit that I have been unsure of this aspect all along - I can follow how to use RK on a 1st Order ODE, and I follow hopw to split the 2nd order ODE into 2 ist orders. I grasp at an intuition level that the 2 split eqtns are dependent so must be advanced together ...but it seems my thinking isn't quite right. Which is that I have f(y)=-4pi^2y, and I use what you have just written to advance that. But because P(y) depends on f(y), I can only advance P once I have advanced y...please correct me where I'm wrong?

 

[DrClaude]

Which is that I have f(y)=-4pi^2y, and I use what you have just written to advance that. But because P(y) depends on f(y), I can only advance P once I have advanced y...please correct me where I'm wrong?

But y depends on p, so how can you advance y? They are both interdependent, which is why we call the differential equations coupled.

As both the wolfman and I said previously, treat the dependent variables as a vector, and update the full vector at each of the intermediate steps in the RK algorithm.

 

[the_wolfman]

I've added comments to your code to show where you go wrong.

k1=temp*yN             
k2=temp*(yN+(k1/2.0)) !k1 estimates how pN changes. Here you treat it as the change in yn!
...
K=(k1+(2.0*k2)+(2.0*k3)+k4)/6.0_dp
pN=pN+K  !Here you correctly use k to update pN

            
j1=temp*pN             
j2=temp*(pN+(j1/2.0)) !J1 estimates how yN changes. Here you treat it as the change in pN!
...
J=(j1+(2.0*j2)+(2.0*j3)+j4)/6.0_dp
yN=yN+J !Here you correctly use j to update yN

The fix is to do something like this

k1=temp*yN               !Update k1 and j1 together
j1=temp*pN             

k2=temp*(yN+(j1/2.0)) !Here j1 updates yN
j2=temp*(pN+(k1/2.0)) !Here k1 updates pN
!k2 depends on j1 and j2 depends of k1. You have to update both together
...
K=(k1+(2.0*k2)+(2.0*k3)+k4)/6.0_dp
J=(j1+(2.0*j2)+(2.0*j3)+j4)/6.0_dp
pN=pN+K
yN=yN+J

Does this make sense?

 

[ognik]

Thank you, I can see clearly how to do the problem now. Could I please impose on either of you further though - I'd really like to understand this a bit better (I did understand that the 2 eqtns were interdependent, and that they could be treated as a vector eqtn, but just couldn't see the interplay between j_i and k_i).

I can now see how that works, and can use that, but I unhappily fall short of understanding the why - is there perhaps a sketch or something that might help me?
Thanks again.

 

[the_wolfman]

In (explicit) Runga-Kutta methods you start by knowing the derivative at the beginning of the time step, and you use this value to estimate the total change in p and y across that time step. Then you use this estimate to make a better estimate for the change in p and y. Which is then used to make a even better estimate... etc.

In doing so, RK methods have to estimate the derivative at intermediate times. However the derrivates in general depend on p(t), y(t), and t. Therefore in-order to accurately approximate the derriavtes at intermediate times you have to know p,y, and t and these times. To do this you have to update p and y simultaneously.

If you google runge-kutta sketch. There are many diagrams illustrating how RK4 evaluates the derivate at multiple points.

 

[ognik]

Hi Wolfman, I can only find sketches for a 1st order ODE application (it may not be possible to sketch what happens with a coupled or 2nd order ODE system?). To rephrase what I am not 'getting':
1) pN is associated with K, so that pN=pN+K
In the code above, you have k1=temp*yN, should that not be k1=pN?
and then k2=(pN+(j1/2.0)) ?
2) You say "k2 depends on j1 and j2 depends on k1" - this is what I am not seeing properly, I mean it makes sense when I read it, but why does it work?

Thanks
Alan

 

[ognik]

Hi again, came across something that makes a little more sense to me, see http://www.phy.davidson.edu/FacHome/dmb/py200/RungeKuttaMethod.htm
It confirms point 1 above, one of the functions (p in mine) is always associated with the k factors.
However it shows that k_i always updates the same variable, as does j_i, so if I had f(t,yn,pn) then k₂ would be =
h*f(t+dt, y_n+k₁/2, p_n + j₁/2)
I think this is the same is what you have been saying, but having the k & j apply to variables within functions (and not to the functions themselves) somehow makes more sense to me. Question is now - have I made sense? :-)

Equally importantly is that my adjusted code still diverges very quickly, I have extracted just the k,j part below - my head is just foggy on this by now, can anyone see what is wrong to make it diverge? I will upload a screenshot of the first cycle in debugger.txt, plus latest pgm.txt.
All help appreciated!

    tN=0.0_dp  ! initial conditions
     yN=1.0_dp
     pN=0.0_dp

Do iter=1,numsteps
      tN=tN+h
       k1=h*pN           ! Update p and y together
       j1=h*const*yN         

       k2=h*(pN+(j1/2.0))
       j2=h*const*(yN+(k1/2.0))   !no t in y(n), so no t+h terms

       k3=h*(pN+(j2/2.0))
       j3=h*const*(yN+(k2/2.0))

       k4=h*(pN+j3)
       j4=h*const*(yN+k3)
   
       K=(k1+(2.0*k2)+(2.0*k3)+k4)/6.0_dp
       J=(j1+(2.0*j2)+(2.0*j3)+j4)/6.0_dp
       pN=pN+K   
       yN=yN+J
...
  end do

 

[the_wolfman]

Sorry if I had a typo in the code I sent you. I think I missed a factor of -4 pi squared.

The problem in your current version of the code is in the last 2 steps. You should add K to yN and J to pN.

I think you're starting to understand how to implement RK. But I'm not sure that you understand how RK works (for one equation or multiple equations). For instances what do jn and kn represent? Can you tell me what each step is doing. I think understanding how RK works will help you understand the implementation.

 

[ognik]

Let me see if I can explain my understanding in my own words: We are working with the slope of a curve. For each point we reach, we try and predict the net point on the curve - using yn+1 = yn + (1/6)(k1 + 2k2 + 2k3 + k4). This is a weighted average of possible next points, with the midpoints being 2 x as important as the end points. By using previous k values to calculate the next k value, there is a built-in correction. Is that close for a 1st order ODE?
For a 2nd order system, we have more than one parameter influencing the slope, so at each step, the slope of each parameter must be predicted separately but simultaneously. So in a general case if we had F(t,x,y) and G(t,x,y), then each step would predict for F(t+dt, x+dx, y+dy) and similar for G. For x we would always k, but the next k would depend on the previous j so take into account the interdependence of the parameters of the 2 single ODE functions. Similarly we would use j for predicting the next y, and j would depend on previous k.
t increases as t+h...
If I have that close, what I don't see is how 'mixing' k and j like that is a satisfactory representation of the interdependence?
In terms of the program, if you are correct about the last 2 steps, then I am still missing some understanding. I start with:
k1=h*pN
j1=h*const*yN
k_i always predicts pN, so why add K to yN?
Finally, I think the program problem is somewhere in the k_i and j_i calculations, they diverge very rapidly, after only a couple of cycles ...
Really appreciate your patience with this, thanks

 

[the_wolfman]

A better name for [itex] k_1 [/itex] is [itex] \Delta y_1 [/itex], a better name for [itex] k_2[/itex] is [itex]\Delta y_2[/itex], etc.
Similarly a better name for [itex] j_1 [/itex] is [itex] \Delta p_1 [/itex], etc.

Here [itex] \Delta y_i = \Delta t \frac{dy}{dt} [/itex] and [itex] \Delta p_i = \Delta t \frac{dp}{dt} [/itex], where the derivatives are evaluated at different locations depending on i.

In order to calculate the derivative, we have to know y, p, and t at an intermediate location. In general [itex] \Delta y_i [/itex] depend [itex] \Delta y_{i-1} [/itex] and [itex] \Delta p_{i-1} [/itex], and similarly [itex] \Delta p_i [/itex] depend [itex] \Delta y_{i-1} [/itex] and [itex] \Delta p_{i-1} [/itex].

In your particular problem, [itex] \frac{dy}{dt} [/itex] only depends on [itex] p [/itex] and [itex] \frac{dp}{dt} [/itex] only depends on [itex] y [/itex]. This is way it looks like we are "mixing" the [itex] k [/itex] and [itex] j [/itex] ( [itex] \Delta y [/itex] and [itex] \Delta p[/itex] ). But this problem is just a special case.

You won't get this impression if you use RK4 to solve a more general problem where [itex] \frac{dy}{dt} [/itex] and [itex] \frac{dp}{dt} [/itex]
both depend on [itex] y [/itex] and [itex] p [/itex].

The problem in your code is the two lines

pN=pN+K  
yN=yN+J

You should have

pN=pN+J
yN=yN+K

If you add that value of K in you output file to y(0) you will get 0.891 in approximate agreement with you exact solution.

 

[ognik]

Thanks again Wolfman, I am frustrated to be taking so long to 'get it'...I am one of those who have to fully understand something before I am happy to move on; my Mother always said that being a perfectionist would cause me problems, but never said what to do about it :-)

Naturally, once I fixed those 2 lines, the program ran well - although I was surprised with a lower than expected accuracy - EG: the error (cf exact value) with step size 0.075 was 0.007135691 at t(n)=2.8500, in general around 0.003. At much smaller step sizes I obviously got much better accuracy, but having already experimented with Adams-Bashforth 4 step (only with 1st order ODEs though), it seems that AB-4 step may be at least as accurate as Runge-Kutta which I wasn't expecting ...browsing I found that both methods expect a local error of O(h^5), and global error of O(h^4)...but there is also stability and efficiency to be taken into account... Question: how in practice does one choose one method over another, without getting sidetracked from the actual problem you want to use the method on?

Yes - and thanks, using Δy1 & Δp1 are much less confusing and I can see how the special dependencies in this case lead to some of my initial confusion, once I grasped that k & j were advancing the variables and not the function directly, it became a lot clearer.

Finally, I found sketches like this one to be useful - http://www.google.co.nz/imgres?imgu...ZVbOQAsermAXPioDwAg&tbm=isch&ved=0CCQQMygFMAU
It would have been nice (for me :-)) to have an equivalent geometric picture for the coupled system, showing both dy and dp points, but couldn't find any. Don't suppose anyone has an idea of how to do something like that, maybe in Mathmatica ... ?

 

[the_wolfman]

If you want to compare AB4 and RK4 you should really use the same equation. Otherwise you're comparing apples and oranges. In general there will be some problems where AB4 excels and other problems where RK4 excels. I recommended solving this system of equations using AB4. The generalization of AB4 to systems of equations is similar to the generalization of RK4. And implementing this code is a good test of how much you truly grasp.

Finding an optimal numerical method to solve a dynamical system of equations is not trivial. Many people dedicate their careers to this task. In general it requires understanding both the dynamical system and many different types of numerical methods. Often one can identify traits important to the dynamical system that are particularly difficult to capture numerically. And then we pick numerical methods that excel at representing said traits. We also model problems to help us quickly compare multiple methods.

 

[ognik]

Thanks again, all much clearer now :-)