Carl Love

Carl Love

24805 Reputation

25 Badges

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

MaplePrimes Activity

These are answers submitted by Carl Love

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’(deg/s)"],
    labeldirections= [horizontal, vertical], font= [TIMES, ROMAN, 20]
    sol, [t, diff(theta(t), t)], t= 1..3, numpoints= 1000,
    labels= ["time(sec)", "Theta’(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.

I can't say for sure without seeing it in context, but I suspect that N1 is an initial condition, for example, the value of the formula for n=1. Every recursive formula will have at least one initial condition, although the value of n where it occurs isn't necessarily 1.

The Maple command rsolve can solve some recursion formulas that can be expressed via algebraic equations (in which case they're usually called recurrence relations or difference equatios). Observe here how the specification of initial conditions affects the results:

#No initial condition specified, so Maple chooses N(0) as the initial condition:
rsolve({N(n) = a+N(n-1)}, N(n));
                      N(0) - a + a (n + 1)

#Initial condition specified at n=0:
rsolve({N(n) = a+N(n-1), N(0)=N0}, N(n));
                       N0 - a + a (n + 1)

#Initial condition specified at n=1:
rsolve({N(n) = a+N(n-1), N(1)=N1}, N(n));
                      N1 + a (n + 1) - 2 a


I suspect that the memory explosion occurs when you do the matrix inverse (KL + a0.ML + a1.CL)^(-1). A matrix multiplication of the form X:= A^(-1).B is mathematically equivalent to X:= LinearAlgebra:-LinearSolve(A, B), but the latter is usually much more efficient.

I think that this does what you asked; let me know:

IsLabeledIsomorphicSubgraph:= proc(G1::Graph, G2::Graph)
uses GT= GraphTheory;
local V:= GT:-Vertices(G1), R;
    if {V[]} subset {GT:-Vertices(G2)[]} then 
        R:= GT:-IsSubgraphIsomorphic(G1, GT:-InducedSubgraph(G2, V), 'isomorphism');
        if [R][1] then subs(R[2], GT:-Edges(G1)) else false fi
end proc


@MaPal93 You asked:

  • Are you implying that when Maple's solve() seems to get stuck is because the system has no solution?

No, I'm implying that when there are more equations than unknowns, even just one more, then the system is almost certainly inconsistent, unless you have some strong theoretical reason for believing otherwise. Imagine drawing 3 curves in a plane. It's almost certain that there won't be a single point where all three intersect.

If we can find any inconsistent subset, under any values of parameters, then the whole system is inconsistent. So, the plan below is

1. Let Eq be the set of 12 equations (assumed here to be expressions implicitly equated to 0); let be the set of 6 variables that you want to solve for.
2. Extract the numerators of the equations; the denominators are irrelevant. Call this EqN.
3. Set any variables in EqN that are not in to 1, or any other simple, exact value, possibly 0. Call this EqP.
4. Find subsets of EqP of size 6 that contain all variables in V. Call this Eq6.

The above steps are very easy. They'll take 5 seconds or less in total, like this:

EqN:= numer~(Eq):
EqP:= eval(Eq, indets(EqN, name)=~ 1):
Eq6:= select(E-> V subset indets(E, name), combinat:-choose(EqP, 6)):

5. Take any member of Eq6 and try solve on it:

SolE:= {solve}(Eq6[1], V);

6. Remove any solution that contains free variables:

SolE2:= remove(s-> ormap(evalb, s), SolE);

7. If SolE2 is empty, go back to step 5 and pick another set of 6.
8. Take any solution in SolE2 and use it to evaluate the other 6 equations:

EqV:= simplify(eval(EqP minus Eq6[1], SolE2[1]));

9. If EqV = {0}, then, miraculously, the system is consistent, and you have a solution, at least for the parameter valuation used in step 3.
10. If EqV contains anything that is obviously not equal to 0, then the entire system is inconsistent; you're done; period.

  • Did you check my script?


I call something like f:= x*cos(x); ff:= x-> f, an attempted indirect parameter reference, the x being the parameter. Those references don't work for procedures created with -> or proc. The most-common way to handle them is ff:= unapply(f, [x]). To be more robust, you should pass x into q.

The quotes that you put around x don't do anything.

It can be done like this:

#Important! The next line has no minus sign!
C:= coeffs(expand(lhs(pde)), e^(gamma*tau), 'n'): 
sort(`[]`~(simplify(log[e^(-gamma*tau)]~([n]), symbolic), [C]));

The returned expression is a list of pairs. The 1st element of each pair is the power to which e^(-gamma*tau) is raised; the 2nd is the corresponding coefficient.

When parse returns an expression containing variables (i.e., names, identifiers, function names), and those variables weren't originally global, they won't be the same as the ones passed in.

The two types of substitution that you want can also be done with a simple loop, of course:

SublistSubs:= proc(S::(list=anything), L::list, {type2::truefalse:= false})
local s:= lhs(S), t:= rhs(S), r:= nops(s), T:= Vector(r--, [t], fill= ()), k, V:= <L>;
    for k to nops(L)-r do 
        if L[k..k+r] = s then 
            if type2 then V[k]:= t; k+= r else V[k..(k+= r)]:= T fi
end proc
    [a,b,c]= x, [2,4,5,1,a,4,5,a,b,434,a,b,c,5,5,5,1], 
    type2=~ [true,false]
    [[2, 4, 5, 1, a, 4, 5, a, b, 434, x, 5, 5, 5, 1], 
      [2, 4, 5, 1, a, 4, 5, a, b, 434, x, b, c, 5, 5, 5, 1]]


For any values of kB, and R, both x=-1 and x=1 are roots. No decimal approximations are needed to verify this; it can be done in exact arithmetic with symbolic k, B, and R:

simplify(eval~(P, x=~ [-1,1]));
             [0, 0]

Your expression initially appears quite lengthy. That's superficial. It can be vastly reduced in size by

P1:= simplify(normal(P)); 

And for the purpose of understanding its roots, it can be further reduced by

P1n:= numer(P1);

The expression is now easy to understand, and you can see why those are the roots for x. It's not totally obvious, but could easily be verified by, for example, pencil and paper and only high-school algebra. 

First 7 8 9 10 11 12 13 Last Page 9 of 363