# Help! Programming Runge Kutta

1. Oct 8, 2009

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:

$$\frac{dp}{dt}$$ = aq - bp

$$\frac{dq}{dt}$$ = -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. Oct 8, 2009

### D H

Staff Emeritus
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?

3. Oct 8, 2009

Not sure exactly whick RK approach it is, we will be incrementing time such that ...

yj+1 = y*j+1 + $$\frac{h}{2}$$(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?

4. Oct 8, 2009

### D H

Staff Emeritus
That is Heun's method.

An alternative is the midpoint method:

\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

5. Oct 8, 2009

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 + $$\Delta$$t?

If this is the case, what would my yj be? Would I have a separate yj for my p and q equations?

6. Oct 8, 2009

### D H

Staff Emeritus
Consider this set of six differential equations:

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

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

\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

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
7. Oct 9, 2009

### TrentFGuidry

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

8. Oct 9, 2009

### D H

Staff Emeritus
That is exactly what he is doing, just written a little differently.