Carl Love

Carl Love

28110 Reputation

25 Badges

13 years, 120 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

The local command can be used at the top level. So, all you need to do is

local gamma;

Your transition matrix P is a patterned tridiagonal or band matrix. It can be constructed with a single-line Matrix command using the scan= band option.

I can see no practical value to storing the positive integer powers of this matrix, so I haven't done that. You might as well simply recompute the powers as needed.
 

restart:

Digits:= 15:

N:= 2:  dimP:= 2*N+1:  m:= dimP-2:

(q,p):= (0.3,0.5):  r:= 1-p-q:  sa:= 0.9:  sb:= 1-sa:

P:= Matrix([[q$m, sa], [sa, r$m, sb], [sb, p$m]], scan= band[1,1]);

Matrix(5, 5, {(1, 1) = .9, (1, 2) = .1, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = .3, (2, 2) = .2, (2, 3) = .5, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = .3, (3, 3) = .2, (3, 4) = .5, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = .3, (4, 4) = .2, (4, 5) = .5, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = .9, (5, 5) = .1})

The matrix is stochastic with no entry equal to 1. The theory of Markov chains says that there is a stochastic vector s such that s.P = s. All rows of P^n approach this s as n approaches infinity. Here are two methods to compute s.

 

Method 1: Eigenvectors:

The equation s.P = s means that s is a left eigenvector of P for eigenvalue 1. So the transpose of s is a right eigenvector of the transpose of P. The transpose of a matrix or vector is denoted in Maple by appending ^+.

(e,V):= LinearAlgebra:-Eigenvectors(P^+):

Find the position of 1 among the eigenvalues and call the position p1.

member(1,e,p1);

true

Extract the corresponding eigenvector and transpose it.

s:= V[..,p1]^+;

Vector[row](5, {(1) = -.627245607666573+0.*I, (2) = -.20908186922219074+0.*I, (3) = -.34846978203698503+0.*I, (4) = -.5807829700616427+0.*I, (5) = -.32265720558980143+0.*I})

The zero imaginary parts can be removed by simplify:

s:= simplify(s);

Vector[row](5, {(1) = -.627245607666573, (2) = -.209081869222191, (3) = -.348469782036985, (4) = -.580782970061643, (5) = -.322657205589801})

Convert to a stochastic vector by normalizing with respect to the 1-norm:

s:= s*csgn(s[1])/LinearAlgebra:-Norm(s,1);

Vector[row](5, {(1) = .30037082818294236, (2) = .10012360939431413, (3) = .16687268232385685, (4) = .2781211372064288, (5) = .15451174289246009})

Compare with a high power of P:

P^999;

Matrix(5, 5, {(1, 1) = .30037082818294697, (1, 2) = .10012360939431567, (1, 3) = .16687268232385943, (1, 4) = .27812113720643244, (1, 5) = .15451174289246247, (2, 1) = .30037082818294686, (2, 2) = .10012360939431565, (2, 3) = .1668726823238594, (2, 4) = .27812113720643233, (2, 5) = .15451174289246242, (3, 1) = .30037082818294686, (3, 2) = .10012360939431565, (3, 3) = .16687268232385938, (3, 4) = .27812113720643233, (3, 5) = .15451174289246242, (4, 1) = .30037082818294686, (4, 2) = .10012360939431565, (4, 3) = .1668726823238594, (4, 4) = .27812113720643233, (4, 5) = .15451174289246242, (5, 1) = .30037082818294686, (5, 2) = .10012360939431564, (5, 3) = .1668726823238594, (5, 4) = .27812113720643233, (5, 5) = .15451174289246242})

Method 2: System of linear equations:

Solving the system s.P = s for s is equivalent to solving (P^+ - I).s^+ = 0 (where I is the identity matrix). This system is of course singular. We fix that by adding the condition that s be a stochastic vector.

LinearAlgebra:-LinearSolve(
  <P^+-Matrix(dimP, shape= identity), <seq(1, 1..dimP)>^+>,
  <seq(0, 1..dimP), 1>
)^+;

Vector[row](5, {(1) = .30037082818294225, (2) = .10012360939431397, (3) = .16687268232385655, (4) = .2781211372064275, (5) = .15451174289245972})

 


 

Download TransitionMatrix.mw

I think that you have enough sophistication and experience with a variety of languages[*1] to realize that any statement of the form "ASCII string X and ASCII string Y mean exactly the same thing in computer language L, which is abstract mathematical object M" is ontologically problematic. The level of evaluation (parser, automatic simplifier, etc.) needs to be considered. For your example, lprint (among numerous possible evaluations) reveals that they're not the same thing at the default level of evaluation. 

The only form of derivative for initial and boundary conditions containing derivatives that's consistently supported by the documentation is D. If some component accepts another form, that's merely grace (to use Preben's word) or a bonus (to use Vv's word). So it's only at the highest level of evaluation---inside the mind of the programmer---that they could be considered equal.

[*1]Judging by your 752 postings on MaplePrimes, every word of which I've read.

If LL is your list of lists, then the list of the averages of the sublists can computed by (add/nops)~(LL). This works for both symbolic and numeric data. 

You can't meaningfully assign values to both a name and an indexed form of the same name. In this case, you've done this with and r[1] and and k[1]. You can avoid this problem while still having the prettyprinted form of the index appear as a subscript by using r__1. If you don't care about whether the index appears as a subscript, you can use and r1.

There is no problem with assigning to multiple indexed forms of the same main name as long as you don't assign to the main stem. For example, you can assign to r[0] and r[1] as long as you don't assign to r.

The command that you're trying to use, ExpandSteps, is only for polynomial expansion, not for matrix arithmetic. The command that applies to your situation is Student:-LinearAlgebra:-GaussJordanEliminationTutor; however, its use is limited to matrices at most 5 x 5.

I see that your matrix has 9-digit exact integers and exact imaginary values. So what's the point of seeing the steps? Do you suspect an error in Maple's computation?

Maple handles regular expressions. See these help pages:

  • ?StringTools,Regular_Expressions
  • ?StringTools,RegMatch
  • ?StringTools,RegSplit
  • ?StringTools,RegSub
  • ?StringTools,RegSubs

All that VV said is true. I'd just like to add that if you want symbolic variables to be treated like would-be integers until they have values assigned, then use irem instead of mod.

@tomleslie Thanks for testing. The results are not surprising given that I used an array half the size of yours but a more complicated indexing scheme, and otherwise our codes are remarkably similar. Here's an update---much terser but not significantly faster---that probably satisfies the OP's titular no-loop request (thus I made this an Answer):

EratosthenesSieve:= proc(N::posint)
description "Returns list of primes <= N";
local p, k, P:= rtable(1..iquo(N-1,2), k-> 2*k+1);
   if N < 5 then return remove(`>`, [2, 3], N) fi;
   P[[seq(seq(iquo(k-1, 2), k= p^2..N, 2*p), p= thisproc(isqrt(N))[2..])]]:= ();
   [2, seq(P)]
end proc:

 

It is useful---and often essential---to avoid repetition when generating combinatorial objects. It can, for example, reduce the time complexity from O(n^2) to O(n), which can be the difference between doability and nondoability.

Here is a generalization of your problem: Let B be the base set from which the integers are drawn, e the exponent, and n the number of integers to draw (with replacement, but regardless of order). Then your problem can be done by:

B:= [$0..5]:  e:= 2:  n:= 2:
[seq(add(B[[seq(C)]]^~e), C= Iterator:-MultiCombination([n$nops(B)], n))];

Note that I used no sets, only lists, yet there are no repetitions. Thus, each desired object was generated only once. The Iterator package is generally very good about that.

Just change the last line to

Maximize(simplify(S), {a^2+b^2+c^2 = 1}) assuming real;

This is needed to get rid of the absolute-value / complex-modulus operators, which are messing up the differentiability of the objective function.

I can confirm using your worksheet that I get the same syntax error as you do. I don't know why. But you can do this:

q1:= (n1, n2)-> `if`(n1 < 0 or n2 < 0 or abs(n1-n2) <> 1, 0, sqrt(max(n1,n2)/2));

That gives me no syntax error.

The Maple equivalent of the Matlab syntax Q(:, 1), where Q is a Matrix, is Q[.., 1].

To answer your overall question---to wit:

  • I would like to ask if there is any command in Maple while solving the equations that are responsible for possible approximations when solving equations

If you use the command solve or any command in its "family", which include the package SolveTools and numerous other commands and packages, and your input does not include any decimal numbers, then there are no approximations used. If the command cannot find an exact solution, then it returns nothing.

What makes you think that the solution that you got is wrong?

In the worksheet attached to your original Question, you give the following formula, which you attribute to VV:

  • The total number of set partitions of n distinct items into e sets of length L, with n = e*L, is n!/e!/(L!)^e.

To be more precise, that is the number of ways of doing it if the subsets are unlabeled. If they are labeled (as in red team, green team, orange team), then the formula is n!/L!^e. Both of these formulas can be easily generalized to the situation where the subsets aren't necessarily equally sized.

Let be a list of positive integers representing the sizes of the blocks of the partition. If the blocks are labeled, then the number of partitions is add(B)!/mul(B!~). This can also be computed in Maple as combinat:-multinomial(add(B), B[]). If the blocks are unlabeled, then the number of partitions is add(B)!/mul(B!~)/mul(rhs~(Statistics:-Tally(B))!~). This can also be computed in Maple as Number(Iterator:-SetPartitionFixedSize(B)).

For B = [5,5,5], the labeled-block formula gives 756756 and the unlabeled 126126. The former leads to the probability 1/1001 and the latter 6/1001. The denominator from my original Answer is binomial(15,5)*binomial(10,5)/3!, which also equals 126126.

 

First 143 144 145 146 147 148 149 Last Page 145 of 395