# How to do multiple integral numerically?

I cannot calculate the following type of integrals numerically :

$\int_0^1 dy \int_0^y f(x) dx$

$f(x)$ can be a complicated function. The problem is due to the fact that the upper limit of one of the integrals is not a number. Can anyone suggest how to compute it. Many thanks!

Edit : Actually my aim is to calculate a residue. I didn't want to bother you with these long expressions. Here is what I wish to compute (numerically at least).

$$\text{Res}\left(\frac{\int_0^x \text{P}(w) S(w) \, dw}{\text{P}(x)},\left\{x,\frac{1}{\pi }\right\}\right)$$ where, $$\text{P}(\text{z}) \text{ := }\frac{1-\pi ^4 z^4}{z^3}$$ $$S(\text{z})\text{ := } \pi ~ A(z) F(z)-\pi ^2 B(z)+\frac{\pi F'(z)}{2 (1-\pi z)}$$ $$F(\text{z})\text{ :=}\int_0^z \frac{\int_0^y \pi A(u) \text{P}(u) \, du}{\text{P}(y)} \, dy-\frac{3 \left(\log \left(\pi ^2 z^2+1\right)+\log (\pi z+1)+\log (1-\pi z)\right)}{64 \pi ^4}$$ $$A(\text{z})\text{ :=}\frac{\pi z (\pi z+2)+3}{4 z (\pi z-1) (\pi z+1) \left(\pi ^2 z^2+1\right)}$$ $$B(\text{z})\text{ :=}\frac{(\pi z (\pi z+2)+3) (\pi z (\pi z (\pi z+1)+1)+5)}{16 (\pi z-1) (\pi z+1)^2 \left(\pi ^2 z^2+1\right)^2}$$

• Can you give us an example of $y$, and possibly a reduced version of $f(x)$ as well? May 13, 2015 at 8:45
• Possible duplicate of mathematica.stackexchange.com/questions/17766/… May 13, 2015 at 9:38
• May 13, 2015 at 9:40
• It's a repeated integral; what you have is equivalent to $$\int_0^1 (1-x)f(x)\mathrm dx$$ That being said, and do forgive me for my bluntness, you've been quite vague as to the nature of your $f(x)$, and this makes it difficult to assist you. Methods that are efficient for one class fail spectacularly for another, so more elaboration is necessary. May 13, 2015 at 9:54
• You see, you should have included those expressions at the outset. It's gnarly now, but it should be amenable to some simplification. May 15, 2015 at 0:41

Something like this should work (assuming that y is not a complicated function, but it is easily generalized to that case):

f[x_?NumberQ]:=(*your function*)

I1[y_?NumberQ]:=NIntegrate[f[x],{x,0,y}]

I2=NIntegrate[I1[y],{y,0,1}]


As an example:

f[x_?NumberQ] := x;
I1[y_?NumberQ] := NIntegrate[f[x], {x, 0, y}];
NIntegrate[I1[y], {y, 0, 1}]


with the result:

0.166667

You can check with Integrate that theresult is correct.

The _?NumberQ in the definition of the functions mean that the function will be evaluated only if the input variable has a numerical value. A link to its use.

And here a link to a problem similar to yours.

• Thanks @mattiav27 ! Can you suggest how to modify it for complicated functions. Actually my functions are really complicated. Your code works fine for simple functions, I have checked. I don't understand why it isn't working for complicated functions. May 13, 2015 at 10:15
• @pinu I need your function: edit you first post May 13, 2015 at 10:34
• @pinu: meanwhile you can take a look at this old question of mine: there are two nested numerical integrations with complicated functions mathematica.stackexchange.com/a/56212/8822 May 13, 2015 at 12:20