## 4135 Reputation

13 years, 281 days

## Solution...

 > restart;
 > with(plots):
 > with(plottools):
 > el := x^2/5^2 + y^2/3^2 + z^2/2^2 = 1; > sp := (x-a)^2 + (y-1/2)^2 + (z-1/2)^2 = 1/2^2; Normals to the ellipsoid and sphere at the point of tangency should be collinear.  This gives us three equations.

 > :  # normal to the ellipoid :  # normal to the sphere %% =~ lambda*%: eqs := convert(%, list)[]; The point of tangency lies on both the ellipsoid and the sphere.  That gives us two more equations.

 > sys := eqs, el, sp; Solve five equations in five unknowns:

 > sol := fsolve({sys}, {x,y,z,lambda,a}, a=0..5); Let's plot the result.  Here is a plot of the ellipsoid.

 > EL := display(scale(sphere(1), 5, 3, 2), style=wireframe):

And here is a plot of the sphere.

 > SP := display(sphere(eval([a,1/2,1/2], sol), 1/2), style=surface, color=red):

And this is the two together.

 > display([EL, SP], scaling=constrained, style=surface); ## Look for equilibria first...

It this, and similar problems, look for the equilibria first.  That tells you where the interesting things happen.  In your case, the equilibria are at y=0 and y=1/9=0.111.  Therefore it makes sense to limit the y direction to the range 0 to 0.2.  Negative y values are not relevant to this problem since populations cannot be negative.  Negative time values are also of no interest since we are interested in predicting the future based on what we have now.  These considerations limit the plotting region to t > 0 and 0 < y < 0.2.  The upper limit for t is determined by trial-and-error through the expectation that the solution should converge to the stable equilibrium point.

I have modified your worksheet to account for those comments.  Here is what we get: Worksheet: mw.mw

## Increasing the Java heap...

In Linux, one normally starts Maple through the command

maple -x &

In cases where the java memory is apt to be exhausted (often due to an extensive animation), one can request an increased  Java heap by changing that to

maple -x -j 4096 &

I think that option is undocumented.   The default value of the Java heap is 512.

That's how it works in Linux.  There may be similar hooks in other operating systems.

If you want insulated ends, you need to supply Neumann, not Dirichlet, conditions.  I have made that change in the worksheet below.

 > restart;
 > pde := diff(u(x,t),t) = a^2*diff(u(x,t),x,x); > ibc := D(u)(0,t)=0,  D(u)(L,t)=0, u(x,0)=x*(L-x); > %pdsolve(pde, {ibc}, numeric): eval(%, {a=1, L=1}):  # set the parameter values pdsol := value(%); > pdsol:-animate(t=0..0.3, frames=40); Accuracy issues

Looking a bit more closely to the solution presented above, I see some disconcerting issues.  We know that insulated boundary conditions imply that the heat equation preserves the average of the solution.  Since the initial condition is x*(L-x), the average is L^2/6.  Therefore the solution should approach the constant value L^2/6.  In our case L=1, therefore that constant should be 1/6 which is approximately 0.1667.  In the animation shown above the solution converges to something like 0.12 which is very far from the target.

We may increase pdsolve's accuracy by reducing the spacestep parameter whose default value is L/20.  To get a limit which is reasonably close to 0.1667, we need to take spacestep=L/1000 (see the demo below) which is a very inefficient way of solving a PDE.  There may be options to pdsolve to obtain the solution more efficiently but I can't tell.  Perhaps someone else who is more knowledgeable about pdsolve's options can chime in.

Something else that I must point out is that the wavy oscillations in the solution curves produced by the default spacestep is a numerical artifact.  Note that those oscillations  are not present in the more accurate solutions.

 > restart;
 > pde := diff(u(x,t),t) = a^2*diff(u(x,t),x,x); > ibc := D(u)(0,t)=0,  D(u)(L,t)=0, u(x,0)=x*(L-x); >

spacestep=1/20 (the default)

 > %pdsolve(pde, {ibc}, numeric, spacestep=1/20): eval(%, {a=1, L=1}):  # define parameter values pdsol := value(%):
 > pdsol:-animate([[1/6, color=blue], [u(x,t), color=red]], t=0..0.3, frames=40); >

spacestep=1/100

 > %pdsolve(pde, {ibc}, numeric, spacestep=1/100): eval(%, {a=1, L=1}):  # define parameter values pdsol := value(%):
 > pdsol:-animate([[1/6, color=blue], [u(x,t), color=red]], t=0..0.3, frames=40); >

spacestep=1/1000

 > %pdsolve(pde, {ibc}, numeric, spacestep=1/1000): eval(%, {a=1, L=1}):  # define parameter values pdsol := value(%):
 > pdsol:-animate([[1/6, color=blue], [u(x,t), color=red]], t=0..0.3, frames=40); ## Solution...

 > restart;
 > de := diff(y(x),x,x) + a*y(x) = 0; > bc_base := y(0)=0, y(L)=0; bc_extra := D(y)(0)=1;    # extra boundary condition to avoid the trivial solution bc := bc_base, bc_extra;   > dsol := dsolve({de,bc}, {y,a}); Originally _Z1, renamed _Z1~:
is assumed to be: integer

Originally _B1, renamed _B1~:

is assumed to be: OrProp(0,1)

A particular solution:

 > eval(dsol, {_Z1=5, _B1=0}): simplify(%)  assuming L>0;  # solution eval([de,bc], %);           # verify solution  Another particular solution:

 > eval(dsol, {_Z1=5, _B1=1}): simplify(%)  assuming L>0;  # solution eval([de,bc], %);           # verify solution  Solving with Neumann boundary condition on the left end

 > bc_base := D(y)(0)=0, y(L)=0; bc_extra := y(0)=1;    # extra boundary condition to avoid the trivial solution bc := bc_base, bc_extra;   > dsol := dsolve({de,bc}, {y,a}); > dsol := simplify(dsol) assuming L > 0; A particular solution:

 > eval(dsol, {_Z2=5}): simplify(%)  assuming L>0;  # solution eval([de,bc], %);           # verify solution  Another particular solution:

 > eval(dsol, {_Z2=-5}): simplify(%)  assuming L>0;  # solution eval([de,bc], %);           # verify solution  The solutions shown above are correct but they can be simplified (by hand, not Maple!)  See the worksheet below for the simplification.

## Converting to autonomous...

A field plot for a nonautonomous system doesn't make much sense.  It is like asking for a picture of a stormy sea.  You may take a snapshot at one moment, but that picture will change at the next moment—there is no one picture.

It is possible, however, to convert a nonautonomous ODE to an autonomous one at the expense of introducing an extra dimension.  Specifically, one introduces a new unknown, let's say z(t), through the equation dz/dt=1, z(0)=0, which is equivalent to z(t)=t.  Then we replace t by z(t) on the right-hand side of the equation,.

Your equation is of the second order.  We convert it to a first order system, and then append the equation dz/dt=1 and thus obtain an autonomous system of 3 equations and then produce the corresponding 3D field plot. See the details in the worksheet.

PS: Your title says "solve" but the code fragment that you have shown attempts to do a field plot which is something quite different.  In what I have written above, I have addressed the field plot issue.  If yiou are interested in a solution, then you will need to supply an initial condition.

## Hints...

I will show the solution of your homework problem's first part which can be done in Maple.  The rest, which calls for a straightforward mathematical argument, will have to be done by hand.  I have provided an outline below.

```restart;
eq := diff(z(t), t) = (a + b*I)*z(t) + G(abs(z(t))^2)*z(t);
Eq := eval(eq, z(t)=r(t)*exp(theta(t)*I));
de1 := evalc(Re(Eq));
de2 := evalc(Im(Eq));
```

Then:

1. Solve {de1, de2} as an algebraic (not differential) system of two linear equations in the two unknowns {dr/dt, dtheta/dt}.  (This step can be done in Maple.)
2. Apply the assumptions on a, b, and G to conclude that the origin is the only the equilibrium.
3. Show that the origin repels orbits.
4. Show that the origin attracts far-away orbits.
5. Apply Poincare-Bendixon to conclude that there exists a limit cycle.

## T*...

Here is how to construct T*.  I don't quite understand the definitions of the Ybar matrices. Perhaps you can explain with an example.

 > restart;
 > with(LinearAlgebra):
 > N := 3; > A := IdentityMatrix(N+1); > B := Matrix([t^k\$k=0..N]); > interface(rtablesize=(N+1)^2):
 > T_star := KroneckerProduct(A,B); PS: The following interprets the notation in the defintion of Ybar in a way that leads to the solution that yiou are looking for.

 > restart;
 > N := 3:
 > Y1 := ; > Y2 := ; > Y1_bar := <(y[2,k]*Y1)\$k=0..N>; > Y2_bar := <(y[1,k]*Y2)\$k=0..N>; ## Cost function...

Suppose the const function is f(x) = x^2.  Then you would plot it with:
plot(x^2, x=0..4);
If you have some other cost function in mind, then you should tell us what it is.

## Banded storage...

If the zero entries of B are not going to change in your calculations, then you should use banded storage for B.  That way, only the nonzero entries of B are stored, thus saving a ton of memory if your N is large.

 > restart;
 > N:=5; > B := LinearAlgebra:-BandMatrix([,[\$1..N]],0,N+1); You can change the entries on the main diagonal or on the superdiagonal if needed:

 > B[1,2] := 100; But you cannot change the entries elsewhere because there is no memory set aside
for them:

 > B[2,1] := 100;

Error, attempt to assign a value outside Matrix bands

 >

## Solution...

Here is one way of doing it.  Make the field "undefined" outside of the desired region.

 > restart;
 > f := map2(piecewise,  x^2+y^2<=1, [y, -sin(x)-(1/10)*y], undefined); > plots:-fieldplot(f, x = -1 .. 1, y = 0 .. 1, scaling=constrained); ## No differentiation necessary...

You have: .

You want to get .

No differentiation is necessary. You may produce v through Here is how we verify this in Maple.

 > restart;
 > u = exp(-lambda*z*sqrt(x^2 + y^2)); v = z*u/(1 - 2*ln(u)/z); simplify(%, {%%}) assuming positive;   ## Another way...

[f(k) \$k=0..5];

## Arbitrary length strings...

Some years ago I wrote this worksheet to make strings of arbitrary length from STL letters.  See if it is useful to you.

3D text

2015-05-31

2015-06-17 revised (see below)

2017-02-05 added a "thickening effect" to produce an "Egyptian" look.

The idea for this worksheet comes from a post by Acer in MaplePrimes:

http://www.mapleprimes.com/questions/204335-Can-Rotate-3d-Text-Like-This-Be-Done-In-Maple

The goal is to produce text with 3D letters.  It requires the letter shapes to be available in the STL format.

I downloaded a set of such letter shapes in a file named All_Alphabet_Letters_A-Z.zip  from

http://www.thingiverse.com/thing:15198

and saved it in the /usr/local/share/char-shapes-stl directory.

This worksheet defines two procs named string_to_plot3d() and string_to_plot3d_cylinder() which take a string as the first argument and a number called "gap" as the second argument.  These procs extract the 3D letters corresponding to the characters in the string, arrange them in the 3D space with "gap" distance between the characters, and return the corresponding PLOT3D structures.  The first proc arranges the characters in a linear fashion.  The second proc wraps the characters around a cylinder.

 > restart;
 > with(plottools): with(plots): with(StringTools):
 > letters_dir := "/usr/local/share/char-shapes-stl/":

This table translates single character names to their corresponding STL file names.

 > STL_shapes_table := table([   seq(c = sprintf("Letter_%s.stl", c), c in "ABCDEFGHIJKLMNOPQRSTUVWXYZ"),   "a" = "Letter_a_lower_case.stl",   "b" = "Letter_b_lower_case.stl",   "c" = "Letter_c_lower_case.stl",   "1" = "one.stl",   "2" = "two.stl",   "3" = "three.stl",   "4" = "four.stl",   "5" = "five.stl",   "6" = "six.stl",   "7" = "seven.stl",   "8" = "eight.stl",   "9" = "nine.stl",   "0" = "zero.stl",   "@" = "Letter_at.stl",   "&" = "and_sign.stl",   #"(" = "bracket_open_close.stl",  # this gives a composite "()" shape, not   #")" = "bracket_open_close.stl",  # very useful   "\$" = "dollar_sign.stl",   "€" = "euro_sign.stl",   "!" = "exclamation_mark.stl",   "-" = "minus_sign.stl",   "#" = "number_sign.stl",   "+" = "plus_sign.stl",   "£" = "pound_sign.stl",   "?" = "question_mark.stl",   "§" = "section_sign.stl",   "¥" = "yen_sign.stl", NULL]):

The Import() command, introduced in Maple2015, reads an STL shape and returns a PLOT3D structure:

 > Import(cat(letters_dir, STL_shapes_table["@"]), title=""); This is a front-end to  Import().  It throws an error if the requested character is unavailable.

 > getletter := proc(c::string)   if not assigned(STL_shapes_table[c]) then     error "shape %1 not available", c   else     Import(cat(letters_dir, STL_shapes_table[c]), title="");   end if; end proc:
 > getletter("A"); > string_to_plot3d := proc(str, gap)   local n, i, shape, a, b, p, x;   n := length(str);   p := NULL;   x := 0;   for i from 1 to n do     if str[i] = " " then       x := x + gap;       next;     end if;     shape := getletter(str[i]);     a := op(1, getdata(shape, 'rangesonly'));     b := op(2, getdata(shape, 'rangesonly'));     translate(shape, x-a, 0, 0);     p := p, select(has, [op(%)], POLYGONS)[];     x := x + b - a + gap;   end do;   return PLOT3D(p); end proc:

Test 1:

 > string_to_plot3d("MAPS", 5): display(%, scaling=constrained, color="Chocolate",     style=patchnogrid, axes=none, orientation=[-80,30,0]); Test 2:

 > string_to_plot3d("THIS IS A TEST", 10): display(%, scaling=constrained, color="Chocolate",    style=patchnogrid, axes=none, orientation=[-110,-30,0]); Here we wrap the string around a cylinder.  Additionally, we apply a transform
that thickens the letters toward their bases, giving them a sort of "Egyptian" look.

 > string_to_plot3d_cylinder := proc(str, gap)   local n, i, shape, a, b, p, x, L, R;   n := length(str);   p := NULL;   x := 0;   for i from 1 to n do     if str[i] = " " then       x := x + gap;       next;     end if;     shape := getletter(str[i]);     a := op(1, getdata(shape, 'rangesonly'));     b := op(2, getdata(shape, 'rangesonly'));     translate(shape, x-a, 0, 0);     p := p, select(has, [op(%)], POLYGONS)[];     x := x + b - a + gap;   end do;   p := PLOT3D(p);   # the "2.5" below is the thickening factor.  Letters are 2.5 times thicker at the base than at the top.   p := transform((x,y,z)->[x,y,z*(1-(2.5-1)*y/40)])(p);   L := x;   R := L/(2*Pi);   return transform((x,y,z) -> [(R+z)*cos(x/R), (R+z)*sin(x/R), y])(p); end proc:

Test 3:

 > string_to_plot3d_cylinder(" THIS IS A TEST ", 10): display(%, color="Chocolate", scaling=constrained, axes=none, style=patchnogrid); >

## Animation...

 > restart;
 > with(plots,animate):

Let the axis be vertical and point up.  A rope of length is attached to the axis

at the point .  The other end at at is loose.

When the rope is set into motion, let be the horizontal displacement
at time of the point whose vertical coordinate is Through applying Newton's law of motion we find that solves the
partial differential equation where (not necessarily a constant) is the rope's linear density, and is

the acceleration due to gravity.  This is derived under the assumption that

the displacement is small.

The details of the derivation and a solution through separation of

variables may be found at:

Here I will go ahead and use that solution in order to produce
an animation in Maple.

The solution is actually expressed as an infinite sum but here I will take

only the first few terms for the sake of simplicity, like this:

 > y := Sum(a[n]*BesselJ(0, 2*omega(n)*sqrt(x/g))*cos(omega(n)*t), n=1..4); where

 > omega := n -> 1/2*sqrt(g/L)*BesselJZeros(0,n); and are arbitrary constant coefficients.  In practice the s are determined

in terms of the initial conditions but I will not go into that.  Instead, I will just

pick arbitrary values for them.  To get realistic solutions, we need the to be

reasonably small so that is small as it is stipulated in the derivation of the

equation of motion.

Here are the parameter values.  Change them if you want.

 > params := g=1, L=1; Experiment: Let's take and the rest of the s equal to zero.  Then

we get the following animation:

 > eval(value(y), {a=0.2, a=0.1, a=0, a=0}): %animate(plot, [[%, x, x=0..L], thickness=5, color="Indigo"],          t=0..36.4, frames=100, scaling=constrained): subs(params, %): value(%); 1 2 3 4 5 6 7 Last Page 1 of 35