As an example it is shown, how one convert a (simple) trigonometric equation into a polynomial problem and use Maple
to find a symbolic answer for the equation. The idea is to use the so called Joukowsky transform, which maps the circle
to the interval [-1, ... , +1].

I would have liked to simplify the result (as it is real in my case), but gave up. May be others have a good idea for that.

Using the procedure permsPosNeg as given in the blog entry Positive and negative permutations, including the improvement obtained in the blog entry Refactoring Maple code, consider the following procedure:

John Fredsted posted some interesting code dealing with permutations, and I suggested a small improvement.  Here, I want to continue the story of that improvement.  First, let us focus one particular line of code:

  posMaps,negMaps := seq(map((perm::list) ->
		(x::list) -> [seq(op(perm[i],x),i=1..nops(perm))]
	,perms),perms in [posPerms,negPerms]);

which uses a lexically scoped procedure to perform the permutations. The first thing to notice is that op(perm[i],x) is really equivalent to x[perm[i]]. Now that we have that perky op gone, we see that the resulting code expression will return unevaluated if x is unknown, unlike op which throws an error message (correctly so!). So now instead of using scoping, we can let Maple actually evaluate the inner perm[i] calls, and use unapply to recover a routine. Putting that together gives us my suggestion:

posMaps,negMaps := seq(map((perm::list) -> unapply([seq(x[perm[i]],i=1..nops(perm))],x),perms), 
perms in [posPerms,negPerms]);

In a recent blog entry, I proposed an easy way to plot the region between two curves. Later I read from an earlier blog entry that “filled=true” in implicitplot can produce amazing effects for plotting regions. Inspired by the blog entry, I’d like to introduce another easy way to plot the region between two curves.

To plot the region between y=f(x) and y=g(x) (x=a..b), we just need the following code:
with(plots):

f:=x->f(x): g:=x->g(x):

implicitplot(y=0, x=a..b, y=f(x)..g(x), filled=true, coloring=[green,green]);

The key to the success of this code is that Maple 8 allows varying range for the second variable y (i.e. y=f(x)..g(x)). However I was sorry to find that this is not allowed in Maple 11 (This will be addressed later.) .

Example 1  The region between y=x and y=x^2.

with(plots):

f:=x->x:g:=x->x^2:

a:=0: b:=1:

region:=implicitplot(y=0,x=a..b,y=f(x)..g(x),filled=true,coloring=[yellow,yellow]):

F:=plot(f(x),x=a-0.2..b+0.2,thickness=3,color=red):

G:=plot(g(x),x=a-0.2..b+0.2,thickness=3,color=blue):

display(F,G,region,scaling=constrained);

Example 2 The region between y=sin(x) and y=cos(x).

with(plots):

f:=x->sin(x):g:=x->cos(x):

a:=0: b:=6:

region:=implicitplot(y=0,x=a..b,y=f(x)..g(x),filled=true,coloring=[grey,grey]):

F:=plot(f(x),x=a-0.2..b+0.2,thickness=3,color=red):

G:=plot(g(x),x=a-0.2..b+0.2,thickness=3,color=blue):

display(F,G,region,scaling=constrained);

Now if we reverse the order of the range options from “x=a..b, y=f(x)..g(x)” to “y=f(x)..g(x), x=a..b”, some strange but interesting thing will happen (See Example 3）.

Example 3

with(plots):

f:=x->sin(x):g:=x->cos(x):

a:=-1: b:=6:

region:=implicitplot(y=0, y=f(x)..g(x),x=a..b, filled=true,coloring=[grey,grey]):

F:=plot(f(x),x=a-0.2..b+0.2,thickness=3,color=red):

G:=plot(g(x),x=a-0.2..b+0.2,thickness=3,color=blue):

display(F,G,region,scaling=constrained);
 

It can be seen that the region in Example 2 has been reflected with respect to the line y=x. But this is not bad because can use this phenomenon to plot regions between curves x=f(y) and x=g(y) (See Example 4).

Example 4  The region between x=y^2/2 and x=y^4/4-y^2/2.

with(plots):

f:=y->y^4/4-y^2/2: g:=y->y^2/2:

region:=implicitplot(y=0,x=f(y)..g(y),y=0..2,filled=true,coloring=[grey,grey]):

F:=plot([f(y),y,y=-1..2.3],thickness=3):

G:=plot([g(y),y,y=-1..2.3],thickness=3,color=blue):

display(region,F,G,scaling=constrained);

Finally some questions to be answered or discussed.

    1. Is “coloring” used in the examples an option in the package Plots? But I cannot find it in the Help (Typing ? coloring produces no results.) .
   2. Why the strange but interesting thing happens in Example 3 ?

   3. Why the above method cannot be realized in Maple 11?
   If we input the following code in Maple 11,
with(plots):

implicitplot(y=0,x=0..1,y=x..x^2,filled=true,coloring=[yellow,yellow]);with(plots);
An error warning will occur:
Error, (in plots/implicitplot) invalid input: invalid range for second variable

This means varying range for the second variable y (eg. y=x..x^2) is not allowed in Maple 11 , but which is allowed in Maple 8. If this is true, then I doubt if Maple 11 is really stronger than its earlier versions in all respects.

In a comment to a Mapleprimes thread, Jacques mentioned an old suggestion of Kahan's that numerical computations should return an estimate of conditioning alongside a result.

I mentioned in this comment an approach for numerical estimation of (all) roots of a univariate polynomial with real or complex numeric coeffficients that is based upon computing eigenvalues of a companion matrix. Here below is some rough code to inplement that idea, but which also returns condition number estimates associated with the eigenvalues.

I include an example of the badly conditioned Wilkinson's polynomial. It is possible that better results could be obtained by using a Lagrange basis representation of that polynomial, but I didn't try to figure out how that would work in an analogous way. The standard Maple utility, fsolve, has no problem with this example.

This forum question led to a discussion of a bitwise magazine review that compared Mathematica 5.2 and Maple 10. In that review the author struggled to get the following numeric integral to compute accurately and quickly in Maple.

evalf(Int(BesselJ(0, 50001*x)*x*exp(I*(355*x^2*1/2)), x = .35 .. 1));

Below, I reproduce an attempt at computing an accurate result quickly in Maple. I'm copying it here because that thread got quite long and messy.

In the blog  Plots of twisted ribbons, the author gave an interesting description of plotting twisted ribbons. In this blog , we give a similar description of twisted ribbons and give the geometrical interpretations of this definition.

Let r(phi)=[a*cos(phi), a*sin(phi), 0] (phi=0..2*Pi) be a circle in the xy-plane, and P be a point on the circle. Let QR be a line segment (with length of 2) passing through the point P and let P be the middle point of QR. Also, QR is coplanar with the z-axis.

Now let P rotate about the z-axis at the angular velocity of phi, where phi is the angle between OP and the x-axis. At the same time, the line segment QR is rotating about its middle point, P, at the angular velocity of theta (where theta is the angle between PQ and the z-axis and theta is dependent on phi, eg, theta=k*phi). In the whole process, QR will remain coplanar with the z-axis.

Apparently, the locus of the line segment QR is a twisted ribbon. When theta=phi/2, we have the Mobius strip.

2. The equation of the twisted ribbon

Now we try to find the equation of the surface. Clearly,

vector(OP)= [a*cos(phi), a*sin(phi), 0].

And with some geometrical manipulations, we have

vector(PQ)=[sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)].

So the vector equation of the twisted ribbon is

V(phi, t)=vector(OP)+t*vector(PQ)

And the parametric equation is
x=a*cos(phi)+t*sin(theta)*cos(phi),
y=a*sin(phi)+t*sin(theta)*sin(phi)

z=t*cos(theta)

(where theta=k*phi (k a constant), phi=0..2*Pi and t=-b..b (b determines the width of the ribbon.))

Or
x=(a+t*sin(k*phi))*cos(phi),
y=(a+t*sin(k*phi))*sin(phi)

z=t*cos(k*phi)

(where k a constant, phi=0..2*Pi and t= -b..b)

When we take k=1/2, 1, 3/2, 2,…, we have different twisted ribbons.

When k=1/2, we have the equation of the Mobius strip:
x=(a+t*sin(phi/2))*cos(phi),
y=(a+t*sin(phi/2))*sin(phi)

z=t*cos(ph)/2

(phi=0..2*Pi and t= -b..b) .

k=1

k=2

John Fredsted suggested using the following procedure (slightly modified) to determine whether an expression was deeply algebraic.

isDeeplyAlgebraic := proc(x)
	if not type(x,'algebraic') then false
	elif type(x,'atomic') then true
	else andmap(procname,x)
	end if
end proc:

TypeTools:-AddType('DeeplyAlgebraic1'
                   , proc(x)
                         if not x :: 'algebraic' then false;
                         elif x :: 'atomic' then true;
                         else andmap(type,x,'DeeplyAlgebraic1');
                         end if;
                     end proc);

The first is a completely flat expression, the second is deeply nested.  The following graphs plot the time required to determine whether each expression is "deeply algebraic" as n increases, with each approach. The graph on the left is the time required to evaluate expr1, the graph on right is for expr2.  The red plot corresponds to the procedure, the green plot corresponds to the type. As you can see, for both flat and nested expressions, the procedure is significantly faster than the type.

I then did some testing to see whether the type matching could be improved.  A more efficient use of the type mechanism is to use a structured type rather than a procedure.  Alas, I don't believe that it is possible to create a purely structured type, with no use of 'satisfies', that is equivalent to what we want.  Here is the best I could come up

TypeTools:-AddType('DeeplyAlgebraic3'
                   , 'And'('algebraic'
                           , 'Or'('atomic'
                                  , 'satisfies'(x->andmap(type,x,'DeeplyAlgebraic3'))
                                 )));

Adding that to each graph gives the following two plots

The yellow line (p3) corresponds to this new type.  For expr1, the flat expression, it is significantly faster than the previous type, and approaches the speed of the standalone procedure.  Alas, for expr2, the nested expression, it is even slower than the previous type.  However, the reason it is slower is that with a nested expression the 'satisfies' part of the type has to be evaluated, which generates a call to a user-level procedure.

This observation suggests that if the type could be expressed as a structured type, with no use of 'satisfies', it might be significantly more efficient.  While I can see no way to do that with the desired predicate, it is possible to construct a type specific to these two examples:

TypeTools:-AddType('nestedF', 'And(algebraic,Or(atomic,function(Or(atomic,nestedF))))');

This isn't equivalent to the original predicate because it only works with functions, not operators (+, *, etc).

Here are the results with that type

Now we are making progress.  The blue plot (p4) is the time for this restricted type.  It is significantly faster than even the standalone procedure.

Note that if there were a structured type, say allops, that returned true if all of the operands of an expression match the given type, then we would be able to construct a purely structured type that matches our original predicate, that is,

TypeTools:-AddType('DeeplyAlgebraic4'
                   , 'And(algebraic
                          , Or(atomic
                               , allops(DeeplyAlgebraic4)))');

I want to plot the normal vector of the Mobius strip represented by the parametric equations:
x=(2+u*cos(v/2))*cos(v):

What could be done with a module whose ModuleLoad routine redefined itself?

Could such a routine do some action, and then cover its tracks effectively by overwriting itself?

Would there be any way to use march() to examine the .mla archive member, in which that ModuleLoad routine is stored, without accessing the name of the module? Presumably any invocation of the actual module name would result in its being accessed from the library and hence trigger its ModuleLoad routine.

Maple's foldl and foldr procedures provide a convenient means to generate an n-ary function from a binary operator. For example, from another thread, one can use foldl to create an n-ary and procedure from Maple's binary `and` function:

There are some routines in Maple's library which, when called the first time, redefine themselves.

One plausible explanation for this is that the new versions are session dependent (external calls, say) while also more efficient to call (repeatedly).

For example, consider StringTools:-Join which seems typical of that package. First, consider it before it's been called at all.

> restart:
> showstat(StringTools:-Join);
 
StringTools:-Join := proc(...

In Maple we can easily plot a space curve if it is given by a parametric form:

x=x(t), y=y(t), z=z(t) (t=a..b).

However, Maple does not give a easy way to plot a space curve if it is given by the intersection of two surfaces:
F(x,y,z)=0, G(x,y,z)=0.

 So, I wonder if there are easy ways to solve this plotting problem?

For example, I'd like to know how to plot the following space curve In Maple :
x^2+y^2+3*z^2=1, 2*x+3*y+z=0.

Example 1  The region between y=x and y=x^2.

with(plots):

f:=x->x: g:=x->x^2:

F:=spacecurve([x,f(x),0],x=-0.2..1.2,color=blue, thickness=3):

G:=spacecurve([x,g(x),0],x=-0.2..1.2,color=red, thickness=3):

region:=plot3d([x,y,0],x=0..1,y=g(x)..f(x),color=grey, style=patchnogrid):

display(region,F,G,axes=normal,orientation=[270,0], scaling=constrained);

Example 2 The region between y=sin(x) and y=cos(x).

with(plots):

f:=x->sin(x): g:=x->cos(x):

F:=spacecurve([x,f(x),0], x=-0.5..6.5,color=blue,thickness=3):

G:=spacecurve([x,g(x),0], x=-0.5..6.5,color=red,thickness=3):

Region:=plot3d([x,y,0], x=0..6, y=g(x)..f(x), color=green,style=contour):

display(Region, F, G, axes=normal, orientation=[270,0],scaling=constrained);

Example 3  The region between x=y^2/2 and x=y^4/4-y^2/2.

with(plots):

f:=y->y^2/2: g:=y->y^4/4-y^2/2:

F:=spacecurve([f(y), y, 0], y=-0.1..2.1, color=blue, thickness=3):

G:=spacecurve([g(y), y, 0], y=-0.2..2.1, color=red, thickness=3):

Region:=plot3d([x, y, 0], y=0..2, x=g(y)..f(y), color=gray, style=patch,grid=[10,10]):

display(Region, F, G, axes=normal, orientation=[270,0], scaling=constrained);

Example 4  Bird’s eye view of Example 3 (Just change the orientation.)
with(plots):

f:=y->y^2/2: g:=y->y^4/4-y^2/2:

F:=spacecurve([f(y), y, 0], y=-0.1..2.1, color=blue, thickness=3):

G:=spacecurve([g(y), y, 0], y=-0.2..2.1, color=red, thickness=3):

Region:=plot3d([x, y, 0], y=0..2, x=g(y)..f(y), color=gray, style=patch,grid=[10,10]):

display(Region, F, G, axes=normal, orientation=[300,55], scaling=constrained);

I use Maple primarily in a Unix text window (TTY or "command-line" Maple), so I am used to seeing common subexpression labeling in the output of my computations. However, in Maple 11, GUI users don't see subexpression labeling by default.  I'll begin by talking about subexpression labeling as it appears in TTY Maple, then I'll talk about it in the GUI.

For starters lets look at an example in which a single subexpression is labeled: