mmcdara

3815 Reputation

17 Badges

6 years, 118 days

MaplePrimes Activity


These are answers submitted by mmcdara

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

Only little adjustments
fixed.mw

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