How can I use Runge Kutta to solve coupled first order DE's in Matlab?

[HoosierDaddy]

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 n_t = p(t) + q(t). Here, n_t = 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!

 

[D H]

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?

 

[HoosierDaddy]

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

y_j+1 = y*_j+1 + \frac{h}{2}(k₂ - k₁)

Where y*_j+1 = y_j + hk₁

and k₁ = f(x_j, y_j) and k₂ = f(x_j+1, y*_j+1)

Does this identify which RK2 method I am attempting to use?

 

[D H]

That is Heun's method.

An alternative is the midpoint method:

\aligned<br /> y_{j+1} &amp;= y_j + hk_{j+1/2} \\<br /> x_{j+1/2} &amp;= x_j +h/2 \\<br /> y_{j+1/2} &amp;= y_j + hk_j/2 \\<br /> k_j &amp;= f(x_j,y_j) \\<br /> k_{j+1/2} &amp;= f(x_{j+1/2},y_{j+1/2})<br /> \endaligned

 

[HoosierDaddy]

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. x_j = t_j, x_j+1 = t_j+1 = t_j + \Deltat?

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

 

[D H]

Consider this set of six differential equations:

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

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

\aligned<br /> \mathbf u &amp;= \left[ r_x \, r_y \,r_z \,\,v_x \,v_y \,v_z \right]^T \\<br /> \frac{\mathbf u}{dt} &amp;= \left[ v_x \,v_y \,v_z \,\,a_x \,a_y \,a_z \right]^T<br /> \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.

 

[TrentFGuidry]

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

y_j+1 = y*_j+1 + \frac{h}{2}(k₂ - k₁)

Where y*_j+1 = y_j + hk₁

and k₁ = f(x_j, y_j) and k₂ = f(x_j+1, y*_j+1)

Does this identify which RK2 method I am attempting to use?

Are you sure it is k₂ - k₁ and not k₂ + k₁?

The full Heun is

k₀ = f(x_i, y_i)
k₁ = f(x_i + h, y_i + hk₀)

y_i+1 = y_i + (1/2 k₀ + 1/2 k₁)h

 

[D H]

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