mmcdara

4029 Reputation

17 Badges

6 years, 183 days

MaplePrimes Activity


These are answers submitted by mmcdara

s this possible?

  • formally : NO YES
  • numerically YES

Here is my first reply where solve gave no results. I keep the file for it contains a plot you might find useful.

Half_Width_mmcdara.mw

And here is the file where the two points are computed formally:

restart

f := (x, x__0, S, C) -> -C*(exp(S*(x - x__0)^2*(-1/2)) - 1)*Heaviside(x - x__0)

proc (x, x__0, S, C) options operator, arrow; -C*(exp(-(1/2)*S*(x-x__0)^2)-1)*Heaviside(x-x__0) end proc

(1)

X0   := 0:
fDOG := unapply(f(x, X0, S__a, 1) - f(x, X0, S__d, 1), (x, S__a, S__d))

proc (x, S__a, S__d) options operator, arrow; -(exp(-(1/2)*S__a*x^2)-1)*Heaviside(x)+(exp(-(1/2)*S__d*x^2)-1)*Heaviside(x) end proc

(2)

DfDOG := diff(fDOG(x, S__a, S__d), x) assuming x > X0

S__a*x*exp(-(1/2)*S__a*x^2)-S__d*x*exp(-(1/2)*S__d*x^2)

(3)

solve({DfDOG=0, x > X0}, x) assuming S__a > S__d;
tpeak := eval(x, %[1])

[{x = 2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)}, x = x]

 

2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)

(4)

hpeak := fDOG(tpeak, S__a, S__d) assuming S__a > S__d, S__d > 0

-exp(-S__a*ln(S__a/S__d)/(S__a-S__d))+exp(-S__d*ln(S__a/S__d)/(S__a-S__d))

(5)

half_hpeak := hpeak/2;

-(1/2)*exp(-S__a*ln(S__a/S__d)/(S__a-S__d))+(1/2)*exp(-S__d*ln(S__a/S__d)/(S__a-S__d))

(6)

eq := fDOG(x, S__a, S__d) = half_hpeak assuming x > X0;

MyChange := ln(S__a/S__d) = (u/sqrt(2))^2*(S__a-S__d); # u is > 0
eval(eq, MyChange);
eval(%, [S__a=2, S__d=1]);
 

-exp(-(1/2)*S__a*x^2)+exp(-(1/2)*S__d*x^2) = -(1/2)*exp(-S__a*ln(S__a/S__d)/(S__a-S__d))+(1/2)*exp(-S__d*ln(S__a/S__d)/(S__a-S__d))

 

ln(S__a/S__d) = (1/2)*u^2*(S__a-S__d)

 

-exp(-(1/2)*S__a*x^2)+exp(-(1/2)*S__d*x^2) = -(1/2)*exp(-(1/2)*S__a*u^2)+(1/2)*exp(-(1/2)*S__d*u^2)

 

-exp(-x^2)+exp(-(1/2)*x^2) = -(1/2)*exp(-u^2)+(1/2)*exp(-(1/2)*u^2)

(7)

xofu := [solve(%, x)];
 

[(-2*ln(1/2+(1/2)*(1+2*(exp(-(1/2)*u^2))^2-2*exp(-(1/2)*u^2))^(1/2)))^(1/2), -(-2*ln(1/2+(1/2)*(1+2*(exp(-(1/2)*u^2))^2-2*exp(-(1/2)*u^2))^(1/2)))^(1/2), (-2*ln(1/2-(1/2)*(1+2*(exp(-(1/2)*u^2))^2-2*exp(-(1/2)*u^2))^(1/2)))^(1/2), -(-2*ln(1/2-(1/2)*(1+2*(exp(-(1/2)*u^2))^2-2*exp(-(1/2)*u^2))^(1/2)))^(1/2)]

 

2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d), -2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)

(8)

solve(MyChange, u);
x_half_peak := eval(xofu, u=%[1]);

2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d), -2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)

 

[(-2*ln(1/2+(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), -(-2*ln(1/2+(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), (-2*ln(1/2-(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), -(-2*ln(1/2-(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2)]

(9)

{abs~(x_half_peak)[]}

{2^(1/2)*abs(ln(1/2-(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), 2^(1/2)*abs(ln(1/2+(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2)}

(10)

sa := 2:
sd := 1:

eval(x_half_peak, [S__a=sa, S__d=sd]):
select(is@`>`, %, 0);
evalf(%);

[(-2*ln(1/2+(1/4)*2^(1/2)))^(1/2), (-2*ln(1/2-(1/4)*2^(1/2)))^(1/2)]

 

[.5627560464, 1.960150176]

(11)

# Applying "select" directly on x_half_peak doesn't work.
# I believe the following is correct but this must be checked

`correct?` := {abs~(x_half_peak)[]};
eval(`correct?`, [S__a=sa, S__d=sd]):
evalf(%);

{2^(1/2)*abs(ln(1/2-(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), 2^(1/2)*abs(ln(1/2+(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2)}

 

{.5627560463, 1.960150176}

(12)

 

Download Formal__Half_Width_mmcdara.mw

This has been made possible after introducing a change of variables

u = sqrt(2)*sqrt((S__a-S__d)*ln(S__a/S__d))/(S__a-S__d)


 

restart:

H_lines := seq(CURVES([[0, i], [1, i]], LINESTYLE(2)), i in [seq](0..1, 0.1)):
V_lines := seq(CURVES([[i, 0], [i, 1]], LINESTYLE(2)), i in [seq](0..1, 0.1)):
H_vec   := plottools:-arrow([0, 0], [0.1, 0], .005, .02, .4, color=black):
V_vec   := plottools:-arrow([0, 0], [0, 0.1], .005, .02, .4, color=black):

p := PLOT(H_lines, V_lines, H_vec, V_vec):
plots:-display(p, axes=none);

 

T := (a, b) -> plottools:-transform((x, y) -> [x+a*y, y+b*x]):

plots:-display(T(2, 3)(p), axes=none)

 

 


 

Download One_method.mw

for some unknown reason Maple (2015... what is the version you use?) is wrong when compared ro reference Statistical Softwares

Result obtained with R:

shapiro.test(S)

	Shapiro-Wilk normality test

data:  S
W = 0.98259, p-value = 0.3481

for fun: SWWT.mw

For your information:
I had already mentionned years ago some problemes with Statistics[ChiSquareGoodnessOfFitTest]

As you seem familiar with dedicated statistical softwares, I can't help but recommend you to use them to verify the results provided by Maple, before using this latter with full confidence.

It is highly likely that the formal integration cannot be done.

Ir remains you the numerical integration, which requires setting the values of p, r and "i" (the imaginary unit?).


 

restart:

local exp;  # In Maple "exp" stands for "exponential",
            # thus you must say explicitely that "exp" is a local name

Is "i" the imaginary unit?
If it is so you must use "I"instead

exp := -sin(alpha)*i*r*(-sin(alpha)*cos(phi)*cos(theta)+sin(theta)*cos(alpha))/(4*sqrt(-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)*Pi*(-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)*(-2+sqrt(2))*Pi)

-(1/4)*sin(alpha)*i*r*(-sin(alpha)*cos(phi)*cos(theta)+sin(theta)*cos(alpha))/((-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)^(3/2)*Pi^2*(-2+2^(1/2)))

(1)

You do not need any assumptions for they are implicitely defined by the integration bounds.
The conditions about theta can or cannot be usefull, we will see this later

Int(exp*p^2*sin(alpha), p = 0 .. 1, alpha = 0 .. (1/4)*Pi, phi = 0 .. 2*Pi)

Int(-(1/4)*sin(alpha)^2*i*r*(-sin(alpha)*cos(phi)*cos(theta)+sin(theta)*cos(alpha))*p^2/((-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)^(3/2)*Pi^2*(-2+2^(1/2))), p = 0 .. 1, alpha = 0 .. (1/4)*Pi, phi = 0 .. 2*Pi)

(2)

Let us to compute this "simple integral

I_p := Int(exp*p^2*sin(alpha), p = 0 .. 1)

Int(-(1/4)*sin(alpha)^2*i*r*(-sin(alpha)*cos(phi)*cos(theta)+sin(theta)*cos(alpha))*p^2/((-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)^(3/2)*Pi^2*(-2+2^(1/2))), p = 0 .. 1)

(3)

# The integrand is of the form

b := add(B[i]*p^i, i=0..2):
c := add(C[i]*p^i, i=0..2):

J_p := A*p^2 / ( sqrt(b)*c )

A*p^2/((p^2*B[2]+p*B[1]+B[0])^(1/2)*(p^2*C[2]+p*C[1]+C[0]))

(4)

# The formal integration alreaqy fayls on this even simpler expression

infolevel[int] := 4:
p^2 / sqrt(b);
int(%, p=0..1)

p^2/(p^2*B[2]+p*B[1]+B[0])^(1/2)

 

Definite Integration:   Integrating expression on p=0..1
Definite Integration:   Using the integrators [distribution, piecewise, series, o, polynomial, ln, lookup, cook, ratpoly, elliptic, elliptictrig, meijergspecial, improper, asymptotic, ftoc, ftocms, meijerg, contour]
Definite Integration:   Trying method distribution.
Definite Integration:   Trying method piecewise.
Definite Integration:   Trying method series.
Definite Integration:   Trying method o.
Definite Integration:   Trying method polynomial.
Definite Integration:   Trying method ln.
Definite Integration:   Trying method lookup.
LookUp Integrator:   unable to find the specified integral in the table
Definite Integration:   Trying method cook.
Definite Integration:   Trying method ratpoly.
Definite Integration:   Trying method elliptic.
int/elliptic: trying elliptic integration
Definite Integration:   Trying method elliptictrig.
Definite Integration:   Trying method meijergspecial.
Definite Integration:   Trying method improper.
Definite Integration:   Trying method asymptotic.
Definite Integration:   Trying method ftoc.

Warning,  computation interrupted

 

Attempt to do numerical integration

infolevel[int] := 0:

I_theta := proc(__theta, __r, __i,  meth)
  eval(exp*p^2*sin(alpha), [theta=__theta, r=__r, i=1]):
  __i*evalf(Int(%, p=0..1, alpha=0..(1/4)*Pi, phi=0..2*Pi, method = meth))
end proc:

I_theta(0, 1, I, _d01ajc);     # fails

I*(Int(Int(Int(-0.4324151988e-1*sin(alpha)^3*cos(phi)*p^2/(1.-2.*cos(alpha)*p+p^2)^(3/2), p = 0. .. 1.), alpha = 0. .. .7853981635), phi = 0. .. 6.283185308))

(5)

I_theta(0, 1, I, _CubaVegas);  # interrupted

Warning,  computation interrupted

 

I_theta(0, 1, I, _MonteCarlo); # fails
 

I*(Int(Int(Int(-0.4324151988e-1*sin(alpha)^3*cos(phi)*p^2/(1.-2.*cos(alpha)*p+p^2)^(3/2), p = 0. .. 1.), alpha = 0. .. .7853981635), phi = 0. .. 6.283185308))

(6)

Monte-Carlo integration can be done "by hand":

with(Statistics):
S_N     := 10^3:   # Change this value to increase the precision
                   # A good practice is to repeat the integration with different input
                   # samples of the same size and assess the integration error.
                   # From a statistical point of view it can be better to do K Monte-Carlo
                   # integrations with S_N points that a single one with K*S_N points.
                   # The Bootstrap method (see help pages) can be use then to refine
                   # the estimation of the integral.

S_p     := Sample(Uniform(0, 1)   , S_N):
S_alpha := Sample(Uniform(0, Pi/4), S_N):
S_phi   := Sample(Uniform(0, 2*Pi), S_N):

S_I := proc(__theta, __i, __r, M)
  local f:
  f := unapply(eval(exp*p^2*sin(alpha), [theta=__theta, r=__r, i=1]), [p, alpha, phi]):
  __i * add(f(S_p[n], S_alpha[n], S_phi[n]), n=1..M) / M:
  return evalf(simplify(%));
end proc:
  

# Example

StringTools:-FormatTime("%H:%M:%S");
DocumentTools:-Tabulate(
  [
    plot3d('Re(S_I(P, R, I, 100))', P=0..1, R=0..1, labels=["p", "r", "Int"], grid=[20, 20], title="Real part"),
    plot3d('Im(S_I(P, R, I, 100))', P=0..1, R=0..1, labels=["p", "r", "Int"], grid=[20, 20], title="Imaginary part")
  ],
  width=60
):
StringTools:-FormatTime("%H%M%S");

"10:00:11"

 

"100014"

(7)

`I'm not in a hurry` := false:

if `I'm not in a hurry` then
  DocumentTools:-Tabulate(
    [
      plot3d('Re(S_I(P, R, I, S__N))', P=0..1, R=0..1, labels=["p", "r", "Int"], title="Real part"),
      plot3d('Im(S_I(P, R, I, S__N))', P=0..1, R=0..1, labels=["p", "r", "Int"], title="Imaginary part")
    ],
    width=60
  ):
end if:

# Sketch of the intergation error
# (just to give you an idea of the variability of the Monte Carlo intagration)



S_I2 := proc(__theta, __i, __r, M1, M2)
  local f:
  f := unapply(eval(exp*p^2*sin(alpha), [theta=__theta, r=__r, i=1]), [p, alpha, phi]):
  __i * add(f(S_p[n], S_alpha[n], S_phi[n]), n=M1..M2) / (M2-M1+1):
  return evalf(simplify(%));
end proc:

for k from 1 to 10 do
  printf("Integration over block %2d = %+1.3e  %+1.3e * I\n", k, (Re, Im)(S_I2(0, 1, I, 1+(k-1)*100, k*100)));
end do;

Integration over block  1 = -7.288e-05  -1.275e-04 * I
Integration over block  2 = -4.107e-05  -7.174e-05 * I
Integration over block  3 = -1.142e-06  +8.639e-06 * I
Integration over block  4 = -1.725e-06  +1.745e-05 * I
Integration over block  5 = +1.180e-04  +7.142e-05 * I
Integration over block  6 = +8.563e-05  +2.254e-04 * I
Integration over block  7 = -1.454e-06  -1.150e-04 * I
Integration over block  8 = +6.503e-05  +1.921e-04 * I
Integration over block  9 = -7.424e-05  -2.054e-04 * I
Integration over block 10 = +2.240e-05  -3.212e-05 * I

 

V := Vector(S_N, n -> evalf(eval(exp*p^2*sin(alpha), [p=S_p[n], alpha=S_alpha[n], phi=S_phi[n], theta=0, r=1]))):
Mean(eval(V, i=1));
Boot := Bootstrap('Mean', eval(V, i=1), replications=1000, output=['value', 'standarderror']);

HFloat(-1.9365481621754508e-4)

 

[HFloat(-1.824804060830804e-4), 0.352362351056818162e-3]

(8)

# The Bootstrap estimation of the integral is val = -0.000186623004800676
# and its error is err = 0.000350218565372046231
#
# A 95% bilateral confidence interval is

Boot[1]+Boot[2]*Quantile(Normal(0, 1), 0.025) .. Boot[1]+Boot[2]*Quantile(Normal(0, 1), 0.975)

HFloat(-8.730979236620875e-4) .. HFloat(5.081371114959267e-4)

(9)

 


 

Download Int_mmcdara.mw

 

A possibility
(it seems that all dthe A[i] but A[2] are null, does this seem right to you?)

restart

with(LinearAlgebra):

with(PDEtools):

Setup(mathematicalnotation = true);

Setup(mathematicalnotation = true)

(1)

assume(x::real);

assume(t::real);

assume(xi::real);

assume(tau::real);

interface(showassumed = 0);

0

(2)

alias(u = u(x, t), q = q(x, t));

u, q

(3)

Eq := diff(u, `$`(x, 2))-u^2-a*u-b*x-c*x^2 = 0;

diff(diff(u, x), x)-u^2-a*u-b*x-c*x^2 = 0

(4)

va := {u = A[0]/(x-x[0])^p};

{u = A[0]/(x-x[0])^p}

 

{diff(diff(u, x), x) = A[0]*(x-x[0])^(-p-2)*p*(p+1)}

 

A[0]*(x-x[0])^(-p-2)*p*(p+1)

(5)

Eq1 := simplify(eval(Eq, va), size);

A[0]*p^2/((x-x[0])^p*(x-x[0])^2)+A[0]*p/((x-x[0])^p*(x-x[0])^2)-A[0]^2/((x-x[0])^p)^2-a*A[0]/(x-x[0])^p-b*x-c*x^2 = 0

(6)

va1 := {x = x[0]+eta, x-x[0] = eta};

{x = x[0]+eta, x-x[0] = eta}

(7)

Eq2 := simplify(eta^(p+2)*simplify(eval(Eq1, va1)));

A[0]*p^2+A[0]*p-eta^(-p+2)*A[0]^2-eta^2*a*A[0]-eta^(p+3)*b-eta^(p+2)*b*x[0]-eta^(p+4)*c-2*eta^(p+3)*c*x[0]-eta^(p+2)*c*x[0]^2 = 0

(8)

``

simplify(eval(Eq2, [p = 2]), size);

-c*eta^6-2*c*eta^5*x[0]-c*eta^4*x[0]^2-b*eta^5-b*eta^4*x[0]-a*eta^2*A[0]-A[0]^2+6*A[0] = 0

(9)

u = expand(sum(A[m]*eta^m, m = 0 .. infinity))/eta^2;

u = (sum(A[m]*eta^m, m = 0 .. infinity))/eta^2

(10)

Eq3 := subs(u = expand(sum(A[m]*eta^m, m = 0 .. infinity))/eta^2, Eq);

diff(diff((sum(A[m]*eta^m, m = 0 .. infinity))/eta^2, x), x)-(sum(A[m]*eta^m, m = 0 .. infinity))^2/eta^4-a*(sum(A[m]*eta^m, m = 0 .. infinity))/eta^2-b*x-c*x^2 = 0

(11)

GenSys := proc(N)
   local expr         := collect(numer(normal(eval(lhs(Eq3), infinity=N))), eta);
   local coefficients := [coeffs(expr, eta, 'k')];
   local unknowns     := [select(has, indets(coefficients, name), A)[]];
   local NbC          := numelems(coefficients);
  
   coefficients := select(has, coefficients=~[0$NbC], A):
   NbC          := numelems(coefficients);

   return [coefficients, unknowns]
end proc:

deg := 1:
g   := GenSys(deg);
sol := solve(g[1], g[2]);

[[-A[0]^2 = 0, -a*A[1] = 0, -a*A[0]-A[1]^2 = 0, -2*A[0]*A[1] = 0], [A[0], A[1]]]

 

[[A[0] = 0, A[1] = 0]]

(12)

deg := 2:
g   := GenSys(deg);
sol := solve(g[1], g[2]);

[[-A[0]^2 = 0, -c*x^2-a*A[2]-b*x-A[2]^2 = 0, -a*A[1]-2*A[1]*A[2] = 0, -a*A[0]-2*A[0]*A[2]-A[1]^2 = 0, -2*A[0]*A[1] = 0], [A[0], A[1], A[2]]]

 

[[A[0] = 0, A[1] = 0, A[2] = RootOf(c*x^2+_Z^2+_Z*a+b*x)]]

(13)

vals := allvalues(eval(A[2], sol[])):
map(u -> [remove(has, sol[], A[2])[], A[2]=u], [vals]):
print~(%):

[A[0] = 0, A[1] = 0, A[2] = -(1/2)*a+(1/2)*(-4*c*x^2+a^2-4*b*x)^(1/2)]

 

[A[0] = 0, A[1] = 0, A[2] = -(1/2)*a-(1/2)*(-4*c*x^2+a^2-4*b*x)^(1/2)]

(14)

 

NULL

 

Download PA_mmcdara.mw

Improvement (automatic treatment of equations of the form A[m]=RootOf(...))

deg := 10: # for instance
g   := GenSys(deg):
sol := solve(g[1], g[2]):

NonZero   := map(u -> if rhs(u) <> 0 then u end if, sol[]);

if NonZero <> [] then 
  NonZero_u := lhs~(NonZero):
  NonZero_v := [allvalues~(rhs~(NonZero))]:
  map(
    i -> 
      seq(
        [
          remove(has, sol[], NonZero_u)[], 
          NonZero_u[i]=NonZero_v[i][j]
        ], 
        j=1..numelems(NonZero_v[i])
      ),
      [$1..numelems(NonZero_u)]
  ):
  print~(%):
else
  print(sol[])
end if:

 


Looking at the integrand it appears it is almost constant wrt y and oscillating wrt x.
Using evalf(Int(integrand, x=..., y=..., method=...)) seems to be challenging for no integration method seems perfecly suited to handle this situation.
Here is a proposal (all depends if the precision you need)

(PS: I use CurveFitting:-Spline only because I'm a lazzy person)

restart:

afa:=0.3:
vh:=3.5:
u:=3.12:
mu:=5.5:
gama:=-4*10^(-29)*(1-cos(6*afa))*(1-1*10^(-8)*I):
d1:=1.78*10^(-9):
d2:=48.22*10^(-9):
HBAR:=1.05457266*10^(-34):
ME:=9.1093897*10^(-31):
ELEC:=1.60217733*10^(-19):
Kh:=2.95*10^10:
kc:=sqrt(2*ME*ELEC/HBAR^2):
k:=kc*sqrt(mu):
k0:=sqrt(k^2-k^2*sin(x)^2):
kh:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2 - k^2 * sin(x)^2 * sin(y)^2):
khg:=sqrt(k^2-(2*Kh*sin(afa/2)*sin(afa/2)-k*sin(x)*cos(y))^2-(2*Kh*sin(afa/2)*cos(afa/2)+k*sin(x)*sin(y))^2):
kg1:=sqrt(k^2-(Kh*cos(Pi/3-afa)-k*sin(x)*cos(y))^2-(Kh*sin(Pi/3-afa)+k*sin(x)*sin(y))^2):
kg2:=sqrt(k^2-(Kh*cos(afa)-k*sin(x)*cos(y))^2-(k*sin(x)*sin(y)-Kh*sin(afa))^2):
k0pl:=sqrt(k^2-k^2*sin(x)^2+kc^2*vh-kc^2*u):
k0mi:=sqrt(k^2-k^2*sin(x)^2-kc^2*vh-kc^2*u):
khpl:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2-k^2*sin(x)^2*sin(y)^2+kc^2*vh-kc^2*u):
khmi:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2-k^2*sin(x)^2*sin(y)^2-kc^2*vh-kc^2*u):
k0plpl:=sqrt(k^2-k^2*sin(x)^2+2*kc^2*vh):
k0mimi:=sqrt(k^2-k^2*sin(x)^2-2*kc^2*vh):
khplpl:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2-k^2*sin(x)^2*sin(y)^2+2*kc^2*vh):
khmimi:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2-k^2*sin(x)^2*sin(y)^2-2*kc^2*vh):
khgplpl:=sqrt(k^2-(2*Kh*sin(afa/2)*sin(afa/2)-k*sin(x)*cos(y))^2-(2*Kh*sin(afa/2)*cos(afa/2)+k*sin(x)*sin(y))^2+2*kc^2*vh):
khgmimi:=sqrt(k^2-(2*Kh*sin(afa/2)*sin(afa/2)-k*sin(x)*cos(y))^2-(2*Kh*sin(afa/2)*cos(afa/2)+k*sin(x)*sin(y))^2-2*kc^2*vh):
kg1plpl:=sqrt(k^2-(Kh*cos(Pi/3-afa)-k*sin(x)*cos(y))^2-(Kh*sin(Pi/3-afa)+k*sin(x)*sin(y))^2+2*kc^2*vh):
kg1mimi:=sqrt(k^2-(Kh*cos(Pi/3-afa)-k*sin(x)*cos(y))^2-(Kh*sin(Pi/3-afa)+k*sin(x)*sin(y))^2-2*kc^2*vh):
kg2plpl:=sqrt(k^2-(Kh*cos(afa)-k*sin(x)*cos(y))^2-(k*sin(x)*sin(y)-Kh*sin(afa))^2+2*kc^2*vh):
kg2mimi:=sqrt(k^2-(Kh*cos(afa)-k*sin(x)*cos(y))^2-(k*sin(x)*sin(y)-Kh*sin(afa))^2-2*kc^2*vh):
A1:=1/(1+I*ME*gama/(HBAR^2*k0pl))*exp(I*k0pl*d1)/2:
B1:=1/(1+I*ME*gama/(HBAR^2*k0pl))*exp(I*k0pl*d1)/2:
A2:=1/(1+I*ME*gama/(HBAR^2*khpl))*exp(I*khpl*d1)/2:
B2:=1/(1+I*ME*gama/(HBAR^2*khpl))*exp(I*khpl*d1)/2:
A3:=1/(1+I*ME*gama/(HBAR^2*k0mi))*exp(I*k0mi*d1)/2:
B3:=1/(1+I*ME*gama/(HBAR^2*k0mi))*exp(I*k0mi*d1)/2:
A4:=1/(1+I*ME*gama/(HBAR^2*khmi))*exp(I*khmi*d1)/2:
B4:=1/(1+I*ME*gama/(HBAR^2*khmi))*exp(I*khmi*d1)/2:

T1:=1/4*Re(abs(A1)^2*k0plpl*exp(I*(k0plpl-conjugate(k0plpl))*d2)+abs(A1)^2*kg1plpl*exp(I*(kg1plpl-conjugate(kg1plpl))*d2)+abs(A2)^2*khplpl*exp(I*(khplpl-conjugate(khplpl))*d2)+abs(A2)^2*khgplpl*exp(I*(khgplpl-conjugate(khgplpl))*d2)+abs(B3)^2*k0mimi*exp(I*(k0mimi-conjugate(k0mimi))*d2)+abs(B3)^2*kg1mimi*exp(I*(kg1mimi-conjugate(kg1mimi))*d2)+abs(B4)^2*khmimi*exp(I*(khmimi-conjugate(khmimi))*d2)+abs(B4)^2*khgmimi*exp(I*(khgmimi-conjugate(khgmimi))*d2)+abs(B1+A3)^2*k0*exp(I*(k0-conjugate(k0))*d2)+abs(B1-A3)^2*kg1*exp(I*(kg1-conjugate(kg1))*d2)+abs(A4-B2)^2*kh*exp(I*(kh-conjugate(kh))*d2)+abs(A4+B2)^2*khg*exp(I*(khg-conjugate(khg))*d2)+conjugate(A1)*B3*k0mimi*exp(I*(k0mimi-conjugate(k0plpl))*d2)+A1*conjugate(B3)*k0plpl*exp(I*(k0plpl-conjugate(k0mimi))*d2)+conjugate(A1)*(B1+A3)*k0*exp(I*(k0-conjugate(k0plpl))*d2)+A1*conjugate(B1+A3)*k0plpl*exp(I*(k0plpl-conjugate(k0))*d2)+conjugate(B3)*(B1+A3)*k0*exp(I*(k0-conjugate(k0mimi))*d2)+B3*conjugate(B1+A3)*k0mimi*exp(I*(k0mimi-conjugate(k0))*d2)-conjugate(A1)*B3*kg1mimi*exp(I*(kg1mimi-conjugate(kg1plpl))*d2)-A1*conjugate(B3)*kg1plpl*exp(I*(kg1plpl-conjugate(kg1mimi))*d2)+conjugate(A1)*(A3-B1)*kg1*exp(I*(kg1-conjugate(kg1plpl))*d2)+A1*conjugate(A3-B1)*kg1plpl*exp(I*(kg1plpl-conjugate(kg1))*d2)+conjugate(B3)*(B1-A3)*kg1*exp(I*(kg1-conjugate(kg1mimi))*d2)+B3*conjugate(B1-A3)*kg1mimi*exp(I*(kg1mimi-conjugate(kg1))*d2)-conjugate(A2)*B4*khmimi*exp(I*(khmimi-conjugate(khplpl))*d2)-A2*conjugate(B4)*khplpl*exp(I*(khplpl-conjugate(khmimi))*d2)+conjugate(A2)*(B2-A4)*kh*exp(I*(kh-conjugate(khplpl))*d2)+A2*conjugate(B2-A4)*khplpl*exp(I*(khplpl-conjugate(kh))*d2)+conjugate(B4)*(A4-B2)*kh*exp(I*(kh-conjugate(khmimi))*d2)+B4*conjugate(A4-B2)*khmimi*exp(I*(khmimi-conjugate(kh))*d2)+conjugate(A2)*B4*khgmimi*exp(I*(khgmimi-conjugate(khgplpl))*d2)+A2*conjugate(B4)*khgplpl*exp(I*(khgplpl-conjugate(khgmimi))*d2)-conjugate(A2)*(A4+B2)*khg*exp(I*(khg-conjugate(khgplpl))*d2)-A2*conjugate(A4+B2)*khgplpl*exp(I*(khgplpl-conjugate(khg))*d2)-conjugate(B4)*(A4+B2)*khg*exp(I*(khg-conjugate(khgmimi))*d2)-B4*conjugate(A4+B2)*khgmimi*exp(I*(khgmimi-conjugate(khg))*d2)):


T2:=1/4*Re(abs(A1)^2*k0plpl*exp(I*(k0plpl-conjugate(k0plpl))*d2)+abs(A1)^2*kg2plpl*exp(I*(kg2plpl-conjugate(kg2plpl))*d2)+abs(A2)^2*khplpl*exp(I*(khplpl-conjugate(khplpl))*d2)+abs(A2)^2*khgplpl*exp(I*(khgplpl-conjugate(khgplpl))*d2)+abs(B3)^2*k0mimi*exp(I*(k0mimi-conjugate(k0mimi))*d2)+abs(B3)^2*kg2mimi*exp(I*(kg2mimi-conjugate(kg2mimi))*d2)+abs(B4)^2*khmimi*exp(I*(khmimi-conjugate(khmimi))*d2)+abs(B4)^2*khgmimi*exp(I*(khgmimi-conjugate(khgmimi))*d2)+abs(B1+A3)^2*k0*exp(I*(k0-conjugate(k0))*d2)+abs(B1-A3)^2*kg2*exp(I*(kg2-conjugate(kg2))*d2)+abs(A4-B2)^2*kh*exp(I*(kh-conjugate(kh))*d2)+abs(A4+B2)^2*khg*exp(I*(khg-conjugate(khg))*d2)+conjugate(A1)*B3*k0mimi*exp(I*(k0mimi-conjugate(k0plpl))*d2)+A1*conjugate(B3)*k0plpl*exp(I*(k0plpl-conjugate(k0mimi))*d2)+conjugate(A1)*(B1+A3)*k0*exp(I*(k0-conjugate(k0plpl))*d2)+A1*conjugate(B1+A3)*k0plpl*exp(I*(k0plpl-conjugate(k0))*d2)+conjugate(B3)*(B1+A3)*k0*exp(I*(k0-conjugate(k0mimi))*d2)+B3*conjugate(B1+A3)*k0mimi*exp(I*(k0mimi-conjugate(k0))*d2)-conjugate(A1)*B3*kg2mimi*exp(I*(kg2mimi-conjugate(kg2plpl))*d2)-A1*conjugate(B3)*kg2plpl*exp(I*(kg2plpl-conjugate(kg2mimi))*d2)+conjugate(A1)*(A3-B1)*kg2*exp(I*(kg2-conjugate(kg2plpl))*d2)+A1*conjugate(A3-B1)*kg2plpl*exp(I*(kg2plpl-conjugate(kg2))*d2)+conjugate(B3)*(B1-A3)*kg2*exp(I*(kg2-conjugate(kg2mimi))*d2)+B3*conjugate(B1-A3)*kg2mimi*exp(I*(kg2mimi-conjugate(kg2))*d2)-conjugate(A2)*B4*khmimi*exp(I*(khmimi-conjugate(khplpl))*d2)-A2*conjugate(B4)*khplpl*exp(I*(khplpl-conjugate(khmimi))*d2)+conjugate(A2)*(B2-A4)*kh*exp(I*(kh-conjugate(khplpl))*d2)+A2*conjugate(B2-A4)*khplpl*exp(I*(khplpl-conjugate(kh))*d2)+conjugate(B4)*(A4-B2)*kh*exp(I*(kh-conjugate(khmimi))*d2)+B4*conjugate(A4-B2)*khmimi*exp(I*(khmimi-conjugate(kh))*d2)+conjugate(A2)*B4*khgmimi*exp(I*(khgmimi-conjugate(khgplpl))*d2)+A2*conjugate(B4)*khgplpl*exp(I*(khgplpl-conjugate(khgmimi))*d2)-conjugate(A2)*(A4+B2)*khg*exp(I*(khg-conjugate(khgplpl))*d2)-A2*conjugate(A4+B2)*khgplpl*exp(I*(khgplpl-conjugate(khg))*d2)-conjugate(B4)*(A4+B2)*khg*exp(I*(khg-conjugate(khgmimi))*d2)-B4*conjugate(A4+B2)*khgmimi*exp(I*(khgmimi-conjugate(khg))*d2)):

#evalf(Int(k*sin(x)*T1,x=0..Pi/2,y=-Pi/6..-Pi/6+afa))

indets(T1, name)

{x, y}

(1)

T1:=1/4*(abs(A1)^2*k0plpl*exp(I*(k0plpl-conjugate(k0plpl))*d2)+abs(A1)^2*kg1plpl*exp(I*(kg1plpl-conjugate(kg1plpl))*d2)+abs(A2)^2*khplpl*exp(I*(khplpl-conjugate(khplpl))*d2)+abs(A2)^2*khgplpl*exp(I*(khgplpl-conjugate(khgplpl))*d2)+abs(B3)^2*k0mimi*exp(I*(k0mimi-conjugate(k0mimi))*d2)+abs(B3)^2*kg1mimi*exp(I*(kg1mimi-conjugate(kg1mimi))*d2)+abs(B4)^2*khmimi*exp(I*(khmimi-conjugate(khmimi))*d2)+abs(B4)^2*khgmimi*exp(I*(khgmimi-conjugate(khgmimi))*d2)+abs(B1+A3)^2*k0*exp(I*(k0-conjugate(k0))*d2)+abs(B1-A3)^2*kg1*exp(I*(kg1-conjugate(kg1))*d2)+abs(A4-B2)^2*kh*exp(I*(kh-conjugate(kh))*d2)+abs(A4+B2)^2*khg*exp(I*(khg-conjugate(khg))*d2)+conjugate(A1)*B3*k0mimi*exp(I*(k0mimi-conjugate(k0plpl))*d2)+A1*conjugate(B3)*k0plpl*exp(I*(k0plpl-conjugate(k0mimi))*d2)+conjugate(A1)*(B1+A3)*k0*exp(I*(k0-conjugate(k0plpl))*d2)+A1*conjugate(B1+A3)*k0plpl*exp(I*(k0plpl-conjugate(k0))*d2)+conjugate(B3)*(B1+A3)*k0*exp(I*(k0-conjugate(k0mimi))*d2)+B3*conjugate(B1+A3)*k0mimi*exp(I*(k0mimi-conjugate(k0))*d2)-conjugate(A1)*B3*kg1mimi*exp(I*(kg1mimi-conjugate(kg1plpl))*d2)-A1*conjugate(B3)*kg1plpl*exp(I*(kg1plpl-conjugate(kg1mimi))*d2)+conjugate(A1)*(A3-B1)*kg1*exp(I*(kg1-conjugate(kg1plpl))*d2)+A1*conjugate(A3-B1)*kg1plpl*exp(I*(kg1plpl-conjugate(kg1))*d2)+conjugate(B3)*(B1-A3)*kg1*exp(I*(kg1-conjugate(kg1mimi))*d2)+B3*conjugate(B1-A3)*kg1mimi*exp(I*(kg1mimi-conjugate(kg1))*d2)-conjugate(A2)*B4*khmimi*exp(I*(khmimi-conjugate(khplpl))*d2)-A2*conjugate(B4)*khplpl*exp(I*(khplpl-conjugate(khmimi))*d2)+conjugate(A2)*(B2-A4)*kh*exp(I*(kh-conjugate(khplpl))*d2)+A2*conjugate(B2-A4)*khplpl*exp(I*(khplpl-conjugate(kh))*d2)+conjugate(B4)*(A4-B2)*kh*exp(I*(kh-conjugate(khmimi))*d2)+B4*conjugate(A4-B2)*khmimi*exp(I*(khmimi-conjugate(kh))*d2)+conjugate(A2)*B4*khgmimi*exp(I*(khgmimi-conjugate(khgplpl))*d2)+A2*conjugate(B4)*khgplpl*exp(I*(khgplpl-conjugate(khgmimi))*d2)-conjugate(A2)*(A4+B2)*khg*exp(I*(khg-conjugate(khgplpl))*d2)-A2*conjugate(A4+B2)*khgplpl*exp(I*(khgplpl-conjugate(khg))*d2)-conjugate(B4)*(A4+B2)*khg*exp(I*(khg-conjugate(khgmimi))*d2)-B4*conjugate(A4+B2)*khgmimi*exp(I*(khgmimi-conjugate(khg))*d2)):

indets(k*sin(x)*T1, name)

{x, y}

(2)

# evalf(Int(k*sin(x)*T1,x=0..Pi/2,y=-Pi/6..-Pi/6+afa))

Warning,  computation interrupted

 

# the integrand seems almost invariant regarding y:

plot3d(Re(k*sin(x)*T1),x=0..Pi/2,y=-Pi/6..-Pi/6+afa, labels=[x, y, z])

 

 

method = _d01akc

--

for finite interval of integration, oscillating integrands; uses adaptive Gauss 30-point and Kronrod 61-point rules.

 

 

 

# as the integrand is x-oscillating and y-constant (or almost), a suggestion is to do this

i1 := evalf(Int(eval(k*sin(x)*Re(T1), y=-Pi/6. ), x=0..Pi/2, method=_d01akc));
i2 := evalf(Int(eval(k*sin(x)*Re(T1), y=-Pi/12.), x=0..Pi/2, method=_d01akc));
i3 := evalf(Int(eval(k*sin(x)*Re(T1), y=     0.), x=0..Pi/2, method=_d01akc));
i4 := evalf(Int(eval(k*sin(x)*Re(T1), y=+Pi/12.), x=0..Pi/2, method=_d01akc));
i5 := evalf(Int(eval(k*sin(x)*Re(T1), y=+Pi/6. ), x=0..Pi/2, method=_d01akc));

0.1179159365e20

 

0.1171354144e20

 

0.1171354144e20

 

0.1171354144e20

 

0.1171354144e20

(3)

# here is an evaluation of the integral

CurveFitting:-Spline(Pi*~[seq(-1/6..1/6, 1/12)], [seq(i||k, k=1..5)], t, degree=1):
approx_1 := int(%, t=-Pi/6..Pi/6);

0.1227660893e20

(4)

# improvements: T1 depends on y only in a thin layer -Pi/6 .. -Pi/6.*(0.94)
evalf(Int(eval(k*sin(x)*Re(T1), y=-Pi/6.*(0.94)), x=0..Pi/2, method=_d01akc));

0.1171354144e20

(5)

# thus this sevond approximation of the integral

CurveFitting:-Spline([-Pi/6, -Pi/6*0.94, Pi/6], [i1, i2, i2], t, degree=1):
approx_2 := int(%, t=-Pi/6..Pi/6);

0.1226761795e20

(6)

 


 

Download text_program_mmcdara.mw

 

As dualxisplot can't do the job (as Carl said it is not intended to this), maybe yu will be interested in that:

KindOf.mw

solve(...) tells you that the solutions verify a fourth order polynomial equation with inderterminate _Z (help page undocumentednames).
To get the explicit expressions of these solutions use allvalues(%)

restart:
A := Vector(3, symbol=a):
B := Vector(3, symbol=b):
LinearAlgebra:-KroneckerProduct(A, B^+)

If you read carefully the ArrayInterpolation help page you see that the plot command is matrixplot.
The object ArrayInterpolation build is an array, not a function (I use MAPLE 2015, maybe it's different in newer versions, so please check this) so this is not possible to use contourplot, contourplot3d nor even plot3d.

Here is a (very costly) workaround based on bilinear interpolation (aka Q1 interpolation in finite elements) on square cells.
If you are only interesting in drawing a single level curve (of "value" C), one can buid a more astute algorithm which interpolates only over square cells whose the range of nodal values contains "C".
 

restart;

with(CurveFitting):

with(plots):

alpha := Array([seq(0..10,1.)]):
beta  := Array([seq(0..10,1.)]);
PA := numelems(alpha);
PB := numelems(beta);

beta := Vector(4, {(1) = ` 1 .. 11 `*Array, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

11

 

11

(1)

NN := Matrix(PA, PB, (i, j) -> cos(alpha[i]/3+beta[j]/3)*sin(alpha[i]*beta[j]/9))

NN := Vector(4, {(1) = ` 11 x 11 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(2)

A  := Matrix(PA, PB, (i, j) -> alpha[i]):
B  := Matrix(PA, PB, (i, j) -> beta [j]):
AB := ArrayTools[Concatenate](3, A, B);
F  := ArrayInterpolation([beta, alpha], NN, AB, method=nearest):

matrixplot(F);

AB := Vector(4, {(1) = ` 1..11 x 1..11 x 1..2 `*Array, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

 

phi := (a, b) -> [(1-a)*(1-b), (1-a)*b, a*b, a*(1-b)];

proc (a, b) options operator, arrow; [(1-a)*(1-b), (1-a)*b, a*b, a*(1-b)] end proc

(3)

FV := (a, b) -> piecewise(
        seq(
          seq(
            op(
              [
                verify(a, alpha[i]..alpha[i+1], interval)
                and
                verify(b, beta[j]..beta[j+1], interval),
                add(phi((a-alpha[i])/hA, (b-beta[j])/hB) *~ [ NN[i, j], NN[i, j+1], NN[i+1, j+1], NN[i+1, j]])
              ]
            ),
            j=1..PB-1
          ),
          i=1..PA-1
        )
      ):

plot3d('FV(a, b)', a=0..10, b=0..10):

contourplot3d('FV(a, b)', a=0..10, b=0..10, contours=[0.2], filledregions=true)

 

contourplot('FV(a, b)', a=0..10, b=0..10, contours=[0.2], view=[0..10, 0..10])

 

 

 

Download Interp_2D.mw

I you use Maple 2019 or a more recent version, I advice you to look the kriging interpolator (If I'm not mistaken it is a member of the CurveFitting package).
If my memory doesnt fool me I belive that something about 2D interpolation has been published in the post section of Mapleprimes about a year ago

 

@ijuptilk 

Well, in this case the things are rather simple.
The idea is to decompose the a priori plotting domain into subdomains where the function is continuous.
For instance
 

restart:

f := y/(x+1) + 1/(y^2-y)

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

(1)

eps := 1e-3:

x_dom := [-4, 3]:
y_dom := [-3, 4]:

s := discont(f, x):
if s <> NULL then
  x_dom := sort(
             [
               x_dom[],
               map(z ->op([z-eps, z+eps]), convert(s, list))[]
             ]
           )
end if;

s := discont(f, y):
if s <> NULL then
  y_dom := sort(
             [
               y_dom[],
               map(z ->op([z-eps, z+eps]), convert(s, list))[]
             ]
           )
end if;

x_cont := [seq(x_dom[k]..x_dom[k+1], k in [seq](1..numelems(x_dom)-1, 2))];
y_cont := [seq(y_dom[k]..y_dom[k+1], k in [seq](1..numelems(y_dom)-1, 2))];


r := rand(0. .. 1.):
plots:-display(
  seq(seq(plot3d(f, x=a, y=b, style=surface, color=ColorTools:-Color([r(), r(), r()])), a in x_cont), b in y_cont),
  view=[default, default, -10..10]
)

[-4, -1.001, -.999, 3]

 

[-3, -0.1e-2, 0.1e-2, .999, 1.001, 4]

 

[-4 .. -1.001, -.999 .. 3]

 

[-3 .. -0.1e-2, 0.1e-2 .. .999, 1.001 .. 4]

 

 

 


Download plot3d_discont.mw


Other types of dicontinuities :
(II)  When the discontinuity is located on a curve of equation f=g(y) (or y=g(x)) one can proceed the same way.
(III) Things are more tetchy if the discontuity line is of the form g(x, y)=0.

A simple procedure can be devisied for "type I" discontinuities (as in your f function).
This same procedure might probably be modifed tio handle "type II" discontinuities (x=g(y) or y=g(x))

... in an infinite number of ways.

(which should not be that surprising given that z depends only on 8 quantities, and even only 4 with their conjugate).

Without any further information about the structure of the matrix, any non singular 3x3 matrix can be modified for its determinant be equal to z.
Here is a trivial example
Download MatricesWhoseDeterminantEqual_Z.mw

Many other simple examples can be found; for instance
MatricesWhoseDeterminantEqual_z_updated.mw
MatricesWhoseDeterminantEqual_z_reupdated.mw

Personal experience.

DeepLearning is poorly documented looks more like an attempt to give a "deep learning" varnish to Maple than to really deliver a package aimed to solve practical problems.
The Classifier is, at best, acceptable, the Regressor is far from classical standards (in particular there is no analyzing tools of the solution you got).
From having compared Maple/DeepLearning to R/Keras* several times, my conclusion is that if you want to solve something other than toy problems, forget Maple.

Of course, this is my own opinion that some may disagree with... so make your own.

* Keras is a high level API to TensorFlow (that you can use in Python and R)
https://www.rstudio.com/blog/keras-for-r/

... which differs from Carl's 

You will find in the attached file:

  • Two procedures named check_1 and check_2 which both return messages of the form "rows number n of M1 and M2 are identical" (or different).
    Not exactly what you ask for but it would be obvious to deduce if M1 = M2 or not.
     
  • A third procedure named check_3 which compares M1 and M2 row by row "uo to a given precision" in case M1 and M2 are matrices of floats.

Where my understanding differs from Carl's is that

  • compare row Rn of M1 to row Rn of M2
  • instead of checking if "for every row R1 of A there exists a row R2 of such that R1 "equals" R2 ...i"

if I'm wrong, forget my answer.
check.mw



The function you want to integrate is nothing but the Probability Density Function (PDF) of a Beta(n-x, x+1) random variable.
So the result (faster and less memory consuming than your approach with Carl's solution)
 

restart:

with(Statistics):

PDF('Beta'(n-x, x+1), t)

piecewise(t < 0, 0, t < 1, t^(-1+n-x)*(1-t)^x/Beta(n-x, x+1), 0)

(1)

f := (x, n, p) -> eval(CDF('Beta'(n-x, x+1), t), t=p):
plot(f(x, 10, 1/2), x=0..10)

 

 


 

Download ProbabilityPointOfView.mw

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