Numerical Integration of 2nd Order DE

strokebow

I am intending to use Runge Kutta 4th order to numerically solve a system of coupled equations:

[itex]\frac{d^{2}x}{dt^{2}}[/itex] = K₁ * x * cos(t) + ( (K₂ * [itex]\frac{dy}{dt}[/itex]) - [itex]\frac{dz}{dt}[/itex] )

[itex]\frac{d^{2}y}{dt^{2}}[/itex]= -K₁ * y * cos(t) + ( (K₂ * [itex]\frac{dz}{dt}[/itex]) - [itex]\frac{dx}{dt}[/itex] )

[itex]\frac{d^{2}z}{dt^{2}}[/itex]= ( (K₂ * [itex]\frac{dx}{dt}[/itex]) - [itex]\frac{dy}{dt}[/itex] )


I'm really a bit stuck to be honest. I've only ever used RK4 on 1st order linear ODEs. I've been reading around alot but not making much progress.

Initial values for [itex]\frac{dx}{dt}[/itex], [itex]\frac{dy}{dt}[/itex], [itex]\frac{dz}{dt}[/itex] are all known. The constants K are known.

Can anyone please help? Thanks

bigfooted

Introduce the derivatives as new variables, say u=dx/dt and solve the system of 6 first order equations, where the second derivatives are replaced by the first derivatives du/dt etc.

so the last equation will be
dw/dt = k2*u - v
but you could also write it as:
dw/dt = k2*dx/dt - dy/dt

They will lead to the same solution, but their convergence behavior will be slightly different.

strokebow

Thanks very much for your input.

So I would define:
u = [itex]\frac{dx}{dt}[/itex]
v = [itex]\frac{dy}{dt}[/itex]
w = [itex]\frac{dz}{dt}[/itex]
and
[itex]\frac{du}{dt}[/itex] = [itex]\frac{d^{2}x}{dt^{2}}[/itex]
[itex]\frac{dv}{dt}[/itex] = [itex]\frac{d^{2}y}{dt^{2}}[/itex]
[itex]\frac{dw}{dt}[/itex] = [itex]\frac{d^{2}z}{dt^{2}}[/itex]

Therefore my equations become:
(1'): [itex]\frac{du}{dt}[/itex] = K[itex]_{1}[/itex] * x * cos(t) + ( (K[itex]_{2}[/itex] * v) - w )
(2'): [itex]\frac{dv}{dt}[/itex] = -K[itex]_{1}[/itex] * y * cos(t) + ( (K[itex]_{2}[/itex] * w) - u )
(3'): [itex]\frac{dw}{dt}[/itex] = ( (K[itex]_{2}[/itex] * u) - v )

Now for applying RK4 . . .
How would I do that simultaneously for both. I've found a couple of exams on the net but its quite difficult to follow since all the equations are inter-related.
How would I go about this for this case?

Any help is much appreciated!

strokebow

Not sure what you mean by this big foot?

DrClaude

You now have a system of 6 1st-order equations to solve, (1')-(3') plus
(4'): [itex]\frac{dx}{dt} = u[/itex]
(5'): [itex]\frac{dy}{dt} = v[/itex]
(6'): [itex]\frac{dz}{dt} = w[/itex]

strokebow

Thanks for your input.

I appreciate that we now have a system of 1st order DE's. Its just that trying to apply RK4 to this system concurrently is quite confusing.
I'm trying to follow the example given on the net (e.g., http://www.nsc.liu.se/~boein/f77to90/rk.html)

But am finding it tough to implement.

Especially since these equations are coupled with more than 2 variables. Like with eqn (1') its depending on both v and w.

DrClaude

Maybe this formulation will be more useful. Set a vector [tex]\mathbf{a} = (u,v,w,x,y,z)[/tex] and a vector function [tex]\mathbf{f}(\mathbf{a},t) = (\frac{du}{dt},\frac{dv}{dt},\frac{dw}{dt},\frac{dx}{dt},\frac{dy}{dt}, \frac{dz}{dt})[/tex] The 4th-order RK algorithm then reads
[tex]
\mathbf{k}_1 = h \mathbf{f}(\mathbf{a}_n, t_n)
[/tex]
[tex]
\mathbf{k}_2 = h \mathbf{f}(\mathbf{a}_n + \frac{\mathbf{k}_1}{2}, t_n + \frac{h}{2})
[/tex]
[tex]
\mathbf{k}_3 = h \mathbf{f}(\mathbf{a}_n + \frac{\mathbf{k}_2}{2}, t_n + \frac{h}{2})
[/tex]
[tex]
\mathbf{k}_4 = h \mathbf{f}(\mathbf{a}_n + \mathbf{k}_3, t_n + h)
[/tex]
[tex]
\mathbf{a}_{n+1} = \mathbf{a}_n + \frac{\mathbf{k}_1}{6} + \frac{\mathbf{k}_2}{3} + \frac{\mathbf{k}_3}{3} + \frac{\mathbf{k}_4}{6}
[/tex]

Hope this helps.

strokebow

Thanks for your help.

So (replacing the 'k's' from the RK4 algorithm with 'm's' and 'n's')

For the first step, would this be correct:
m1u = h x (1')
m1v = h x (2')
m1w = h x (3')

n1u = h x (4')
n1v = h x (5')
n1w = h x (6')

?

DrClaude

I'm not sure I follow you.

The "first step" is the calculation of [itex]\mathbf{k}_1[/itex]. The first element of [itex]\mathbf{f}(\mathbf{a}_n,t_n)[/itex] is
[tex]
\left.\frac{du}{dt}\right|_{t_n} = K_1 x_n \cos(t_n) + ( K_2 v_n - w_n )
[/tex]
so the first element of [itex]\mathbf{k}_1[/itex] is
[tex]
h \times \left[ K_1 x_n \cos(t_n) + ( K_2 v_n - w_n ) \right]
[/tex] and so on. For instance, the 4th element of [itex]\mathbf{k}_1[/itex] is
[tex]
u_n
[/tex]
Once you have calculated all six elements of [itex]\mathbf{k}_1[/itex], the "second step" is do the the same for [itex]\mathbf{k}_2[/itex], which will depend on the calculation of [itex]\mathbf{k}_1[/itex]. Once you have calculated all [itex]\mathbf{k}[/itex]'s, you can advance the solution one time step with the formula for [itex]\mathbf{a}_{n+1}[/itex].

The easiest way to implement this numerically is not to program it explicitly, bu to write a function that given an array of six values [itex](u,v,w,x,y,z)[/itex] and a value of [itex]t[/itex], returns an array of the six values [itex](du/dt,dv/dt,dw/dt,dx/dt,dy/dt,dz/dt)[/itex]. The RK integrator then just makes 4 calls to this function to calculate [itex]\mathbf{k}_1[/itex]-[itex]\mathbf{k}_4[/itex], and then the new values of [itex](u,v,w,x,y,z)[/itex].

strokebow

Thanks Dr C. Appreciate you taking time out to help me with this.

So, there are 6 equations for each step. i.e., k[itex]_{1}[/itex] has 6 equations. Writing them explicitly (I take your point that a generic function is best for programming but for now I am just trying to understand this solution):
[itex]k_{1}u = h (K_1 x \cos(t) + ( K_2 v - w ))
[/itex]

[itex]k_{1}v = h (-K_1 y \cos(t) + ( K_2 v - w ))
[/itex]

[itex]k_{1}w = h (K_2 u - v)
[/itex]

[itex]k_{1}u' = ?
[/itex]

[itex]k_{1}v' = ?
[/itex]

[itex]k_{1}w' = ?
[/itex]

When you said:
Do you mean that we just consider these (4th, 5th and 6th elements of the vector) to be constants?

Because we can derive expressions for u, v and w from (1'), (2') and (3')

DrClaude

Lets take it one step back. You converted your system of 3 2nd-order differential equations into a system of 6 first-order differential equations. You have increased the number of variables to 6, but that is a small price to pay compared with the gain 2nd to 1st order.

The equations now read:
[tex]\frac{du}{dt} = K_{1} x \cos(t) + ( K_{2} v - w )[/tex]
[tex]\frac{dv}{dt} = -K_{1} y \cos(t) + ( K_{2} w - u )[/tex]
[tex]\frac{dw}{dt} = (K_{2} u - v )[/tex]
[tex]\frac{dx}{dt} = u[/tex]
[tex]\frac{dy}{dt} = v[/tex]
[tex]\frac{dz}{dt} = w[/tex]

The first step is to calculate the elements of [itex]\mathbf{k}_1[/itex], which I will call [itex](u',v',w',x',y',z')[/itex].
[tex] u' = h \left[ K_{1} x_n \cos(t_n) + ( K_{2} v_n - w_n ) \right] [/tex]
[tex] v' = h \left[ -K_{1} y_n \cos(t_n) + ( K_{2} w_n - u_n ) \right] [/tex]
[tex]w' = h (K_{2} u_n - v_n )[/tex]
[tex]x' = h u_n[/tex]
[tex]y' = h v_n[/tex]
[tex]z' = h w_n[/tex]

Now, you need to calculate [itex]\mathbf{k}_2 = (u'',v'',w'',x'',y'',z'')[/itex].
[tex] u'' = h \left[ K_{1} \left( x_n + \frac{x'}{2}\right) \cos(t_n + \frac{h}{2}) + ( K_{2} \left(v_n + \frac{v'}{2} \right) - \left(w_n + \frac{w'}{2} \right) ) \right] [/tex]
and so on for the other variables. Finally, once you have all values [itex]u',u'',u''',u''''[/itex], you find
[tex]
u_{n+1} = u_n + \frac{u'}{6} + \frac{u''}{3} + \frac{u'''}{3} + \frac{u''''}{6}
[/tex]
and so on for the other variables.

strokebow

Dr C. Thanks very much. That is very clear.

I'm going to have a go at writing it all out explicitly and then I'll post it on here.

Cheers!!