# Multivariable numerical-symbolic integration

I need to find the symbolic expression of a multivariable integral. The real integrand contains 11 variables (I have to integrate "only" four of them) so example code is simplified here:

x, y, z = var('x','y','z')
integrand = sin(x)*cos(y)*z
result = integrand.integral(x, 0, pi).integral(z, 0, 1)


It works just fine generally but given the complexity of the integral no algorithm gives me a result and it returns just:

integrate(integrate(sin(x)*cos(y), x, 0, pi), z, 0, 1)


Is there a built in way to numerically evaluate an integral on some variables that returns symbolic result with the numerical integration found coefficients (or constants)? Namely, in this case:

cos(y)


EDIT1: I managed to get a very very simple example of the function I need:

def numsym_integral( f, variable, a, b, points = 1000):
dx = (b-a)/points
result = 0
for i in range(0,points):
result = result + f.substitute({variable:a + i*dx}) * dx
return result


It still doesn't work though in my specific case because of the complexity of the expression. I get the number of terms gets multiplied by "points" per integration, leading to something too big. Could something built in be more efficient?

edit retag close merge delete

Don't you mean :

integrate(integrate(sin(x)*cos(y), y, c, d), x, a, b)


?

In so, the numerical equivalent would be

numerical_integrate(lambda u:numerical_integrate(lambda v:sin(u)*cos(v), c, d), a, b)


Illustration :

sage: var("x, y, a, b, c, d")
(x, y, a, b, c, d)
sage: integrate(integrate(sin(x)*cos(y), y, c, d), x, a, b)
-(cos(a) - cos(b))*(sin(c) - sin(d))
sage: integrate(integrate(sin(x)*cos(y), y, c, d), x, a, b).subs({a:0, b:pi/2, c:0, d:pi/2})
1
sage: numerical_integral(lambda u:numerical_integral(lambda v:sin(u)*cos(v), 0, pi/2), 0, pi/2)
(1.0, 1.1102230246251565e-14)


No, I didn't mean that. My bad though, I've not explained clearly enough, I just edited the example code to make it clearer and added a (not working in my case) not built in solution to better explain what I mean.

About your second edit : What you have programmed is a simple version (constant step) of trapeze integration. numerical_integral has lots of further sophistication allowing integration in cases where the trapeze rule fails miserably (singularities, infinite bounds, etc...).

I fail to understand the problem you are trying to solve...

From what I've understood on the docs numerical_integral forces you to set the value of every other non integrating variable in your integrand. I don't need a number, I possibly need an expression.

From what I've understood on the docs numerical_integral forces you to set the value of every other non integrating variable in your integrand. I don't need a number, I possibly need an expression.

That's why you have to nest numerical_integral calls : the first argument must be numerically evaluable at run time (i. e. no symbolic variables).

In other words, the value of a numerical_integral call is a numeric value, not an expression...

Sort by » oldest newest most voted

It may happen that integrate(integrate(f(x,y), x, a, b), y, c, d)fails but that integrate(integrate(f(x, y),y, c, d), x, a, b) succeeds.

IF your integrand $f(x_1,\cdots,x_{11})$ complies with the preconditions of the Fubini-Tonelli theorem, then $\displaystyle\int f(x_1,\cdots,x_{11})\, dx_i\cdots dx_j$ does not depend of the order of the integration variables order $x_1,\cdots,x_j$.

You might try to search the ordered subsets $x_i, \cdots, x_j$ of your integration variables $x_1,\cdots,x_4$ such as $\displaystyle\int f(x_1,\cdots,x_{11})\, dx_i\cdots dx_j$ has an explicit and convenient form $F(x_1,\cdots x_{11})$. There are $\displaystyle\sum_{i=0}^4 i!\,=\,34$ such subsets.

Programatically using the results of such a "brute force" is made possible by the fact that a "non-integrated" result R will have R.operator() being integrate.

Similarly, you may try to loop over the various available integration algorithms.

Using this form, the numerical integration of $F$ over the remaining integration variables is your solution.

more

I didn't know the integration variables order mattered in sage, thanks! My integrand complies with Fubini's theorem. I'm giving your suggestion a try now but it's still really slow and it could solve my use case but not similar ones in the future. If it won't work I'll try to simplify my integral by substitution to avoid composite functions in the denominator which is the reason why the number of terms got multiplied by the number of points after every integration with earlier method.