Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Change the last line to:

plot( [seq(seq(eval(lambda, Nb=j), j in [0.1,0.2,0.3]), F in [0.1,0.2,0.5])], delta2=0.02..0.1);

 

The answer is 42.    (cf Douglas Adams)

 

To see that, let:

f := x -> -197+(6212/15+(-3571/12+(195/2+(-179/12+13*x*(1/15))*x)*x)*x)*x;

proc (x) options operator, arrow; -197+(6212/15+(-3571/12+(195/2+(-179/12+(13/15)*x)*x)*x)*x)*x end proc

Then we have:

f(1), f(2), f(3), f(4), f(5), f(6);

3, 10, 2, 7, 7, 42

Regarding your remrk "the movement seems to be impossible in real life", yes, that is impossible.  If you inspect closely the animation that you have posted, you will note that at certain point the tetrahedra penetrate each other.

However, it can be done with 8 tetrahedra without penetration.  Here it is:

Here is the worksheet that produced it: kaleidocycle.mw

And here is a variant where each tetrahedron is painted in four colors:

and here is the worksheet: kaleidocycle8-alt.mw

——————————————————

 

For whatever it's worth, here is the self-penetrating six-tetrahedron version:

and here is the corresponding worksheet: kaleidocycle6.mw


 

A kaleidocycle of 6 tetrahedra

We construct a loop of 6 regular tetrahedra forming a kaleidocycle, and animate it.
The construction requires the tetrahedra to penetrate each other, so the result is not "realistic".
To have a kaleidocycle with non-penetrating members, we need at last 8 tetrahedra.

2018-02-17

restart;

with(plots): with(plottools):

Each of the regular tetrahedra has edge length 1.
The first tetrahedron has vertices T[i], i=1..4.  As the kaleidocycle twists, the first edge, T[1]--T[2], remains in the plane y=0 in the Cartesian coordinates.  The opposite edge, T[3]--T[4], remains in the plane y=tan(Pi/3)*x.  C[1] and C[2] are the midpoints of those edges.

We describe the overall configuration of the kaleidocycle as a function of the independent variable t which is the angle of the first edge relative to the y axis.  The opposite edge makes an angle of s relative to the z axis.  The angle s depends on t.

We express the T's and C's in terms of the geometric parameters a, b, d which also depend on t.

C[1], C[2] := <a,0,d>, <b*cos(Pi/3), b*sin(Pi/3), -d>;
T[1], T[2], T[3], T[4] :=
  C[1] + <cos(t),0,sin(t)>/2,
  C[1] - <cos(t),0,sin(t)>/2,
  C[2] + <sin(s)*cos(Pi/3),sin(s)*sin(Pi/3),cos(s)>/2,
  C[2] - <sin(s)*cos(Pi/3),sin(s)*sin(Pi/3),cos(s)>/2;

C[1], C[2] := Vector(3, {(1) = a, (2) = 0, (3) = d}), Vector(3, {(1) = (1/2)*b, (2) = (1/2)*b*sqrt(3), (3) = -d})

 

Vector[column](%id = 18446884273665409022), Vector[column](%id = 18446884273665409262), Vector[column](%id = 18446884273665409502), Vector[column](%id = 18446884273665409742)

(1)

Sanity check: The length of edge T[1] -- T[2] should be 1:

simplify((T[1]-T[2])^+ . (T[1]-T[2]));

1

(2)

Sanity check: The length of edge T[3] -- T[4] should be 1:

simplify((T[3]-T[4])^+ . (T[3]-T[4]));

1

(3)

Each of the remaining four edges should be of length 1.  This gives us four equations:

edge_lengths :=
       (T[1]-T[3])^+ . (T[1]-T[3]) = 1,
       (T[1]-T[4])^+ . (T[1]-T[4]) = 1,
       (T[2]-T[3])^+ . (T[2]-T[3]) = 1,
       (T[2]-T[4])^+ . (T[2]-T[4]) = 1;

(a+(1/2)*cos(t)-(1/2)*b-(1/4)*sin(s))^2+(-(1/2)*b*3^(1/2)-(1/4)*sin(s)*3^(1/2))^2+(2*d+(1/2)*sin(t)-(1/2)*cos(s))^2 = 1, (a+(1/2)*cos(t)-(1/2)*b+(1/4)*sin(s))^2+(-(1/2)*b*3^(1/2)+(1/4)*sin(s)*3^(1/2))^2+(2*d+(1/2)*sin(t)+(1/2)*cos(s))^2 = 1, (a-(1/2)*cos(t)-(1/2)*b-(1/4)*sin(s))^2+(-(1/2)*b*3^(1/2)-(1/4)*sin(s)*3^(1/2))^2+(2*d-(1/2)*sin(t)-(1/2)*cos(s))^2 = 1, (a-(1/2)*cos(t)-(1/2)*b+(1/4)*sin(s))^2+(-(1/2)*b*3^(1/2)+(1/4)*sin(s)*3^(1/2))^2+(2*d-(1/2)*sin(t)+(1/2)*cos(s))^2 = 1

(4)

The edges T[1] -- T[2] and T[3] -- T[4] should be orthogonal to each other.  This yields a relationship between s and t:

simplify((T[1]-T[2])^+ . (T[3]-T[4])) = 0; isolate(%, s);

(1/2)*cos(t)*sin(s)+sin(t)*cos(s) = 0

 

s = -arctan(2*tan(t))

(5)

We substitute (5) in the equations (4) and thus obtain a a system of four equations in the three unknowns a, b, d, with t as a parameter.  We have one too many equations but there is redundancy in them, so the system is still solvable:

eval({edge_lengths}, (5)):
sol := solve(%, {a,b,d}):
allvalues(%);

{a = (2/3)*(-3/(24*cos(t)^2-32))^(1/2)*(3*cos(t)^2-4), b = -(4/3)*(-3/(24*cos(t)^2-32))^(1/2), d = (-3/(24*cos(t)^2-32))^(1/2)*cos(t)*sin(t)}, {a = -(2/3)*(-3/(24*cos(t)^2-32))^(1/2)*(3*cos(t)^2-4), b = (4/3)*(-3/(24*cos(t)^2-32))^(1/2), d = -(-3/(24*cos(t)^2-32))^(1/2)*cos(t)*sin(t)}

(6)

We pick the solution with positive a:

params := (6)[1][], (5);

a = (2/3)*(-3/(24*cos(t)^2-32))^(1/2)*(3*cos(t)^2-4), b = -(4/3)*(-3/(24*cos(t)^2-32))^(1/2), d = (-3/(24*cos(t)^2-32))^(1/2)*cos(t)*sin(t), s = -arctan(2*tan(t))

(7)

Here are the vertices of the first tetrahedron as a function of t:

subs(params, convert~([T[1],T[2],T[3],T[4]], list)):
vertices := unapply(%, t);

proc (t) options operator, arrow; [[(2/3)*(-3/(24*cos(t)^2-32))^(1/2)*(3*cos(t)^2-4)+(1/2)*cos(t), 0, (-3/(24*cos(t)^2-32))^(1/2)*cos(t)*sin(t)+(1/2)*sin(t)], [(2/3)*(-3/(24*cos(t)^2-32))^(1/2)*(3*cos(t)^2-4)-(1/2)*cos(t), 0, (-3/(24*cos(t)^2-32))^(1/2)*cos(t)*sin(t)-(1/2)*sin(t)], [-(2/3)*(-3/(24*cos(t)^2-32))^(1/2)-(1/2)*tan(t)/(1+4*tan(t)^2)^(1/2), -(2/3)*(-3/(24*cos(t)^2-32))^(1/2)*3^(1/2)-(1/2)*tan(t)*3^(1/2)/(1+4*tan(t)^2)^(1/2), -(-3/(24*cos(t)^2-32))^(1/2)*cos(t)*sin(t)+(1/2)/(1+4*tan(t)^2)^(1/2)], [-(2/3)*(-3/(24*cos(t)^2-32))^(1/2)+(1/2)*tan(t)/(1+4*tan(t)^2)^(1/2), -(2/3)*(-3/(24*cos(t)^2-32))^(1/2)*3^(1/2)+(1/2)*tan(t)*3^(1/2)/(1+4*tan(t)^2)^(1/2), -(-3/(24*cos(t)^2-32))^(1/2)*cos(t)*sin(t)-(1/2)/(1+4*tan(t)^2)^(1/2)]] end proc

(8)

Proc to plot the kaleidocycle at a give t:

kaleidocycle6 := proc(t)
  local T;
  T[1] := tetrahedron(vertices(t));
  T[2] := reflect(T[1], [[0,0,0], [cos(Pi/3),sin(Pi/3),0], [cos(Pi/3),sin(Pi/3),1]]);
  T[3] := rotate(T[1], 0, 0, 2*Pi/3);
  T[4] := rotate(T[2], 0, 0, 2*Pi/3);
  T[5] := rotate(T[1], 0, 0, 4*Pi/3);
  T[6] := rotate(T[2], 0, 0, 4*Pi/3);
  display(T[i] $i=1..6,
    color=[red,green,blue,yellow,magenta,cyan]);
end proc:

Animate the kaleidocycle.  Note the self-penetration!

animate(kaleidocycle6, [t], t=0..Pi,
   scaling=constrained, axes=none, paraminfo=false, frames=50);

 

 

 

 


 


 


 

 

 

 

Posting several screenfuls of Maple code onto a web page is not particularly helpful and can discourage people from reading. Why not just upload your Maple worksheet instead?  In the window where you write your message for posting, note the big fat green arrow.  Click on it to upload.

I waded through what you have posted and have only a vague idea of what you are attempting to  do.  It appears that you wish to plot an uppercase letter "A" and then rotate it by some angle.  If so, then don't build the angle in the definition of A.  Just make a plain "A" in the normal upright position and name it something, let's say A.  To rotate A about a point [x,y] by angle alpha, do
plottools:-rotate(A, alpha, [x,y]);
That's all.

 

Let t be the angle that the vector makes relative to the positive x axis.  Let L be the vector's length.  Then the vector is given by

v := L*<cos(t), sin(t)>;

 

In answering your previous questions, Kitonum has shown you the proper syntax for defining position vectors.  Let's call them r1 and r2 in this case.  You need to do:
r1 := t -> 2*t^2*_i+16*_j+(10*t-12)*_k;
r2 := t -> (20-6*t)*_i+4*t^2*_j+2*t^2*_k;

Then you may verify that r1(2) and r2(2) are equal, that is, the airplanes collide at t=2.

The velocities are the derivatives of the position vectors.  They are computed via:

v1 := D(r1);
v2 := D(r2);

The velocities at the time of collision are given by v1(2) and v2(2).  You
should be able to calculate the angle between those vectors.

Additionally:

  1.  Heed rlopez's advice on the proper use of "evaluate" versus "solve".
  2.  You will receive better help if instead of posting code fragments, you upload a complete worksheet.  In the window where you edit your message before sending, observe the big fat green arrow.  Click on it to upload worksheet.

 

Consider the functions "`w__i`(x,y), i=1,2," defined for all x, y, and let
    "`D__i`={(x,y) : `w__i`(x,y)>0 }".

We say that the function w__i characterizes the domain D__i.  One can then verify that

• 

the function w__1+w__2+sqrt(w__1^2+w__2^2) characterizes the domain `union`(D__1, D__2)

• 

the function w__1+w__2-sqrt(w__1^2+w__2^2)characterizes the domain `intersect`(D__1, `D__2.`)

This idea was introduced/popularized by Rvachev; see, e.g.,

http://appliedmechanicsreviews.asmedigitalcollection.asme.org/article.aspx?articleid=1395415

 

This enables us to build a characterizing function for a complicated domain by splitting the domain
into simpler domains. Here I illustrate the method by creating the 8-shaped domain seen at the end of
this worksheet.

 

restart;

A function which is positive on .25 < abs(r) and abs(r) < 1 and nonpositive otherwise:

R := -(r^2-1^2)*(r^2-0.25^2);

-(r^2-1)*(r^2-0.625e-1)

plot(R, r=-2..2, view=-0.5..0.5);

Function of x, yobtained by rotating the graph of R about the vertical axis:

w1 := subs(r=x^2+y^2, R);

-((x^2+y^2)^2-1)*((x^2+y^2)^2-0.625e-1)

w1 characterizes one of the two loops of the eight-shaped figure:

plot3d(w1, x=-2..2, y=-2..2, view=0..0.2);

Translate w1 in the y direction to obtain the characterizing function of the eight-shape's other loop:

w2 := subs(y=y-1.8, w1);

-((x^2+(y-1.8)^2)^2-1)*((x^2+(y-1.8)^2)^2-0.625e-1)

The function w1 + w2 + sqrt(w1^2 + w2^2) is positive where one or both of w1 and w2 are positive,
and nonpositive otherwise.

By taking the max(0, ...) we replace the negative parts by zero:

f := max(0, w1 + w2 + sqrt(w1^2+w2^2));

max(0, -((x^2+y^2)^2-1)*((x^2+y^2)^2-0.625e-1)-((x^2+(y-1.8)^2)^2-1)*((x^2+(y-1.8)^2)^2-0.625e-1)+(((x^2+y^2)^2-1)^2*((x^2+y^2)^2-0.625e-1)^2+((x^2+(y-1.8)^2)^2-1)^2*((x^2+(y-1.8)^2)^2-0.625e-1)^2)^(1/2))

Here is what f  looks like:

plot3d(f, x=-2..2, y=-2..4);

The function g is 1 where f is positive, and zero otherwise:

g := piecewise(f > 0, 1, 0);

g := piecewise(0 < max(0, -((x^2+y^2)^2-1)*((x^2+y^2)^2-0.625e-1)-((x^2+(y-1.8)^2)^2-1)*((x^2+(y-1.8)^2)^2-0.625e-1)+sqrt(((x^2+y^2)^2-1)^2*((x^2+y^2)^2-0.625e-1)^2+((x^2+(y-1.8)^2)^2-1)^2*((x^2+(y-1.8)^2)^2-0.625e-1)^2)), 1, 0)

Here is what g looks like:

plot3d(g, x=-1.4..1.4, y=-1.4..3.1, grid=[200,200],
  style=surface, color=gold, scaling=constrained,
  orientation=[-135,-60,180], axes=none);

Download mw.mw

I have no idea what your code is meant to plot, but if you wish to shift the origin of an existing plot to x=1, y=2, here is how:
newplot := plottools:-translate(oldplot,1,2);

 

I don't know how to do this in Maple, but it is very easy to solve by hand.


Integrating the differential equation diff(u(x, z, t), z)+C = 0 we obtain u(x, z, t) = -C*z+f(x, t),

where f(x, t) is an arbitrary function.  In view of this, the condition u(x, h(x, t), t) = 0 takes the form

-C*h(x, t)+f(x, t) = 0, and therefore f(x, t) = C*h(x, t). Plugging this back into the previous

expression for u we arrive at the solution "u(x,z,t)=-C*z+C*h(x,t)."

In line 13 from the bottom, delete the extra "end if".

My interpretation of the problem's statement is different from acer's and kitonum's.  I view the first equation as defining γ as the root xx of the equation on the right-hand side. Similarly, the second equation defines ψ as the root xx of the equation on the right-hand side.  With that in mind, here is how to proceed.

restart;

local gamma;

eq1 := gamma*sinh(gamma)+xi*psi*cosh(gamma) = 0;

gamma*sinh(gamma)+xi*psi*cosh(gamma) = 0

eq2 := -gamma^2+beta+sqrt(psi/w)*coth(sqrt(psi/w))-psi = 0;

-gamma^2+beta+(psi/w)^(1/2)*coth((psi/w)^(1/2))-psi = 0

In eq2 we see that we need psi/w >= 0if we are to have real solutions.

From eq1 we see that xi*psi = -gamma*tanh(gamma),  that is,"psi/(w)= -1/(xi*w)*gamma*tanh(gamma)." 

But gamma*tanh(gamma) >= 0, therefore psi/w >= 0 if and only if  xi and w have opposite signs.

params := xi=-1, w=1, beta=2;

xi = -1, w = 1, beta = 2

fsolve(eval({eq1,eq2}, {params}));

{gamma = 1.449541816, psi = 1.298212894}

Download mw.mw

You have a set of equation SysEq which are linear in the variables dq.  The coefficient matrix is obtained through
J := LinearAlgebra:-GenerateMatrix(SysEq, dq)[1];

In your case J is an 8x11 matrix.  To display it, do:
evalm(J);

 

@assma Here is the code to do that calculation.  I assume that your 6.28 is actually meant to be 2π.  Change it to anything else as needed.

The 2D version

err2 := abs(wr(x,y) - w3(x,y));

range2 := x=0..2*Pi, y=0..2*Pi;

First, be sure that this one works correctly:

maximize(err2, range2);

Only then measure the usage statistics:

CodeTools:-Usage(maximize(err2, range2));

 

The 3D version

err3 := abs(wr(x,y,z) - w3(x,y,z));

range3 := x=0..2*Pi, y=0..2*Pi, z=0..2*Pi;

maximize(err3, range3);

CodeTools:-Usage(maximize(err3, range3));

 

Download mw.mw

 

restart;

with(Statistics):

pts := [[0, 0], [1, 13], [2, 21], [6, 45], [12, 54], [15, 77]];

[[0, 0], [1, 13], [2, 21], [6, 45], [12, 54], [15, 77]]

The equation of the straight line that goes through the point (2,21) is  y = 21+a*(x-2),

that is, y = a*x-2*a+21. The LinearFit function does not accept a model involving an additive
constant
such as the 21 above, but that's alright.  We subtract 21 from the data, fit with the model

y = a*x-2*a, and then add 21 to the result.

pts1 := (x->x[1])~(pts);

[0, 1, 2, 6, 12, 15]

pts2 := (x->x[2]-21)~(pts);  # note the subtracted 21

[-21, -8, 0, 24, 33, 56]

L := LinearFit(a*x-2*a, pts1, pts2, x) + 21;  # note the added 21

HFloat(4.151724137931035)*x+HFloat(12.69655172413793)

Let's verify that L passes through (2,21):

eval(L, x=2);

HFloat(21.0)

plot([pts, L], x=0..15);

Download linearfit.mw

An object of the form [u1, u2, ..., un] is called a list in Maple. If that list is named X, then X[i] is the ith entry of that list. Each of the ui can itself be a list. If every ui is a list, then we say X is a list of lists. In that case the notation X[i][j] indicates the jth element of the ith list in X.

The argument J of your proc is expected to be a list of list.

  • The first line of the proc creates a list L from the list J as follows. Let's say the ith entry of L is [a,b,c]. Then the ith entry of L is the two-element list
        [a, (b+c*I)/a + 50*(1+I)/abs(J[1][2])]
    where I is Maple's notation for sqrt(-1), and abs() is the absolute value.
    Note: You would want J[1][2] to be nonzero, otherwise you will be dividing by zero.  In the example in the last line of your post is J[1][2] zero.  That's bad.
  • The second line of the proc defines a vector R whose ith entry equals the ith entry in L.
  • The next line says that:
        T := abs(L[1][3]);
        r := L[1][4];
    But L[1][4]  makes no sense since it attempts to pick the 4th entry of the list L[1] which is a list of two entries; there is no 4th entry.  Check for typos.

The rest of the code should be self-evident.

 

First 33 34 35 36 37 38 39 Last Page 35 of 53