Carl Love

## 24805 Reputation

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

## copy, ArrayTools, and ?rtable_indexing...

First off, please do not use the old lowercase constructor commands matrix, array, vector. The new constructors are capitalized---Matrix, Array, Vector---and they are much more efficient and versatile. I will not help you anymore if you continue to use the lowercase constructors, as I'd rather just forget about them. The rest of this post will deal strictly with the newer versions. Also, anything that I say below about Arrays could equally well be said about Matrices or Vectors.

You wrote:

for l to 5 do: for i to 6 do: for j to 6 do: L[i,j] := P[i,j,l]: od: od: V[l]:=L: od:
The code fails in that V[1]=V[2]=...v[5] = P[i,j,5].

The immediate problem is that V[l]:= L needs to be V[l]:= copy(L); however, assigning individual entries in a loop is an inefficient way to do this (see below for a more efficient way). The reason the copy is needed is a somewhat arcane topic. You can read some about it at ?copy, but the topic requires a moderate amount of computer programming knowledge to understand. The take-away is Don't directly assign one Array to another (unless you understand the consequences); it does not make an independent copy of the Array.

There are many commands in the ?ArrayTools package for efficiently handling things like changing dimensions, copying parts of Arrays, etc.

You wrote:

Also, given a matrix, is there a way to extract a column vector from that matrix? I.e. in FORTAN is M is a matrix, M(:,3) would return the third column of M. Is there a similar command in Maple?

Yes, and the syntax is almost identical: M(.., 3) or M[.., 3]. The available combinations for selecting and rearranging subparts of Arrays is vast, and discussed at ?rtable_indexing.

So, applying this subpart selecting to your original loop, you could do the whole job as

V:= Array(1..5, i-> P[.., .., i]);

## Post example...

My guess (and it's only a guess until I see the example) is that you have a parameter with the same name as a global variable and that you're unintentionally using both in the main function. I'd be able to tell in 5 seconds if I saw an example. Making it a procedure will not help. You should never have to cut-and-paste to fix these things; almost everything can be handled programmatically/algebraically in Maple.

## Generalization makes this much simpler...

Sometimes---more often than I would expect---solving a more general problem is easier yet still gives a complete solution to the original problem. Such is the case here. The problem at hand is equivalent to finding all directed paths of a given length in a graph. What you call routes are formally called paths: a sequence of distinct vertices, each adjacent to the preceding.

Once the data structure is defined as I have done in the procedure below, finding the paths is essentially a one-liner. There's no need to consider special cases like vertices on corners, sides, etc.

AllPaths:= proc(
G::GraphTheory:-Graph
,maxlength::nonnegint:= GraphTheory:-NumberOfVertices(G)
)
#Returns a table R such that R[k] is a list of all directed paths of length k in G,
#where k ranges from 0 to maxlength. G can be directed or undirected.
#Each path is represented as a list of vertices.
uses GT= GraphTheory;
local
k  #current length
,p  #current path
,v  #current vertex
,V:= GT:-Vertices(G)  #constant
,E:= table(V =~ (`{}`@op)~ (GT:-Departures(G)))  #constant table of edges
,R:= table([0= `[]`~ (V)])  #Initialize with zero-length paths.
;
for k to maxlength do
#For each path of length k-1, make all possible 1-vx extensions.
R[k]:= [seq](seq([p[], v], v in E[p[-1]] minus {p[]}), p in R[k-1])
end do;
eval(R)
end proc:

Now apply it to the problem at hand:

squares:= [seq](seq([i,j], i= 1..4), j= 1..4):

Name the vertices conveniently.

V:= (v-> cat(['a','b','c','d'][v[1]], v[2]))~ (squares);

V := [a1, b1, c1, d1, a2, b2, c2, d2, a3, b3, c3, d3, a4, b4, c4, d4]

Vsq:= table(V =~ squares):

Define the edges of the graph. Horizontal, vertical, and diagonal adjacency are all covered by simply saying that the infinity norm of the difference of the points is 1.

Sq44:= GraphTheory:-Graph(
select(e-> LinearAlgebra:-Norm(<Vsq[e[1]] - Vsq[e[2]]>) = 1, combinat:-choose({V[]}, 2))
):

Generate the paths and time it.

CodeTools:-Usage(assign(Pths, AllPaths(Sq44, 6))):

memory used=11.31MiB, alloc change=4.27MiB, cpu time=219.00ms, real time=210.00ms

Explore the results a bit.

nops(Pths[6]);

68272

Get a well-distributed sample of the paths.

Samp:= [seq](Pths[6][2^k], k= 1..16);

Samp := [[a1, a2, a3, a4, b3, b2, c1], [a1, a2, a3, a4, b3, b2, c3], [a1, a2, a3, a4, b3, c2, b2], [a1, a2, a3, a4, b3, c3, c2], [a1, a2, a3, a4, b4, c3, c4], [a1, a2, a3, b2, b3, c3, d4], [a1, a2, a3, b2, c3, d3, c4], [a1, a2, a3, b4, c3, b3, c4], [a1, a2, b1, c2, c3, b4, a3], [a1, a2, b3, b4, a3, b2, c2], [a1, b1, c1, c2, d3, c3, b3], [a2, a1, b2, a3, b3, c2, d1], [a2, b3, c3, d2, c1, c2, b2], [a4, b4, c3, c2, b1, c1, b2], [b4, c3, d2, c1, b1, a2, a1], [d4, c3, d2, c2, b3, b4, a3]]

Develop a few tools for visualizing the paths.

To center a vertex in its square:

pt:= v-> Vsq[v] -~ 1/2:

Background to plot on:

grid:= plot(0, view= [0..4, 0..4], gridlines, scaling= constrained):

pathplot:= p->
plots:-display(
grid
,plot(
[[pt](p[1]), [pt](p[-1]), pt~(p)]
,style= [point, point, line], color= [green, red, yellow]
,symbol= diamond, symbolsize= 30
)
)
:

And finally here are the 16 sample paths:

plots:-display(Matrix(4,4, pathplot~ (Samp)));

(Sorry, I don't know how to upload an array of plots. If you execute the above code yourself, you will get a 4x4 array of plots---one for each of the sample paths.)

## Correct. Now what is radius?...

Correct, the shell method is easier in this case.

Give Maple the following command to plot the region in question:

plots:-implicitplot([y= x^2+4*x+6, x=1, y=2, x=3], x= 0..4, y= 0..30);

Suppose you took just one arbitrary point (x,y) inside the region and rotated that point in the manner suggested by the problem. What would be the radius of the circle made by that rotation?

## Do you think it is shell, or disc?...

Which method do you think will work better for this problem: the shell method, or the disc method?

## display an array...

I don't know if you mean two plots side-by-side, each composed of three plots; or if you mean a 3 x 2 array of plots. The former can be done with

display(< display(A) | display(B) >);

The latter can be done with

display(< < A > | < B > >);

## flexp needs to be a function...

I found two problems. The first one I don't really understand, but anyway, I needed to retype Zpratt in the plot command. I could tell there was something wrong with it because the "pratt" was not in italics.

The second problem is more mathematical. You need to make flexp a function of k, i.e. you need a "k ->" in front of the expression. Then when you refer to flexp in the two final Z function definitions, make it flexp(k).

As for putting the plots together, the command is plots:-display, not simply display. So,

plots:-display([A,B, etc.]);

Finally, there are linestyles other than dash. They are solid, dot, dash, dashdot, longdash, spacedash, or spacedot.

## Can't be plotted, it's not a function...

The integral is simply infinity. It doesn't depend on any variable. The x is gone after you integrate.

## piecewise...

Preben's answer is great, and it will work for procedures in general. I just want to direct you to a great command that is made especially for your kind of function: It's called ?piecewise.

## Return <> return...

Maple is case sensitive, and Return is not the same thing as return, which is what Maple uses. (Indeed, all of Maple's ?keywords are in lowercase.) One of the procedures must use the capitalized version.

Just guessing here, but was one of the sets of procedures translated automatically from Mathematica?

## expand, but with "frontend"...

You need to expand (i.e., distribute sin(z) over the rest) the first expression; otherwise, it won't contain the subexpression that you told applyrule to look for. But simply using the command expand will do too much: It will also apply a trig identity expand cos(2*z). So do this

frontend(expand, [exprsn1]);

Then use applyrule on the result of this.

## How to plot...

First, get rid of the line k:= 0, 1..npts/2. Do a restart to eliminate the effect of that command. Then, at the end, do

plot(Zuncomp(k), k= 0..npts/2, labels= ["k", "Zuncomp"])

## Another workaround...

Another workaround is to set Digits:= 20.

## Assumptions are not cumulative...

Assumptions via assume on one variable are not cumulative. Your second assume command erases the effect of the first. The second one is unnecessary anyway: It is implied by the first. If you actually want to make multiple assumptions on a single variable, put them together in one assume command or use the additionally command. See the examples at ?assume for details. If you remove the second assume, you will get your expected result.

By the way, you misspelled alpha as alph, but that doesn't change the result.

I agree with Preben about assuming. I think that it makes for cleaner, easier-to-read, less-error-prone, and more-flexible worksheets. Although there is still a place for assume if you're trying to get Maple to do a symbolic proof.

## Bug in Standard GUI...

This is a bug in the Standard GUI. Your commands run fine in the command-line Maple. If you have access to the Classic GUI (I don't), try that also.

A warning to anyone trying to investigate this: If you try to probe this, say with printlevel:= 100, the GUI will lose connection with the kernel (without any message), and you'll be stuck. Setting interface(prettyprint= 0) does not help. Another curiosity: The output, i.e. the returned value, of the OP's commands is literally the symbol `Non-fatal error while reading data from kernel.`; that is what soltn is set to; it is not an error message in the usual sense. Therefore, it must the kernel that is sending the message!

 First 356 357 358 359 360 361 362 Page 358 of 363
﻿