Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

I see that you have been struggling with animation of rigid body motion.  I think a good deal of the trouble is due your dependence on the ugly/messy Euler angles and Euler rotation matrices.  If you move forward from Euler's time by about a century, you will discover the beautiful algebra of quaternions and their use in describing rotations in 3D.

Unfortunately these have not found their way yet into the standard curriculum.  At one point I started writing a textbook on the subject based on my lecture notes, but have not had the time/energy/motivation to finish.  Maybe some day…

In the mean time, you may find the attached Maple worksheet of some use.  It illustrates the use of quaternions in rigid body dynamics.  You will see that there are no Euler angles or rotation matrices in it.  I think the contents should be easy to grasp and modify to suit your needs although it gives no explanation of the theory behind it.  Here is a sample animation produced in that worksheet:

rigid-body-dynamics-with-quaternions.mw

Edit:
But wait!  Looking over this now I see that in my attemp to simplify things I may have gone too far and written nonsense.  Don't bother with it for now.  I will double-check and upload a corrected version later.

Edit:
OK, false alarm. What I had posted was correct after all. Here I am posting a revised version which expands some of the documentation.  The calculations remain the same.

rigid-body-dynamics-with-quaternions-v2.mw

Note that the pair of algebraic equations 2A + B = 0, 2B + A = 0 implies that A = 0 and B = 0.

Comparing these with your eq2 and eq3, we conclude that the second derivatives of y(x) and z(x) are zero, and therefore their first derivatives are constants.

But that contradicts your boundary conditions D(y)(0)=0, D(y)(1)=1, hence no solution.

 

The boundary condition D[1](u)(0,t)=1000 should be D[1](u)(0,t)=-1000.  Remember that heat flows from hot to cold, that is, against the gradient.

Method 1:

pointplot([seq([n/L, g[n]], n = -L .. L)],  color = blue, connect);

Method 2:

plot([seq([n/L, g[n]], n = -L .. L)],  color = blue);

The equation is solvable if and only if

Int(f(x)*sin(x), x = 0 .. Pi) = 0  

To see this, forget about the boundary conditions for the moment.  Find (by hand, not Maple) the general solution of the differential equation which will involve two arbitrary constants, let's say C1 and C2.   This is done through what is known as the method of variation of constants which you will find in any elementary differential equations textbook.

Then applying the boundary conditions leads to two linear algebraic equations in the two unknowns C1 and C2.  You will immediately see that the solvability of the system is equivalent to the integral condition stated above.

Advice: It is best to do this in paper and pencil first if you want to understand what is going on. Afterward you can do the same thing in Maple (it's very easy).  If you get stuck, show your work then I or someone else will help you.

 

<(f(t) $t=-3..3)>;

Vector[column]([[f(-3)], [f(-2)], [f(-1)], [f(0)], [f(1)], [f(2)], [f(3)]])

 

Your file consists of ASCII NUL characters—exactly 2180223 of them.  There's nothing else in it.  Sorry.

Details:

On a Unix system we verify this by doing:

 tr -d '\0' < Noter_til_matematik_1.mw > zz.mw
 ls -l zz.mw

which deletes all NUL characters in the file and writes what remains to the file zz.mw.  We see that zz.mw has size zero, therefore there is nothing in it.

You didn't give the value of C_A0.  Try setting it to 1 and then doing:

plots:-odeplot(dsol, [[W,C_A], [W,X(W)]], W=0..10000);

Maple's Import() command reads an STL object into a PLOT3D structure:
        P := Import("file.stl");

Suppose that your pendulum hinges about a line L in the xyz space.
We identify L by specifying two points on it, for example:
        pt1 := [0,0,0];  pt2 := [0,0,1];
Change these to whatever corresponds to the pendulum in your STL.

Preben has shown how to obtain the function TH(t) that gives
the angle of rotation of the pendulum about the hinge at time t.
Therefore the pendulum's rotated position is given by:
        rotate(P, TH(t), [pt1,pt2]);
where rotate() comes from Maple's plottools package.

Let's say you want to animate the motion over the time interval 0..T
through n frames.  We generate the frames through
        frames := seq(rotate(P, TH((k-1)/(n-1)*T), [pt1,pt2]), k=1..n):
then display the animation:
        display([frames], insequence);

Your line:

For[i=1, i<m+2, i++, w[i,1]=0; w[i,n+1]=0];

looks vaguely similar to the following C code:

for (i=1; i<m+2; i++) {
    w[i,1] = 0;
    w[i, n+1] = 0;
}

If that's what you really mean, then the Maple equivalent is:

for i from 1 to m+1 do   # note m+1, not m+2
    w[i, 1] := 0;        # note :=
    w[i, n+1] := 0;
end do;

 

Is this what you mean?

with(plots):  with(plottools):
Cone := cone([0,0,1], 1, -1, grid=[51,2], color=gold, transparency=0.4):
Cyls := display([seq(cylinder([0,0,0], 1-h, h, color=red, strips=50),
    h=0.01..1, 0.02)], insequence):
display([Cone, Cyls], style=surface, lightmodel=light1, axes=none);

 

I think what you are asking is how to solve the equation Av = 0, where A is a 4x4 matrix and v is a vector.

If so, then

LinearAlgebra:-NullSpace(A);

will give you what you want.

If you meant something else, then you will have to explain what your notation means.

 

It's easy to calculate the modulus.  It is also easy to verify that it is not always less than 1.  Why do you expect it to be less than 1?

z := coth(x+I*y);

Compute the (square of the) modulus:

z*conjugate(z):
expand(%) assuming x::real, y::real:
m := simplify(%);

(cosh(x)^2+cos(y)^2-1)/(cosh(x)^2-cos(y)^2)

Now m is the square of the modulus of z at any x and y.  Let's evaluate it at a point on the boundary of your square domain with n=0:

subs(x=(n+1/2)*Pi, y=0, n=0, m):
evalf(%);

The result is greater than 1.

Try with n=1:

subs(x=(n+1/2)*Pi, y=0, n=1, m):
evalf(%);

Again, the result is greater than 1.

 

At any given location the daylight fraction varies with the time of the year.  Days are longer in the summer and shorter in the winter.  The calculations in your worksheet do not account for that.

I have attached a worksheet which introduces a function Z(tau, lambda, t), where tau is the tilt of Earth's axis, lambda is the latitude of a location, and t is the time of the year.  In a full year t goes from 0 to 2*Pi.  The function Z is the fraction of daylight at a latitude lambda at time t of the year.

Here are the graphs of Z(Pi/6, Pi/4, t) and Z(Pi/6, lambda, Pi/2):

 

daylight2.mw

Note added on Dec 28:

I wasn't happy with the ugly calculations in the worksheet which I posted above.  I examined the details and saw that the complications were due to a non-optimal choice of basis vectors.  With a better choice of basis vectors we get a cleaner calculation.  See:

daylight3.mw

The differences between the two worksheets are cosmetic.  The resulting graphs are identical.

 

Perhaps this is what you are looking for:

http://www.mapleprimes.com/questions/208795-Periodic-Piecewise-Function

 

 

First 47 48 49 50 51 52 53 Last Page 49 of 58