## 6071 Reputation

15 years, 209 days

## A couple of fixes...

Your PDE has two independent variables, x and t. In your formulation, you are changing x but not t, but PDEtools:-dchange expects that you change both independent variables, let's say from (x,t) to (xi,tau).  You want tau = t which is okay; you just have to tell dchange that that's what you want. Therefore we change your code to:

```restart;
tr := { x = xi + mu*tau, t = tau, u(x,t) = U(xi) };
pde := diff(u(x,t),t) +p*u(x,t)*diff(u(x,t),x) + q* diff(u(x,t),x\$3) = 0;
PDEtools:-dchange(tr, pde, [xi, tau, U]);```

This yields the expected result

## Phase portrait...

```restart;
with(plots):
local x, y, L;
plottools:-polygon([[x, y-L/3], [x+L,y], [x,y+L/3], [x+L/6,y]], _rest);
end proc:
display(%);
end proc:
display(
plot(x^2/4, x=-4..4, thickness=3, color=blue),
plot([0,lambda, lambda=-2..4], thickness=3, color=blue),
plot(
[[[-4,-1],[4,-1]], [[-4,0],[4,0]], [[-4,1],[4,1]], [[-4,2],[4,2]]],
linestyle=dash, color="Green"),
pointplot([[0,-1], [0,0], [-2,1], [0,1], [2,1],
[-2*sqrt(2),2], [0,2], [2*sqrt(2),2]],
symbol=solidcircle, symbolsize=20, color=black),
[-2,-1,-0.3], [2,-1,-0.3],
[-2.7,0,-0.3], [2.7,0,-0.3],
[-3.3,1,-0.3], [3.3,1,-0.3],
[-1.0,1,0.3], [1.0,1,0.3],
[-3.8,2,-0.3], [3.8,2,-0.3],
[-1.5,2,0.3], [1.5,2,0.3]
], color=red),
axes=boxed, labels=[x,lambda], view=-2..3, scaling=constrained);
```

## Solution...

```restart;
with(plots):
eq1 := 2*x^2 + 14*y^2 = 7;
eq2 := 2*x^2 - 2*y^2 = 3;

#Intersection point in the first quadrant
pt := fsolve({eq1,eq2}, {x=0..2, y=0..1});

# Calculate the normal vector to eq1, then turn it by 90 degrees
# to make it into a tangent vector.
# Additionally, normalize the tangent vector to unit length.
n1 := < diff(lhs(eq1),x), diff(lhs(eq1),y) >:
< -%[2], %[1] >:
t1 := % / sqrt(%^+ . %);

# Repeat the calculations for eq2:
n2 := < diff(lhs(eq2),x), diff(lhs(eq2),y) >:
< -%[2], %[1] >:
t2 := % / sqrt(%^+ . %);

display(
implicitplot([eq1, eq2], x=-2..2, y=-1..1, color=[red,blue]),
arrow(eval([<x,y>, 0.8*t1], pt)[], color=red),
arrow(eval([<x,y>, 0.8*t2], pt)[], color=blue),
pointplot(eval(<x,y>, pt), symbol=solidcircle, symbolsize=24, color=white),
pointplot(eval(<x,y>, pt), symbol=circle, symbolsize=24, color=black),
scaling=constrained);
```

## Make two changes...

Change
S1 := plot3d(F, x = -Pi/2 .. Pi/2, t = 0 .. 2*Pi, color = "Green");
to
S1 := plot3d(F, x = -Pi/2 .. Pi/2, t = 0 .. 2*Pi, color = "Green", transparency=0.4);

Change
S2 := plot3d(G, x = 0 .. 100, t = 0 .. 2*Pi, color = "Cyan");
to
S2 := plot3d(G, x = 0 .. 1, t = 0 .. 2*Pi, color = "Cyan");

## Cannot be done...

I assume that the first r2 in your equation is meant to be r1.

If so, then the data is inconsistent since for all x and y, the left-hand side of the equation is always greater than the right-hand side, and therefore the equality cannot hold.

You can see that in the picture below where I plot the equation's left- and right-hand sides for a range and x and y.

```restart;
eq := (m1 + m2)*(x^2+y^2) + 2*m1/r1 + 2*m2/r2 = c;
EQ := subs(
r1 = sqrt((x-d1)^2 + y^2 + z^2),
r2 = sqrt((x-d2)^2 + y^2 + z^2),
d1 = 0.6, d2 = 0.4,
m1 = 10, m2 = 2, c = 20,
z = 0,
eq);
plot3d([rhs(EQ),lhs(EQ)], x=-1..2, y=-1..1,
numpoints=1000, view=0..300, color=[gold, cyan], style=patchcontour);```

## The op command...

Depending on the way you want to use the indices, this alternative to Carl Love's suggestion may also be useful:

op(node[2,3,5]);

which yields 2, 3, 5.

## Plotting contour lines...

It is possible to plot those curves through solving differential equations.but that's too much work for little gain.  Maple's contourplot function is much more effective for that purpose.  Here is how to produce the curves at the bottom of that website:

```restart;
with(plots):
# Some judicious selection of contour levels:
levels1 := [ i^2 \$i=0..6 ] /~ 5 +~ 1;
levels2 := i/10 \$i=1..9;
levels := [levels1[], levels2[]];
display(
contourplot(exp(y^2)*cos(x)^2, x=-Pi..Pi, y=-2..2, contours=levels),
contourplot(y*sin(x), x=-Pi..Pi, y=-2..2, contours=11,
filledregions=true, coloring=[red,yellow]),
scaling=constrained, axes=none);
```

## Perhaps this would be satisfactory...

```restart;
with(plots):
with(LinearAlgebra):
f := x -> x^2/2 - 1;
g := x -> 2 - x^2/4;

F := <f(x), x*cos(t), x*sin(t) >;

S1 := plot3d(F, x=-2..2, t=0..2*Pi, color="Orange"):

F, 1.5*Normalize(diff(F,x), 2):
eval(%, {x=2,t=Pi}):
A1 := arrow(%, color=red):

F, 1.5*Normalize(diff(F,t), 2):
eval(%, {x=2,t=Pi}):
A2 := arrow(%, color=red):

G := < g(x), x*cos(t), x*sin(t) >;

S2 := plot3d(G, x=-2..2, t=0..2*Pi, color="Cyan"):

G, 1.5*Normalize(diff(G,x), 2):
eval(%, {x=2,t=Pi}):
A3 := arrow(%, color=red):

eval(F, t=0):
convert(%, list):
T1 := tubeplot(%, x=-2..2, radius=0.05, style=surface, color=yellow):

eval(G, t=0):
convert(%, list):
T2 := tubeplot(%, x=-2..2, radius=0.05, style=surface, color=blue):

display(S1, S2, A1, A2, A3, T1, T2, scaling=constrained,
axes=framed, labels=[x,y,z]);
```

## plots:-display(LH,LG);  ...

`plots:-display(LH,LG);`

## plot3d(y^2, x=0..1, y=-1..1, filled, st...

`plot3d(y^2, x=0..1, y=-1..1, filled, style=surface);`

## Two different ways...

Here is the proper way of defining your Lagrangian:

```L := J/2*(diff(x(t),t)/R)^2
+ m__r/2*diff(x(t),t)^2
+ m__p/2*((diff(x(t),t) + l*diff(phi(t),t)*cos(phi(t)))^2 + (l*diff(phi(t),t)*sin(phi(t)))^2)
+ m__p*g*l*cos(phi(t));```

but be sure to verify that I have not introduced errors in transcribing the image that you have provided.

If you are interested in the final form of the Euler-Lagrange equations and don't care about the intermediate steps, then Maple's EulerLagrange function in the VariationalCalculus package will produce the result directly, as in Method 1 below.  Otherwise look at Method 2.

Method 1:

`tmp1 := VariationalCalculus:-EulerLagrange(L, t, [x(t),phi(t)]);`

tmp1 consists of a set of four expressions—two equations of the form something=K1 and something=K2, which are the integrals of motion (conservation of energy and momentum).  The other two are the differential equations of motion.  In the next step I delete the integrals of motion, and keep the differential equations:

`tmp2 := remove(type, tmp1, equation);`

tmp2 consists of two differential equations as I noted above.  To verify that we indeed have two items, we do:

`nops(tmp2);`

and indeed it returns "2".

Next, we name those equations de1 and de2:

```de1 := tmp2[1] = 0;
de2 := tmp2[2] = 0;```

Then you do whatever you need to do with de1 and de2

Method 2:

If you insist on doing the calculations one-step-at-a-time as you have outlined in your statement, you need to load the Physics package.  The package redefines Maple's "diff" operator which then enables you to calculate derivatives with respect to derivatives.

```with(Physics):
dL_dxdot := diff(L, diff(x(t),t));         # derivative of L with respect to xdot
dL_dphidot := diff(L, diff(phi(t),t));     # derivative of L with respect to phidot
dL_dphi := diff(L, phi(t));                # derivative of L with respect to phi
diff(dL_dxdot, t);```

## Try this...

```p1 := plot(x^2 - 1, x = -2 .. 2, color=blue):
p2 := plots:-pointplot([3/2,5/4], symbol=solidcircle, symbolsize=20, color=white):
p3 := plots:-pointplot([3/2,5/4], symbol=     circle, symbolsize=20, color=black):
plots[display](p1, p2, p3);
```

## Pedal surface...

In the attached worksheet I show how to generalize the definition of the pedal curve to pedal surface and calculate the equation of the pedal surface in the general case.  Then I illustrate the result by calculating the pedal surface of an ellipsoid with respect to its center.  Here is what it looks like.

## No need for extra packages...

No need for extra packages.  Just do this:

```restart;
k := x -> piecewise(x < -6, 1/2*x + 10, x < -3, 1, x < 6*Pi, sin(1/2*x) + 2, x < 25, x - 15);
plot3d(k(x), t=0..2*Pi, x=-20..20, coords=cylindrical);
```

## A better parametrization...

The parametrization of a surface (any surface) is not unque.  The way a surface is parametrized can have a drastic effect in how it is rendered on screen.

The standard parametrization of the oloid does not produce good computer-generated graphics as you have observed. But an oloid has a great deal of symmetries.  It is not difficult to see that it may be decomposed into four congruent patches, each of which has a quite "tame" parametrization for which Maple's default plotting grid is more than adequate.  Here I construct the oloid by putting together the four patches.

```restart;
alpha[1] := sin(t):
alpha[2] := -1/2 - cos(t):
alpha[3] := 0:
beta[1] := 0:
beta[2] := 1/2 - cos(t)/(1+cos(t)):
beta[3] := sqrt(1 + 2*cos(t))/(1+cos(t)):
x := (1-m)*alpha[1] + m*beta[1];
y := (1-m)*alpha[2] + m*beta[2];
z := (1-m)*alpha[3] + m*beta[3];
plot3d([[x,y,z], [x,y,-z], [z,-y,x], [-z,-y,x]], m=0..1, t=-Pi/2..Pi/2,
scaling=constrained, color=[red,green,blue,yellow]);
```

I have colored the patches in four different colors in order to delineate them.  Change surface style and color as it suits your needs.

 3 4 5 6 7 8 9 Last Page 5 of 44
﻿