Help! Programming Runge Kutta

Main Question or Discussion Point

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!
 

Answers and Replies

D H
Staff Emeritus
Science Advisor
Insights Author
15,329
681
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?
 
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?
 
D H
Staff Emeritus
Science Advisor
Insights Author
15,329
681
That is Heun's method.

An alternative is the midpoint method:

[tex]\aligned
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})
\endaligned[/tex]
 
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?
 
D H
Staff Emeritus
Science Advisor
Insights Author
15,329
681
Consider this set of six differential equations:

[tex]\aligned
\frac {d\mathbf r}{dt} &= \mathbf v \\
\frac {d\mathbf v}{dt} &= \mathbf a_{\text{ext}} \equiv \mathbf F_{\text{ext}}/m
\endaligned[/tex]

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

[tex]\aligned
\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
\endaligned[/tex]

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:
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?
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
 
D H
Staff Emeritus
Science Advisor
Insights Author
15,329
681
That is exactly what he is doing, just written a little differently.
 

Related Threads for: Help! Programming Runge Kutta

  • Last Post
Replies
1
Views
6K
  • Last Post
Replies
7
Views
2K
  • Last Post
Replies
4
Views
2K
  • Last Post
Replies
16
Views
1K
  • Last Post
Replies
5
Views
6K
  • Last Post
Replies
6
Views
19K
  • Last Post
Replies
2
Views
2K
  • Last Post
Replies
7
Views
754
Top