Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Although you wrote F:=(x,y)->(x+y,2*x^2+y^2), I suspect that you would rather have F as a vector-valued function.  Additionally, since F is a function from R2 to R2, the proper name of the derivative is the Jacobian, not the gradient.  With these in mind, here is what we do:

restart;

sys := [x+y, 2*x+y^2];

[x+y, y^2+2*x]

F := unapply(Vector(sys), [x,y]):

J := unapply(VectorCalculus:-Jacobian(F(x,y), [x,y]), [x,y]):

Thus we have:

F(x,y);

_rtable[18446884420085923838]

J(x,y);

_rtable[18446884420085924678]

Write a proc f which receives a guess for T and produces a new guess:

f := proc(T)
  local Tnew;
  blah blah ... compute the new T
  return Tnew;
end proc;

Then apply the proc in a loop, as in:

T1 := starting value of T;
T := f(T1);
while abs(T-T1)> 1e-3 do
  T1 := T;
  T := f(T1);
end do;

Upon exit from that loop, the value of T will be approximately what you expect.  The accuracy iis controlled by the value 1e-3 in the code above.  You may increase or decrease that as needed.

PS: I saw rlopez's note after I posted this.  And he is right; the way you have worded makes no sense.  In what I have shown I have assumed that we want two successive values of T to be the same, which agrees rlopez's guess.

PPS: In general there is no guarantee that such an iteration will converge.  You may want to add a counter which is incremented in each pass of the loop, and break out of the loop if the number of iterations exceeds a prescribed maximum.

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

 

First 38 39 40 41 42 43 44 Last Page 40 of 59