mmcdara

3815 Reputation

17 Badges

6 years, 122 days

MaplePrimes Activity


These are answers submitted by mmcdara

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

You should begin by read the help pages and try to do something by yourself.

restart:

with(LinearAlgebra):

A := Matrix(2, 2, [0, a, b, 0]);
lambda, V := Eigenvectors(A);

A := Matrix(2, 2, {(1, 1) = 0, (1, 2) = a, (2, 1) = b, (2, 2) = 0})

 

lambda, V := Vector(2, {(1) = (a*b)^(1/2), (2) = -(a*b)^(1/2)}), Matrix(2, 2, {(1, 1) = a/(a*b)^(1/2), (1, 2) = -a/(a*b)^(1/2), (2, 1) = 1, (2, 2) = 1})

(1)

__lambda, __V := Eigenvectors(Transpose(A));

__lambda, __V := Vector(2, {(1) = (a*b)^(1/2), (2) = -(a*b)^(1/2)}), Matrix(2, 2, {(1, 1) = b/(a*b)^(1/2), (1, 2) = -b/(a*b)^(1/2), (2, 1) = 1, (2, 2) = 1})

(2)

 

Download EigenVectors.mw

The attached file focuses only of "fitting" through Linear Regression.
As you see Statistics:-LinearFit doesn't accept Units, but a simple workaround can be found.

Hope it will help

restart:

 

Fit with Statistics:-LinearFit

 

with(Statistics):

# Example from LinearFit help page

X := Vector([1, 2, 3, 4, 5, 6]):
Y := Vector([2, 3, 4, 3.5, 5.8, 7]):
LinearFit([1, t, t^2], X, Y, t);

HFloat(1.960000000000003)+HFloat(0.16499999999999876)*t+HFloat(0.11071428571428583)*t^2

(1)

# Naive use of LinearFit for dimension quantities doen't work

N  := numelems(X):
XU := X *~ Units:-Unit('s'):
YU := Y *~ Units:-Unit('m'):
LinearFit([1, t, t^2], XU, YU, t);

Error, (in Statistics:-LinearFit) unable to store 'Units:-Unit(s)' when datatype=float[8]

 

# A workaround: write and solve the normal equations to obtain dimension coefficients
#
# Note that "t" is a number (dimensionless quantity).

B := t -> `<|>`(1, t, t^2):

H := Matrix(N, numelems(B(t)), (i, j) -> B(XU[i])[j]):
C := (H^+ . H)^(-1) . (H^+ . YU);  # dimension coefficients
Model := simplify(B(t*Units:-Unit('s')) . C);                 # dimension model

C := Vector(3, {(1) = 1.9600000*Units:-Unit(m), (2) = .1650000*Units:-Unit(m)/Units:-Unit(s), (3) = .11071429*Units:-Unit(m)/Units:-Unit(s)^2})

 

(1.96+.16500*t+.1107142900*t^2)*Units:-Unit('m')

(2)

 

Download LinearFitWithUnits.mw

First definition of S is nothing but a variant of Tomleslie's solution

Second definition seems closer to what you asked (the xi's are not elements of a set but of a list, which seems safer given that a set is automatically ordered... see at the end of my code).

restart

S := (x, h, N) -> x +~ h *~ [$0..N]:

S(x__0, b, 4);
S(x__0, b, 0);

[x__0, x__0+b, x__0+2*b, x__0+3*b, x__0+4*b]

 

[x__0]

(1)

A := N -> cat~('x__', [$0..N]):
S := N -> ListTools:-PartialSums(A(N)) +~ b *~ [$0..N]:

'A(4)' = A(4);
'S(4)' = S(4);

A(4) = [x__0, x__1, x__2, x__3, x__4]

 

S(4) = [x__0, x__0+x__1+b, x__0+x__1+x__2+2*b, x__0+x__1+x__2+x__3+3*b, x__0+x__1+x__2+x__3+x__4+4*b]

(2)

# compare this

A := cat~('x__', [$0..20]); # A is a list

# to this

A :=cat~('x__', {$0..20}); # A is a set

[x__0, x__1, x__2, x__3, x__4, x__5, x__6, x__7, x__8, x__9, x__10, x__11, x__12, x__13, x__14, x__15, x__16, x__17, x__18, x__19, x__20]

 

{x__0, x__1, x__10, x__11, x__12, x__13, x__14, x__15, x__16, x__17, x__18, x__19, x__2, x__20, x__3, x__4, x__5, x__6, x__7, x__8, x__9}

(3)

 

Download MaybeThis.mw

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