Problems with Runge-Kutta method

[matteo86bo]

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:

\frac{d\beta(r,t)}{dt}=A\beta(r,t)^N

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 ...

P.S.:

beta(i,j) mean beta at r_i and t_j
so i and j referes to radius and time grid respectively

 

[Lord Crc]

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.

 

[matteo86bo]

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?

 

[Lord Crc]

Like I said I got a bit confused by your notation. Let's step back a bit. The ODE is given by
\frac{dy(t)}{dt} = f(y(t), t)

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

The midpoint method is then

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

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

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

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 :)

 

[matteo86bo]

thank you man!
I messed up with the grid ... the sign is my fault, but it doesn't matter ...

 

[Lord Crc]

Glad it helped :)