@9009134 

Your constants are not assigned. After giving them numeric values, you may try to compute the integrals, but the integrands being very "odd", you may need to fine-tune int as above.

@9009134 

Try

Lps:=expand(L[p, s](r, theta, phi)*F[p,1](r, theta));
kk:=int(Lps, [theta=0..2*Pi,r=0.5..1], numeric, epsilon=1e-7);
# kk:=int( L[p, s](r, theta, phi)*F[p,1](r, theta), [theta=0..2*Pi,r=0.5..1], numeric);

@9009134 

You keep modifying the worksheet and have eliminated my code.

Now you have functions of 3 variables and want to integrate numerically wrt 2 of them.
It wil not work of course.

@Pascal4QM 

Yes, I know how it works, but I wanted a reference (I was asked by a colleague). It seems to be not documented.

@9009134 

Replace the last part with

for p from 1 to 3 do;
for s from 1 to 3 do;

K1:=int(L[p, s](r, theta, phi)*f[p,1](r, theta, phi), phi=0..2*Pi):
K1e:=expand(K1);
K1e := `if`(type(K1e,`+`), [op(K1e)], [K1e]);
kk:=add(op(1,J)*int( J/op(1,J), [theta=0..2*Pi,r=0.4..1], numeric), J=K1e);

kkk:=int((f[p,1](r, theta, phi))^2, theta=0..2*Pi, phi=0..2*Pi, r=0.4..1);

k[p,s]:=kk/kkk:

evalf(%);
print([p,s]=%)
od;
od;

Please note that this is only a quick workaround for a Maple weakness. It may fail if you change the expressions too much.

@9009134 

You have ignored my suggestion.

Float(undefined)ok.mw

I don't understand why not simply plot the solutions of the system. In your case it is a finite set, but of course it could be an algebraic variety. I don't find that P1 + a*P2 is very intuitive and it's far from the ideal <P1,P2>.

For P1,P2,P3 the ideal is <1> (you can use Groebner:-Basis), so, no solutions for the system.

@acer 

Nice post. Do you have such a tricky int example at hand?

@ecterrab 

Thank you for the answer. I know about op(procname). A concrete example can be seen with showstat(log).
The question is about documentation. It should be under ?functions but it is probably elsewhere. Or, is it not documented?

@Abdoulaye 

Yes, it's possible to use a Riemann sum (or a numeric integration method such as Simpson's) but I don't see why, because an interpolation is anyway needed after that (otherwise the complexity will grow exponentially and it will be impossible to compute e.g. F4 or even F3).

@Abdoulaye 

You may try the following:

- Fix a (uniform) grid X = [x0,...,xM], Y=[y0,...,yN] 

- After computing Fn at the points of the grid (xi,yj), replace Fn with its SplineInterpolation SFn. This way SFn can be integrated, producing (probably) useful approximations for the next iteration. Do the same for Un.

N.B. I am sure that there are better methods for your (undisclosed) original problem.

@MDD 

So, you want "disjoint" monomials.

restart;
vars:=[x,y,z]:
m1,m2,m3:= (a-b)*x*y*z^2, c*x^3, 3*b*y*z^2:
DisjMon:= (m1,m2,vars::list(name)) -> not depends(gcd(m1,m2), vars):
DisjMon(m1,m2,vars), DisjMon(m2,m3,vars), DisjMon(m1,m3,vars);

                                  false, true, false

@MDD 

EqMon tests if two terms have the same monomials, with possible distinct coefficients.
Define exactly what you want to test. Your example is not relevant and you have a typo.
 

@9009134 

Here is the the [2,3] integral.

Digits:=15:
ti:=time[real]():
F23:=(3/4*(501230769273.847*cos(theta)^2*r^2+5012307692738.48*r*cos(theta)+12530769231846.2)*sin(theta)*r^4*cos(theta)*cos(3*phi)/(5+r*cos(theta))^4-3/4*(.835384615384617e-1*cos(theta)^2*r+278461538504.616*cos(theta)^2*r^3+.139230769230770*cos(theta)+2784615385046.16*cos(theta)*r^2+6961538462615.36*r)*r^2*sin(theta)*cos(3*phi)/(5+r*cos(theta))^3-.710076923076924*sin(theta)*r^2*cos(theta)*cos(3*phi)/(5+r*cos(theta))^2-.229730769230770*r*sin(theta)*cos(3*phi)/(5+r*cos(theta))+.187961538461539*r^3*sin(theta)*cos(3*phi)/(5+r*cos(theta))^3+.187961538461539*sin(theta)*r^4*cos(theta)*cos(3*phi)/(5+r*cos(theta))^4-.417692307692310*r*sin(theta)*cos(3*phi)/(5+r*cos(theta))^2-.835384615384620*sin(theta)*r^2*cos(theta)*cos(3*phi)/(5+r*cos(theta))^3)*r^4*cos(theta)*cos(3*phi):
temp:=int(F23, phi=0..2*Pi):
L:=[op(expand(temp))]:
add(op(1,J)*int( J/op(1,J), [theta=0..2*Pi,r=2/10..1], numeric), J=L);
time=time[real]()-ti;
#                        1.09206001690410*10^(-6)
#                          time = 0.427

N.B. It would be nice to isolate yourself the problem before posting.

@acer 

I thought that these local variables have another status.
So, they differ from the regular ones, only by being "quasi" inaccessible.

@Carl Love  My question was not related to your original answer.

Thank you both for the clarification.