cancel
Showing results for
Did you mean:
cancel
Showing results for
Did you mean:

Community Tip - Learn all about PTC Community Badges. Engage with PTC and see how many you can earn! X

## How to generate a semipositive indecomposible matrix of nxn  8-Gravel

## How to generate a semipositive indecomposible matrix of nxn

I am still using MC 11.

I was able to make a random square matrix called A of n dimension, with 1<n<1000 (infinite would have been better but my pc horse power is limited). A needed to be positive, with Max Eigenvalue  lower than 1.

I used the following combination of formulas, which probably is not the optimum, but it works.

1): Z(z,z)=:rnd(1/n)    # so to have numbers that fulfill the condition MaxEval<1

2) A=: matrix(n,n,Z)

The constrains that I would like to add  in generating A are the following.

a) Still a nxn square Matrix with Max Eigenvalue lower than 1.

b) But A should be semi-positive, with random coefficients, including zeros and not just positive as above.

c) However A should be indecomposable (called also irreducible), i.e. A^x  (with x a positive number) should turn a positive matrix.

Thanks for any suggestion.

1 ACCEPTED SOLUTION

Accepted Solutions  23-Emerald III
(To:anthonyQueen)

With the rnd() random number generator you get random numbers from a uniform sitribution. The chance of hitting a specific number (such as 0) is virtually zero; exactly it is 1 over the number of possible random values, since Mathcad's range of real number is pretty large, you chance of hitting 0 therefor is pretty low: approximatly 0. If you want a specific probability of occurrence of any number, you have to restrict the number of possible random values. I admit that I restricted it fiercely. You may want to use:

Z(z,z)=: floor(rnd(n^2))/n^3

or any higher power

I considered another option to include 0, which is

Z(z,z)=: ((rnd(2*n)-n)^2)/n^2

This brings 0 to the middle of the range, but does not increase the probability of hitting it exactly. On first run with n=4, I got 3.33e-6 as smallest value.

Yet another possibility is to additionally declare all values within a certain distance from 0 as 0 itself:

Z(z,z)=:

|r<-rnd(n)/n^2

|return if(r<0.0001,0,r)

But this ruins your uniform distribution.

Another construct could be:

Z(z,z)=: if(rnd(n)<0.0001,0,rnd(r))/n^2

Which could result in non-zero numbers below 0.0001/n^2, because if on first call to rnd() the result is >= 0.0001, a new random number is generated by the second call to rnd(), but I don't exactly know what it does to your (uniform) distribution.

Success!
Luc

4 REPLIES 4  23-Emerald III
(To:anthonyQueen)
Use Z (z,z)=:floor (rnd (n))/n^2.
That will cause a 0 every now and then.

Success!
Luc  8-Gravel
(To:LucMeekes)

Many Thanks Luc!

Seems simple and working. Yet I noticed that it generates only n numbers, zeros included. So for a matrix 3x3 there are only 3 numbers generated, which are put in the 9 cells.  And this may be a problem, since each cell does not have a random semipositive number (I should probably mentioned that in my constraints). I fear that  matrix A could become singular, which is not what I want.  23-Emerald III
(To:anthonyQueen)

With the rnd() random number generator you get random numbers from a uniform sitribution. The chance of hitting a specific number (such as 0) is virtually zero; exactly it is 1 over the number of possible random values, since Mathcad's range of real number is pretty large, you chance of hitting 0 therefor is pretty low: approximatly 0. If you want a specific probability of occurrence of any number, you have to restrict the number of possible random values. I admit that I restricted it fiercely. You may want to use:

Z(z,z)=: floor(rnd(n^2))/n^3

or any higher power

I considered another option to include 0, which is

Z(z,z)=: ((rnd(2*n)-n)^2)/n^2

This brings 0 to the middle of the range, but does not increase the probability of hitting it exactly. On first run with n=4, I got 3.33e-6 as smallest value.

Yet another possibility is to additionally declare all values within a certain distance from 0 as 0 itself:

Z(z,z)=:

|r<-rnd(n)/n^2

|return if(r<0.0001,0,r)

But this ruins your uniform distribution.

Another construct could be:

Z(z,z)=: if(rnd(n)<0.0001,0,rnd(r))/n^2

Which could result in non-zero numbers below 0.0001/n^2, because if on first call to rnd() the result is >= 0.0001, a new random number is generated by the second call to rnd(), but I don't exactly know what it does to your (uniform) distribution.

Success!
Luc  8-Gravel
(To:LucMeekes)
Thanks Luc again. I certainly try them too!! 