Processing math: 100%

First time here? Check out the FAQ!

Ask Your Question
0

Multivariable numerical-symbolic integration

asked 4 years ago

MateCheck gravatar image

updated 4 years ago

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?

Preview: (hide)

Comments

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)[0], a, b)[0]

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], 0, pi/2)
(1.0, 1.1102230246251565e-14)
Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 4 years ago )

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.

MateCheck gravatar imageMateCheck ( 4 years ago )

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

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 4 years ago )

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.

MateCheck gravatar imageMateCheck ( 4 years ago )

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

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 4 years ago )

1 Answer

Sort by » oldest newest most voted
0

answered 4 years ago

Emmanuel Charpentier gravatar image

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(x1,,x11) complies with the preconditions of the Fubini-Tonelli theorem, then f(x1,,x11)dxidxj does not depend of the order of the integration variables order x1,,xj.

You might try to search the ordered subsets xi,,xj of your integration variables x1,,x4 such as f(x1,,x11)dxidxj has an explicit and convenient form F(x1,x11). There are 4i=0i!=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.

Preview: (hide)
link

Comments

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.

MateCheck gravatar imageMateCheck ( 4 years ago )

Nope, still not working.

MateCheck gravatar imageMateCheck ( 4 years ago )

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 4 years ago

Seen: 1,196 times

Last updated: Feb 14 '21