## 3815 Reputation

6 years, 122 days

## It is highly likely that the formal inte...

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)

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) (2)

Let us to compute this "simple integral

 > I_p := Int(exp*p^2*sin(alpha), 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 ) (4)
 > # The formal integration alreaqy fayls on this even simpler expression infolevel[int] := 4: p^2 / sqrt(b); int(%, p=0..1) 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 (5)
 > I_theta(0, 1, I, _CubaVegas);  # interrupted
 > I_theta(0, 1, I, _MonteCarlo); # fails (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");  (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']);  (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+Boot*Quantile(Normal(0, 1), 0.025) .. Boot+Boot*Quantile(Normal(0, 1), 0.975) (9)
 >

## A possibility >...

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

 > > > >  (1)
 > > > > >  (2)
 >  (3)
 >  (4)
 >    (5)
 >  (6)
 >  (7)
 >  (8)
 > >  (9)
 > >  (10)
 > >  (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, g);  (12)
 > deg := 2: g   := GenSys(deg); sol := solve(g, g);  (13)
 > vals := allvalues(eval(A, sol[])): map(u -> [remove(has, sol[], A)[], A=u], [vals]): print~(%):  (14)
 > Improvement (automatic treatment of equations of the form A[m]=RootOf(...))

```deg := 10: # for instance
g   := GenSys(deg):
sol := solve(g, g):

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 i...

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) (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) (2)
 > # evalf(Int(k*sin(x)*T1,x=0..Pi/2,y=-Pi/6..-Pi/6+afa))
 > # 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));     (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); (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)); (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); (6)
 >

## For what it worths.......

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 ...

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(%)

## Like that...

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

## If you read carefully the ArrayInterpola...

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);   (1)
 > NN := Matrix(PA, PB, (i, j) -> cos(alpha[i]/3+beta[j]/3)*sin(alpha[i]*beta[j]/9)) (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);  > phi := (a, b) -> [(1-a)*(1-b), (1-a)*b, a*b, a*(1-b)]; (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]) >

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 th...

@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) (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] )     >

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))

## YES......

... 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

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

## Personal experience.DeepLearning is...

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/

## My understanding......

... 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

## A probability point of view...

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

## You should begin by read the help pages ...

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);  (1)
 > __lambda, __V := Eigenvectors(Transpose(A)); (2)
 >

## The attached file focuses only of "...

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); (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);
 > # 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  (2)
 >

## First definition of S is nothing but a v...

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);  (1)
 > A := N -> cat~('x__', [\$0..N]): S := N -> ListTools:-PartialSums(A(N)) +~ b *~ [\$0..N]: 'A(4)' = A(4); 'S(4)' = S(4);  (2)
 > # compare this A := cat~('x__', [\$0..20]); # A is a list # to this A :=cat~('x__', {\$0..20}); # A is a set  (3)
 >

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