3728 Reputation

17 Badges

6 years, 26 days

MaplePrimes Activity

These are replies submitted by mmcdara




is smarter than what I did to generate the density plot.

I vote up!


When this kind of thing happens ( which is almost the rule when you use 2D math [browse the related questions on mapleprimes]), a simple thing is to lprint the suspicious command.
So, after K__xva[1]+K__xva[1]; type lprint(%)
What you will see is nothing but whatKitonum obtained bycopying/pasting yout code in a 1D math worksheet.

Note that you would have obtained the same thing had you copied the suspicious line into a code snippet.

How to get rid of this?
Two solutions: stop using the 2D mathematical mode, or be even more careful when typing text (any removal of a character can potentially generate this type of problem)... which often means "rewrite the line"


You have removed a hell of a thorn from my side!


I like your second solution.
I vote up


"Is it because we need the denominator out before starting the elementary operations?"
In the present case all the denominators are equal and thus if SA is the SmithForm  of matrix A then the SmithForm SB of this matrix A with all these terms multiplied by the "something" is proportional to "something" times SA.
If this "something" is a polynomial whose leading term as a coefficient equal to C, than SB =(1/C) *~ SA :

A := Matrix([[1,2*x,2*x^2+2*x],[1,6*x,6*x^2+6*x],[1,3,x]]): # help pages' example
S := factor~(SmithForm(A));

p := randpoly(x, degree=3, coeffs = rand(-5..5));
B  := p*~A:
SB := factor~(SmithForm(B, x)):
simplify(SB /~ p);

Things are more complicated if not all the denominatores of the elements of A are equal.
Let N the matrix defined by N[i,j] = numer(A[i,j])
In this case the structure of Smith forms of A and N are the same (for instance zeros are at the same places), but the previous proportionnality no longer holds.
See for Instance


There are indeed two safer ways to use the aliases and I am very eager to receive this kind of advices.

About your final question "Why do you use concatenated names instead of indexed names?"
I found this (IMO) a safe way to avoid conflicts, for instance:

T := 1:  # here t represents a variable
plot( ..., t=0..T); # here I use the same T defined above
# fine length of the tube named T and of elements 1 and 2

L__T := 28 :
L[1] := 10: L[2]:=20

Another situation where I found using "__" useful is when the "main name" is protected; for instance:

# Temperatures of elements 1, 2, ...
 Temperature[1] := 10
Error, attempting to assign to array/table `Temperature` which is protected.  Try declaring `local Temperature`; see ?protect for details.

Temperature__1 := 10;

# accelerations

gamma[r] := 10;
Error, attempting to assign to array/table `gamma` which is protected.  Try declaring `local gamma`; see ?protect for details.

gamma__r := 10;

As soon as I have decided to work with variable names like L__T instead of L[..], it seemed "natural"  to use the concatenation symbol || to construct them programmatically. 


Thanks for this explanation.

About my "slightly suspicious objective":
I study a many degreees of freedom mechanical system. Each of these dof (a mass) is conventionnally refered to by a bame (not a number). 
For instance one can have names [A, B, C].

TO enhance readibility of the differential system whivh governs the dynamics of these masses, I thought to avoid representing the dependent variable.
Here is an example


mass_names := [A, B, C]

[A, B, C]


masses        := seq(M__||m, m in mass_names);
displacements := seq(x__||m, m in mass_names);

stiffnesses   := Matrix(3$2, (i,j) -> `if`(i=j, 0, K__||(cat(mass_names[i],mass_names[j]))));
dampings      := seq(xi__||m, m in mass_names);

masses := `#msub(mi("M"),mi("A"))`, `#msub(mi("M"),mi("B"))`, `#msub(mi("M"),mi("C"))`


displacements := `#msub(mi("x"),mi("A"))`, `#msub(mi("x"),mi("B"))`, `#msub(mi("x"),mi("C"))`


stiffnesses := Matrix(3, 3, {(1, 1) = 0, (1, 2) = `#msub(mi("K"),mi("AB"))`, (1, 3) = `#msub(mi("K"),mi("AC"))`, (2, 1) = `#msub(mi("K"),mi("BA"))`, (2, 2) = 0, (2, 3) = `#msub(mi("K"),mi("BC"))`, (3, 1) = `#msub(mi("K"),mi("CA"))`, (3, 2) = `#msub(mi("K"),mi("CB"))`, (3, 3) = 0})


xi__A, xi__B, xi__C


[ seq(alias(x__||m = (x__||m)(t)), m in mass_names) ]:

a := sort([indices(alias:-ContentToGlobal,':-nolist')]);

[x__A, x__B, x__C]


velocities := diff(a, t);

[diff(x__A, t), diff(x__B, t), diff(x__C, t)]


accelerations := diff(a, t$2);

[diff(diff(x__A, t), t), diff(diff(x__B, t), t), diff(diff(x__C, t), t)]


masses *~ accelerations =~ stiffnesses . `<,>`(displacements) +~ dampings *~ velocities^~2

Vector[column]([[M__A*(diff(x__A(t), `$`(t, 2))) = K__AB*x__B+K__AC*x__C+xi__A*(diff(x__A(t), t))^2], [M__B*(diff(x__B(t), `$`(t, 2))) = K__BA*x__A+K__BC*x__C+xi__B*(diff(x__B(t), t))^2], [M__C*(diff(x__C(t), `$`(t, 2))) = K__CA*x__A+K__CB*x__B+xi__C*(diff(x__C(t), t))^2]])




@Carl Love 

Thank you for this technical detail about your method

Regarding my Answer: I need to have the same bounds for the dimensions to be transposed, even if transposition of rectangular matrices  is well defined (which makes your answer all the more clever)

A := Array(1..3, 1..2, 1..2, (i, j, k) -> i+2*j+3*k):
TA := Array(1..3, 1..2, 1..2, (i, j, k) -> A[j, i, k]):

Error, Array index out of range
      TA[() .. (), () .. (), 1], TA[() .. (), () .. (), 2]

@Carl Love 

I vote up


I'm afraid you're right.
Unless one of a, b or c is 0, in such case solutions are easy to find, I've spent a lot of time trying to find a trick.

sys := [x*y*z + y*z + y = a, x*y*z + x*z + z = b, x*y*z + x*y + x = c]:
                      x y z + y z + y = a
                      x y z + x z + z = b
                      x y z + x y + x = c

# an eample of a particular case

solve(eval(sys, a=0), [x, y, z])
       [[                 2                           ]  
       [[      c        -b  + b c - b + c          b  ]  
       [[x = - -, y = - -----------------, z = - -----], 
       [[      b                c                b - c]  

         [                    b  ]]
         [x = c, y = 0, z = -----]]
         [                  c + 1]]

# only the second solution can lead to (x, y, z) integers (provided c+1 divides b)
# but not the first one:
#   as b divides c, let's rewrite it while setting c=K*b (K non null integer)
#   to see that z cannot be an integer

simplify(eval(%[1], c=K*b))
           [              K b + K - b - 1        1  ]
           [x = -K, y = - ---------------, z = -----]
           [                     K             K - 1]

One can also prove some facts such that

  • a=b=c implies x=y=z, but it remains to prove x is an integer (unless a=b=c=0 which gives x=y=z=0)
  • a=b implies x=c-a, y=z=0  (the other solutions can't be real)

I'd put some hope in observing that a+b+c involves the symmetric functions of the root of a 3rd degree equation


       3 x y z + x y + x z + y z + x + y + z = a + b + c

# let sigma__p the sum of the products of p roots of a 3rd degree equation

rel := a+b+c = 3*sigma__1+sigma__2-sigma__3

          a + b + c = 3 sigma__1 + sigma__2 - sigma__3

# thus x, y, and z must be the integer roots of function F below

F := (U-x)*(U-y)*(U-z);
F := collect(expand(F), U)
                    (U - x) (U - y) (U - z)
        3                 2                              
       U  + (-x - y - z) U  + (x y + x z + y z) U - x y z

... but it seemed to be a dead end.

As I know you are (were) involved in IMO, have you not ever heard of this problem?


At first glance I have had the same reaction as yours.
It was so obvious that I made myself to wonder if the question of the OP was not this one "How to choose a triple (a, b, c) of integers such that x, y and z would be integers too?"

Or stated differently: given a triple (a, b, c) of integers, and three functions f, g and h of (x, y, z), how can I check if the solution of the system {f(x,y,z)=a, g(x,y,z)=b, h(x,y,z)=c} wrt (x,y,z) is a integer?

Thiinvolve solving third degree equations in integer numbers, which is far from simple.
Unless there is some trick (the problem looks like some asked in mathematics problem books).



I've been lured by the way you did the computations.
I hadn't realized that 
In fact if H(n, y) is the nth Hermite polynomial relatively to the measure mu(y)dy, you have written 
HN(n, y) = sqrt(H(n, y)*mu(y)) which makes < HN(n, y), HN(m, y)>=0 wrt the measure dy.
So the HN(n, y)  are the Hermite-Gauss functions... which I hadn't realized first sight


I've just edit my answer

The code you give reminds me of classical problems in probability:

  • u7 is what is called the (type ) Gram-Charlier expansion of the density of probability of some RV X.
    In this case the coefficients a0...a6 are determined from the moments of X.
  • Sometimes X obey an algebraic/ODE/PDE relation (on aODE in your case ?) and a method to build the pdf/cdf of X is to use "chaos polynomials", for instance Hermite's.
    Here again we obtaine something close to what you wrote.

Any connection with this ?

@Mariusz Iwaniuk 

How do you deduce from your first table that the formula

sum(m!*x^(m - n)*LaguerreL(m, n - m, x)^2/n!, n = 0 .. infinity) = exp(z)

given by the OP is right?


Excellent, I vote up!

2 3 4 5 6 7 8 Last Page 4 of 86