Carl Love

Carl Love

24648 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

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.

)

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

 

Download plotcircle.mw

Note that the gridlines and subticks do not occur in the worksheet; that is just an effect of uploading to MaplePrimes.

In the usual case, names that occur inside procedure definitions are not evaluated until the procedure is invoked, even if the name has already been assigned at the time the procedure is defined. Example:

a:= 1:
f:= ()-> a;

                        f:=  () -> a
f();
                               1
a:= 2:
f();
                               2

One can force the names to be evaluated at the time the procedure is defined by using subs. Example:

restart;
a:= 1:
f:= subs('a'= a, ()-> a);

                          f:= () -> 1
a:= 2:
f();
                               1

The first case is what one wants in the vast majority of cases.

@acer It's not automatic simplfication. You can verify that by entering

'1/sqrt(2)';

and observing that the operation is not performed.

I think that the operation can be avoided by overloading the operator `^` for the case where the left operand is an integer and the right operand is a negative fraction. See ?overload. So something like

`save^`:= eval(`^`);
unprotect(`^`);
`^`:= overload([
   proc(a::integer, b::And(negative,fraction))
      option overload;
      'procname'(args) #return unevaluated
   end proc,
   eval(`save^`)
]);

I think that I'm missing some details there, but that's the basic idea.

 

Christopher wrote:

And I'm going to have to disagree with the rationalization for ease of divisibility because what about sqrt(7)/7 - it is not simplified to 1/sqrt(7).  It just so happens that any simplification of numerical values are reduced down to values with sqrts in the numerator.  Perhaps a debate will ensue about which is better.  Personally I think two numbers (one of them being a 1) is better than two numbers where neither is a 1.

Find an elementary school child who can do long division with pencil and paper and have that child do these two problems:

  1. 1.0000 / 2.646
  2. 2.646 / 7.

carrying out each to four decimal places. I guarantee you that that child will find the second problem easier. They'll probably even have the second problem finished in less time that it takes to get the second digit of the first problem.

I'm not saying that this is a good reason for the "simplification" these days; I'm saying that this is why the tradition of "rationalizing the denominator" exists.


Often a for loop can be replaced with seq and some small savings can be achieved that way. And certainly you can do that here. But, in this case, much greater savings can be made by performing the integral under assumptions on A, and then using unapply on the result to make F an operator (procedure).

restart;

F:= unapply(int(1/(1+A*cos(2*Pi*x))^3, x= 0..1), A)

   assuming A > 0, 1 >= A;

proc (A) options operator, arrow; (1/2)*(A^2+2)/((-A^2+1)^(1/2)*(A^4-2*A^2+1)) end proc

limit(F(A), A= 1, left);

infinity

F(1):= infinity;

infinity

plot(F, 0..0.9);

   

Download cos_cube_int.mw

Also, why do you have with(linalg)? You should never use that. The only reason it still exists is so that old code that uses it will still run.

The command that you want is ?zip.

This is rather tricky to do and even to describe: The Matrix that you pass to listcontplot should only contain the z values. The x and y values are determined implicitly by the position in the Matrix. I see that you have 1116 data points. Since 1116 = 36 x 31, I will suppose (just for the sake of exposition, here) that you have 36 x values and 31 y values. In that case your Matrix should have dimensions 31 x 36. 

If you post your matrix of (x,y,z) data, I will work on putting it into the required format.

There are many ways to do that. Which way you use may depend on

  1. what you want to do with the saved data (read it by eye, re-use it later in the program, etc.)
  2. what order you want to access the saved data (randomly, last-in-first-out, first-in-last-out, all at once, etc.)
  3. how much memory you have (can the main memory hold the data? do you need to write to disk? etc.)

The fundamental data structure for this---the one that is most versatile in most of the above cases---is called a table. It is so fundamental that it often does not need to explicitly declared.

Let's say you have this loop:

for k from 1 to N do
    ... some computations to create or modify a vector V ...
end do;

Add one line:

for k from 1 to N do
   ... some computations to create or modify a vector V ...;
   MySavedData[k]:= copy(V)
end do;

At any later time in the program, if you want to retrieve the vector for, say, k=50, do

W:= MySavedData[50];

or, if it is safe and convenient to overwrite V at this point, do

V:= MySavedData[50];

You may need to do the retrieval with copy(MySavedData[50]), depending on whether you intend to modify the vector again and whether you wish to also preserve the original.

Tensor and TensorSort are locals of the Physics module, not exports. You can access all module locals with the :- syntax by giving the command

kernelopts(opaquemodules= false);

at the top level. Many people include this command in their initialization file. After this command, you need to access them as Physics:-Tensor and Physics:-TensorSort, even though the Physics prefix may not appear in the original code.

The problem is too easy to use a computer for. There are 6 columns and only 5 dimensions, therefore, they are not a basis and the matrix is not invertible. There is no computation required to do the problem other than simply counting to 6.

The reason that it is taking so long is that it is first computing Eigenvalues(A) in exact arithmetic, and then applying evalf to that; rather than working in floating point from the start. With the all-in-floating-point approach, it is possible to work at Digits = 25 and get the eigenvalues of a 200x200 matrix in reasonable time. And, as acer suggests, it will be much faster still at Digits = 15.

I hope these examples explain everything. If not, please ask for clarification.

restart;
Digits:= 25:
with(LinearAlgebra):
A:= RandomMatrix(50,50):
Af:= Matrix(A, datatype= float):
CodeTools:-Usage(assign('E1', Eigenvalues(Af)));
memory used=349.99MiB, alloc change=0 bytes, cpu time=2.03s, real time=1.98s
CodeTools:-Usage(assign('E2', evalf(Eigenvalues(A))));
memory used=0.93GiB, alloc change=0 bytes, cpu time=7.38s, real time=7.20s
CodeTools:-Usage(assign('E3', Eigenvalues(A)));
memory used=341.92MiB, alloc change=0 bytes, cpu time=2.20s, real time=2.14s
CodeTools:-Usage(assign('E3', evalf(E3))):
memory used=0.59GiB, alloc change=0 bytes, cpu time=5.06s, real time=4.92s

So the time for E2 is the sum of the times for the two steps of E3.

Norm(sort(Re(E1)) - sort(Re(E2)));
                                    -21
                           3.4 10   
Norm(sort(Im(E1)) - sort(Im(E2)));
                                    -21
                           2.5 10   
Norm(E2 - E3);
                               0.

So E1 and E2 are approximately equal in values; whereas E2 and E3 are identical.


A:= RandomMatrix(200,200, datatype= float):
CodeTools:-Usage(Eigenvalues(A)):
memory used=16.58GiB, alloc change=24.00MiB, cpu time=108.56s, real time=104.13s
Digits:= 15:
A:= RandomMatrix(200,200, datatype=float[8]):
CodeTools:-Usage(Eigenvalues(A)):
memory used=401.12KiB, alloc change=0 bytes, cpu time=219.00ms, real time=484.00ms
A:= RandomMatrix(2^11, 2^11, datatype= float[8]):
CodeTools:-Usage(Eigenvalues(A)):
memory used=32.23MiB, alloc change=32.00MiB, cpu time=26.20s, real time=16.11s


It is obvious that there is linear dependency because the number of vectors (6) is greater than their dimension (5). It remains to find a dependency relation.

A:= Matrix([
   [3,1,4,1,5,9],
   [2,6,5,3,5,8],
   [9,7,9,3,2,3],
   [8,4,6,2,6,4],
   [3,3,8,3,2,7]
]):

A dependency relation means that there is a nonzero solution to A.x = 0 or equivalently that one column can be expressed as a linear combination of the others.

R:= LinearAlgebra:-ReducedRowEchelonForm(A);

R := Matrix(5, 6, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = -23/11, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 4/3, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 0, (3, 5) = 0, (3, 6) = 106/33, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1, (4, 5) = 0, (4, 6) = -214/33, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 1, (5, 6) = 50/33})

The significance of the reduced row echelon form is that R.x = 0 is essentially the same system of equations as A.x = 0, but the R form makes dependency relations obvious. A dependency relation is as follows, where the subscript notation shows how individual columns of a Matrix are specified in Maple.

add(R[k,6]*'A'[..,k], k= 1..5) = 'A'[..,6];

-(23/11)*A[() .. (), 1]+(4/3)*A[() .. (), 2]+(106/33)*A[() .. (), 3]-(214/33)*A[() .. (), 4]+(50/33)*A[() .. (), 5] = A[() .. (), 6]

Explicitly verify the dependency:

eval(%);

(Vector(5, {(1) = 9, (2) = 8, (3) = 3, (4) = 4, (5) = 7})) = (Vector(5, {(1) = 9, (2) = 8, (3) = 3, (4) = 4, (5) = 7}))

 

 

Download RREF.mw

An alternative to the explicit option is to apply allvalues to a result returned by solve which contains RootOfs.

solve({x=k/4,y=-k/3,z=3*k/8,2*x^2+3*y^2+4*z^2=9},{x,y,z,k});
allvalues(%);
 

You wrote:

I will have to draw some tick marks on for r=0 and r=1, in Inkscape. I don't know why there isn't an (obvious) option to have tickmarks but not lines, or to move the lines so that they look nice.

No need to use another program. The tickmarks and the gridlines can be specified independently. So you can get your r=0 and r=1 tickmarks like this:

polarplot([(1+cos(t)^2)*(1/2), cos(t)], t= 0..Pi/2,
   axis[radial]= [gridlines= .25*~[$1..3], tickmarks= 5],
   axis[angular]= [tickmarks= 5, gridlines],
   coordinateview= [DEFAULT, 0..Pi/2],
   angulardirection= clockwise, angularorigin= top
);

Download polarplot.mw

The plot options are not obvious because there are so many of them. I found this solution on page ?plot,axis.

First 344 345 346 347 348 349 350 Last Page 346 of 360