Problems with Runge-Kutta method

  • Thread starter matteo86bo
  • Start date
  • #1
60
0
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:

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

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
 

Answers and Replies

  • #2
335
44
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.
 
  • #3
60
0
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?
 
  • #4
335
44
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 :)
 
  • #5
60
0
thank you man!
I messed up with the grid ... the sign is my fault, but it doesn't matter ...
 
  • #6
335
44
Glad it helped :)
 

Related Threads on Problems with Runge-Kutta method

  • Last Post
Replies
1
Views
6K
  • Last Post
Replies
3
Views
871
  • Last Post
Replies
1
Views
4K
  • Last Post
Replies
6
Views
13K
Replies
4
Views
2K
  • Last Post
Replies
0
Views
8K
  • Last Post
Replies
2
Views
2K
  • Last Post
Replies
4
Views
2K
  • Last Post
Replies
1
Views
2K
  • Last Post
Replies
5
Views
6K
Top