Carl Love

25796 Reputation

10 years, 350 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

A solution is to ?overload the procedure `diff/Int`.

I coded the overload of `diff/Int` solution. It was actually much easier than I was anticipating, just 11 lines of code. This solution illustrates a very generally applicable technique, one that can be used to modify the behavior of almost any Maple library or even built-in procedure. This technique is not often discussed on MaplePrimes, so I suggest that most readers save this Answer in their Favourites.

restart;

`diff/Int/orig`:= eval(`diff/Int`):
unprotect(`diff/Int`);
proc(f::algebraic, x::{name, name=range}, M::(identical(method)=name))
subsindets(
diff(Int(f,x), args[-1]),
specfunc(anything, 'Int'),
subs([_DUMMY= args[3..-2]], J-> 'Int'(op(J), _DUMMY))
)
end proc,

`diff/Int/orig`
]):

Test

 > S:= Int(sin(a*x)/x, x= 0..infinity, method= _DEFAULT);

 > diff(S,a);

 > diff(Int(sin(a*x)/x, x), a);

I have only overloaded the case where the first argument to Int is type algebraic, the second argument is a name or name= range, and the third argument is method= name. All other cases should behave as they did originally. If you need the overload to work for more general cases of Int, please let me know.

Sets not ordered. You have no theta(t)....

You wrote:

Notice that the order of the output from the EulerLagrange command when three functions are input does not correspond to the order of the functions listed, and when they are requested separately, the last function is not included:

The two phenomena are unrelated. As to the order: The EulerLagrange command always gives its result as a set, in curly braces. There is no way to control the order of elements in a Maple set: An internal lexicographic ordering is used so that set-membership checks can be done with binary search and set order is preserved across sessions. (Older Maple used an address-based ordering and set order varied by session.)

You wrote:

T3:=(1/2)*m*(diff(r(t),t)^2+r(t)^2*diff(theta(t),t)^2+diff(z(t),t)^2):
U3:=m*g*z(t):
L3:=T3-U3;

As for the second issue: Your equations above use diff(theta(t), t) but they never use theta(t) undifferentiated. This confuses EulerLagrange, and perhaps it should be rejected as invalid input. If you change diff(theta(t),t) to theta_prime(t) and do not include this in the solution variables, then the union of the individual results for z(t) and r(t) will be the same as the result when both of those are requested.

Extract the data points...

It can be done like this (I just picked a simple curve for the implicitplot).

P:= plots:-implicitplot(y= exp(x), x= -10..10, y= 0..10^4):
data:= op([1,1], P):
plots:-logplot(convert(data, listlist));

viewpoint option...

Change your final display command to

display([P, f(C)], viewpoint= ["circleright", frames= 32]);

See ?plot3d,viewpoint. The path of the "camera" can be arbitrary; I just picked an easy option above.

style= patchnogrid or style= wireframe...

Add the option style= patchnogrid or style= wireframe to your implicitplot3d command.

Arrays with an arbitrary number of dimen...

It can be done, for example, like this:

# Returns an Array with n dimensions each with extent 1..m
DynArray:= (n,m)-> Array((1..m)\$n, args[3..-1]):
A:= DynArray(3,4);

Be careful with this! A has m^n elements. 2^33 = 8 Gig of elements.

unapply...

plot3d(unapply(S,x,y), 0..1, 0..1, labels= [x,y,``]);

Browse...

To look at the contents of a large matrix, right-click on it and select Browse.

As to exporting to Excel: Are you saying that some entries in your Matrix are themselves lists? Do the lists have more than one element? If yes, do they all have the same number of elements?

Worksheet package...

This may be too tedious to be worth it: Use package ?Worksheet. This will read in a worksheet and convert it to a very large function tree. For the first case, find, by whatever ad hoc method, the function name that corresponds to the headings that you want. Then find, by experimenting with op, which operand or suboperand of that function is the heading string that you want.  Now you can extract all occurences that function name from any worksheet using indets(..., specfunc(anything, ...)), and then extract the heading string from those with op.

Two versions of the lattice polynomial...

Here are two versions of contour plots of the lattice polynomial. To deal with the extremely large values, I took the 25th root of the polynomial (since the polynomial is degree 25). One problem with the Mathematica plot is the both extremely positive and extremely negative values are shown in white.

Contourplots.mw

P:= mul(mul((x+I*y)-(a+b*I), a= -2..2), b= -2..2):
plots:-contourplot(
surd(Re(P),25), x= -3..3, y= -3..3,
grid= [100,100], contours= 20,
filled= true, coloring= [brown, green],
axes= boxed
);

plot3d(
surd(Re(P), 25), x= -3..3, y= -3..3,
grid= [100,100],
orientation= [-90,0,0],
);

Some more-elementary solutions...

LinearAlgebra:-IntersectionBasis is a great one-step solution to this problem. But for the sake of pedagogy, which may satisfy the Asker, and to reinforce my own learning, I present here several more solutions of this problem expressed in terms of more elementary operations.

 > restart;
 > with(LinearAlgebra):
 > a1:= < 3,1,4,1 >:  a2:= < 5,9,2,6 >:  a3:= < 5,3,5,8 >: b1:= < 9,7,9,3 >:  b2:= < 2,3,6,4 >:  b3:= < 6,2,8,4 >: A:= < a1 | a2 | a3 >:  B:= < b1 | b2 | b3 >:

Solution 0: (Technique already presented by Markiyan;  included here as a correctness check.)

 > Basis0:= IntersectionBasis([[a1,a2,a3], [b1,b2,b3]]);

Solution1: Use the reduced row-echelon form.

 > R:= ReducedRowEchelonForm(< A | -B >);

 > Rows, Cols:= op(1, R);

 > Basis1:= [Column](A . R[1..3,5..6], 1..Cols-Rows);

Verify that Basis1 and Basis2 are equivalent, i.e. that they have the same span. If the intersection of the spans of these bases is two-dimensional (i.e. the basis of the intersection has two members), then the spans are the same.

 > nops(IntersectionBasis([Basis1, Basis0]));

Solution 2: Use a parametrized solution from LinearSolve.

 > S:= LinearSolve(< A | -B >, ZeroVector(Rows), free= t);

 > P:= indets(S, name);

Get independent solution vectors by sequentially setting each parameter in P to 1 while the others are set to 0.

 > V:= `<|>`(seq(eval(S, {p= 1} union ((P minus {p})=~ 0)), p in P));

Solution 2A. Using columns of A.

 > Basis2A:= [Column](A . V[1..3,..], 1..Cols-Rows);

 > nops(IntersectionBasis([Basis2A, Basis1]));

Solution 2B. Using columns of B.

 > Basis2B:= [Column](B . V[4..6,..], 1..Cols-Rows);

 > nops(IntersectionBasis([Basis2A, Basis2B]));

Solution 3.

 > N:= NullSpace(< A | -B >);

 > Basis3:= (v-> B.v[4..6]) ~ (N);

 > nops(IntersectionBasis([Basis1, Basis3]));

The smoothest that I can get...

Here is the smoothest plot that I can get without using image-processing tools (expanding on an idea of Markiyan's presented in the thread he refers to above).

P3:= plots:-implicitplot(
y= sin(x), x= 0..2*Pi, y= -1.05..1.05,
grid= [32,32], gridrefine= 4, crossingrefine= 4, resolution= 2^12,
rangeasview, thickness= 2, color= brown,
axis[1]= [tickmarks= piticks]
):
op([1,1,2,1,2], P3); #Check number of points used.
84
print(P3);

eval(x(t), ...)...

Assign the solution given by dsolve to a variable:

soln:= dsolve(..., x(t));

where you have to substitute your actual dsolve command. Then, to plot it, do

plot(eval(x(t), soln), t= 0..10);

Less complicated way...

eval(V[dis], Units:-Unit= 1);

(A comment irrelevant to using the above command but essential to understanding why it works: The global Unit is initially assigned (i.e. before any packages are are loaded) to

proc() Units:-Unit(args) end proc,

and `print/Unit` is initially assigned to

proc()
local Unit;
if IsWorksheetInterface('Standard') or _EnvInRecursiveUnitsModulePrint = true  then Unit(args)
else Unit([args])
end if
end proc.

)

convert(..., `local`) and ``(...)...

You need to make a copy of the names x, y, and r with convert(..., `local`). You store the copy under a different name, and when that name is invoked, it produces the copy of the original name, which has no value assigned, even if the original original was assigned. Whew! That's probably very confusing, so here's an example:

restart;
B:= convert(a, `local`);
a
a:= 1;
1
B;
a
a;
1

The second thing that's needed for your plot is that we want the point to display as (x0, y0) rather than as x0y0. So we invoke that with the empty-name function as ``(x[0], y[0]), and the quote marks do not appear in the typeset display. One other change that I made to your code is to put the alignments inside the brackets [] of the individual text elements, because I think that the point should be aligned {below, left} and that the radius should be aligned {above, right}.

 > restart;
 > with(plots):
 > R:= convert(r, `local`): X:= convert(x, `local`): Y:= convert(y, `local`):   x[0]:= 1: y[0]:= 1: r[0]:= 1: plotCircule:= plot(    [x[0]+r[0]*sin(t), y[0]+r[0]*cos(t), t= -Pi .. Pi] ): plotRadius:= arrow(    {[[x[0], y[0]], [r[0]*cos((1/4)*Pi), r[0]*sin((1/4)*Pi)]]},    shape= arrow ): plotTags:= textplot([    [x[0], y[0],       'typeset'(``(X[0],Y[0])), 'align'= {'below', 'left'}    ],    [x[0] + r[0]*cos((1/4)*Pi), y[0] + r[0]*sin((1/4)*Pi),       'typeset'(R[0]), 'align'= {'above', 'right'}    ] ]): display([plotCircule, plotRadius, plotTags]);

 >