mmcdara

620 Reputation

2 Badges

3 years, 14 days

MaplePrimes Activity


These are answers submitted by mmcdara

I do not think you will be able to obtain an analytic expression.
Here is an approximate one:


 

restart:

z := sqrt(G*(2-G))+(1-G)*(arccos(G-1)+(1/5)*Pi)+k*sqrt(G*(2-G))*cos(sqrt((1+k)/k)*arccos(G-1)+(1/5)*Pi)+sqrt(k/(1+k))*(1-G)*k*sin(sqrt((1+k)/k)*arccos(G-1)+(1/5)*Pi)

(G*(2-G))^(1/2)+(1-G)*(arccos(G-1)+(1/5)*Pi)+k*(G*(2-G))^(1/2)*cos(((1+k)/k)^(1/2)*arccos(G-1)+(1/5)*Pi)+(k/(1+k))^(1/2)*(1-G)*k*sin(((1+k)/k)^(1/2)*arccos(G-1)+(1/5)*Pi)

(1)

p := plots:-contourplot(z, G=-1..2, k=1..10, contours=[0], grid=[400, 400]):

op(-1, p);
c := op(1..-2, op(1, p)):
M := convert(op~(1, [c]), Matrix);

"AXESLABELS[Typesetting:-mi("G",italic = "true",mathvariant = "italic"), Typesetting:-mi("k",italic = "true",mathvariant = "italic")]"

 

Matrix(%id = 18446744078394417014)

(2)

plots:-display(p, gridlines=true)

 

# here is the approximate relation k=f(G) such that z(k=f(G), G) = 0)

fit := Statistics:-LinearFit([1, G, G^2, G^3], M, [G]);

plots:-display(p, plot(fit, G=min(M[..,1])..max(M[..,1]), color=blue, thickness=7, transparency=0.8, gridlines=true))

HFloat(5.167125933526261)-HFloat(6.356106757782967)*G+HFloat(3.0043651079986984)*G^2-HFloat(0.767811803757935)*G^3

 

 

zg := eval(z, k=fit):  # = z(k=fit(G), G):

fit_error := evalf[4](map(u -> eval(zg, G=u), [seq(0..1, 0.01)]))

[0.15e-1, 0.3e-2, 0.5e-2, 0.2e-2, 0.2e-2, 0., 0.2e-2, -0.4e-2, -0.3e-2, -0.3e-2, -0.2e-2, 0.1e-2, -0.4e-2, -0.4e-2, -0.1e-2, -0.4e-2, -0.2e-2, 0.5e-3, -0.27e-2, -0.16e-2, -0.35e-2, -0.6e-3, 0., 0.2e-3, -0.16e-2, 0.8e-3, 0.5e-3, -0.22e-2, -0.6e-3, -0.20e-2, 0.8e-3, 0.10e-2, -0.5e-3, 0.16e-2, -0.6e-3, 0.42e-2, 0.42e-2, 0.30e-2, 0.95e-3, 0.309e-2, 0.246e-2, 0.280e-2, 0.157e-2, 0.433e-2, 0.3612e-2, 0.191e-2, 0.247e-2, 0.188e-2, 0.351e-2, 0.24e-3, 0.256e-2, 0.295e-2, 0.275e-2, 0.299e-2, 0.65e-3, 0.271e-2, 0.229e-2, 0.58e-3, 0.40e-3, 0.24e-2, 0.12e-2, 0.19e-2, 0.5e-3, 0.6e-3, 0.10e-2, -0.2e-3, -0.6e-3, -0.10e-3, -0.92e-3, -0.98e-3, -0.19e-3, -0.199e-2, -0.288e-2, -0.349e-2, -0.258e-2, -0.166e-2, -0.225e-2, -0.358e-2, -0.265e-2, -0.176e-2, -0.248e-2, -0.318e-2, -0.339e-2, -0.304e-2, -0.336e-2, -0.198e-2, -0.170e-2, -0.224e-2, -0.183e-2, -0.236e-2, -0.79e-3, -0.107e-2, -0.133e-2, -0.57e-3, 0.30e-3, 0.29e-3, 0.140e-2, 0.1598e-2, 0.2962e-2, 0.3399e-2, 0.4400e-2]

(3)

 


 

Download LevelCurve2.mw

Here is a solution where each level curve is a colred differently

Maybe you would also have a legend for each level curve?
In this case you could ne interested by what acer did here 2D contour plot and legend

PS I removed the ouput for the uploaded display doesn't have colored level curves (but they are)
 

restart:

with(plots):

f := sin(x*y):
r := x = -2 .. 2, y = -2 .. 2:

cols := [seq(ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]),i=1..19)]:
plots:-display(
  plot3d(f, r, axes = framed, style=surface, color=gray),
  seq(contourplot3d(f, r, contours=[(j-10)/10], color=cols[j]), j=1..19)
);  

 

 


 

Download contourplot.mw

Luckily I had already done it, I just wasted a little time finding the code.
From memory I believe I recoded in Maple the code from R (not completely sure).

Here it is. 
I remember having done some tests and it seemed ok by the time , but feel free to inform me if it's not.


 

restart:

LSD_Sampling := proc(alpha, N)
  local c, V, X, n, u, q:

  uses Statistics:

  c     := log(1-alpha):
  V     := Sample(Uniform(0, 1), N):
  X     := Vector[column](N):
  for n from 1 to N do
     if V[n] > alpha then
        X[n] := 1:
     else
        u := Sample(Uniform(0, 1), 1)[1]:
        q := 1 - exp(c*u):
        if V[n] < q^2 then
           X[n] := ceil(1+log(V[n])/log(q)):
        elif V[n] < q then
           X[n] := 1:
        else
           X[n] := 2:
        end if:
     end if:
  end do:

  return X:
end proc:

Statistics:-Histogram(LSD_Sampling(0.9, 10000))

 

 


 

Download LogarithmicSeriesSampling.mw


BTW: I'm curious, what application of the Logarithmic Series distribution are you interested in?
 

Even Re(code) doesn't work. The reason is: Maple doesn't know if the variables are real or complex !

Assuming all the quanties are real a simple way is

subs(J=I, collect(subs(I=J, code), J));

Other possibility: use assume, or assuming (see help pages)

I already ask this question a few years ago, you can trace the answers here

Question: I'd like to have some tips about Version Control (source control) process

The solution I use:

  1. I assume you have 2 mw files, let's say File1.mw and File2.mw whoch correspond to different versions of the same code.
     
  2. Open File1.mw and Export it to maple input format (I choose the same name as the mw file)
     
  3. Open File2.mw and Export it to maple input format; two files File1.mpl and File2.mpl are created
     
  4. If not already installed, download NotePad++ (free on Windows ; Notepad++ has contextual coloring for Maple syntax)
     
  5. Open NotePad++ and load File1.mpl and File2.mpl
     
  6. In the menu of NotePad++ choose Run > Compare (this from memory because I do not have NotePad++ on this laptop)
    The differences between the two versions are then easily visble

Once done, you have the choice to keep developping with the Java interface of Maple (as usual), or to develop "within" NotePad++.
In this latter case you will have to load an mpl file in a Maple session to verify that its content gives the expected result.
It's up to you, I personally prefer developping "within" a worksheet, saving in a new mw file, and exporting it in a new mpl file.

PS : there exist probably tools apart NotePad++ that can "understand" Maple's syntax and propose contextual coloring.

I would like to complete acer's answer. Personally I like to use notations like theta__i instead of theta[i], because i is not only limited to integer values.
But using theta__i instead of theta[i] requires some precautions.
You will find below a few examples.

Be careful, when used in a procedure, variables of the form theta__||n (last example) are not implicitely considered as local but as global: this may cause unwanted behaviours.

 

# differences :

a := theta[2]:
b := theta__2:

whattype(a), whattype(b);

printf("%a, %a", theta[2], theta__2);
print(a,b);

indexed, symbol

 

theta[2], theta__2

 

theta[2], theta__2

(1)

# use of theta__  ; first case: wrong way

restart:

vars := seq(theta__i, i=1..4);

mean := Array(1..4, i -> diff(w(vars), theta__i));

 

vars := `#msub(mi("&theta;",fontstyle = "normal"),mi("i"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("i"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("i"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("i"))`

 

Array(%id = 18446744078232662014)

(2)

# use of theta__  ; second case: good way

restart:

vars := seq(theta__||i, i=1..4);

mean := Array(1..4, i -> diff(w(vars), theta__||i));

 

vars := `#msub(mi("&theta;",fontstyle = "normal"),mi("1"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("2"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("3"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("4"))`

 

Array(%id = 18446744078232662014)

(3)

# Why using theta__ instead of theta[..] can be useful?
# Because the index n in theta__n can be of non numeric type:

restart:
MyIndices := [a, b, c, d];
vars := seq(theta__||n, n in MyIndices);

mean := Array(1..4, i -> diff(w(vars), vars[i]));

MyIndices := [a, b, c, d]

 

vars := `#msub(mi("&theta;",fontstyle = "normal"),mi("a"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("b"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("c"))`, `#msub(mi("&theta;",fontstyle = "normal"),mi("d"))`

 

Array(%id = 18446744078232662014)

(4)

restart:

f := proc()
local n:
for n from 1 to 4 do
  theta__||n := n
end do:
end proc:

f():
seq(theta__||i, i=1..4);

1, 2, 3, 4

(5)

 


 

Download indexed_theta.mw

 

Two ways (probably among others)

restart:
vars := seq(theta[i], i=1..4);
mean := Array(1..4, i -> diff(w(vars), theta[i]));

# or
mean := Array(1..4, i -> diff(w(vars), vars[i]));

Download mean.mw

In case you would be interested, give a lookk to this site :

https://oeis.org/?language=english

Next just copy-paste 2,3,4,6,8,9,10,12,14,15,16,18 in the search field and enjoy

 

I personally prefer to define f, f2 and f3 defore instanciate x and y (even if the result is the same).
The main point is in the definition of f3 which should had been f3 := (x, y) -> f(x, y)+f2(x, y)


 

restart:

f  := (x,y) -> x+y;
f2 := (x,y) -> x^2+y^3;
f3 := (x,y) -> f(x,y)+f2(x,y);

proc (x, y) options operator, arrow; x+y end proc

 

proc (x, y) options operator, arrow; x^2+y^3 end proc

 

proc (x, y) options operator, arrow; f(x, y)+f2(x, y) end proc

(1)

f(x,y) ; f2(x,y) ; f3(x,y)

x+y

 

y^3+x^2

 

y^3+x^2+x+y

(2)

x := 5;
y := 10;
f(x, y) ; f2(x, y) ; f3(x, y)

5

 

10

 

15

 

1025

 

1040

(3)

x := 'x':
y := 'y':
x;
y;
f(x, y) ; f2(x, y) ; f3(x, y)

x

 

y

 

x+y

 

y^3+x^2

 

y^3+x^2+x+y

(4)

restart:

x := 5;
y := 10;

5

 

10

(5)

f  := (x,y) -> x+y;
f2 := (x,y) -> x^2+y^3;
f3 := (x,y) -> f(x,y)+f2(x,y);

proc (x, y) options operator, arrow; x+y end proc

 

proc (x, y) options operator, arrow; x^2+y^3 end proc

 

proc (x, y) options operator, arrow; f(x, y)+f2(x, y) end proc

(6)

f(x, y) ; f2(x, y) ; f3(x, y)

15

 

1025

 

1040

(7)

 


 

Download f.mw

Maybe this?


 

restart

u := sin(x+y):

MyPlot := plot3d(u, x = 0 .. 2*Pi, y = 0 .. 2*Pi):

plottools:-getdata(MyPlot);

["grid", [0. .. 6.28318530717958, 0. .. 6.28318530717958, -1. .. 1.], Vector(4, {(1) = ` 1..49 x 1..49 `*Array, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*C_order})]

(1)

M := plottools:-getdata(MyPlot)[-1];

M := Vector(4, {(1) = ` 1..49 x 1..49 `*Array, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*C_order})

(2)

# search Export in the help pages for other formats

MyFile := FileTools:-JoinPath(["Desktop/Maple/M.csv"], base=homedir):
Export(MyFile, M):

56423

(3)

 


 

Download Export.mw

I believe Christopher has put the finger on the problem when saying "...  is it the y portion? "

The structure of the pde is rather simple with the form F . grad(Phi) = 0 (I use Phi instead of w).
For these equations the key is in finding the characteristic curves. The problem seems to be: Maple cannot find a closed form expression for the second component of these curves.

(Nota: Maple find the characteristic curves in closed form if alpha=0 ... and successes too in solving the pde)


 

restart:

pde :=  diff(Phi(x,y,z),x)
        +
        (y^2- a*exp(alpha*x)*(x*y-1))*diff(Phi(x,y,z),y)
        +
        (c*exp(beta*x)*z^2+b*exp(-beta*x))*diff(Phi(x,y,z),z)
        =
        0;

diff(Phi(x, y, z), x)+(y^2-a*exp(alpha*x)*(x*y-1))*(diff(Phi(x, y, z), y))+(c*exp(beta*x)*z^2+b*exp(-beta*x))*(diff(Phi(x, y, z), z)) = 0

(1)

# Write the pde as F . grad(Phi) = 0 where F is the vector (A, B, C)
# Next parameterize F as a function of s : F(s) = (A(s), B(s), C(s))

A := s -> 1;
B := s -> (y(s)^2- a*exp(alpha*x(s))*(x(s)*y(s)-1));
C := s -> (c*exp(beta*x(s))*z(s)^2+b*exp(-beta*x(s)));

proc (s) options operator, arrow; 1 end proc

 

proc (s) options operator, arrow; y(s)^2-a*exp(alpha*x(s))*(x(s)*y(s)-1) end proc

 

proc (s) options operator, arrow; c*exp(beta*x(s))*z(s)^2+b*exp(-beta*x(s)) end proc

(2)

# Define the characteristic curve associated to the pde as the curve defined by (U(s), V(s), W(s))
# where :
#   diff(U(s), s) = A(s)
#   diff(V(s), s) = B(s)
#   diff(W(s), s) = C(s)
#
# Then Phi(s) is constant along each characteristic curve
#
# Characteristic curve, component 1

eq1 := diff(U(s), s) = A(s);
dsolve(eq1, U(s)):
U := unapply(rhs(%), s);

diff(U(s), s) = 1

 

proc (s) options operator, arrow; s+_C1 end proc

(3)

# Characteristic curve, component 3


eq3 := subs({x(s)=U(s), z(s)=W(s)}, diff(W(s), s) = C(s));

dsolve(eq3, W(s))

diff(W(s), s) = c*exp(beta*(s+_C1))*W(s)^2+b*exp(-beta*(s+_C1))

 

W(s) = -(1/2)*(exp(beta*_C1)*exp(beta*s)*beta+tan((1/2)*((exp(beta*_C1))^2*(exp(beta*s))^2*(4*b*c-beta^2))^(1/2)*(_C2-s)/(exp(beta*_C1)*exp(beta*s)))*((exp(beta*_C1))^2*(exp(beta*s))^2*(4*b*c-beta^2))^(1/2))/((exp(beta*_C1))^2*(exp(beta*s))^2*c)

(4)

# Characteristic curve, component 2


eq2 := subs({x(s)=U(s), y(s)=V(s)}, diff(V(s), s) = expand(B(s)));

infolevel[dsolve] := 4:
dsolve(eq2, V(s))

diff(V(s), s) = V(s)^2-a*exp(alpha*(s+_C1))*(s+_C1)*V(s)+a*exp(alpha*(s+_C1))

 

Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature
trying 1st order linear
trying Bernoulli
trying separable
trying inverse linear
trying homogeneous types:
trying Chini
differential order: 1; looking for linear symmetries
trying exact
Looking for potential symmetries
trying Riccati
trying inverse_Riccati
trying 1st order ODE linearizable_by_differentiation
--- Trying Lie symmetry methods, 1st order ---
 -> Computing symmetries using: way = 4
 -> Computing symmetries using: way = 2
 -> Computing symmetries using: way = 6
trying symmetry patterns for 1st order ODEs
-> trying a symmetry pattern of the form [F(x)*G(y), 0]
-> trying a symmetry pattern of the form [0, F(x)*G(y)]
-> trying symmetry patterns of the forms [F(x),G(y)] and [G(y),F(x)]

 -> Computing symmetries using: way = HINT
   -> Calling odsolve with the ODE diff(y(x) x) = 2*y(x)/x y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear

      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x) = y(x)/x y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)+y(x)*alpha y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful

   -> Calling odsolve with the ODE diff(y(x) x)+exp(alpha*(x+_C1))*a*K[1]*(x+_C1) y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature

      <- quadrature successful
   -> Calling odsolve with the ODE diff(y(x) x)+(y(x)*_C1*alpha+y(x)*alpha*x+y(x)+2*K[1])/(x+_C1) y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear

      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)+K[1] y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)+y(x)*alpha-K[1] y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)+y(x)*(_C1*alpha+alpha*x+1)/(x+_C1) y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x) = 0 y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x) = -y(x)*alpha y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)-2*y(x)/x y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)-(x*alpha*K[1]+y(x))/x y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)-(x*_C1*alpha*K[1]+y(x)*_C1+x*K[1]-K[1]*alpha)/(_C1*x-1) y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful

   -> Calling odsolve with the ODE diff(y(x) x)-y(x)/x y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)-y(x)*_C1/(_C1*x-1) y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
   -> Calling odsolve with the ODE diff(y(x) x)+(exp(alpha*_C1)*a*K[1]-2*y(x))/x y(x)
      *** Sublevel 2 ***
      Methods for first order ODEs:
      --- Trying classification methods ---
      trying a quadrature
      trying 1st order linear
      <- 1st order linear successful
-> trying a symmetry pattern of the form [F(x),G(x)]
-> trying a symmetry pattern of the form [F(y),G(y)]
-> trying a symmetry pattern of the form [F(x)+G(y), 0]

-> trying a symmetry pattern of the form [0, F(x)+G(y)]
-> trying a symmetry pattern of the form [F(x),G(x)*y+H(x)]

-> trying a symmetry pattern of conformal type

 

# Note there is no explicit solution found: I suppose this is the reason why Maple cannot solve the pde

eq2 := subs({x(s)=U(s), y(s)=V(s), alpha=0}, diff(V(s), s) = expand(B(s)));

infolevel[dsolve] := 4:
dsolve(eq2, V(s))

diff(V(s), s) = V(s)^2-a*exp(0)*(s+_C1)*V(s)+a*exp(0)

 

Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature
trying 1st order linear
trying Bernoulli
trying separable
trying inverse linear
trying homogeneous types:
trying Chini
Chini's absolute invariant is: a*(s+_C1)^2
differential order: 1; looking for linear symmetries
trying exact
Looking for potential symmetries

Found potential symmetries: [0 1] [exp(-V)/a exp(-V)*(s+_C1)]
Proceeding with integration step.

 

V(s) = -(-_C1*exp(-(1/2)*a*s^2)/exp(s*_C1*a)-s*exp(-(1/2)*a*s^2)/exp(s*_C1*a)+1/(-(1/2)*exp(-(1/2)*_C1^2*a)*erf((1/2)*(-2*a)^(1/2)*(s+_C1))*(-2*Pi*a)^(1/2)+_C2))*exp(s*_C1*a)*a/exp(-(1/2)*a*s^2)

(5)

pdsolve(subs(alpha=0, pde), Phi(x,y,z))

Methods for first order ODEs:

--- Trying classification methods ---
trying a quadrature
trying 1st order linear
trying Bernoulli
trying separable
trying inverse linear
trying homogeneous types:
trying Chini
Chini's absolute invariant is: a*x^2
differential order: 1; looking for linear symmetries
trying exact
Looking for potential symmetries
Found potential symmetries: [0 1] [exp(-y) exp(-y)*a*x]
Proceeding with integration step.

Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature
trying 1st order linear
trying Bernoulli
trying separable
trying inverse linear
trying homogeneous types:
trying Chini
Chini's absolute invariant is: beta^2/(c*b)
<- Chini successful

Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature
trying 1st order linear
<- 1st order linear successful

 

Phi(x, y, z) = _F1(-(-erf((1/2)*(-2*a)^(1/2)*x)*a*x+y*erf((1/2)*(-2*a)^(1/2)*x)+(-2*a/Pi)^(1/2)*exp((1/2)*a*x^2))/((-2*a/Pi)^(1/2)*(-a*x+y)), beta*(2*beta*arctan(beta*(2*c*exp(beta*x)*z+beta)/(4*b*beta^2*c-beta^4)^(1/2))-(beta^2*(4*b*c-beta^2))^(1/2)*x)/(beta^2*(4*b*c-beta^2))^(1/2))

(6)

 


 

Download CharacteristicCurves.mw


 

Maybe this ?

Watchout: you did not say what to do if A inter B is the empty set
 


 

restart:

p := proc(A::set, B::set, C::set)
   local AB, ABC:
   AB := A intersect B:
   if AB <> {} then
      ABC := AB intersect C:
      if ABC <> {} then
         return ABC
      else
         return `union`(seq(AB *~ C[n], n=1..nops(c)))
      end if;
   else
      error  "It's not specified what to do when A inter B = {} "
   end if:
end proc:

a := {x6,x4,x2,x7,x8,x9,x10}:  
b := {x2,x3,x5,x8}:
c := {x4,x9,x11,x12,x13}:

p(a, b, c)

{x2*x11, x2*x12, x2*x13, x2*x4, x2*x9, x8*x11, x8*x12, x8*x13, x8*x4, x8*x9}

(1)

a := {x6,x4}:  
b := {x5,x8}:
c := {x12,x13}:

p(a, b, c)

Error, (in p) It's not specified what to do when A inter B = {}

 

a := {x6,x4}:  
b := {x4, x5,x8}:
c := {x4, x12,x13}:

p(a, b, c)

{x4}

(2)

 


 

Download sets.mw

Sorry, I do not know about this game and I can give an answer because this later depends on the strategy used to form the new teams.
I consider here 2 strategies (one could probably imagine a few others) . One gives the result proposed by CarlLove after Kitonum's correction, the second a rather different result.

You will find the correct values and a numerical simulation for each strategy


 

restart:

with(Statistics):
with(combinat):

# How is the game managed (I do not know about it)?
#
# Strategy 1:
# One decides to form the red team first by randomly picking 5 out of 15 members
# Next one decides to form the green team by randomly picking 5 out of the 10 remaining members
# The 5 remaining ones form the last (orange) team
#
#
# Strategy 2
# One randomly form 3 new teams of 5 members and randomly decide to color thes teams as red, green, orange
 

 

STRATEGY 1

 

ProbabilityToPutFiveBluesInTheRedTeam := mul((6-n)/(15-n), n=0..4);
ProbabilityToPutTheLastBlueInTheOrangeTeam := 1/2;

ProbabilityStrategy1 := ProbabilityToPutFiveBluesInTheRedTeam * ProbabilityToPutTheLastBlueInTheOrangeTeam ;

2/1001

 

1/2

 

1/1001

(1)

# simulation strategy 1

OriginalBlueTeam   := [B$6]:
OriginalYellowTeam := [Y$9]:
candidates         := [OriginalBlueTeam[], OriginalYellowTeam[]]:


N := 10^6:


NewTeams := Matrix(N, 3):
for n from 1 to N do
   x := randperm(candidates):
   NewTeams[n, 1] := numelems(select(has, x[ 1.. 5], B)):
   NewTeams[n, 2] := numelems(select(has, x[ 6..10], B)):
   NewTeams[n, 3] := numelems(select(has, x[11..15], B)):
end do:


DesiredNumberOfBluesInNewTeams := [5, 0, 1]; # teams ordered red, green, orange

FrequencyStrategy1 := add( map(n -> if [entries(NewTeams[n,..], nolist)] = [5, 0, 1] then 1 end if, [$1..N]) ) / N

[5, 0, 1]

 

201/200000

(2)

evalf([ProbabilityStrategy1, FrequencyStrategy1]);

[0.9990009990e-3, 0.1005000000e-2]

(3)

 

STRATEGY 2

 

ProbabilityToPutFiveBluesInTheSameTeam   := mul((6-n)/(15-n), n=0..4) * 3;
ProbabilityToPutTheLastBlueInAnotherTeam := 1;

ProbabilityStrategy2 := ProbabilityToPutFiveBluesInTheSameTeam * ProbabilityToPutTheLastBlueInAnotherTeam ;

6/1001

 

1

 

6/1001

(4)

# simulation strategy 2

OriginalBlueTeam   := [B$6]:
OriginalYellowTeam := [Y$9]:
candidates         := [OriginalBlueTeam[], OriginalYellowTeam[]]:


N := 10^6:


NewTeams := Matrix(N, 3):
for n from 1 to N do
   x := randperm(candidates):
   NewTeams[n, 1] := numelems(select(has, x[ 1.. 5], B)):
   NewTeams[n, 2] := numelems(select(has, x[ 6..10], B)):
   NewTeams[n, 3] := numelems(select(has, x[11..15], B)):
end do:


DesiredPartitions := { permute([5, 0, 1])[] };

FrequencyStrategy2 := add( map(n -> if member([entries(NewTeams[n,..], nolist)], DesiredPartitions) then 1 end if, [$1..N]) ) / N;

{[0, 1, 5], [0, 5, 1], [1, 0, 5], [1, 5, 0], [5, 0, 1], [5, 1, 0]}

 

3003/500000

(5)

evalf([ProbabilityStrategy2, FrequencyStrategy2]);

[0.5994005994e-2, 0.6006000000e-2]

(6)

 


 

Download Teams.mw

plots:-implicitplot(
     X^3-24.478115*X^2/Y-(0.2038793409e-2+19.08455282*Y/(Y-97.539))*X-.2550630228/Y,
     X=0..50,
     Y=0..1
);

You can add some options (see implicitplot help page), for onstance  gridrefine=2 to get a smoother plot.
 

The Grid Package and the Grid Computing Toolbox are two differents things :

  • the former is included in your Maple license
  • the latter has to be purchased separatetly

In the Maple 2015.2 help pages it's written 
Before installing the Grid Computing Toolbox, you must install and activate Maple 18. For details and installation instructions, see the Install.html file on the product disc.  
which seems to suggest that the Grid Computing Toolbox was introduced with Maple 18 (???)

1 2 3 4 5 6 7 Last Page 1 of 9