847 Reputation

12 Badges

3 years, 183 days

MaplePrimes Activity

These are answers submitted by mmcdara

Several people have already answered your question in different ways

To complete their propositions, here is a time and memory comparison of different ways to answer yout problem


with(Statistics):   # the most natural (?) and one of the fastest if you want 0 and 1 outputs
B := RandomVariable(Bernoulli(1/2)):

numThrows:= 10^5:

CodeTools:-Usage( Sample(B, numThrows) ):

memory used=0.87MiB, alloc change=0 bytes, cpu time=7.00ms, real time=7.00ms, gc time=0ns



with(StringTools) : # a simple way to get "H" and "T" ouputs

numThrows:= 10^5:
CodeTools:-Usage( Random(numThrows, "HT" ) ):                     # a string of the form "HHT..."
CodeTools:-Usage( Explode(Random(numThrows, "HT" )) ):            # a list of the form [ "H", "H", "T", ...
CodeTools:-Usage( Vector(Explode(Random(numThrows, "HT" ))) ):    # a vector (here fot comparison with tomleslie's solution below)

memory used=154.93KiB, alloc change=0 bytes, cpu time=193.00ms, real time=193.00ms, gc time=0ns
memory used=0.86MiB, alloc change=0 bytes, cpu time=192.00ms, real time=193.00ms, gc time=0ns

memory used=22.27MiB, alloc change=0 bytes, cpu time=383.00ms, real time=384.00ms, gc time=41.38ms



# tomleslie  ... probably the fastest (?) to deliver "H" and "T" ouputs
numThrows:= 10^5:

r:= rand(1..2):

CodeTools:-Usage(  Vector[row](numThrows, i->r()) ):
CodeTools:-Usage(  Vector[row](numThrows, i->["H","T"][r()]) ):

memory used=0.81MiB, alloc change=0 bytes, cpu time=33.00ms, real time=33.00ms, gc time=0ns
memory used=0.77MiB, alloc change=0 bytes, cpu time=38.00ms, real time=38.00ms, gc time=0ns



with(LinearAlgebra): # advantage: there exist a very fast way to tansform 0 & 1 into "H" & "T"

numThrows:= 10^5:

V := CodeTools:-Usage(  RandomVector[row](numThrows,generator=rand(0..1)) ):  # outputs are 0 or 1
W := CodeTools:-Usage(  subs({0="H", 1="T"}, V) ):                            # outputs are "H" or "T"

memory used=4.70MiB, alloc change=2.00MiB, cpu time=50.00ms, real time=51.00ms, gc time=5.63ms
memory used=0.76MiB, alloc change=0 bytes, cpu time=3.00ms, real time=3.00ms, gc time=0ns



with(SignalProcessing):  # comparable to Statistics

numThrows:= 10^5:

V := CodeTools:-Usage( floor~(GenerateUniform(numThrows, 0, 2)) ):  # the fastest (?), use Hfloats (default)

UseHardwareFloats := false:
CodeTools:-Usage( floor~(GenerateUniform(numThrows, 0, 2)) ):  # don't use Hfloats

memory used=1.79MiB, alloc change=0 bytes, cpu time=7.00ms, real time=8.00ms, gc time=0ns

memory used=203.74MiB, alloc change=47.02MiB, cpu time=1.58s, real time=1.59s, gc time=88.87ms



with(combinat):  # just for fun

numThrows:= 10^5:

CodeTools:-Usage( randperm(["T"$numThrows, "H"$numThrows])[1..numThrows] ):

memory used=20.67MiB, alloc change=4.58MiB, cpu time=161.00ms, real time=162.00ms, gc time=61.65ms






After some digressions unrelated to your initial question, here is my definitive contribution



Seed := randomize():


# Number of dice

N := 3;



# These two lines can be replaced by any algorithm amid those presented in this
# same thread (just choose the one you want)

S := add(RV(DiscreteUniform(1, 6)), k=1..N):
P := Probability~(S =~ [$N..6*N])

[1/216, 1/72, 1/36, 5/108, 5/72, 7/72, 25/216, 1/8, 1/8, 25/216, 7/72, 5/72, 5/108, 1/36, 1/72, 1/216]


# Define a discrete RV with outcomes of probabilities given by P
# Watch out: by default (I do not know if and how this could be changed), this RV has outcomes 1..5*N+1
# This raises an issue that will be fixed when the resuklts are ploted (the outcomes are N..6*N)

Dice := RV(ProbabilityTable(P))



# Sample "Dice" as any other RV

M := 100:
X := Sample(Dice, M):

h := Histogram(X, minbins=5*N+1):

# The next operation is aimed to superimpose the exact mss function (aka P) and its
# empirical estimation h.
# I transform h into a list that can be plot with ColumnGraph.
# Watch out:
#      Maple 2018: use op~(2, op~(1, [op(1..-2, op(1, h))] ));
#      Maple 2015: use map(u -> print(u[-1][2]), [op(1..-2, op(1, h))]);


H := map(u -> u[-1][2], [op(1..-2, op(1, h))]):

    [H, P],
    legend=["Sampled", "Exact"],
    color=[red, yellow],
    title=cat("Sample size = ", M),
    tickmarks=[[seq(i=i+N, i=0..5*N)], default],
    size=[600, 400]

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`









A few remarks:

  1. Use Quantile(X, k+i+j, numeric) for a faster generation:
    CodeTools:-Usage(Matrix(100, 10, (i,j) -> Quantile(X, (i-1)/100+(j-1)/1000, numeric))):
    memory used=14.25MiB, alloc change=0 bytes, cpu time=158.00ms, real time=159.00ms, gc time=0ns
    CodeTools:-Usage(Matrix(100, 10, (i,j) -> Quantile(X, (i-1)/100+(j-1)/1000))):
    memory used=430.05MiB, alloc change=0 bytes, cpu time=6.28s, real time=5.90s, gc time=555.98ms

    Note the CPU time (in this particular case) can be halved: Quantile(X, p) = -Quantile(X, 1-p)
  2. Avoid computing Quantile(X, 0) which is obviously -infinity and not -15.681...

Here is a version which uses DocumentTools for a smarter rendering


As you wrote it, your maplet can't work: Evaluate needs a function as argument, not a sequence of instructions.
Here is the correct way to code it in a simple case (no "for a in .... end do" loop ... but you will be able to generalize with from Carl's answer).

PS: the code I wrote conforms to your's, but it could be improved (remark first that there is no need to click on the "Display function" function button to visualize the MathML code in the MMLV field). It's also better to define "Input function  f(x)=" within a Label component, ...



f := proc()
local g, TL2:

g := proc(func, slp)
  local a;
  a := Tangent(func,x=slp):
end proc:

    "Input function  f(x)=", TextField['func'](20)
    Button("Display function", Action(Evaluate('MMLV'='MathML[Export](func)')))
    "Input slope", TextField['slp'](10)
    Button("Tangent line", Action( Evaluate('function'=g, Argument(func), Argument(slp))))


end proc:






You had forgotten that the projected dx*dy surface wasn't the surface itself:




SetCoordinates(cartesian[x, y, z]):

R := sqrt(2):
s := x^2+y^2+z^2-R^2;
u := VectorField([x, 1, z]);

N := Gradient(s);



u := Vector(3, {(1) = x, (2) = 1, (3) = z}, attributes = [vectorfield, coords = cartesian[x, y, z]])


N := Vector(3, {(1) = 2*x, (2) = 2*y, (3) = 2*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])


n := N/sqrt(add(N[k]^2, k = 1 .. 3));
u . n;

un := algsubs(x^2+y^2+z^2=R^2, %):
un := simplify(un)

n := x*`#mover(mi("e"),mo("_"))`[x]/sqrt(x^2+y^2+z^2)+y*`#mover(mi("e"),mo("_"))`[y]/sqrt(x^2+y^2+z^2)+z*`#mover(mi("e"),mo("_"))`[z]/sqrt(x^2+y^2+z^2)






intr1 := solve(subs(z = 0, s), y);

intr2 := solve(subs(z = 0, y = 0, s), x);

(-x^2+2)^(1/2), -(-x^2+2)^(1/2)


2^(1/2), -2^(1/2)


# forgot that the projected dx*dy surface is not the surface irself

l     := sqrt(x^2+y^2);
h     := sqrt(2-x^2-y^2);
phi   := arctan(h/l);
scale := simplify(1/sin(phi)) assuming x^2+y^2 > 0;

# surface of the sphere

R1 := 2*int(int(scale, y = intr1[2] .. intr1[1]), x = intr2[2] .. intr2[1], numeric);














R1 := 2*Student:-MultivariateCalculus:-MultiInt(scale*un, y = intr1[2] .. intr1[1], x = intr2[2] .. intr2[1])






du := Divergence(u)



intr1 := solve(s, z);

intr2 := solve(subs(z = 0, s), y);

intr3 := solve(subs(z = 0, y = 0, s));

R2 := int(int(int(du, z = intr1[2] .. intr1[1]), y = intr2[2] .. intr2[1]), x = intr3[2] .. intr3[1]);


(-x^2-y^2+2)^(1/2), -(-x^2-y^2+2)^(1/2)


(-x^2+2)^(1/2), -(-x^2+2)^(1/2)


2^(1/2), -2^(1/2)









In your particular case it's far more easier to verify the Gass theorem in spherical coordinates

I'll give a closer look to your code later


SetCoordinates(cartesian[x, y, z]):

MySphere := x^2+y^2+z^2=2;
u := VectorField([x, 1, z]);

N := Gradient(lhs(MySphere));

x^2+y^2+z^2 = 2


u := Vector(3, {(1) = x, (2) = 1, (3) = z}, attributes = [vectorfield, coords = cartesian[x, y, z]])


N := Vector(3, {(1) = 2*x, (2) = 2*y, (3) = 2*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])


n := N/sqrt(add(N[k]^2, k = 1 .. 3));
n := subs(MySphere, n);

n := x*`#mover(mi("e"),mo("_"))`[x]/sqrt(x^2+y^2+z^2)+y*`#mover(mi("e"),mo("_"))`[y]/sqrt(x^2+y^2+z^2)+z*`#mover(mi("e"),mo("_"))`[z]/sqrt(x^2+y^2+z^2)


n := Vector(3, {(1) = (1/2)*2^(1/2)*x, (2) = (1/2)*2^(1/2)*y, (3) = (1/2)*2^(1/2)*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])


un := DotProduct(u, n);
du := Divergence(u);





SphereRadius := sqrt(rhs(MySphere));



J := (r, phi) -> r^2*sin(phi);

int(int(int(du*J(r, phi), theta=0..2*Pi), phi=0..Pi), r=0..SphereRadius);

proc (r, phi) options operator, arrow; VectorCalculus:-`*`(r^2, sin(phi)) end proc






int(int(subs(r=SphereRadius, unsph)*J(SphereRadius, phi), theta=0..2*Pi), phi=0..Pi);






Not the solution, just a way to derive the equations by applying differential operators.
I think you have an ambiguity about the pressure (you declare p bar as a constant but it has a priori different values at z=-1 and z=+1 (unless pa is null)

For the solution: I'm not familiar with pdsolve and PDETools, but I'm afraid that you will have to do the things by hand. 
If I'm not mistaken this is a classical poiseuille flow with uniform velocity at the inlet (parabolic profile expected at the outlet) and I remember solving it by hand long ago at school. Maybe you will have to reproduce explicitely the same modus operandi in Maple?

Wait for a more enlightened reply about pdsolve/PDEtools




SetCoordinates(cylindrical[r,theta, z]):

alias(u=u(r, 'theta', z));
alias(v=v(r, 'theta', z));
alias(p=p(r, 'theta', z));



u, v


u, v, p


V := VectorField(< v, 0, u >);

V := Vector(3, {(1) = v(r, theta, z), (2) = 0, (3) = u(r, theta, z)}, attributes = [vectorfield, coords = cylindrical[r, theta, z]])


(v+r*(diff(v, r))+r*(diff(u, z)))/r


Continuity_equation := eval(expand(%));

v/r+diff(v, r)+diff(u, z)


Diffusion_term := beta^2 *~ map(expand, eval(%));

Diffusion_term := Vector(3, {(1) = beta^2*(-v(r, theta, z)/r^2+(diff(v(r, theta, z), r))/r+diff(diff(v(r, theta, z), r), r)+(diff(diff(v(r, theta, z), theta), theta))/r^2+diff(diff(v(r, theta, z), z), z)), (2) = 2*beta^2*(diff(v(r, theta, z), theta))/r^2, (3) = beta^2*(diff(diff(u(r, theta, z), z), z)+(diff(u(r, theta, z), r))/r+diff(diff(u(r, theta, z), r), r)+(diff(diff(u(r, theta, z), theta), theta))/r^2)})


Pressure_term := 1 *~ Gradient(p)

Pressure_term := Vector(3, {(1) = diff(p(r, theta, z), r), (2) = (diff(p(r, theta, z), theta))/r, (3) = diff(p(r, theta, z), z)})


Advection_term := 1 *~ Vector[column](3, i -> DotProduct(V, Gradient~(V[i])))

Advection_term := Vector(3, {(1) = v(r, theta, z)*(diff(v(r, theta, z), r))+u(r, theta, z)*(diff(v(r, theta, z), z)), (2) = 0, (3) = v(r, theta, z)*(diff(u(r, theta, z), r))+u(r, theta, z)*(diff(u(r, theta, z), z))})


C := :-Vector[column](3, i -> Advection_term[i] = - Pressure_term[i] + Diffusion_term[i])[[1, 3]]

C := Vector(2, {(1) = v(r, theta, z)*(diff(v(r, theta, z), r))+u(r, theta, z)*(diff(v(r, theta, z), z)) = -(diff(p(r, theta, z), r))+beta^2*(-v(r, theta, z)/r^2+(diff(v(r, theta, z), r))/r+diff(diff(v(r, theta, z), r), r)+(diff(diff(v(r, theta, z), theta), theta))/r^2+diff(diff(v(r, theta, z), z), z)), (2) = v(r, theta, z)*(diff(u(r, theta, z), r))+u(r, theta, z)*(diff(u(r, theta, z), z)) = -(diff(p(r, theta, z), z))+beta^2*(diff(diff(u(r, theta, z), z), z)+(diff(u(r, theta, z), r))/r+diff(diff(u(r, theta, z), r), r)+(diff(diff(u(r, theta, z), theta), theta))/r^2)})


NS := [Continuity_equation, entries(C, nolist)]

[v/r+diff(v, r)+diff(u, z), v*(diff(v, r))+u*(diff(v, z)) = -(diff(p, r))+beta^2*(-v/r^2+(diff(v, r))/r+diff(diff(v, r), r)+(diff(diff(v, theta), theta))/r^2+diff(diff(v, z), z)), v*(diff(u, r))+u*(diff(u, z)) = -(diff(p, z))+beta^2*(diff(diff(u, z), z)+(diff(u, r))/r+diff(diff(u, r), r)+(diff(diff(u, theta), theta))/r^2)]


# You can say "p bar is a constant and { p(z=-1)=0 & p(z=-1)=0 & p(z=-+1)=p__a } unless p__a=0 too  

BC := [
:-D[1](u)(0, 'theta', z) = 0,
v(0, 'theta', z) = 0,
u(1, 'theta', z) = 0,
v(1, 'theta', z) = p+alpha

seq(print(BC[i]), i=1..numelems(BC))

(D[1](u))(0, theta, z) = 0


v(0, theta, z) = 0


u(1, theta, z) = 0


v(1, theta, z) = p+alpha





Replace your line   Diam;
by                           Diam := simplify(Diam);

Nevertheless I would be tempted to say that Maple should simplify automatically.

Your line evalf(Diam) is corrupted.
Just rewrite it

PS:Personnaly I never use this document mode input that I find very capricious and subjected to hide characters mistakenly typed

I'm not sure I understand you perfectly.
What I think is: you have a discrete RV "X" with mass function fX(x)=1/(x+1) for x = 0, ..., k, where k is some given number (>0).

If It is that, a way to proceed is to use a ProbabilityTable Distribution




# Discrete distribution from a probability table
# The sum of all masses must be equal to 1.

P := k -> [seq(1/(x+1), x=0..k)] /~ add(1/(x+1), x=0..k) ;

# example:

K := 5:
N := 10^3:


X := RandomVariable(ProbabilityTable(P(K))):

S := Sample(X, N):


P := proc (k) options operator, arrow; `~`[:-`/`]([seq(1/(x+1), x = 0 .. k)], ` $`, add(1/(x+1), x = 0 .. k)) end proc


"[[Typesetting:-mfrac(Typesetting:-mn("20",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("10",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("20",mathvariant = "normal"),Typesetting:-mn("147",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("5",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("4",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("10",mathvariant = "normal"),Typesetting:-mn("147",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false")]]"


Vector[row]([3., 2., 6., 2., 1., 1., 2., 1., 2., 1.])







If your only concern is to compare the model and the experiment, and if this experiment is a "not too large" collection of triples (t, x, y) , then the option output=Array[ ...], where ... is the list of the experimental t values, is probably the simpler solution.

Regarding now your "comparison metric".
The answer is very simple: there is unfortunately no good way to proceed!
Computing "the Root Mean Square of all (xexperiment(t)-xmodel(t), yexperiment(t)-ymodel(t))" seems very reasonable.
But it can be a very poor solution for some problems, specifically if x or y are of different natures, vary on very different scales and so on.
Luckily it doesn't seem to be the case here.

Another point you must be aware of: if the (t, x, y) experimental values come from measurements, there are probably entailed by "errors". Suppose that x and y are perfectly measured, but that t is not.
For instance the "true" t is T=1, and the measured value is t=1.1: then comparing the experiment at time t=1.1 to the model at time t=1.1 doesn't make sense.
A (the?) correct approach in this case is to use a bayesian framework where one infer the true value of T given its measure t.

Maybe you could be interested in "fitting" your model on the experiments? For instance to assess the value of A ...?
You also mention that "The problem is however that even if the (x,y) curves are rather close, if may not be a good model".
Are you concerned with aapossible "model error"?
Here again, for these two questions, the correct answer relies to the "bayesian calibration of computer model" issue.

Let me know if you are interested by taht.

To conclude: try your  "Root Mean Square of all (xexperiment(t)-xmodel(t), yexperiment(t)-ymodel(t)", examine carefully the residues (the model-to-experiment distance at each time t) and if all looks good, then it's done 

MAT := < header1, header2, <Column(t1[2,1],1..2) | Column(t2[2,1],2) > >:



Assuming the DEFINITION below suits you, you will find some elements in the attached file.


  • Let E the ellipse defined by the equation (x/a)^2+(y/b)^2=1
  • For each point P on E one define s(P) as the measure of the angle (counter clockwise counted) between Ox and OM
  • Let A and B two points on E such that
    • s(A) > =  0
    • s(B) > s(A)
    • s(B) <= 2*Pi
  • Let M any point such that s(A) <= s(M) <= s(B)
  • Let x(A) the abscissa of A and x(B) the one of B

Then the mean value of the distance d(O, M) from O=(0, 0) to M as M moves from A to B is defined by the integral of d(O, M) between x(A) and x(B), divided by x(B)-x(A).

In file you will find a procedure 'md' which computes the mean value of d(O, M). It takes 4 arguments:

  1. b
  2. alpha = s(A)
  3. beta = s(B)

'md' handles (if I'm not mistaken) the 3 cases:

  1. alpha < Pi and beta < Pi
  2. alpha < Pi and beta > Pi
  3. alpha > Pi and beta > Pi 

Before that you will find a few lines that treat more formally the first situation alpha < Pi and beta < Pi

As claimed in this file, I was not capable to obtain a general formal solution for any values of a and b (that's why I arbitrarily wrote a:=2 and b:=1 at the top at the file).
The search of a general formal solution requires using assumptions about a, b, p, q (see the commands in the file) and the 'int' procedure of Maple is not very comfortable with 'assume' or 'assuming' statements (someone here will probably correct me on this point)

Here is a possible (NOT THE) definition of the mean radius.

I replace "arithmetic mean" by mean in the sense (one sense) mentioned by Rouben:

  • I consider the canonical ellipse with principal axis Ox and Oy (natural cartesian coordinates)
  • let M1, M2<> M1 two points  on the ellipse
  • Let M a third point on the ellipse.
     Assuming M moves from M1 to M2 counterclockwise in such a way that the angular sector swept by M is less than 2*Pi,
    one can define the mean radius by the mean value of OM when M moves from M1 to M2.



Maybe the shortest solution?
(just for fun, no practical interest when compared to the previous answers)


splitChange := proc(w)
  local bw:
  bw := map(u -> if w[u+1]<>w[u] then cat(w[u], " ") else w[u] end if, [$1..length(w)-1]);
  StringTools:-Words(cat(op(bw), w[-1]));
end proc:


I guess you took as a reference the article by Faith Chelimo Kosgei?
(  practically the one I was able to find on the subject)

If it so I understand you asked for some help (I 've rarely read such a muddled, obscure and confusing paper!)

Nevertheless, after a lot of work spent to read between the lines of the article, you will find in the attached file the algorithm applied to the 3 toy examples this paper contains.

I leave you the pleasure to apply the method to your own test cases

4 5 6 7 8 9 10 Page 6 of 11