Carl Love

Carl Love

24663 Reputation

25 Badges

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

MaplePrimes Activity

These are answers submitted by Carl Love

Yes, it can be easily done with code that is much shorter and more efficient. Let be that set, or any other structure, to be extracted from. Do

indets(S, 'suffixed'("t", nonnegint) &under String)

The textbook example that you're following starts with two cosine functions of the same period, but different amplitudes and phase shifts. The problem that you're currently trying to solve has a sine and a cosine function of the same period (and different amplitudes and phase shifts). Fortunately, by a trivial identity, sin(x) = cos(x - Pi/2) for any x. Applying that, our two functions are now 

f:= t-> 3*cos(4*t - Pi/2);  g:= t-> 5*cos(4*t - 1);

Applying the textbook's algorithm now, we need the modulus (or abs) and argument (that thing with arctan) of

sc:= 3*exp(-I*Pi/2) + 5*exp(-I*1);

The command polar gives the modulus and argument of any complex-number expression together in one step:

             polar(7.697020822, -1.212177299)

And that's the amplitude and phase shift of f + g given as the textbook's answer.

I'm not sure exactly what you mean by "symbol", which is a formally defined type in Maple. Let e be the expression. You could do 

indets(e, symbol)

The problem with that is that it might leave out some things that you want to see. For example, the formal symbols of f(x[n], y) are fxy, and n, but the above will only return {y}. The type rules at work here may seem mysterious at first, but it's easy to understand when you know that

  1. x[n] is considered a name but not a symbol;
  2. symbol is a subset of name;
  3. indets doesn't dig down into indices (it does dig all the way to the bottom of expressions, but for some reason indices aren't included), so the n isn't found. 

So, instead, you could do

indets(e, {name, function})

The problem with that is that it may return structures that you think should be further broken down into symbols. For example, if e:= f(x[n], y), it'll return {y, x[n], f(x[n], y)}. You might want that broken down into {f, y, x, n}. That would be fairly easy to do if you tell me that that's what you want.

I don't think that it'd ever be possible to do that in Maple for several reasons. The foremost reason that comes to mind is that in Maple procedures are what some computer scientists call "first-class data structures": They can be manipulated like any other variable. I don't think that that's true in Java. So having a simple variable and a procedure with the same name would be essentially the same as having an integer and a float variable with the same name---and I'll assume that you understand why that couldn't happen in any language.

If by degree you mean the degree with respect to a particular variable, such as the x of your example, then it can be done with simplify by giving an extra argument (or arguments) called side relations, which are polynomial expressions assumed to be equal to 0:

p:= x^6 + 2*x^3 + x^2;
simplify(p, {x^5}); #must use { }
simplify(p, {x^3});

Unlike Christian's Answer, this works even if p is not itself a polynomial but rather is a more-general function of polynomials in x. And unlike Tom's Answer, it also works regardless of whether the polynomials within p are expanded into terms. (However, their Answers, unlike mine, would be fairly easy to modify to handle the multivariate case with degree meaning total degree.)

The side relations facility is much more general than these examples show; see help page ?simplify,siderels.

In your Question, you present the following vector of equations (or expressions implicitly equated to 0) to be solved for the w variables. The equations are linear in those variables.

equ:= < 
    0.85*((phi__cs*beta__s*f__con) . (D__pile*w__1)) + <-5/6*H - 5/8*P>, 
    0.85*((phi__cs*beta__s*f__con) . (D__pile*w__2)) + <5/6*H - 5/8*P>, 
    0.85*((phi__cs*beta__s*f__con) . (D__pile*w__3)) + <1/2*H - 3/8*P>

The form of these expressions is quite awkward. Before proceeding, it's mandatory to correct these 2 things:
1. `.` multiplication must be changed to `*`.
2. Internal vectors must be converted to scalars.

There's another optional correction that I'd like to make, because it makes the output nicer:
3. Convert decimal numbers to exact fractions.

All three corrections can be done automatically, regardless of the number of entries:

eqs:= [evalindets](
        convert(equ, rational), #Convert 0.85 to 17/20.
        specop(`.`), `*`@op #Correct improper multiplication.
    rtable, seq #Convert internal vectors to scalars; convert outer vector/matrix to list

Your equations are of course trivial to solve individually. However, I suspect that in general you will want to have equations that contain more than one of the w variables each. The code below will also handle that more general case (as long as they're still linear in those variables).

LA:= LinearAlgebra:
W:= indets(eqs, suffixed(w__))[]; #Extract decision variables.
<W> =~ LA:-LinearSolve(LA:-GenerateMatrix(eqs, [W]));

You say that the awkward format of your original equ is due to it being the output of some other program (or equation). But there is no good reason why that program should produce its output in that form. You should change that other program or get its author to change it. 

Constructing matrices from "blocks" (or submatrices) and permuting the rows and columns of matrices is much easier than anything shown in this thread so far, and it can be done using no packages (no LinearAlgebra, no ListTools, no linalg), no loops, and no index variables (such as the ij, and k that you used).

B1:= 24*<9, -12; -12, 16>; #2x2 block
Z:= Matrix(2,2); #2x2 zero matrix
#3x3 arrangement of 2x2 blocks to produce 6x6 matrix:
K1:= <B1, -B1, Z; -B1, B1, Z; Z, Z, Z>;
K2:= K1[[5.., ..4]$2]; #Apply permutation [5,6,1,2,3,4] to rows and columns.
B3:= <500, 0; 0, 0>;
K3:= <B3, Z, -B3; Z, Z, Z; -B3, Z, B3>;

#This creates a flattened form of your KG:
KG__flat:= <K||(1..3)>;

#If for some strange reason you actually want the non-flat
#KG that you originally had:
KG:= rtable(1..3, 1..1, [K||(1..3)], subtype= Matrix);

Here are some answers to some of your 3 additional questions:

1. How to flatten matrices?
Rhetorical answer: It's much easier to create a flat matrix in the first place (for example, as shown above) than it is to flatten a matrix whose entries are themselves matrices or vectors. So, is there some reason why you can't create them flat in the first place? I'm not saying that flattening is impossible, just that it'd be better to avoid it.

2. How to count the rows or columns of a matrix?
Answer: There are many ways to do it. Here are three ways to count the rows of a matrix M:

[op](1, M)[1];  [upperbound](M)[1];  rhs([rtable_dims](M)[1]);

To count the columns, replace [1] by [2] in any of the above.

3. How to delete rows and columns from a displayed matrix to facilitate the entering of new entries in 2D Input?
Non-answer: I don't know how to do things in 2D Input. However, it's trivial to delete rows or columns, or make them blank, or just re-enter them in 1D input. 


See ?ProgrammingGuide,Appendix1. In particular, go to "Internal Representation" => "Internal Representations of Data Types" => "BINARY".

BINARY is one of the 62 fundamental data types in Maple (63 if you count "temporary"). Many of those data types cannot have a pointer pointing at them that can exist at the user level. In other words, they cannot appear as the right side of an assignment statement in Maple code, nor can they be otherwise manipulated in Maple expressions; they can only be processed in internal kernel code. BINARY is one of those types.

kernelopts(dagtag = BINARY);

Error, object at address is binary

The call to p1 inside p2 should be 

p1(a, b, ':-debug'= debug)

The ':-debug' refers to the global keyword; the debug refers to p2's parameter. Putting a parameter in quotes, as your attempted 'debug', never does anything.

As far as I can tell, your main problem is mostly caused by the Maple GUI's "document mode". Anyway, I couldn't get it to work. I converted it to "worksheet mode" and 1-D input. I also removed a bunch of unnecessary junk. Both of these were fairly easy to do, and now it works fine.


(x[1],y[1]):= l[1]*~(sin,-cos)(theta(t));
(x[2],y[2]):= l[2]*~(-sin,cos)(theta(t));
(x[4],y[4]):= (x[1],y[1])+~l[4]*~(-sin,cos)((theta+phi)(t));
(x[3],y[3]):= (x[2],y[2])+~l[3]*~(sin,-cos)((theta-psi)(t));
R:= table([seq](k= [x[k],y[k]], k= 1..4));
dR:= diff~(R, t);

l[1]*sin(theta(t)), -l[1]*cos(theta(t))

-l[2]*sin(theta(t)), l[2]*cos(theta(t))

l[1]*sin(theta(t))-l[4]*sin(theta(t)+phi(t)), -l[1]*cos(theta(t))+l[4]*cos(theta(t)+phi(t))

-l[2]*sin(theta(t))-l[3]*sin(-theta(t)+psi(t)), l[2]*cos(theta(t))-l[3]*cos(-theta(t)+psi(t))

table( [( 1 ) = [l[1]*sin(theta(t)), -l[1]*cos(theta(t))], ( 2 ) = [-l[2]*sin(theta(t)), l[2]*cos(theta(t))], ( 3 ) = [-l[2]*sin(theta(t))-l[3]*sin(-theta(t)+psi(t)), l[2]*cos(theta(t))-l[3]*cos(-theta(t)+psi(t))], ( 4 ) = [l[1]*sin(theta(t))-l[4]*sin(theta(t)+phi(t)), -l[1]*cos(theta(t))+l[4]*cos(theta(t)+phi(t))] ] )

table( [( 1 ) = [l[1]*(diff(theta(t), t))*cos(theta(t)), l[1]*(diff(theta(t), t))*sin(theta(t))], ( 2 ) = [-l[2]*(diff(theta(t), t))*cos(theta(t)), -l[2]*(diff(theta(t), t))*sin(theta(t))], ( 3 ) = [-l[2]*(diff(theta(t), t))*cos(theta(t))-l[3]*(-(diff(theta(t), t))+diff(psi(t), t))*cos(-theta(t)+psi(t)), -l[2]*(diff(theta(t), t))*sin(theta(t))+l[3]*(-(diff(theta(t), t))+diff(psi(t), t))*sin(-theta(t)+psi(t))], ( 4 ) = [l[1]*(diff(theta(t), t))*cos(theta(t))-l[4]*(diff(theta(t), t)+diff(phi(t), t))*cos(theta(t)+phi(t)), l[1]*(diff(theta(t), t))*sin(theta(t))-l[4]*(diff(theta(t), t)+diff(phi(t), t))*sin(theta(t)+phi(t))] ] )

F:= [theta, phi, psi](t):
V:= [l[1], l[2], l[3], l[4], m[1], m[2], diff(F, t)[], (cos,sin)(psi(t)), g]:

T:= combine(simplify(collect(m[1]/2*inner(dR[4]$2)+m[2]/2*inner(dR[3]$2), V)), trig);

-(diff(theta(t), t))^2*l[1]*l[4]*m[1]*cos(phi(t))-(diff(theta(t), t))^2*l[2]*l[3]*m[2]*cos(psi(t))+(1/2)*(diff(theta(t), t))^2*m[1]*l[1]^2+(1/2)*(diff(theta(t), t))^2*m[2]*l[2]^2+(1/2)*m[1]*l[4]^2*(diff(phi(t), t))^2+(1/2)*m[1]*l[4]^2*(diff(theta(t), t))^2+(1/2)*m[2]*l[3]^2*(diff(psi(t), t))^2+(1/2)*m[2]*l[3]^2*(diff(theta(t), t))^2+m[1]*l[4]^2*(diff(theta(t), t))*(diff(phi(t), t))-m[2]*l[3]^2*(diff(theta(t), t))*(diff(psi(t), t))-(diff(theta(t), t))*(diff(phi(t), t))*l[1]*l[4]*m[1]*cos(phi(t))+(diff(theta(t), t))*(diff(psi(t), t))*l[2]*l[3]*m[2]*cos(psi(t))

U:= collect(simplify(expand(g*m[1]*y[4] + g*m[2]*y[3])), V);


L:= simplify(collect(T-U, V));

(1/2)*(-2*cos(phi(t))*l[1]*l[4]*m[1]-2*cos(psi(t))*l[2]*l[3]*m[2]+(l[1]^2+l[4]^2)*m[1]+m[2]*(l[2]^2+l[3]^2))*(diff(theta(t), t))^2+(1/2)*(-2*l[4]*m[1]*(cos(phi(t))*l[1]-l[4])*(diff(phi(t), t))+2*(diff(psi(t), t))*l[3]*m[2]*(cos(psi(t))*l[2]-l[3]))*(diff(theta(t), t))+(1/2)*m[1]*l[4]^2*(diff(phi(t), t))^2+(1/2)*m[2]*l[3]^2*(diff(psi(t), t))^2-g*((cos(phi(t))*l[4]*m[1]-cos(psi(t))*l[3]*m[2]-l[1]*m[1]+l[2]*m[2])*cos(theta(t))-sin(theta(t))*(sin(phi(t))*l[4]*m[1]+sin(psi(t))*l[3]*m[2]))

Eul:= remove(has, VariationalCalculus:-EulerLagrange(L, t, F), K[1]);

{sin(phi(t))*l[1]*l[4]*m[1]*(diff(theta(t), t))^2-g*(-sin(phi(t))*l[4]*m[1]*cos(theta(t))-sin(theta(t))*cos(phi(t))*l[4]*m[1])+l[4]*m[1]*(cos(phi(t))*l[1]-l[4])*(diff(diff(theta(t), t), t))-m[1]*l[4]^2*(diff(diff(phi(t), t), t)), sin(psi(t))*l[2]*l[3]*m[2]*(diff(theta(t), t))^2-g*(sin(psi(t))*l[3]*m[2]*cos(theta(t))-sin(theta(t))*cos(psi(t))*l[3]*m[2])-l[3]*m[2]*(cos(psi(t))*l[2]-l[3])*(diff(diff(theta(t), t), t))-m[2]*l[3]^2*(diff(diff(psi(t), t), t)), -g*(-(cos(phi(t))*l[4]*m[1]-cos(psi(t))*l[3]*m[2]-l[1]*m[1]+l[2]*m[2])*sin(theta(t))-cos(theta(t))*(sin(phi(t))*l[4]*m[1]+sin(psi(t))*l[3]*m[2]))-(2*(diff(phi(t), t))*sin(phi(t))*l[1]*l[4]*m[1]+2*(diff(psi(t), t))*sin(psi(t))*l[2]*l[3]*m[2])*(diff(theta(t), t))-(-2*cos(phi(t))*l[1]*l[4]*m[1]-2*cos(psi(t))*l[2]*l[3]*m[2]+(l[1]^2+l[4]^2)*m[1]+m[2]*(l[2]^2+l[3]^2))*(diff(diff(theta(t), t), t))-l[4]*m[1]*(diff(phi(t), t))^2*sin(phi(t))*l[1]+l[4]*m[1]*(cos(phi(t))*l[1]-l[4])*(diff(diff(phi(t), t), t))-(diff(diff(psi(t), t), t))*l[3]*m[2]*(cos(psi(t))*l[2]-l[3])+(diff(psi(t), t))^2*l[3]*m[2]*sin(psi(t))*l[2]}

    phi(0) = Pi/4, psi(0) = Pi/4, theta(0) = (3*Pi)/4,
    D(phi)(0) = 0, D(psi)(0) = 0, D(theta)(0) = 0

PARAM:= [g = 9.8, l[1] = 10, l[2] = 100, l[3] = 100, l[4] = 21, m[1] = 1000, m[2] = 1]:

sys:= eval(Eul, PARAM)[]:

sol:= dsolve([sys, INITS], numeric):

    sol, [t, diff(psi(t), t)], 1..3, numpoints= 1000,
    labels= ["times(sec)", "Psi&rsquo;(deg/s)"],
    labeldirections= [horizontal, vertical], font= [TIMES, ROMAN, 20]
    sol, [t, diff(theta(t), t)], t= 1..3, numpoints= 1000,
    labels= ["time(sec)", "Theta&rsquo;(deg/s)"],
    labeldirections= [horizontal, vertical], font = [TIMES, ROMAN, 20]






The correct expression is

`R[1]'`:= ...

Note that the paired outer quotes are different than the solitary inner quote. On my keyboard, the outer quotes are on the upper left, under Esc, and the inner quote is to the left of Enter.

Regarding your followup problem: From a set of sets S you want to extract all L-subsets such that 

  1. all pairwise intersections are size 1 (a.k.a., singleton),
  2. every element occurs in exactly 2 of the sets.

A procedure for that is

SingletonCliques:= (S::set(set), L::And(posint, Not({1,2})))->
local L2:= L*(L-1)/2, L3:= 2*L-3;
        s-> local p;
            nops(op~(s)) = L2 and andseq(nops(op~(p)) = L3, p= combinat:-choose(s,2)),
        combinat:-choose(select(nops = L-1, S), L)

Again, this is a purely set operation rather than a graph operation, and the edges are atomic with respect to it. You risk causing confusion by mentioning edges or other graph-theoretic aspects.

You wrote "I know this only affects the display," but that's not true. When you see a constant multiplier or divisor distributed over added terms, that's how the expression is actually stored, not just how it displays.

You could use an inert operator to avoid that:

(x^2 - 2*x - 1) %/ 4;

My interpretation of what you want is this:

d:= 2e-4:
theta:= (z,t)->
        * exp(-1.603200636*t)
T:= [seq](Pi/2+0.1*t, t= 0..7):
    theta~(z,T), z= 0..d,
    color= COLOR~(HUE, 0.87*T/(max-min)(T)), #0.87: see #1 below
    legend= (t=~ evalf[3](T))

#1: 0.87 is a fudge factor to prevent the HUE scale from wrapping 
# around from red at 0 back to red at 1. 


The operation that you asked for is a purely set operation---not a graph operation---and thus the edges are "atomic" with respect to it. In other words, it's irrelevant that the edges are themselves represented as sets; they might as well just be labels. 

This performs your requested operation:

IsPairwiseDisjoint:= (S::set(set))-> local s; evalb(add(nops(s), s= S) = nops(op~(S)))
FixedSizeBlocks:= (S::set(set), L::nonnegint)->
    select(IsPairwiseDisjoint, combinat:-choose(S,L))

Although that's the operation that you requested, I suspect that there's another set operation that'll actually work better for your graph-theoretic purpose. I'll write that up in a separate entry.

4 5 6 7 8 9 10 Last Page 6 of 361