Carl Love

Carl Love

27599 Reputation

25 Badges

12 years, 66 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

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.

You need to make it so that there is no tickmark on the radial axis that is equal to the full length of the axis. In this case, we need to remove the tickmark at r = 1. We do this with the axis[radial]= [tickmarks= ...] option. The other options below are to make the plot otherwise exactly like in your post, which I assume you did with the menus.

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


Download polarplot.mw


cos(3*x) = 13/14:

expand(%);

4*cos(x)^3-3*cos(x) = 13/14

sol:= [solve](%, cos(x)):

simplify(fnormal(evalf(sol -~ cos(arccos(13/14)/3))), zero);

[0., -1.59744226005789, -1.37849148916767]

So the first root returned by solve is the one we want. So the final answer is

cos(arccos(13/14)/3) = sol[1];

cos((1/3)*arccos(13/14)) = (1/28)*(2548+(588*I)*3^(1/2))^(1/3)+7/(2548+(588*I)*3^(1/2))^(1/3)

 


Download trigrad1314.mw

I made a minor change to deal with some of the long display labels. Here's a repost of the worksheet with the new code and some new examples. As for colored labels, I think something may be possible along those lines: colors (and other formatting information) can be encoded into names. I'll think about it.

FunctionToTree

Carl Love 23-Mar-2013

 

A module to convert a function tree representation (such as is returned by ToInert or by various commands in XMLtools) into a tree represented in the GraphTheory package

restart;

 

(* Written by Carl Love, 23-Mar-2013.

module () description "Convert a function tree to a GraphTheory tree."; local Vertices::table, VertexLabels::table, Vertex_index::nonnegint, Edges::table, AddEdge, AddVertex, AddSubTree, Prefix::nonnegint, Labels::table, StripName, NameVertices, ModuleInits, ModuleApply; export Legend::table; end module

Example 1

F:= ToInert(x+2*y*z^2);

_Inert_SUM(_Inert_PROD(_Inert_PROD(_Inert_NAME("y"), _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2))), _Inert_INTPOS(2)), _Inert_NAME("x"))

I came up with a scheme to abbreviate the vertex labels. In the line below, you need to map the second part of the _Inert_ function names to displayable labels. Note that the display label can be the empty symbol ``, which is useful for the leaf nodes, whose first operand becomes part of the label. Also note that this table must map strings (things in double quotes) to symbols (things in single backquotes).

Abbrs:= table([
   "SUM"= `+`, "PROD"= `*`, "POWER"= `^`,
   "EQUATION"= `=`, "EXPSEQ"= `,`, "RANGE"= `..`,
   "INTPOS"= ``, "NAME"= ``, "ASSIGNEDNAME"= ``,
   "FUNCTION"= Func
]):

The Prefix= 7 argument tells it to strip off the first 7 characters of each function name.

G:= FunctionToTree(F, Labels= Abbrs, Prefix= 7);

GRAPHLN(directed, unweighted, [`*`, `*`, `+`, `2`, `2`, `^`, x, y, z], Array(%id = 18446744073901591718), `GRAPHLN/table/48`, 0)

For any given tree, it's hard to say which command of DrawNetwork or DrawGraph will display it better. Neither is great.

GraphTheory:-DrawNetwork(G);

Note that the edge from `*` to x is partially obscured by the edge from `*` to  `^`, which makes it look like there is an edge from `^` to x. You can tell that there isn't by looking for the arrows.

 

The circular presentation of  DrawGraph ensures against obscures edges, but it looks less tree-like.

GraphTheory:-DrawGraph(G);

The module exports a table Legend that shows the correspondence between the vertex labels and the original functional representation.

eval(FunctionToTree:-Legend);

table( [( `+` ) = _Inert_SUM(_Inert_PROD(_Inert_PROD(_Inert_NAME("y"), _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2))), _Inert_INTPOS(2)), _Inert_NAME("x")), ( z ) = _Inert_NAME("z"), ( `*` ) = _Inert_PROD(_Inert_PROD(_Inert_NAME("y"), _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2))), _Inert_INTPOS(2)), ( y ) = _Inert_NAME("y"), ( x ) = _Inert_NAME("x"), ( `*` ) = _Inert_PROD(_Inert_NAME("y"), _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2))), ( `2` ) = _Inert_INTPOS(2), ( `2` ) = _Inert_INTPOS(2), ( `^` ) = _Inert_POWER(_Inert_NAME("z"), _Inert_INTPOS(2)) ] )

These can be converted to their original expressional forms:

FromInert ~ (%);

table( [( `+` ) = 2*y*z^2+x, ( z ) = z, ( `*` ) = 2*y*z^2, ( y ) = y, ( x ) = x, ( `*` ) = y*z^2, ( `2` ) = 2, ( `2` ) = 2, ( `^` ) = z^2 ] )

Example 2

F:= ToInert([1,2,3,4]);

_Inert_LIST(_Inert_EXPSEQ(_Inert_INTPOS(1), _Inert_INTPOS(2), _Inert_INTPOS(3), _Inert_INTPOS(4)))

G2:= FunctionToTree(F, Labels= Abbrs, Prefix= 7);

GRAPHLN(directed, unweighted, [`,`, `1`, `2`, `3`, `4`, LIST], Array(%id = 18446744073901597270), `GRAPHLN/table/53`, 0)

GraphTheory:-DrawNetwork(G2);

eval(FunctionToTree:-Legend);

table( [( `,` ) = _Inert_EXPSEQ(_Inert_INTPOS(1), _Inert_INTPOS(2), _Inert_INTPOS(3), _Inert_INTPOS(4)), ( `4` ) = _Inert_INTPOS(4), ( `1` ) = _Inert_INTPOS(1), ( `3` ) = _Inert_INTPOS(3), ( `2` ) = _Inert_INTPOS(2), ( LIST ) = _Inert_LIST(_Inert_EXPSEQ(_Inert_INTPOS(1), _Inert_INTPOS(2), _Inert_INTPOS(3), _Inert_INTPOS(4))) ] )

FromInert ~ (%);

table( [( `,` ) = 1, 2, 3, 4, ( `4` ) = 4, ( `1` ) = 1, ( `3` ) = 3, ( `2` ) = 2, ( LIST ) = [1, 2, 3, 4] ] )

Example 3

J:= Int(x, x= 0..1);

Int(x, x = 0 .. 1)

F:= ToInert(J);

_Inert_FUNCTION(_Inert_NAME("Int", _Inert_ATTRIBUTE(_Inert_EXPSEQ(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))), _Inert_NAME("_syslib")))), _Inert_EXPSEQ(_Inert_NAME("x"), _Inert_EQUATION(_Inert_NAME("x"), _Inert_RANGE(_Inert_INTPOS(0), _Inert_INTPOS(1)))))

G3:= FunctionToTree(F, Labels= Abbrs, Prefix= 7);

GRAPHLN(directed, unweighted, [`,`, `..`, `0`, `1`, `=`, Func, Int, x, x], Array(%id = 18446744073901600270), `GRAPHLN/table/56`, 0)

GraphTheory:-DrawGraph(G3);

GraphTheory:-DrawNetwork(G3);

eval(FunctionToTree:-Legend);

table( [( `,` ) = _Inert_EXPSEQ(_Inert_NAME("x"), _Inert_EQUATION(_Inert_NAME("x"), _Inert_RANGE(_Inert_INTPOS(0), _Inert_INTPOS(1)))), ( `..` ) = _Inert_RANGE(_Inert_INTPOS(0), _Inert_INTPOS(1)), ( x ) = _Inert_NAME("x"), ( `1` ) = _Inert_INTPOS(1), ( `=` ) = _Inert_EQUATION(_Inert_NAME("x"), _Inert_RANGE(_Inert_INTPOS(0), _Inert_INTPOS(1))), ( x ) = _Inert_NAME("x"), ( Int ) = _Inert_NAME("Int", _Inert_ATTRIBUTE(_Inert_EXPSEQ(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))), _Inert_NAME("_syslib")))), ( `0` ) = _Inert_INTPOS(0), ( Func ) = _Inert_FUNCTION(_Inert_NAME("Int", _Inert_ATTRIBUTE(_Inert_EXPSEQ(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))), _Inert_NAME("_syslib")))), _Inert_EXPSEQ(_Inert_NAME("x"), _Inert_EQUATION(_Inert_NAME("x"), _Inert_RANGE(_Inert_INTPOS(0), _Inert_INTPOS(1))))) ] )

FromInert ~ (%);

table( [( `,` ) = x, x = 0 .. 1, ( `..` ) = 0 .. 1, ( x ) = x, ( `1` ) = 1, ( `=` ) = x = 0 .. 1, ( x ) = x, ( Int ) = Int, ( `0` ) = 0, ( Func ) = Int(x, x = 0 .. 1) ] )

Example 4

F:= cos(x+2*y)^2+1;

cos(x+2*y)^2+1

F:= ToInert(F);

_Inert_SUM(_Inert_POWER(_Inert_FUNCTION(_Inert_ASSIGNEDNAME("cos", "PROC", _Inert_ATTRIBUTE(_Inert_EXPSEQ(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))), _Inert_NAME("_syslib")))), _Inert_EXPSEQ(_Inert_SUM(_Inert_NAME("x"), _Inert_PROD(_Inert_NAME("y"), _Inert_INTPOS(2))))), _Inert_INTPOS(2)), _Inert_INTPOS(1))

G4:= FunctionToTree(F, Labels= Abbrs, Prefix= 7);

GRAPHLN(directed, unweighted, [`*`, `+`, `+`, `,`, `1`, `2`, `2`, Func, `^`, cos, x, y], Array(%id = 18446744073901605462), `GRAPHLN/table/61`, 0)

GraphTheory:-DrawGraph(G4);

GraphTheory:-DrawNetwork(G4);

FromInert ~ (FunctionToTree:-Legend);

table( [( cos ) = cos, ( `+` ) = cos(x+2*y)^2+1, ( `,` ) = x+2*y, ( `*` ) = 2*y, ( y ) = y, ( `1` ) = 1, ( x ) = x, ( `+` ) = x+2*y, ( `2` ) = 2, ( `^` ) = cos(x+2*y)^2, ( `2` ) = 2, ( Func ) = cos(x+2*y) ] )

 


Download FunctionTree.mw

Your expression for the integrand has unbalanced parentheses. I think you're missing a left parenthesis "(" at the very beginning. Once you repost a correction, I might be able to help you more. But I think that it's unlikely that you will get a symbolic answer.

The command is plots:-listcontplot(M, contours= [$1..20]). But you need to lay out your points in a grid form, which involves some sorting and approximating. I vaguely recall seeing a tool for this a many years ago. I will look for it.

Does the error occur in CalculateMatrix or after? Put a print statement after the CalculateMatrix call to find out.

It is a system of linear equations. It would be much more efficient to solve it with matrix operations (i.e. LinearAlgebra:-LinearSolve). Would a floating-point solution be acceptable? That would use much less memory.

You can convert the arrows, or anything else, to Cartesian coordinates, and put them on a polar plot like this:

 

N:= 12:

magnitude:= RandomTools:-Generate(list(float(range= 0..1, method= uniform), N)):

phase:= RandomTools:-Generate(list(float(range= 0..360, method= uniform), N)):

A:= < < magnitude[] > | < phase[] > >:

ExcelTools:-Export(A, "C:/Users/Carl/desktop/data.xls"):

A:= convert(ExcelTools:-Import("C:/Users/Carl/desktop/data.xls"), Matrix):

Convert to radians:

A[.., 2]:= evalf(Pi)/180 *~ A[.., 2]:

Polar to Cartesian converter:

xy:= P-> < P[1]*cos(P[2]), P[1]*sin(P[2]) >:

p1:= plots:-arrow([seq](xy(A[k,..]), k= 1..N), shape= arrow, color= red):

Dummy plot to give us the polar background:

p2:= plots:-polarplot(1.1*max(A[.., 1]), color= white, transparency= 1):

plots:-display(p1,p2);


Download polarplot.mw

Use Sum (with a capital S), instead of sum, and you will get the answer that you expect.

I have no idea what that escaped variable _O and the rest of that huge mess represents. Surely the result of differentiating the lowercase sum is a bug. However, I will note that the fault lies with sum, not diff, because you get a nearly identical mess if you just do the derivative in your head and enter

sum(sin(theta[j](t))*cos(theta[j](t))*diff(theta[j](t), t), j= 1..i-1);

First 376 377 378 379 380 381 382 Last Page 378 of 392