mmcdara

7891 Reputation

22 Badges

9 years, 63 days

MaplePrimes Activity


These are replies submitted by mmcdara

@Preben Alsholm 

Thank you Preben, I'm gonna thest this.

@jalal  @Kitonum

Kitonum's methods will surely return different numbers but if you saying "generate different numbers" means that "all places in the vector are to be occupied by random numbers drawn from independent identically distributed random variables with discrete uniform distribution over 10, ..., 20, this is not a good solution.

A simple proof is that if the vector has length 11, then the 11th number to place is not at all random.
More of this you will not be capable to generate a random vector whose length is larger than the number of "random" numbers (11 in the case 10..20)

So,  @jalal if you want to generate random numbers you must accept some numbers to appear more than once!
Your question "How to have only different integers in the first column ?" is not consistent with the "random" requirement.

PS : you must have been very unlucky to obtain third times the same number on the first column, for the probability of such an event is 20 / 161051 ~  0.000124...

@Preben Alsholm 

No harm, I'm not even sure it is possible.
I had put some hope on the  "pre-variables" (for instance using pre(v(t)) or something like that; results are in the attached file) but it seems to be a dead end.

I keep digging
 

restart:

with(plots):

# x(t)        body displacement
# v(t)        body velocity
# M__p     body mass
# V__p     body volume
# P__0      pre-stress of the spring
# K           spring stiffness
# C           solid-solid friction coefficient
# rho__h   density of the fluid the body is immersed in
#
# gamma__longi(t) longitudinal acceleration

gamma__longi := t -> A*t;

proc (t) options operator, arrow; A*t end proc

(1)

 

Preben's proposal (yellow highlighted)

equations := {
               diff(x(t), t) = v(t),
               M__p * diff(v(t), t) = (M__p - V__p*rho__h) * gamma__longi(t)
                                      - (P__0 + K*x(t))
                                      - (b(t)*(-C)+(1-b(t))*(C))
            };


# take-off time
# (obtained by equating to 0 the rhs of the 2nd edo after having set x(t)=v(t)=0 and removed the friction term)


rest_state := eval(equations[1], [x(t)=0, v(t)=0]);
T__off     := solve(eval(rest_state, C=0), t);


# initial conditions

IC := { x(T__off) = 0, v(T__off) = 0, b(T__off) = 0 };


# full system

sys := equations union IC:


# solution

P_SOL := dsolve(
                    sys,
                    numeric,
                    parameters=[A, M__p, V__p, rho__h, P__0, K, C],
                    discrete_variables=[b(t)::boolean],
               #    events=[[[v(t)=0, And(fallingedge(v(t)))],b(t)=1-b(t)]]
               #    events=[[[v(t)=0, And(risingedge(v(t)))],b(t)=1-b(t)]]
                    events=[[[v(t)=0, And(change(v(t)))],b(t)=1-b(t)]]
                ):

{M__p*(diff(v(t), t)) = (-V__p*rho__h+M__p)*A*t-P__0-K*x(t)+b(t)*C-(1-b(t))*C, diff(x(t), t) = v(t)}

 

0 = (-V__p*rho__h+M__p)*A*t-P__0+b(t)*C-(1-b(t))*C

 

P__0/((-V__p*rho__h+M__p)*A)

 

{b(P__0/((-V__p*rho__h+M__p)*A)) = 0, v(P__0/((-V__p*rho__h+M__p)*A)) = 0, x(P__0/((-V__p*rho__h+M__p)*A)) = 0}

(2)

P_DATA := [
            A      = 100,
            M__p   = 10,
            V__p   = 1,
            rho__h = 1,
            P__0   = 1000,
            K      = 10000,
            C      = 0
          ]:

T__0 := eval(T__off, P_DATA);


P_SOL(parameters=P_DATA):
P_SOL(10):
Frictionless_P_SOL := odeplot(P_SOL, [t, x(t)], t=T__0..2, color=blue):



P_DATA := [
            A      = 100,
            M__p   = 10,
            V__p   = 1,
            rho__h = 1,
            P__0   = 1000,
            K      = 10000,
            C      = 100
          ]:

T__0 := eval(T__off, P_DATA);


P_SOL(parameters=P_DATA):
P_SOL(10):
StrongFriction_P_SOL := odeplot(P_SOL, [t, x(t)], t=T__0..2, color=red):

display(Frictionless_P_SOL, StrongFriction_P_SOL);

10/9

 

10/9

 

 

 


 

Download bad_event_2.mw

@acer 

I discovered this specifity of Maine 2 days ago when following the elections on the web (I'm not American).
So my answer is simply: I don't know if this splitting matters.

A few days ago my daughter asked me how many electors there were.  
I answered  in a professorial manner "of course 2*269+1"... until my certitude felt down as I heard that it could be a 269-269 tie. Looking to Wikipedia I found the complete electoral college had 538 members and that a tie wasn't impossible, just unlikely ("Still, an overall 269-269 tie remains decidedly unlikely, but it is possiblein-the-presidential-race-what-happens-in-an-electoral-college-tie)

So i asked to my myself how unlikely this situation could be and if what could be "extreme" configurations that could lead to such ties. 9% is not really an unlikely situation from a statistical point of view; the fact is that it seems unlikely because it never (?) happened before.
I guess a lot of people in the US must have studied this situation already?
At the very end the tranger thing is the complexity of the "loophole" the US constitution has put in place (given that an electoral college with an odd number of members seems so natural a solution).

A related question could be "Given the historical color of a state, what is the probability for a tie to occur?"
To show how abstract my simulations are, there are about 20% of chances for a tie to happen if California and Texas share the same color. It is still on this date an unrealistic situation.
If we discard it the probability of a tie falls down to about 7%.

states := map(convert, states, set):
t := 0:
for s in states do
  if numelems(`intersect`(s, {California, Texas})) = 2 then t:=t+1 end if:
  if numelems(`intersect`(s, {California, Texas})) = {} then t:=t+1 end if:
end do:
evalf(t/nbties);
                          0.2032459955





 

from compile help page:

There are additional restrictions imposed by the numeric code compiler on top of those imposed by the CodeGeneration package.

  •   No nested procedures can be translated.
  •   No modules can be translated.
  •   Exception handling (other than merely raising an exception) cannot be translated.
  •   Procedures that return arrays or allocate new arrays cannot be translated.
  •   Procedures that call other procedures about which insufficient type information can be obtained at compile time cannot be translated.
  •   Integer data are limited to word-sized integers on the current hardware. Computations that produce larger integers will raise an overflow exception.
  •   Floating-point computations, including complex floating-point operations, are carried out entirely in hardware precision (like evalhf).
  • Since all supported hardware uses radix 2 floating-point arithmetic, while Maple's software floating-point system is a radix 10 system, computations that depend crucially on the radix will produce different results.

     

Having already been concerned by this problem I'm very interested by an authorized answer.
For instance is there a way to bypass the limitation concerning arrays, matrices, vectors, lists, ... ?

 

@Preben Alsholm 

if it's not too much asking I would like you to help me to improve your model base on these two elements:

friction := - (b(t)*f+(1-b(t))*(-f))
events=[[v(t)=0,b(t)=1-b(t)]]

This model works perfectly well on the problem I submitted, but fails as soon as I add a spring on the system.

The full model is given in the attached file:

  • first part is the solution obtained with your proposal;
  • the second is the one a regularizad Coulomb's model gives;
  • the last is the explanation of the wrong (I do not want to be rude) solution your proposal gives (I do not pretend the regularized solution is good but it is less disturbing than the first one).

I have identified the problem with "your" model. It comes from the expression of the event wich flip from 0 (1) to 1 (0) the discrete variable b(t) each time v(t) = 0.
In fact this sould be done each time v(t) changes its sign.
I have tried, without any success, to express this in terms of event.

Could you help me to fix this?

TIA

bad_event.mw

Years ago I remember my having asked if it was possible to export the output of DocumenTools[Tabulate] just the way it can be done for a plot.
If I don't mistake the answer was no (I searched for the Q/A but I failed).

When I face this issue I just do a screen capture of the region of interest and convert the result in a jpg, pdf or other format.

@Rouben Rostamian  

Good idea indeed!
I've just checked out your idea and it is an elegant one.
Unfortunately it remains quite slow (nevertheless faster than what I did) and I'm not sure at all this can be dome faster.

Please have a look to the two last plots in the updated file and the different computational times I got.
odeplot_differences.mw

Meanwhile Preben has proposed a way to solve the true problem (the one without any regularization) that cuts short all the questions I was asking myself.
So I think he probably gave me a definitive solution.
In any cas thnks for your involvement
 

@Preben Alsholm 

The Prior is: this is a very good idea
Actually I had already seen the possibility to use boolean variables in the help page dedicated to events but it is so abstruse and mysterious  (IMA) that I always spend a lot of time to understand how to use them (without even being sure to get the expected result :-) ):

And the Posterior is: this is a very excellent idea to obtain the exact (even if numerically approximated) solution.
It is "infinitely" faster than the method based on the regularization by the hyperbolic tangent and doesn't require adjusting maxfun.
At the very end it avoids seeking for a value of c which would realize a balance between assuracy and low computational time.
Really a very useful way to tackle the problem.

 

@Carl Love  @salma951

Maybe the interest of the exercise is to froce the student to use many different things: sampling (and maple posseses a lot of things to do that), conditional structure, interactive input ("the zero case" could be tackled using Maplets for instance), random picking of points from a matrix (not si immediate, personally I would use combinat:-randperm) .....
Probably a more structured exercise invoking the same themes would have been more attractive.

The "0 case" seems indeed strange ; if the numbers are floats ine the range -100..100 the probability to get 0 is exactly null. If "in the range from -100 to 100" means "all the integers from -100 do 100 included end equally likely" this probability becomes (1/201)^2.
Thus so small a value that the the student will probably never face this "0 case".
I would say here that the problem is

  1. not completely specified (integers or floats, what does random mean?),
    more of this the phrase "If Maple returns negative figures for r and h" is strange, it should be 
    "If Maple returns negative figures for r OR h"
  2. not cleverly set because the small probability of occurrence  of the "0 case" will not allow the student to verify of if the way it handles it is correct or not


 

For you Salma there are two pieces of code you could take inspiration from
This one is for interactively choosing a value of r (it's not an snwer to your point 1, just a helping hand so that you can go further (for your problem change -2..2 to -100..100)

with(Maplets):
with(Maplets[Elements]):

maplet := Maplet([
   ["Enter the value of r: ", TextField['IB1'](10)],
   [Button("OK", Shutdown(['IB1'])), Button("Cancel", Shutdown())]
]):

roll := rand(-2..2);
if roll() = 0 then 
  result := parse(Display(maplet)[]);
end if;

This one is about the random selection of r and h from the matrix.
I have supposed here that you have to pick two different numbers, thus [1..2] in the second line; if you are allowed to pick twice the same number, try to imagine how you can do this.
Two hints: do twice the selection of a single number by using randperm, or use rand(...) or Sample with an ad hoc distribution.

M := Matrix(3,3,[-65, 1, 32, 1, -12, 3, 12, -2, 2]):
r, h := combinat:-randperm([entries(M, nolist)])[1..2][];
                             1, 32


 

Maybe this link could interest you  An example with orientable and non-orientable...

@Carl Love 

Hi, 
I totally agree.
But, as a friend of mine used to say "even in the worst movie there is always a handfull of seconds which save it from the catastrophe".
So maybe there is a nugget somewhere?
In fact I think his exercice has the interest to demonstrate the power of Maple for this latter can found the PDF of VolCone in a closed form (without no help).
More of this one can think to this exercice as an example about how to create a new distribution and sample it?


 

restart:

with(Statistics):

r := RandomVariable(Uniform(-100, 100)):
h := RandomVariable(Uniform(-100, 100)):

V := 1/3*Pi*abs(r)^2*abs(h)

(1/3)*Pi*abs(_R)^2*abs(_R0)

(1)

pdf := PDF(V, t)

pdf := 3*piecewise(3*t/Pi < 0, 0, 3*t/Pi <= 1000000, -(1/3000000)*(-1000+sqrt(3)*sqrt(t/Pi))*sqrt(3)/sqrt(t/Pi), 1000000 < 3*t/Pi, 0)/Pi

(2)

s := Support(V)

RealRange(0, Open((1000000/3)*Pi))

(3)

f := unapply(pdf, t):

U := Distribution(PDF = (t -> f(t)), Support=op(1, s)..op(op(2, s)));

Vol := RandomVariable(U)

_m4514979488

 

_R1

(4)

Histogram(Sample(Vol, 10^4, method=envelope))

 

 


 

Download aesthetic.mw

@Momo 

The transpose Q^+ of Q is stochastic matrix.
Then its stationnary state of Q^+ is a 2x2 matrix whose all the rows are identical.
Each of these rows is the left eigenvector of Q^+ associated to the eigenvalue 1.

 

restart:

Q := Matrix(2, 2, [0.9, 0.2, 0.1, 0.8]);

# define a row vector of length 2
V    := Vector[row](2, symbol=v);

# here are its entries
vars := [entries(V, nolist)];

# the left eigenvector of Q^+ associated to the eigenvalue 1 verifies
# Vp . Q^+ - Vp = 0.

Vp   := eval(V, solve([entries(V.Q^+ - V, nolist)], [entries(V, nolist)])[1] )^+;


# If you are looking for the stationnary matrix you must "normalize" Vp in 
# order that Vp[1] + Vp[2] = 1

Vp   := eval(V, solve([entries(V.Q^+-V, nolist), add(vars)=1], [entries(V, nolist)])[1] )^+;

# the stationary state of Q then is
Q__Stationnary := < Vp | Vp >;  # = [[2/3, 2/3], [1/3, 1/3]]

 

@acer 
Thank you very much acer!

Im really impressed by all the "simple" tricks you have used in f4 and f4a (inverting p and n, evalhf, inplace, the way you write "add" in place of "my" Heaviside/add operations).
At the end I'm surprised that your "shazam" procedure, a "add" with an explicit loop, be more efficient than the native "add" once compiled. This is the kind of thing I would never have thought of.

You wrote "If you have to repeat the whole computation several times, then the "memory used" number can have more effect..", I agree with that: it is often an issue for me, probably even harder to manage that the one of time efficiency (on the one hand, you can afford to wait for a solution, but you can't cope with the saturation of your PC's memory).

I'm going right now to apply all this stuff to a slightly more complex situation.

If it's not too much to ask I'd like to ask you question related to my initial one.
How would you do this effectively (sampling L Binomial(N, P) each with a different value of P)

# a in [0, 1) and  b in [a, 1] are known
N := 10:   # for instance
L := 1000: # for instance
p := Sample(Uniform(a, b), L)
k := Vector(L, i -> Sample(Binomial(N, p[i]), 1)[1])

Here is the R's code and a sketch of its performance

f <- function(L, N) { 
  +     p <- runif(L, 0, 1) 
  +     rbinom(L, N, p)      # p being of length L, R understands "one single sample per value of p"
  + }
system.time(f(10^6, 10^4))
utilisateur     système      écoulé 
0.140       0.006       0.145 

Feel free to tell me if this requires too much work from your side and if you prefer that I open a specific new thread.

Thanks again

Already an easy-to-find  improvement 

pi0 := (n, p) -> (1-p)^n;
pi1 := (n, p) -> n*(1-p)^(n-1) * p;

f3 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi0~(N, Q): 
  P1    := pi1~(N, Q): 
  P01   := P0 + P1: 
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := Heaviside~(P01-U):
  Got0  := Heaviside~(P0-U):
  N01   := add(Got01):
  N0    := add(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f3(10^5, 0, 5e-4, 25000));

memory used=38.20MiB, alloc change=0 bytes, cpu time=380.00ms, real time=309.00ms, gc time=132.12ms

 



PS : Obviously the values of N0 and N1 can be obtained analytically (wich would be the fastest way to answer the issue), but for many reasons this is not what I am interested in.

restart:

pi0 := (n, p) -> (1-p)^n;
pi1 := (n, p) -> n*(1-p)^(n-1) * p;


pi  := (n, p, k) -> binomial(n, k) * (1-p)^(n-k) * p^k;

proc (n, p) options operator, arrow; (1-p)^n end proc

 

proc (n, p) options operator, arrow; n*(1-p)^(n-1)*p end proc

 

proc (n, p, k) options operator, arrow; binomial(n, k)*(1-p)^(n-k)*p^k end proc

(1)

assumptions := n > 0, p > 0, p < 1, a >=0, b > a, b <= 1;

ek0 := ( int(pi0(n, p), p=a..b) / (b-a) ) * R assuming assumptions;
ek1 := ( int(pi1(n, p), p=a..b) / (b-a) ) * R assuming assumptions;

0 < n, 0 < p, p < 1, 0 <= a, a < b, b <= 1

 

((1-b)^n*b+(1-a)^(n+1)-(1-b)^n)*R/((n+1)*(b-a))

 

((1-a)^n*a*n-(1-b)^n*b*n+(1-a)^n-(1-b)^n)*R/((n+1)*(b-a))

(2)

data := [a=0, b=5e-4, n=25000, R=10^6];
eval([ek0, ek1], data)

[a = 0, b = 0.5e-3, n = 25000, R = 1000000]

 

[79996.50310, 79992.78806]

(3)

 


 

Download MakeItFaster_Exact.mw

First 99 100 101 102 103 104 105 Last Page 101 of 154