Register to reply

Numerical Integration of 2nd Order DE

by strokebow
Tags: integration, numerical, order
Share this thread:
strokebow
#1
Feb13-13, 06:40 AM
P: 101
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] = K1 * x * cos(t) + ( (K2 * [itex]\frac{dy}{dt}[/itex]) - [itex]\frac{dz}{dt}[/itex] )

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

[itex]\frac{d^{2}z}{dt^{2}}[/itex]= ( (K2 * [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
Phys.Org News Partner Science news on Phys.org
'Office life' of bacteria may be their weak spot
Lunar explorers will walk at higher speeds than thought
Philips introduces BlueTouch, PulseRelief control for pain relief
bigfooted
#2
Feb13-13, 11:30 AM
P: 298
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
#3
Feb13-13, 02:56 PM
P: 101
Quote Quote by bigfooted View Post
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.
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
#4
Feb14-13, 04:06 AM
P: 101
Numerical Integration of 2nd Order DE

Quote Quote by bigfooted View Post
They will lead to the same solution, but their convergence behavior will be slightly different.
Not sure what you mean by this big foot?
DrClaude
#5
Feb14-13, 04:35 AM
Sci Advisor
PF Gold
DrClaude's Avatar
P: 1,427
Quote Quote by strokebow View Post
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?
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
#6
Feb14-13, 04:54 AM
P: 101
Quote Quote by DrClaude View Post
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]
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
#7
Feb14-13, 07:03 AM
Sci Advisor
PF Gold
DrClaude's Avatar
P: 1,427
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
#8
Feb14-13, 07:30 AM
P: 101
Quote Quote by DrClaude View Post
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.
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
#9
Feb14-13, 07:50 AM
Sci Advisor
PF Gold
DrClaude's Avatar
P: 1,427
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
#10
Feb14-13, 08:17 AM
P: 101
Quote Quote by DrClaude View Post
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].
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:
Quote Quote by DrClaude View Post
For instance, the 4th element of [itex]\mathbf{k}_1[/itex] is
[tex]
u_n
[/tex]
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
#11
Feb14-13, 08:47 AM
Sci Advisor
PF Gold
DrClaude's Avatar
P: 1,427
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
#12
Feb14-13, 08:56 AM
P: 101
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!!


Register to reply

Related Discussions
Numerical integration, use given random distribution as integration points Calculus 12
Second Order Numerical Integration w/ Neumann Boundary Conditions Differential Equations 1
Numerical analysis (composite numerical integration) Engineering, Comp Sci, & Technology Homework 9
Numerical solution of 2nd order ODE Calculus & Beyond Homework 3
Second order ODE (numerical) Differential Equations 10