Carl Love

## 24648 Reputation

10 years, 47 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

## Working with assumed variables...

When you make assumptions on a variable, a copy of that variable is made. If you then assign to that same variable, then you'll have the original and the copy being used in different ways, and things get confusing. The takeaway advice for the inexperienced user is Don't assign to assumed variables.

• a:=solve(z1,a);
• z2;

You should change those to

• soln:= solve(z1, {a});
• eval(z2, soln);

When used in this way, eval has most of the desired effect of an assignment statement without the side-effect of changing the global state. You can of course assign the result of eval to another variable.

## blackbox...

Better than what you asked for---the solutions for a static list---we can make a "blackbox" solution that provides the x for any given y. This is what symbolic computation systems excel at.

 > restart;
 > eqn:= x = y*(25+2*x)^(2/3)/(25*(.566*(25*x)^(2/3)));

 > X:= Y-> fsolve(eval(eqn, y= Y), x);

 > X(1);

 > X(365);

 > plot(X, 1..365, labels= [y,x]);

 >

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);