Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help! Programming Runge Kutta

  1. Oct 8, 2009 #1
    Hello All,

    I need to utilize a Runge Kutta second order approach to solve two coupled first order DE's simultaneously given some initial conditions and a conservation relationship.

    The DE's are as follows:

    [tex]\frac{dp}{dt}[/tex] = aq - bp

    [tex]\frac{dq}{dt}[/tex] = -aq + bp

    Where a and b are constants, a=3, b=4, and the conservation relationship is nt = p(t) + q(t). Here, nt = 3.

    I'm having trouble getting the ball rolling with writing the R.K. algorithm for this, will be programming in matlab. I have solved the equations analytically by hand and have found the following:

    p(t) = -0.055972e-3t + 3.055972e-4
    q(t) = -0.055972e-7t + 3.055972

    Any help on how to get started with the Runge Kutta algorithm to solve numerically would be appreciated!
  2. jcsd
  3. Oct 8, 2009 #2

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    You have to do a little work before we can help you. What do you know about Runge Kutta methods? About RK2 in particular? Which particular variant of RK2 (there are more than one) are you supposed to use?
  4. Oct 8, 2009 #3
    Not sure exactly whick RK approach it is, we will be incrementing time such that ...

    yj+1 = y*j+1 + [tex]\frac{h}{2}[/tex](k2 - k1)

    Where y*j+1 = yj + hk1

    and k1 = f(xj, yj) and k2 = f(xj+1, y*j+1)

    Does this identify which RK2 method I am attempting to use?
  5. Oct 8, 2009 #4

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    That is Heun's method.

    An alternative is the midpoint method:

    y_{j+1} &= y_j + hk_{j+1/2} \\
    x_{j+1/2} &= x_j +h/2 \\
    y_{j+1/2} &= y_j + hk_j/2 \\
    k_j &= f(x_j,y_j) \\
    k_{j+1/2} &= f(x_{j+1/2},y_{j+1/2})
  6. Oct 8, 2009 #5
    Thanks D H, I believe we are required to use Heun's method. Can I just draw a parallel to Heun's method with my values? i.e. xj = tj, xj+1 = tj+1 = tj + [tex]\Delta[/tex]t?

    If this is the case, what would my yj be? Would I have a separate yj for my p and q equations?
  7. Oct 8, 2009 #6

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Consider this set of six differential equations:

    \frac {d\mathbf r}{dt} &= \mathbf v \\
    \frac {d\mathbf v}{dt} &= \mathbf a_{\text{ext}} \equiv \mathbf F_{\text{ext}}/m

    In other words, Newton's laws of motion. Another way to look at this is as a six vector:

    \mathbf u &= \left[ r_x \, r_y \,r_z \,\,v_x \,v_y \,v_z \right]^T \\
    \frac{\mathbf u}{dt} &= \left[ v_x \,v_y \,v_z \,\,a_x \,a_y \,a_z \right]^T

    Integration techniques such as Runge Kutta can be used on vectors as well as on scalars. Molecular simulations, for example, essentially work on vectors with thousands of elements or more.
    Last edited: Oct 8, 2009
  8. Oct 9, 2009 #7
    Are you sure it is k2 - k1 and not k2 + k1?

    The full Heun is

    k0 = f(xi, yi)
    k1 = f(xi + h, yi + hk0)

    yi+1 = yi + (1/2 k0 + 1/2 k1)h
  9. Oct 9, 2009 #8

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    That is exactly what he is doing, just written a little differently.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook