
#1
Feb1313, 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] = K_{1} * x * cos(t) + ( (K_{2} * [itex]\frac{dy}{dt}[/itex])  [itex]\frac{dz}{dt}[/itex] ) [itex]\frac{d^{2}y}{dt^{2}}[/itex]= K_{1} * y * cos(t) + ( (K_{2} * [itex]\frac{dz}{dt}[/itex])  [itex]\frac{dx}{dt}[/itex] ) [itex]\frac{d^{2}z}{dt^{2}}[/itex]= ( (K_{2} * [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 



#2
Feb1313, 11:30 AM

P: 262

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. 



#3
Feb1313, 02:56 PM

P: 101

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 interrelated. How would I go about this for this case? Any help is much appreciated! 



#4
Feb1413, 04:06 AM

P: 101

Numerical Integration of 2nd Order DE 



#5
Feb1413, 04:35 AM

Sci Advisor
PF Gold
P: 1,111

(4'): [itex]\frac{dx}{dt} = u[/itex] (5'): [itex]\frac{dy}{dt} = v[/itex] (6'): [itex]\frac{dz}{dt} = w[/itex] 



#6
Feb1413, 04:54 AM

P: 101

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. 



#7
Feb1413, 07:03 AM

Sci Advisor
PF Gold
P: 1,111

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 4thorder 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. 



#8
Feb1413, 07:30 AM

P: 101

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') ? 



#9
Feb1413, 07:50 AM

Sci Advisor
PF Gold
P: 1,111

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]. 



#10
Feb1413, 08:17 AM

P: 101

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: Because we can derive expressions for u, v and w from (1'), (2') and (3') 



#11
Feb1413, 08:47 AM

Sci Advisor
PF Gold
P: 1,111

Lets take it one step back. You converted your system of 3 2ndorder differential equations into a system of 6 firstorder 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. 



#12
Feb1413, 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 