Solving Systems of Vector Differential Equations

[Arrowstar]

Hello,

I have a question regarding solving differential equations, in particular systems of differential equations, that involve vectors. For example, my particular problem is attached to this post (solution variables are alpha and Ac). While I'm no stranger to either vectors or solving ODEs, I'm not particularly sure how to go about solving a system like this when it involves cross -products, differentiation, and the like. Any tips or good places to start?

 

[agentgonzo]

It's been a while since I've done anything this high-level in maths, but if you let b = d(alpha)/dt, then the latter part of the first equation is
... + b x (b x r1).
Due to the way this works out, (b x r1) must be orthogonal to the plane containing b and r1, which means that (b x (b x r1)) must be in the plane containing b and r1, so you can simplify down to (b x (b x r1)) ≡ mb + nr1 for scalar n and m.
This gives
a1 = ac + alpha x r1 + (mb + nr1)
You can then dot with r1 to get rid of the cross term
a1.r1 = ac.r1 + r1.(alpha x r1) + m(r1.b) + n(r1.r1)
= ac.r1 + 0 + m(r1.b) + n(len(r1)^2)).

Not sure if this helps, but may be a way forward. As I said, it's been a while since I did vectors or DE.

 

[Arrowstar]

Interesting, I'll have to look into that. What does len() do? Thank you for a place to start.

 

[agentgonzo]

Interesting, I'll have to look into that. What does len() do? Thank you for a place to start.

len would be the length of the vector (or modulus).

ie, len(x) = |x|

Sorry, been coding for too long. I've lost the ability to get the correct mathematical terms.

 

[Kozak]

You can use Lagrange's identity:

Then, after rewriting them per component, you've got 3 differential equations with known quotients.

 

[Arrowstar]

Darkwalker: what is the significance of the [] and ()? I'm afraid I've never seen that before....

 

[Quick_Nick]

EDIT: Alright. I seem to have been more or less right. Hopefully this image will look better.

http://upload.wikimedia.org/math/c/a/9/ca921a2ba29d5c8ab3e2db70d4b64a27.png

 

[Arrowstar]

Okay. A, B, and C are vectors, no? So you end up with a scalar times a vector, and another scalar times a vector... okay. So then if I break those up into X Y Z components and solve, I suppose I'm all set. To clarify:

d(vector)/dt = dx/dt i + dy/dt j + dz/dt k

Where i, j, and k are the cartesian components of vector, and i, j, and k represent the unit vectors in the x, y, and z axes, respectively. Correct?

 

[Kozak]

Yeah, that's right.

 

[Arrowstar]

While we're talking about vector rules and properties, let me ask another question. When I take the cross of two vectors, I end up with another vector. Presuming my vectors are in [x y z] components, will the result of cross product end up in [x y z] format as well? Explained another way, does

[x1 y1 z1] X [x2 y2 z2] = [x3 y3 z3]

make sense? I believe it does, I just want to makes sure when I go to split components up after taking a cross product, I still end up with something like [x y z].

 

[agentgonzo]

Yes, the result of a cross product is a vector of the same dimension as the two input vectors:
http://en.wikipedia.org/wiki/Cross_product#Matrix_notation

 

[Kozak]

Yes, the result of a cross product is a vector of the same dimension as the two input vectors:
http://en.wikipedia.org/wiki/Cross_product#Matrix_notation

Yeah, specifically n=3. Cross product is defined for 3-component vectors only.

 

[Arrowstar]

You can use Lagrange's identity:

Then, after rewriting them per component, you've got 3 differential equations with known quotients.

So I've done this, I have my system of 6 ODEs. Refering to the original system I posted in the OP, I have equations that refer to Ac_x, Ac_y, Ac_z, alpha_x, alpha_y, and alpha_z (and d(alpha_x)/dt, etc etc etc). Now... who's got a good way of solving this beast? :rofl: I've run it through TI's Derive, but that comes up with just "[]" and my MATLAB install (and thus the symbolic_toolbox) seems to be having issues. I don't suppose there exists trial versions of Maple or Mathematica?

---------- Post added at 01:00 PM ---------- Previous post was at 11:31 AM ----------

After getting MATLAB working, I have attempted to solve the system of equations, which appears as:

a1x = acx - r1x*(omega_y^2 + omega_z^2) + r1y*(omega_x*omega_y - alpha_z) + r1z*(alpha_y + omega_x*omega_z)
a1y = acy + r1x*(alpha_z + omega_x*omega_y) - r1y*(omega_x^2 + omega_z^2) + r1z*(omega_y*omega_z - alpha_x)
a1z = acz + r1x*(omega_x*omega_z - alpha_y) + r1y*(alpha_x + omega_y*omega_z) - r1z*(omega_x^2 + omega_y^2) 
a2x = acx - r2x*(omega_y^2 + omega_z^2) + r2y*(omega_x*omega_y - alpha_z) + r2z*(alpha_y + omega_x*omega_z)
a2y = acy + r2x*(alpha_z + omega_x*omega_y) - r2y*(omega_x^2 + omega_z^2) + r2z*(omega_y*omega_z - alpha_x)
a2z = acz + r2x*(omega_x*omega_z - alpha_y) + r2y*(alpha_x + omega_y*omega_z) - r2z*(omega_x^2 + omega_y^2)

Solving for the Ac terms (x y z) and the omega terms (x y z) is what I'm after. The reason I solve for omega is because omega is d(alpha)/dt, which means I immediately have a solution I can plug into the numerical integrator I've written. The Ac terms are also unknowns, and not related to alpha or omega directly, so I might as well solve for it as well.

Now, when I try to solve that system the way I described, it appears as if the solver algorithm cannot find a solution (runs and runs and runs...). Admittedly, I'm keeping omega a distinct term, meaning that the symbolic software in MATLAB doesn't understand the d(alpha)/dt = omega relationship. I suppose this could be the source of the issue, but since MATLABs differential equation solver doesn't appear to accept some variables as known constants, I'm having issues using that function. Anyone have any hints?

 

[Arrowstar]

So I found a paper that supposedly addresses this question. The downside is that it appears to be slightly bogus... see the attachment. The author goes ahead and calculates the centroid of a triangle as normal (sums the coordinates and divides by as many points as there are). However, he seems to do the exact same thing for the velocity and acceleration of the centroid! That is, he does:

c_dot = (1/3)*(v1 + v2 + v3)

where v1, v2, v3 are the velocity vectors of the three points of the triangle and c_dot is the velocity vector of the centroid of the triangle. He does the same thing for the acceleration of the centroid.

Does this seem odd/incorrect to anyone else? It'd be a shame if he was wrong, the rest of the paper is incredibly clear and it'd be an excellent resource to have to throw out.

 

[agentgonzo]

The author goes ahead and calculates the centroid of a triangle as normal (sums the coordinates and divides by as many points as there are). However, he seems to do the exact same thing for the velocity and acceleration of the centroid! That is, he does:

c_dot = (1/3)*(v1 + v2 + v3)

Yes, that seems to be correct to me. Why wouldn't it be?