From: rusin@washington.math.niu.edu (Dave Rusin) Subject: Re: system of first order odes Date: 20 Apr 1995 00:00:00 GMT Newsgroups: sci.math In article <3mnsh8$qs6@alpha.epas.utoronto.ca>, Colin Lynch wrote: >I need to solve a system of first order odes: > >xdot(t) = (a/2) y(t) - (b/2) sin(a t) z(t), >ydot(t) = -(a/2) x(t) + (b/2) cos(a t) z(t), >zdot(t) = (b/2) sin(a t) x(t) - (b/2) cos(a t) y(t) Well, no one else has come forward so I'll put in my two cents. Simple manipulations show first of all that x^2+y^2+z^2 has derivative of zero, meaning the particle is moving around on a sphere. I'll bet you knew that when you set up the equations. I also calculate that 3z+2b(sin(at)y+cos(at)x) has a derivative of zero, so that you can easily eliminate z from your set of equations. This gives you a way to reduce the problem to a 1-variable ODE. dave ============================================================================== From: rusin@washington.math.niu.edu (Dave Rusin) Subject: Re: system of first order odes Date: 20 Apr 1995 00:00:00 GMT Newsgroups: sci.math In article <3mnsh8$qs6@alpha.epas.utoronto.ca>, Colin Lynch wrote: >I need to solve a system of first order odes: I thought a little more about your system of differential equations. (I don't know why I do this; I guess I'm incorrigible.) The answer is not all that hard to write down, and I'm sure slicker methods exist for getting there, but perhaps someone reading this can learn a trick or two. The system of ODEs can be presented in matrix form as (x) ( 0 a/2 -(b/2)sin(at) ) (x) diff (y) = (-a/2 0 +(b/2)cos(at) ) * (y) (z) ((b/2)sin(at) -(b/2)cos(at) 0 ) (z) Let me remark before I begin that you can easily assume a=1 by scaling time t appropriately. I won't need to do that though; on the other hand, I will find it very convenient to switch to another -- moving -- orthonormal coordinate system. As I remarked in a previous post, the two expressions x^2+y^2+z^2 and (b cos(at))x + (b sin(at))y + (3a)z have derivatives equal to zero, meaning they are constant over the particle's journey. In particular, at any time t the particle will lie in a certain sphere (which will not change over time) and in a certain plane (which will change, but in a very specific way, as we'll see below). Thus the particle is at that moment in the intersection of the two, which is a circle. So you can imagine a circle sweeping around the sphere, like a searchlight against the dome of the sky, and the particle moving around the side of the circle as the circle moves. Let me be more specific about the plane. If C is the constant such that (b cos(at))x + (b sin(at))y + (3a)z = C for all times t, then at any fixed time t, the set of points (x,y,z) satisfying this equation lie on a plane perpendicular to the vector v1 = ( b cos(at) , b sin(at) , 3a ) Now, it's easy to write down two vectors perpendicular to this one: v2 = ( -sin(at) , cos(at) , 0 ) and v3 = ( 3a cos(at), 3a sin(at), -b ) will do. So the plane through the origin which is perpendicular to v1 is the span of these two vectors. We can get all the other planes which are perpendicular to v1 by translating this one in the direction of v1. In particular, the plane in question may be described as the set of points { C1 v1 + C2 v2 + C3 v3 : C2, C3 arbitrary} where however C1 is fixed -- different choices of C1 will give all the different planes described by the different choices of C. The geometry is a little easier if we shrink our vectors. Let r = sqrt( b^2 + 9 a^2 ) and let w1 = v1/r w2=v2 w3 = v3/r. Then the wi are of length one and still orthogonal, so that a typical linear combination C1 w1 + C2 w2 + C3 w3 lies on the sphere of radius R iff C1^2 + C2^2 + C3^2 = R^2. Also, taking dot products with v1 = r w1, we see this linear combination lies on the plane iff C1 r = C. Thus, once the sphere radius and the plane displacement are fixed, so is C1 = C/r and the C2 and C3 are arbitrary except that they must satisfy C2^2 + C3^2 = (R^2 - C^2/r^2). Let C4 be the square root of this last term; then we can parameterize the pairs C2 and C3 which are needed as the coordinates of a circle of radius C4. Thus, finally, we have a description of the intersection of the plane and the sphere: it's the set of points C1 w1 + C4 cos(u) w2 + C4 sin(u) w3 As time goes on, the vectors w1 w2 and w3 are swinging around the unit sphere. Think of a radar antenna: w1 is the vector sticking out front (where the thing-a-ma-jig of the antenna is), w2 is painted on the radar dish parallel to the ground, and w3 is painted on the radar dish perpendicular to that. Your particle is being swung round as if on a lasso attached at the end of w1 and rotating parallel to the radar dish. There is only the "minor" issue of determining _where_ on that moving circle the particle is at each time t: in other words, we have to find u in terms of t. Our description of the thing's position makes it possible to compute its velocity as C1 w1' + C4 cos(u) w2' + C4 sin(u) w3' + ( - C4 sin(u) w2 + C4 cos(u) w3 ) (du/dt) It's easy to compute v1'= ab v2, v3' = 3a^2 v2, and to note that (3a v1 - b v3 ) = r^2.(0,0,1), so that v2' = (-a/b) ( v1 - 3a [3a v1 - b v3] / r^2) Thus we get w1' = (ab/r) w2 w2' = (-ab/r w1 - 3a^2/r w3 ) w3' = (3a^2/r) w2 so that the velocity of the particle is (-ab/r)(C4 cos(u)) w1 + (C1 ab/r + C4 3a^2/r)w2 + (-C4 sin(u) ) (du/dt) w2 (-3a^2/r)(C4 cos(u)) w3 + ( C4 cos(u) ) (du/dt) w3 On the other hand, the velocity of the particle is _supposed_ to be (a certain matrix A) times (its position). If I've done it properly the matrix products are A w1 = (ab/r) w2 A w2 = (-ab/r) w1 + (3a^2+b^2)/(2r) w3 A w3 = -(3a^2+b^2)/(2r) w2 so that the differential equations ask that the velocity be (-ab/r) C4 cos(u) w1 + (ab/r) C1 - (3a^2+b^2)/(2r) C4 sin(u) w2 + (3a^2+b^2)/(2r) C4 cos(u) w3 Comparing the two formulas for the velocity will tell us how the particle is to move around in its circle: after cancellation, we obtain the condition that (C4=0 or) du/dt = r/2. This is _constant_, so that the particle is moving evenly around its circle (the person swinging the lasso is doing so smoothly). Tossing in the constant of integration here, we get u = (r/2)t + C5, so that the general solution of the system is C1 w1 + C4 cos(rt/2+ C5) w2 + C4 sin(rt/2+C5) w3 where w1 = ( b cos(at) , b sin(at) , 3a )/r w2 = ( -sin(at) , cos(at) , 0 ) and w3 = ( 3a cos(at), 3a sin(at), -b )/r and r = sqrt( b^2 + 9 a^2 ) There are, of course, three arbitrary constants C1, C4, and C5. Using the angle-addition formulas for sine and cosine, the general solution is a linear combination of these three specific solutions: p1 = v1 = ( b cos(at) , b sin(at) , 3a ) p2 = cos(rt/2) w2 + sin(rt/2) w3 p3 = sin(rt/2) w2 - cos(rt/2) w3 I leave it to the reader to express these last two fundamental solutions in the original coordinate system. dave PS -- as a rule of thumb, it's a lot easier to solve a system of equations if I have a sense of what phenomenon they're supposed to model.