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

Problems with Runge-Kutta method

  1. May 16, 2010 #1
    Hi everyone,
    I've a little problem with my algorithm.
    I have two grids: one for time and the other one for the radius.
    I need to evaluate this equation:


    I try solving this with a simple Runge Kutta II method, but I'm not convinced at all that this works properly.

    Let's check my code:

    k1 = deltaT*(-A*beta(i,j-1)**N
    k2 = deltaT*(-A*(beta(i,j-1)+k1/2.)**N)

    beta(i,j) = beta(i,j-1) + k2

    What do you think about that?
    This is not the simplest case, just the one you can read in Numerical Recipes ...


    beta(i,j) mean beta at r_i and t_j
    so i and j referes to radius and time grid respectively
  2. jcsd
  3. May 17, 2010 #2
    I'm having difficulties making sense of your code. It looks like you're trying to the midpoint method, but in that case the evaluation of beta when computing k2 is all wrong. When you take the half-step you will, in general, not end up on a grid vertex, so you can't use "beta(i,j-1)" there.
  4. May 17, 2010 #3
    Should it look like that?

    k2 = deltaT*(-A*( (beta(i,j-1)+beta(i,j))/2. +k1/2.)**N)

    i know i need to evaluate beta in t+deltaT/2.. can this work?
  5. May 18, 2010 #4
    Like I said I got a bit confused by your notation. Let's step back a bit. The ODE is given by
    [tex]\frac{dy(t)}{dt} = f(y(t), t)[/tex]

    In your case we have [itex]y(t) = \beta(r, t)[/itex] and [itex]f(y(t), t) = Ay(t)^N[/itex]. Note the use of y(t) in the right hand of f.

    The midpoint method is then

    [tex]\tilde{y}_{j+1} = y_j + \frac{h}{2}f(y_j, t_j)[/tex]
    [tex]y_{j+1} = y_j + hf(\tilde{y}_{j+1}, t_j + \frac{h}{2})[/tex]

    where h is the step length ([itex]\Delta t[/itex]). Inserting the expression for f we get

    [tex]\tilde{y}_{j+1} = y_j + \frac{h}{2}f(y_j, t_j) = y_j + \frac{h}{2}Ay_j^N[/tex]
    [tex]y_{j+1} = y_j + hf(\tilde{y}_{j+1}, t_j + \frac{h}{2}) = y_j + hA\tilde{y}_{j+1}^N[/tex]

    Here y_j = beta(i,j-1) using your notation. So in the second line we don't access the original beta, instead we use the midpoint we found. Also note that you seem to have flipped signs on A between the ODE and the code, while I haven't (not sure which sign is the correct one).

    Hope this helps, and I hope I didn't mess up somewhere :)
  6. May 18, 2010 #5
    thank you man!
    I messed up with the grid ... the sign is my fault, but it doesn't matter ...
  7. May 19, 2010 #6
    Glad it helped :)
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook