Carl Love

## 25005 Reputation

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

## With < < > >...

seq(Matrix([[seq(a[k], k = x+m .. y+m)]]), m = 0 .. 2);

^^^^^^^

Obviously, the result of such a command is a sequence of Matrices, not a Matrix.

Here's my shortest code:

< seq(< a[x+m..y+m] >^%T, m= 0..2 ) >;

In particular, note the for any list L and valid indices a and b,

[seq(L[k], k= a..b)]

is the same thing as simply

L[a..b]

## principal branches and surd or RealDomai...

In Maple, when a negative real number is raised to a fractional power with an odd denominator, the result is not real. As it says at ?arithop or ?^,

For non-integer (complex) numeric powers y and (complex) numeric x, the expression x^y is evaluated as exp(y*ln(x)), where ln(x) is evaluated using the principal branch of the logarithm.

The function surd can be used to compute real odd roots of negative real numbers.  See ?surd for more information.

Since ln(-1) = I*Pi, the upshot of that first sentence is that for integer n, (-1)^(1/n) = cos(Pi/n) + I*sin(Pi/n):

1   1    (1/2)
- + - I 3
2   2

cos(Pi/3)+I*sin(Pi/3);
1   1    (1/2)
- + - I 3
2   2
But

surd(-1, 3);
-1

or

use RealDomain in (-1)^(1/3) end use;
-1

(see ?RealDomain).

## coeff(b, x, 0)...

In this case, you need to circumvent the automatic simplification of x^0 to 1. You do this with

coeff(b, x, 0).

## plots:-surfdata...

We make Points1 into an (M+1) x (N+1) x 3 listlistlist:

A:= [seq](Points1[(M+1)*k+1..(M+1)*(k+1)], k= 0..N);

plots:-surfdata(A);

And you can add the same options that you used for the pointplot3d command.

Please remove the with(linalg) from your worksheet. It works without it.

## map Re...

If they can be considered the same thing and also be considered the same as 1, then it is trivial:

Re ~ (sol);

Then proceed with whatever else you were doing with the repetitions.

## Easier to remove repetitions than to fin...

If your goal is simply to remove the repetitions from a list L, and you don't care about the order, then it is trivial: convert to a set:

{L[]};

If you want to keep it a list, but you still don't care about the order, it is still trivial: convert to a set, then convert back to a list:

[{L[]}[]];

When order doesn't matter, removing repetitions, whether or not they exist, is much easier than detecting or finding them: The former can be done entirely in Maple's kernel.

## A list is not an rtable...

An rtable is a Matrix, Vector, or Array. A list is not an rtable. FindRepetitions returns a list (so converting it to a list doesn't do anything). Why don't you use the numelems command that acer recommended? If your Maple version is less than 15, use nops for a list.

## Another three simple ways...

If the implicit-form solution can be solved for the independent variable, as it can in this case, then there are several more very easy ways to get it.  Here's three: (Sorry, explicit uploading of worksheets seems broken right now, so I am uploading this in plaintext (prettyprint=1).)

`restart;ode:= x - y(x)/diff(y(x), x) = y(x);                            y(x)                               x - -------- = y(x)                           d                                       --- y(x)                                  dx            (1)dsolve(ode);                       /        /     x    \      \             y(x) = exp|LambertW|- --------| + _C1|                       \        \  exp(_C1)/      /isolate(subs(y(x)= y, %), x);                     x = -y (-_C1 + ln(y))(2) The name of the conversion below is y_x regardless of the names actually used in the ODE.convert(ode, y_x);                       d              x(y)                      --- x(y) = -1 + ----                       dy              y  dsolve(%);                    x(y) = (-ln(y) + _C1) y(3) The order of  the substitutions below is significant.subs(diff(y(x),x)= 1/diff(X(y),y), y(x)= y, x= x(y), X= x,  ode);                             / d      \                        x(y) - y |--- x(y)| = y                             \ dy     /    dsolve(%);                    x(y) = (-ln(y) + _C1) y`

## y*conjugate(y) = abs(y)^2...

For any complex y, y*conjugate(y) = |y|^2 (although I would've guessed that you already knew this identity). The identity can be applied automatically via

simplify(z);

The resulting equation can be solved for abs(y) via

solve(%, abs(y));

indicating that there is a circle of solutions in the complex plane.

## while loop and Norm...

How about this? With the data in your worksheet, this took exactly 200 iterations. The resulting plot is a straight line (as well as can be seen with the naked eye). Note the the default for Norm is the infinity norm, which is equivalent to checking the absolute values of all of the entries. In other words, the convergence criterion is that the distance for both coordinates of all n-1 points is less than epsilon.

ErrX:= Vector(n, [0,1]):
ErrY:= Vector(n, [0,1]):
rho:= 1.2:  #relaxation factor
epsilon:= 10.^(-5):
while LinearAlgebra:-Norm(< ErrX, ErrY >) >= epsilon do
for k from 2 to n do
Xk:= X[k] - rho*GX(X, Y, k)/GXX(X, Y, k);
ErrX[k]:= Xk - X[k];
X[k] := Xk;
Yk:= Y[k] - rho*GY(X, Y, k)/GYY(X, Y, k);
ErrY[k]:= Yk - Y[k];
Y[k]:= Yk
end do
end do;

## Matrix filled by recursive proc...

restart;
m:= proc(r::nonnegint, c::nonnegint)
option remember;
if r*c=0 then 0 else c*thisproc(r-1,c)+a*thisproc(r-1,c-1) end if
end proc:
m(0,0):= 1:

M:= Matrix(5,5, m);

Matrix indices start at 1, not 0, so the (trivial) 0th row and column are not shown. If you want indices that start somewhere other thn 1, you need to use an Array. But then you won't get the neat tabular layout in the GUI.

## map, Threads, curry...

Your procedure myp can be reduced to

myp:= (a,b,c)-> zip(`+`, a, b)*c;

Is it possible to use something like
#map2(myp,[\$1..4],{[x,y,z,w],[x2,y2,z2,w2],[x3,y3,z3,w3]},[4,10,0.5]);
to achieve this?

seq(map2(myp, [\$1..4], {[x,y,z,w],[x2,y2,z2,w2],[x3,y3,z3,w3]}, s), s= [4,10,0.5]);

(I don't know why the line above breaks bad. It appears as one line in the editor, and I've tried several times to fix it.)

Is there a document I can read apart from the ?map to look for some general tips to 'simutanous' computing both symbolically and numerically in Maple?

See ?operators,elementwise and ?curry.

And is there any documentation on using multicore CPUs?

## rtable_num_elems(..., All)...

When applied to an rtable (Matrix, Vector, Array), nops is not the number of elements; rather, it is the number of operands in the internal structure. You can see the internal structure of A via op(A). For the number of elements, use rtable_num_elems(A, 'All').

## Some plots...

I'd like to show you the relevant commands, but I also want you to do your homework. So I'll compromise by posting the following without comments, most output elided, and without an attached file.

restart;
X:= t-> 3*cos(t):  Y:= t-> 2*sin(t):  t0:= 0.5:
kappa:= eval(VectorCalculus:-Curvature(<X(t),Y(t)>), t= t0):
dist:= (X1,X2)-> sqrt((X1[1]-X2[1])^2+(X1[2]-X2[2])^2):
slope:= (D(Y)/D(X))(t0):
x0:= X(t0):  y0:= Y(t0):
eq1:= (Cy-y0)/(Cx-x0) = -1/slope:
eq2:= dist([Cx,Cy],[x0,y0]) = 1/kappa:
fsolve({eq1,eq2}, {Cy,Cx});
{Cx = 4.13904333495107, Cy = 2.19319067267216}
fsolve({eq1,eq2}, {Cx,Cy}, avoid= %);
{Cx = 1.1264520363911743, Cy = -0.27548851825535093}
assign(%);
tan_eq:= slope*(x-x0)+y0:
Ellipse:= plot([X,Y,0..2*Pi], color= blue, scaling= constrained, thickness= 2):
Circle:= plot([Cx+cos(t)/kappa, Cy+sin(t)/kappa, t= 0..2*Pi], color= red):
Center:= plot([[Cx,Cy]], style= point, symbolsize= 25, symbol= cross, color= red):
Tangent:= plot(tan_eq, x= 0..4, color= green, thickness= 2):
plots:-display([Ellipse, Circle, Center, Tangent]);

## Missing a factor of V...

There are some problems with your terminology. You say gg is a vector. It doesn't make sense to raise a vector to a power. Perhaps g is a vector and gg represents the vector's dot product with itself, which would be a scalar. What you call the "gradient of psi" is the derivative of V with respect to psi. You can get this as

diff(V, psi);

But that will not express the derivative in terms of V.

You expression GVpsi is off by a factor of V. In other words, the numerator of the coefficient in front of the expression should be V instead of 1. It is difficult to get Maple to produce this expression as the derivative of V. However, you can verify it by subtracting diff(V, psi) and simplifying, obtaining 0. The expression can be derived "by hand" by a process called "logarithmic differentiation", which you can look up in a calculus textbook, which is often used to take the derivative of complicated expressions of the form f(x)^g(x).

 > restart;
 > V1:= (1-alpha+alpha*g^(1-1/psi))^(1/(1-1/psi));

 > dV:= diff(V1, psi);

 > dV2:= V*(alpha*V^(1/psi-1)*g^(1-1/psi)*ln(g)-ln(V))/psi^2/(1-1/psi);

 > eval(dV2, V= V1);

 > simplify(% - dV, symbolic);

 >