For this you can use GenerateMatrix in the LinearAlgebra package.

For this you can use GenerateMatrix in the LinearAlgebra package.

Proving that a number is irrational is generally very difficult (with some exceptions, notably algebraic numbers).  This is not an area where an amateur is likely to succeed. 

You wrote:

For instance, if you had an infinite number of irrational numbers and half of them were the sqrt(2) and the other half were –sqrt(2), their sum would not have to be irrational, would it?

In that case, their sum would be neither rational nor irrational, it would simply not exist, i.e. it would be a divergent series.

Something like this (where R[i,j] is the rate for car i passing car j when car j is in front of car i, and these are given numbers):


> G:= Matrix(N!, N!, storage=sparse, datatype=float[8]);
   perms:= combinat[permute](N):
   for i from 1 to N! do
     permi:= perms[i];
     for k from 1 to N-1 do
       permj:= subsop(k=permi[k+1],k+1=permi[k],permi);
       member(permj,perms,'j');
       G[i,j]:= R[permi[k],permi[k+1]];
       G[i,i]:= G[i,i] - G[i,j];
     end do;
    end do:

To find the equilibrium vector:

> pi:= op(1,LinearAlgebra:-NullSpace(G^%T));
   pi:= pi/add(pi[i],i=1..N!);

Something like this (where R[i,j] is the rate for car i passing car j when car j is in front of car i, and these are given numbers):


> G:= Matrix(N!, N!, storage=sparse, datatype=float[8]);
   perms:= combinat[permute](N):
   for i from 1 to N! do
     permi:= perms[i];
     for k from 1 to N-1 do
       permj:= subsop(k=permi[k+1],k+1=permi[k],permi);
       member(permj,perms,'j');
       G[i,j]:= R[permi[k],permi[k+1]];
       G[i,i]:= G[i,i] - G[i,j];
     end do;
    end do:

To find the equilibrium vector:

> pi:= op(1,LinearAlgebra:-NullSpace(G^%T));
   pi:= pi/add(pi[i],i=1..N!);

For the N-car case, the N! states are all permutations of the N cars.  You have a possible transition from one state to another if the second is obtained from the first by interchanging two adjacent elements in the permutation; if those elements are i and j, the matrix element is R[i,j].  It's a big matrix, but sparse: only N nonzero entries in each row or column, including the diagonal entries which make the sum over each row 0.  The sparseness would allow Maple to store the Matrix in a reasonable amount of memory, if you use the option storage=sparse.  The MatrixExponential will not be feasible to compute, because that would not be a sparse matrix.  Perhaps a particular row of the MatrixExponential, corresponding to the probabilities at time t after starting in a particular state, could be done with dsolve/numeric, though I'm not sure it would be up to solving such a huge system.  But implementing, say, a classical Runge-Kutta method to solve it wouldn't be too hard.  The equilibrium vector should be OK, I think, using LinearSolve.

For the N-car case, the N! states are all permutations of the N cars.  You have a possible transition from one state to another if the second is obtained from the first by interchanging two adjacent elements in the permutation; if those elements are i and j, the matrix element is R[i,j].  It's a big matrix, but sparse: only N nonzero entries in each row or column, including the diagonal entries which make the sum over each row 0.  The sparseness would allow Maple to store the Matrix in a reasonable amount of memory, if you use the option storage=sparse.  The MatrixExponential will not be feasible to compute, because that would not be a sparse matrix.  Perhaps a particular row of the MatrixExponential, corresponding to the probabilities at time t after starting in a particular state, could be done with dsolve/numeric, though I'm not sure it would be up to solving such a huge system.  But implementing, say, a classical Runge-Kutta method to solve it wouldn't be too hard.  The equilibrium vector should be OK, I think, using LinearSolve.

Another possibility is to produce separate plots and combine them with display in the plots package.

> eps:= .0001;
   f:= piecewise(x < 0, sin(x), x < 2*Pi, cos(x), 2);
   plots[display](plot(f, x=-5..0-eps, colour=red,linestyle=1,thickness=3),
                       plot(f, x=0..2*Pi-eps, colour=blue, linestyle=2, thickness=3),
                       plot(f, x=2*Pi .. 10, colour=black, linestyle=3, thickness=3));

Another possibility is to produce separate plots and combine them with display in the plots package.

> eps:= .0001;
   f:= piecewise(x < 0, sin(x), x < 2*Pi, cos(x), 2);
   plots[display](plot(f, x=-5..0-eps, colour=red,linestyle=1,thickness=3),
                       plot(f, x=0..2*Pi-eps, colour=blue, linestyle=2, thickness=3),
                       plot(f, x=2*Pi .. 10, colour=black, linestyle=3, thickness=3));

The problem with MatrixExponential in case 3, I think, is that the characteristic polynomial has an irreducible cubic factor, and MatrixExponential insists on using the complicated explicit formula for the roots of that factor rather than a more manageable RootOf.  We can do better "by hand". 

> with(LinearAlgebra):    
   Q3:= eval(Q,parameters3);
   CP:= factor(CharacteristicPolynomial(Q3,s));

> e1:= 0:
alias(e2=RootOf(op(2,CP),s,index=1),
  e3=RootOf(op(2,CP),s,index=2),
  e4=RootOf(op(3,CP),s,index=1),
  e5=RootOf(op(3,CP),s,index=2),
  e6=RootOf(op(3,CP),s,index=3));
 E:=map(simplify,Matrix([seq(op(NullSpace(Q3-e)),e=[e1,e2,e3,e4,e5,e6])]));
 DD:= DiagonalMatrix(map(e -> exp(t*e),[e1,e2,e3,e4,e5,e6]));
 Einv:= map(simplify,E^(-1)):
 P:= E . DD . Einv; # this is MatrixExponential(Q3,t)

The problem with MatrixExponential in case 3, I think, is that the characteristic polynomial has an irreducible cubic factor, and MatrixExponential insists on using the complicated explicit formula for the roots of that factor rather than a more manageable RootOf.  We can do better "by hand". 

> with(LinearAlgebra):    
   Q3:= eval(Q,parameters3);
   CP:= factor(CharacteristicPolynomial(Q3,s));

> e1:= 0:
alias(e2=RootOf(op(2,CP),s,index=1),
  e3=RootOf(op(2,CP),s,index=2),
  e4=RootOf(op(3,CP),s,index=1),
  e5=RootOf(op(3,CP),s,index=2),
  e6=RootOf(op(3,CP),s,index=3));
 E:=map(simplify,Matrix([seq(op(NullSpace(Q3-e)),e=[e1,e2,e3,e4,e5,e6])]));
 DD:= DiagonalMatrix(map(e -> exp(t*e),[e1,e2,e3,e4,e5,e6]));
 Einv:= map(simplify,E^(-1)):
 P:= E . DD . Einv; # this is MatrixExponential(Q3,t)

The point is that you're not evaluating it at j=0 and j=1, you're simplifying an inert Sum.
By using Sum, you're specifically telling Maple not to just evaluate the summand at j=0 and j=1 and add them.  So it performs what simplifications it can on the summand, oblivious to the fact that these simplifications are invalid when j=0 or j=1.

The "specialization problem" in general is basically the fact that computer algebra systems consider symbolic variables to have "generic" values; the simplifications they perform should be valid "generically", but may run into trouble when you plug in specific values of the variables.  For example, x/x will be simplified to 1, but this would not work for x=0.  For a less trivial example, consider

int(cos(n*x), x=0 .. Pi);

Maple evaluates this as sin(Pi*n)/n.  If you further assume n is an integer, you get 0.  But of course this is not true for n=0. 

The point is that you're not evaluating it at j=0 and j=1, you're simplifying an inert Sum.
By using Sum, you're specifically telling Maple not to just evaluate the summand at j=0 and j=1 and add them.  So it performs what simplifications it can on the summand, oblivious to the fact that these simplifications are invalid when j=0 or j=1.

The "specialization problem" in general is basically the fact that computer algebra systems consider symbolic variables to have "generic" values; the simplifications they perform should be valid "generically", but may run into trouble when you plug in specific values of the variables.  For example, x/x will be simplified to 1, but this would not work for x=0.  For a less trivial example, consider

int(cos(n*x), x=0 .. Pi);

Maple evaluates this as sin(Pi*n)/n.  If you further assume n is an integer, you get 0.  But of course this is not true for n=0. 

@elango8 : if you do this in Worksheet mode rather than Document mode you'll see
the P[1] := ...,, P[2] := ..., etc. 
If you don't want to assign values but just see the results, then you could use = rather than :=.

@elango8 : if you do this in Worksheet mode rather than Document mode you'll see
the P[1] := ...,, P[2] := ..., etc. 
If you don't want to assign values but just see the results, then you could use = rather than :=.