Matemáticas con Mathematica

Práctica 9. Integración de funciones de varias variables

Francisco Javier Pérez González
Departamento de Análisis Matemático
Universidad de Granada

Inicialización

<< Graphics`Master`

\(\(\(Off[General::spell1]; 
    \)
 (*para\ que\ Mathematica\ no\ nos\ moleste\ con\ avisos\ de\ error\ 
        en\ los\ nombres\ de\ los\ comandos*) \)\)

(* ThisMathematica code was developed by Major Glenn Sjoden at the US Air Force 
Academy's Department of Mathematica l Sciences and later adapted to it's current format by Lt Col Steve Hadfield *)
\(Riemann2D[f_, xvars_, yvars_, \ nn_, mm_, \ opts___] := 
        CompoundExpression[
          
cell[xo_, yo_, zo_, x1_, y1_, z1_, colr_] := 
            
Graphics3D[
              Table[{FaceForm[colr], 
                  
Cuboid[{xo, yo, zo}, {x1, y1, z1}]}, {i, 
                  0, 1, 2}, {j, 1, 1, 2}, {k, 1, 1, 2}]], 
	
          cellw[xo_, yo_, zo_, x1_, y1_, z1_, colr_] := 
            
WireFrame[
              Graphics3D[
                Table[{colr, Cuboid[{xo, yo, zo}, {x1, y1, z1}]}, 
                  
{i, 0, 1, 2}, {j, 1, 1, 2}, {k, 1, 1, 
                    2}]]], 
xlo = xvars[\([2]\)], 
          xhi = xvars[\([3]\)], 
ylo = yvars[\([2]\)], 
          yhi = yvars[\([3]\)], 
          
dex = N[\((xhi - xlo)\)\/nn], 
          dey = N[\((yhi - ylo)\)\/mm], 
          
p3d = 
            Plot3D[f, xvars, yvars, 
              
BoxRatios -> {1, 1, 1.5}, 
              AxesLabel -> {x, y, z}, DisplayFunction -> Identity], 
          
RVol = 0, 
          
b3d = 
            Table[{xb = xa + dex; yb = ya + dey; 
                
xavg = \(xa + xb\)\/2; 
                yavg = \(ya + yb\)\/2; 
                
zb = 
                  f /. {xvars[\([1]\)] -> xavg, yvars[\([1]\)] -> yavg}; 
                RVol = RVol + \((xb - xa)\) \((yb - ya)\)*zb; 
                
cell[xa, ya, 0, xb, yb, zb, Blue]}, 
              
{xa, xlo, xhi - dex, dex}, {ya, ylo, 
                yhi - dey, dey}], \ 
          Show[{p3d, b3d}, \ opts, 
            PlotLabel -> 
              StyleForm[
                "Vol≅"\ TraditionalForm[\ RVol // N], 
                FontSize -> 12]]]; \)

fig1:=Module[{f},f[x_, y_] = (Sqrt[36 - 3  (x^2) - 6  (y^2)]);
G1=Plot3D[If[-2 ≤ x ≤ 2  &&  
                        -2 ≤ y ≤ 2, f[x, y], 0],{x,  -2.5,  2.5}, {y,  -2.5, 2.5},BoxRatios -> {1, 1, 1.5},BoxStyle  ->  RGBColor[1, 0, 0],Background -> RGBColor[1, 1, 0.75],AxesEdge 
                      -> {{-1, -1}, {1, -1}, {-1, 
                          -1}},Ticks 
                      -> {{{-2, "-2"}, {0, "0"}, {2, 
                            "2"}}, {{2, "2"}, {0, "0"}}, {{2, 
                            "2"}, {4, "4"}, {6, "6"}}},PlotRange 
                      -> {{-2.5, 2.5}, {-2.5, 2.5}, {0, 6}},AxesLabel -> {"X", "Y", "Z"},DisplayFunction -> Identity];G2=ParametricPlot3D[{0, y, t*f[0, y]}, {y,  -2, 2},{t, 0, 1},BoxRatios -> {1, 1, 1.5},BoxStyle  ->  RGBColor[1, 0, 0],Background -> RGBColor[1, 1, 0.75],AxesEdge 
                      -> {{-1, -1}, {1, -1}, {-1, 
                          -1}},Ticks -> {{{0, "X=0"}}, {{0, "Y=0"}, {2, 
                            "2"}, {-2, "-2"}}, {{2, "2"}, {4, 
                            "4"}, {6, "6"}}},PlotRange 
                      -> {{-2.5, 2.5}, {-2.5, 2.5}, {0, 6}},DisplayFunction -> Identity];
G3=FilledPlot[f[0, y],{y, -2, 2},PlotRange -> {{-3, 3}, {0, 6.5}},AspectRatio -> Automatic,Fills -> RGBColor[0, 1, 1],DisplayFunction -> Identity,Background -> RGBColor[1, 1, 0.75],AxesLabel -> {"y", "z=f(0,y)"},AxesStyle  ->  RGBColor[1, 0, 0],Ticks 
                      -> {{{-2, "-2",  .01}, {0, 
                            "0",  .01}, {2, "2",  .01}}, {{2, 
                            "2",  .01}, {4, "4",  .01}, {6, 
                            "6",  .01}}}];
Show[GraphicsArray[{G1, G2, G3}], 
                ImageSize -> {528, 220}, 
                ImageMargins -> {{0, 0}, {0, 0}}, 
                ImageRegion -> {{0, 1}, {0, 1}}, 
                DisplayFunction -> $DisplayFunction]];

polares := 
      Show[Graphics[{Line[{{- .3, 0}, {2, 0}}], 
            Line[{{0, - .3}, {0, 2}}], Text["Y", {0.1, 1.9}], 
            Text["X", {2,  .1}], 
            Circle[{0, 0},  .4, {0, π/4}], {Hue[ .6], 
              Text["(x,y)", {1.1, 1.15}]}, {PointSize[ .035], Hue[ .6], 
              Point[{1, 1}], Hue[0], Point[{0, 1}], Point[{1, 0}], 
              Text["ρcosθ", {1, - .15}], 
              Text["ρsenθ", { .33, 1.15}], 
              Text["θ", { .5,  .2}]}, {Line[{{0, 0}, {1, 1}}], 
              Dashing[{ .02,  .02}], Line[{{0, 1}, {1, 1}, {1, 0}}]}}], 
        Background -> RGBColor[1, 1, 0.75], 
        AspectRatio -> Automatic, 
        TextStyle -> {FontSize -> 12, FontWeight -> "Bold"}, 
        ImageSize -> {390, 162}];

esfericas := 
      Show[Graphics[{
            Line[{{0, 6}, {0, 0}, {2, 4}, {0, 0}, {6, 0}, {0, 0}, {-4, 
                  -4}}], {Dashing[{ .012,  .012}], 
              Line[{{0, 5.5}, {2, 4}, {2, -1.5}, {0, 0}}]}, {
              Dashing[{ .012,  .012}], Line[{{2, -1.5}, {3.5, 0}}]}, {
              Dashing[{ .012,  .012}], 
              Line[{{2, -1.5}, {-1.5, -1.5}}]}, 
            Circle[{0, 0}, 2, {ArcTan[2], π/2}], 
            Circle[{0, 0},  .6, {-3 π/4, ArcTan[- .75]}], {
              Hue[0], PointSize[ .02],  Point[{0.01, 5.5}], 
              Text["ρcosφ", {1, 5.6}], 
              Text["ρ", {1.3, 2.12}], 
              Text["ρsenφ", {2.2, -1.8}], 
              Text["ρsenφcosθ", {-1.7, 
                  -1.8}], 
              Text["ρsenφsenθ", {3.7,  .35}], 
              Text["φ", { .485, 2.47}], 
              Text["θ", {0, -1}]}, {Hue[0], PointSize[ .02],  
              Point[{3.5, 0}]}, {Hue[0], PointSize[ .02],  
              Point[{-1.5, -1.5}]}, {Hue[0.6], PointSize[ .02],  
              Point[{2, 4}], Text["(x,y,z)", {2.8, 4.4}]}, 
            Text["X", {-3.4, -3.9}], Text["Z", {0.3, 5.95}], 
            Text["Y", {5.95, 0.25}]}], AspectRatio -> Automatic, 
        Background -> RGBColor[1, 1, 0.75], 
        TextStyle -> {FontSize -> 14, FontWeight -> "Bold"}];

Integración de funciones de dos y tres variables

Para tener una idea intuitiva de lo que se hace cuando se calcula una integral doble, consideremos una función (todas las funciones se suponen continuas) positiva de dos variables f : A⟶ℝ , f(x, y) ≥ 0 para (x, y) ∈ A donde A ⊂ (ℝ^2). En tal caso, la integral doble de la función f en el conjunto A es el volumen del "cilindro" en (ℝ^3) cuya base es el conjunto A y que tiene como "tapadera" la gráfica de f.

∫(∫_A) f(x, y) d(x, 
                y) = 
        volumen {(x, y, z) :  (x, y) ∈  A,  
            0 ≤ z ≤  f(x, y)}

Para calcular el volumen del conjunto  {(x, y, z) :  (x, y) ∈  
            A,  0 ≤ z ≤  f(x, y)} lo que se hace es obtener aproximaciones y pasar al límite. Consideremos el caso en que  A =  [a, b]×[c, d]  y sean a = (Subscript[x, 0]) < (Subscript[x, 1]) < (Subscript[x, 2]) <   ...   < 
            (Subscript[x, n - 1]) < (Subscript[x, n]) = b y  c = (Subscript[y, 0]) < (Subscript[y, 1]) < (Subscript[y, 2]) <   ...   < 
            (Subscript[y, m - 1]) < (Subscript[y, m]) = d particiones de los intervalos [a, b] y  [c, d] respectivamente.  Elijamos puntos ((Subscript[s, k]) , 
          (Subscript[t, j])) ∈ [(Subscript[x, k - 1]), (Subscript[x, k])]×[(Subscript[y, j - 1]), 
            (Subscript[y, j])], y formemos el ortoedro cuya base es el rectángulo [(Subscript[x, k - 1]), (Subscript[x, k])]×[(Subscript[y, j - 1]), 
          (Subscript[y, j])] y altura igual a f((Subscript[s, k]), (Subscript[t, j])). El volumen de dicho ortoedro es f((Subscript[s, k]), (Subscript[t, j]))((Subscript[x, k]) - (Subscript[x, k - 1]))((Subscript[y, j]) - 
              (Subscript[y, j - 1])). La suma

(∑_(k = 1)^n)(∑_(j = 1)^m)f((Subscript[s, k]), (Subscript[t, j]))((Subscript[x, k]) - (Subscript[x, k - 1]))((Subscript[y, j]) - (Subscript[y, j - 1]))

se llama una suma de Riemann de la función f y su valor es una aproximación al volumen buscado. La integral  de la función f en el rectángulo A = [a, b]×[c, d] se define como el límite de las sumas de Riemann cuando las longitudes de los intervalos de las particiones tienden a cero.

Consideremos la función f(x, y) = 1 + (x^2) + (y^2). Como ya sabes, el comando "Plot3D[ ]" permite representar la gráfica de dicha función sobre un rectángulo.

\(f[x_, y_] = 1 + x^2 + y^2; 

    Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}, 
      Background -> RGBColor[1, 1, 0.75], 
      PlotRange -> {{\(- .2\), 1.2}, {\(- .2\), 1.2}, Automatic}, 
      AspectRatio -> 1, BoxRatios -> {1, 1, 1}, 
      PlotPoints -> 30]; \)

Vamos a calcular la integral de f sobre el cuadrado unidad [0,1]x[0,1]. Tenemos que calcular el volumen del cilindro en \(TraditionalForm\`\(\(ℝ\^3\)\(\ \)\(cuya\)\(\ \)\(base\)\(\ \)\(es\)\(\ \)\(dicho\)\(\ \)\(cuadrado\)\(\ \)\(y\)\(\ \)\(cuya\)\(\ \)\(tapadera\)\(\ \)\(es\)\(\ \)\(la\)\(\ \)\(gráfica\)\(\ \)\(de\)\(\ \)\(la\)\(\ \)\(función . \ 
          Para\)\(\ \)\(verlo\)\(\ \)\(mejor\)\(\ \)\(representamos\)\(\ \)\(la\)\(\ \)\(función\)\(\ \)\)\)

If[
      0 ≤ x ≤ 1  && 0 ≤ y ≤ 1, 
      1 + (x^2) + (y^2), 0]

Plot3D[ 
      If[0 ≤ x ≤ 1 && 0 ≤ y ≤ 1, 
        1 + x^2 + y^2, 0],  {x, - .2, 1.2},  {y, - .2,   1.2}, 
      Background -> RGBColor[1, 1, 0.75], AspectRatio -> 1, 
      BoxRatios -> {1, 1, 1}, PlotPoints -> 60];

De forma análoga a como se hace para funciones de una variable, consideramos particiones en cada uno de los intervalos [0,1]  lo que nos proporciona una cuadrícula que divide el cuadrado [0,1]x[0,1] en pequeños rectángulos. Aquí puedes ver la cuadrícula correspondiente a dividir en 10 partes de igual longitud cada intervalo. En cada uno de los cuadraditos de la cuadrícula se ha representado su centro.

Module[{
        averg1 = 
          (Drop[Range[0, 1,  .1], 1] + Drop[Range[0, 1,  .1], -1])/
            2}, 
CartesianMap[Identity, 
         {0, 1, 0.1}, {0, 1, 0.1}, 
        
Ticks 
          -> {Range[0, 1,  .1], Range[0, 1,  .1]}, 
        AxesLabel -> {"x", "y"}, PlotStyle -> {Blue, Red}, 
        Epilog -> {PointSize[0.015], Hue[1], 
            Flatten[Outer[Point, averg1, averg1] /. 
                Point[x_, y_] -> Point[{x, y}]]}]];

Si multiplicamos el área de un cuadradito por el valor de la función en su centro, obtenemos el volumen de un ortoedro cuya base es el cuadradito y cuya altura es el valor de la función en su centro. Si ahora sumamos los volúmenes de todos esos ortoedros obtenemos una suma de Riemann que es un valor aproximado del volumen exacto del cilindro. Aquí puedes ver representado gráficamente este proceso. El comando "Riemann2D[f[x,y], {x,a,b}, {y,c,d}, n, m, Options]" hace una partición del rectángulo  [a, b]×[c, d] en  n×m subrectángulos, calcula el valor de f en el punto medio de cada uno de ellos y representa la gráfica de f y la suma de Riemann correspondiente con el valor aproximado del volumen. Admite opciones como "Plot3D".

Show[
      Riemann2D[1 + (x^2) + (y^2), {x, 0, 1}, {y, 0, 1}, 10, 10, 
        Background -> RGBColor[1, 1, 0.75], AspectRatio -> 1, 
        PlotRange -> {{- .2, 1.2}, {- .2, 1.2}, Automatic}, 
        BoxRatios -> {1, 1, 1}], 
      DisplayFunction -> $DisplayFunction];

Puedes comprobar cómo la aproximación mejora conforme aumenta el número de puntos de división

Show[
      GraphicsArray[
        Table[Riemann2D[1 + (x^2) + (y^2), {x, 0, 1}, {y, 0, 1}, i, i, 
            DisplayFunction -> Identity, Axes -> False, 
            PlotRange -> {{- .2, 1.2}, {- .2, 1.2}, Automatic}, 
            Background -> RGBColor[1, 1, 0.75]], {i, 8, 14, 3}]], 
      DisplayFunction -> $DisplayFunction, ImageSize -> {528, 220}, 
      ImageMargins -> {{0, 0}, {0, 0}}, 
      ImageRegion -> {{0, 1}, {0, 1}}];

El valor de la integral es

Integrate[1 + x^2 + y^2, {x, 0, 1}, {y, 0, 1}]
N[%]

Como puedes ver la aproximación obtenida no es mala.

Para integrar una función f[x, y] sobre un conjunto acotado A que no sea un rectángulo, lo que se hace es encerrar el conjunto en un rectángulo A ⊂ [a, b] x[c, d] e integrar en dicho rectángulo la función g[x, y] que coincide con f[x, y]cuando  (x, y) ∈ A y g[x, y] = 0en otro caso. Dicha función puedes definirla en Mathematica en la forma g[x, y] = If[cjto[x, y], f[x, y], 0] donde cjto[x,y] debe ser una función lógica cuyo valor sea True cuando (x, y) ∈ A y False en otro caso. Con la ayuda de esta misma función puedes representar la gráfica de la función f cuando las variables (x, y) están en el conjunto A. Aquí puedes ver el cilindro cuya base es el círculo de centro (1/2,1/2) y radio 1/2 y cuya tapadera es la gráfica de la funcón.

cjto[x_, y_] := 
      0 ≤ x ≤ 1 && 
        (1/2) - (Sqrt[(1/4) - (((x - (1/2)))^2)]) ≤ y ≤  
          (1/2) + (Sqrt[(1/4) - (((x - (1/2)))^2)])   g[x_, y_] := If[cjto[x, y], 1 + (x^2) + (y^2), 0]Plot3D[g[x, y], {x, 0, 1}, {y, 0, 1}, 
      Background -> RGBColor[1, 1, 0.75], 
      PlotRange -> {{- .2, 1.2}, {- .2, 1.2}, Automatic}, 
      AspectRatio -> 1, BoxRatios -> {1, 1, 1}, 
      PlotRange -> {0, 3}, PlotPoints -> 50];

Nota: A partir de la versión 4 también puedes usar para esto la función "Boole[ ]".

Aquí puedes ver la aproximación de las sumas de Riemann a la integral cuando integramos la función en el círculo de centro (1/2,1/2) y radio 1/2.

Show[
      Riemann2D[g[x, y], {x, 0, 1}, {y, 0, 1}, 20, 20, 
        Background -> RGBColor[1, 1, 0.75], AspectRatio -> 1, 
        PlotRange -> {{- .2, 1.2}, {- .2, 1.2}, Automatic}, 
        BoxRatios -> {1, 1, 1}], 
      DisplayFunction -> $DisplayFunction];

El valor de la integral es

Integrate[
      1 + x^2 + y^2, {x, 0, 1}, {y, 1/2 - (Sqrt[1/4 - (x - 1/2)^2]), 
        1/2 + (Sqrt[1/4 - (x - 1/2)^2])}]
N[%]Clear["@"]

La integral de una función de tres variables se define de forma parecida. Dicha integral, para el caso de una función positiva f  definida en un conjunto A de (ℝ^3)se interpreta como el "volumen" del "cilindro" en  (ℝ^4) que tiene como base el conjunto A y como "tapadera" la gráfica de f.

∫∫(∫_A) f(x, y, 
                  z) d(x, y, z) = 
        volumen {(x, y, z, w) :  (x, y, z) ∈ A ,  
            0 ≤ w ≤  f(x, y, z)}

Observa que la integral de una función positiva de una variable mide el área de un conjunto de (ℝ^2), la integral de una función positiva de dos variables mide el volumen de un conjunto de (ℝ^3)y la integral de una función positiva de tres variables mide el "volumen" de un conjunto de (ℝ^4). No obstante, en el caso particular de que integremos la función constante igual a 1 en un conjunto A de (ℝ^2), resulta que
						∫(∫_A) 
                d(x, y)=área(A) 
sin más que tener en cuenta que el volumen de un cilindro recto de altura constante igual a 1 coincide con el área de la base. Análogamente si integramos la función constante igual a 1 en un conjunto A de (ℝ^3), resulta que
						∫∫(∫_A) 
              d(x, y, z) = volumen (A)

Ejercicio 1

Teorema de Fubini

El teorema de Fubini permite calcular una integral doble haciendo dos integrales simples. No es difícil comprender lo que dice el teorema ni tampoco lo es entender por qué se cumple. De hecho, no es la primera vez que en este curso aparece dicho teorema, ya lo usamos para calcular volúmenes de sólidos de revolución cuando estudiamos las aplicaciones de las Integrales simples.  Un caso particular del teorema de Fubini afirma que el volumen de una región en (ℝ^3)es igual a la integral del área de sus secciones por planos paralelos a uno dado. Por tanto, el teorema de Fubini permite calcular volúmenes calculando áreas de secciones planas. 
Consideremos una función positiva,  f , definida en el rectángulo  A = [a, b]×[c, d]. Para calcular el volumen del conjunto  Ω = {(x, y, 
              z) :  (x, y) ∈  A,  
          0 ≤ z ≤  f(x, y)} podemos proceder como sigue. Para cada (Subscript[x, 0]) ∈ [a, b] fijo calculamos el área de la sección, Ω((Subscript[x, 0])),  que se obtiene cortando Ω  con el plano de ecuación X = (Subscript[x, 0]). Fíjate que  Ω((Subscript[x, 0])) es una sección de Ω perpendicular al eje OX  y, por tanto, paralela al plano YZ. Como

Ω(
          (Subscript[x, 0])) = {((Subscript[x, 0]), y, 
                z) : (x, y, z) ∈ Ω} = {((Subscript[x, 0]), 
                y, z) : y ∈ [c, d] ,  
            0 ≤ z ≤ f((Subscript[x, 0]), y)}

se tiene que Ω((Subscript[x, 0])) es la región del plano X = (Subscript[x, 0]) comprendida entre la curva  z=f((Subscript[x, 0]), y), el eje OY y las rectas y = a, y = b. Como sabes, el área de dicha región viene dada por λ(Ω(
            (Subscript[x, 0]))) = (∫_c^d) f((Subscript[x, 0]), y) dy (si tienes alguna duda, observa que la proyección de Ω((Subscript[x, 0])) sobre el plano YZ  es el conjunto {(y, z) : y ∈ [c, d] ,  
        0 ≤ z ≤ f((Subscript[x, 0]), y)}).  Ahora hay que integrar las áreas de las secciones Ω(x) cuando x ∈ [a, b] y obtenemos finalmente que

\(TraditionalForm\`∫\(∫\_A\( f(x, y)\) \(d(x, 
                y)\)\) = \(volumen\ \((Ω)\) = \(∫\_a\%b\( λ(Ω(
                    x))\) dx = ∫\_a\%b\([∫\_c\%d\( f(x, y)\) dy]\) dx\)\)\)

Razonando de forma análoga, considerando secciones Ω(y) de Ω paralelas al plano XZ, se obtiene la igualdad

\(∫\(∫\_A\( f(x, y)\) \(d(x, 
                    y)\)\) = \(volumen\ \((Ω)\) = \(∫\_c\%d\( λ(Ω(
                        y))\) dy = ∫\_c\%d\([∫\_a\%b\( f(x, y)\) dx]\) dy\)\)\)

Las integrales (∫_a^b)[(∫_c^d) f(x, 
                      y) dy] dx y (∫_c^d)[(∫_a^b) f(x, 
                      y) dx] dy se llaman integrales iteradas y, como consecuencia del teorema de Fubini, son iguales y su valor común es igual a la integral ∫(∫_A) f(x, y) d(x, 
              y). Observa que las integrales iteradas son dos integrales simples. Para calcular (∫_c^d) f(x, y) dy lo que se hace es integrar respecto a la variable y considerando x fija. Para ello lo que se hace es obtener una primitiva de la función y⟶f(x, y) y usar la regla de Barrow. Fíjate que una primitiva de la función y⟶f(x, y) puede describirse como una primitiva parcial de f(x, y) con respecto a y. ¿Te recuerda esto a la derivación parcial?
La representación gráfica siguiente puede ayudarte a entender lo que se hace. La función representada es f(x, y) = (Sqrt[36 - 3  (x^2) - 6  (y^2)]) en el rectángulo  [-2, 2]×[-2, 
          2]. Puedes ver el "cilindro" bajo la gráfica de la función, la sección del mismo por el plano X=0 y la proyección de dicha sección sobre el plano XZ.

fig1;

Para calcular una integral  
      ∫(∫_A) f(x, y) d(x, y) cuando el recinto de integración A no es rectángulo se procede de la misma forma. La única diferencia es que ahora tenemos que empezar por determinar los valores de x tales que el plano X=x corta al "cilindro" bajo la gráfica de f , es decir, tenemos que determinar la proyección de A sobre el eje OX . Supongamos que dicha proyección sea un intervalo [a, b]. Ahora, para cada x ∈ [a, b] hay que calcular el área de la sección Ω(x), o lo que es igual, el área de la región en el plano YZ comprendida entre el eje OY y la curva z = f(x, y) donde la variable y está en el conjunto A(x) = {y : (x, y) ∈ A}. Supongamos que A(x) sea un intervalo (tampoco pasa nada si es unión de varios intervalos). Entonces tenemos que

\(TraditionalForm
      \`∫\(∫\_A\( f(x, y)\) \(d(x, y)\)\) = 
        ∫\_a\%b
            \([∫\_\(A(x)\)\(f(x, y)\) dy]\) 
            dx\).

Análogamente se obtiene que

∫(∫_A)f(x, y) d(x, y) = 
            (∫_c^d)
                [(∫_(A(y)))f(x, y) dx] 
                dy

donde [c, d] es la proyección de A sobre el eje OY, y para cada y ∈ [c, d], A(y) = {x : (x, y) ∈ A}.

En los casos más corrientes el conjunto A suele ser un conjunto de tipo I o de tipo II (recuerda que los estudiamos en el cuaderno de Aplicaciones de la Integral). Esto es A = {(x, y) : a ≤ x ≤ b,  
          g(x) ≤ y ≤ h(x)} (tipo I) o bien A = {(x, y) : 
            c ≤ y ≤ d,  ϕ(y) ≤ 
            x ≤ σ(y)} (tipo II). En tales casos tenemos que

\(TraditionalForm\`∫\(∫\_A\( f(x, y)\) \(d(x, 
                y)\)\) = \(∫\_a\%b\([∫\_\(A(x)\)\(f(x, 
                      y)\) dy]\) dx = ∫\_a\%b\([∫\_\(g(x)\)\%\(h(x)\)\(f(x, 
                      y)\) dy]\) dx\)\)   (A de tipo I)

\(∫\(∫\_A\( f(x, y)\) \(d(x, 
                    y)\)\) = \(∫\_c\%d\([∫\_\(A(y)\)\(f(
                          x, y)\) dx]\) dy = ∫\_c\%d\([∫\_\(ϕ(y)\)\%\(σ(y)\)\(f(x, 
                          y)\) dx]\) dy\)\)  (A de tipo II)

De forma análoga a lo antes visto, el teorema de Fubini permite calcular integrales triples sin más que calcular tres integrales simples. Para el caso de una función  f  definida en el rectángulo de (ℝ^3) \(TraditionalForm
      \`A = \([a, b]\)×\([c, d]\)×\([u, v]\)\) se tiene que

\(TraditionalForm
      \`∫\(∫\(∫\_A\( f(x, y, z)\) 
                \(d(x, y, z)\)\)\) = 
        ∫\_a\%b
            \([∫\_c\%d
                  \([∫\_u\%v\( f(x, y, z)\) dz]\) 
                  dy]\) dx\)

Observa que ahora hay seis integrales iteradas pero el valor de todas ellas es el mismo. Naturalmente, cuando A es un subconjunto de (ℝ^3) hay más posibilidades, pero la idea es siempre la misma: se obtiene primero la proyección de A sobre uno de los ejes o sobre uno de los planos coordenados, y para cada punto fijado en dicha proyección se obtiene el conjunto de los puntos de A que lo proyectan. 
Si, por ejemplo, la proyección de A sobre el eje OZ es un intervalo J = [s, t], y para cada z ∈ J es \(TraditionalForm
      \`A(z) = {\((x, y)\) : \((x, y, z)\) ∈ A}\), entonces:
			∫∫(∫_A) f(x, y, z) 
                    d(x, y, z)=(∫_s^(t 
              ))[∫(∫_(A(z)))f(x, y, z) 
                      d(x, y)] dz
En el caso en que A sea un conjunto de tipo I en (ℝ^3), es decir, A puede representarse en la forma \(TraditionalForm
      \`A = {\((x, y, z)\) : \((x, y)\) ∈ Ω, \ 
          g(x, y) ≤ z ≤ \ h(x, y)}\), donde Ω es la proyección de A sobre el plano XY, y  g , h  son funciones reales definidas en Ω, entonces:
			\(TraditionalForm
      \`∫\(∫\(∫\_A\( f(x, y, z)\) 
                d\((x, y, z)\)\)\) = 
        ∫\(∫\_\(Ω\ 
              \)\([∫\_\(g(x, y)\)\%\(h(x, y)\)\(f(x, y, z)\) 
                    dz]\) d\((x, y)\)\)\)

Los comandos Integrate[ ], NIntegrate[ ]

El comando "Integrate[f[x,y],{x,a,b},{y,c,d}]" calcula simbólicamente (cuando puede) la integral de f[x,y] cuando la variable x ∈ [a, b], y ∈ [c, d]. La integración se realiza primero respecto de la variable y.

\(TraditionalForm
    \`Integrate[f[x, y], {x, a, b}, {y, c, d}] = 
      ∫\_a\%b\([∫\_c\%d f[x, y] dy]\) 
          dx\)

Los límites de integración para la variable y pueden depender de la variable x. El comando "NIntegrate[f[x,y],{x,a,b},{y,c,d}]" calcula numéricamente la integral. También puedes usar, con las modificaciones previsibles, estos comandos para funciones de tres variables. A partir de la versión 5 de Mathematica estos comandos se han mejorado bastante y pueden integrar funciones del tipo "If[cjto[x,y],f[x,y],0]" para proporcionar directamente la integral de f[x, y] en el conjunto definido por "cjto[x,y]".

Ejemplo 1

Calcular la integral de la función f(x, y) = ((x y)/(1 + (y^2))) en el conjunto \(TraditionalForm
      \`A = {\((x, y)\) : x\^2 + y\^2 ≤ 1, \ x ≥ 0, 
          \ y ≥ 0}\).

Integrate[((x*y)/(1 + (y^2))), {x, 0, 1}, {y, 0, (Sqrt[1 - (x^2)])}]

Mathematica ha calculado primero la integral

Integrate[((x*y)/(1 + (y^2))), {y, 0, (Sqrt[1 - (x^2)])}]

y la función obtenida la ha integrado en [0,1].

Integrate[(1/2) x Log[2 - (x^2)], {x, 0, 1}]

Podríamos haberlo hecho también invirtiendo el orden de integración.

Integrate[((x*y)/(1 + (y^2))), {y, 0, 1}, {x, 0, (Sqrt[1 - (y^2)])}]

\(Integrate[\(x*y\)\/\(1 + y\^2\), {x, 0, \@\(1 - y\^2\)}]
\)Integrate[%, {y, 0, 1}]

Ejercicio 2

Ejercicio 3

Teorema del cambio de variables

Para funciones de una variable sabemos que \(TraditionalForm
      \`∫\_a\%b\( f(x)\) dx = 
        ∫\_c\%d\( f(g(t))\) g' \((t)\) dt\) donde se supone que a = g(c) y b = g(d). Supongamos que la función g es inyectiva, entonces g debe ser creciente o decreciente. Si es decreciente se tiene que d < c y g' (x) ≤ 0, por lo que \(TraditionalForm
      \`∫\_c\%d\( f(g(t))\) g' \((t)\) dt = 
        \(\(-\(∫\_d\%c\( f(g(t))\) g' \((t)\) 
                dt\)\) = 
          ∫\_d\%c f(g(t)) | g' \((t)\) | dt\)\). Podemos, por tanto, cuando g es inyectiva, escribir en todos los casos

\(TraditionalForm
      \`∫\_a\%b\( f(x)\) dx = 
        ∫\_α\%β f(g(t)) | g' \((t)\) | 
          dt\)

donde g es una biyección del intervalo [α, β] sobre [a, b]. Esta fórmula se generaliza para funciones de varias variables dando lugar al teorema del cambio de variables.

El teorema del cambio de variable para integrales dobles afirma que

\(TraditionalForm
      \`∫\(∫\_A\( f(x, y)\) \(d(x, y)\)\) = 
        ∫\(∫\_B\ f(g(u, v))\) | det\ \(\(J\_g\)(u, v)\) | 
          d(u, v)\)

donde se supone que la función g es una biyección de B sobre A de clase (C^1)(sus funciones componentes tienen derivadas parciales continuas) con determinante jacobiano distinto de cero, esto es, det (Subscript[J, g])(u, v) ≠  0 para todo (u, v) ∈ B. En esta fórmula se interpreta que la función g hace un cambio de coordenadas pues permite asignar a cada punto (x, y) ∈ A el único punto (u, v) ∈ B tal que g(u, v) = (x, y). 
Coordenadas polares. La función \(TraditionalForm
      \`g(ρ, θ) = 
        \((ρ\ cos\ θ, \ ρ\ sen\ θ)\)\)es una biyección de Ω = \(\(ℝ\^+\)×\)]-π, π[ sobre \(TraditionalForm
      \`ℝ\^2\{\((x, 0)\) : x ≤ 0}\) (el plano al que se ha quitado el semieje de abscisas negativo). Las componentes de g tienen derivadas parciales continuas y fácilmente se comprueba que \(TraditionalForm
      \`det\ \(\(J\_g\)(ρ, θ)\) = ρ > 0\). Los números ρ y θ dados por x = ρ cos θ,  
      y = ρ sen θ donde ρ > 0, -π<θ≤π se llaman las coordenadas polares del punto (x, y).

polares;

Como aplicación del teorema del cambio de variables obtenemos la igualdad
		\(TraditionalForm\`∫\(∫\_A\( f(x, y)\) \(d(x, 
                y)\)\) = ∫\(∫\_B\ \(f(ρ\ cos\ θ, \ ρ\ sen\ θ)\) ρ\ \(d(ρ, θ)\)\)\) 
donde B={(ρ, θ):\((ρ\ cos\ θ, \ ρ\ sen\ θ)\)∈ A}. Si, por ejemplo, el conjunto A es de tipo I,  A = {(x, y) : a ≤ x ≤ b,  
          g(x) ≤ y ≤ h(x)}, entonces B={(ρ, θ):\(a ≤ ρ\ cos\ θ ≤ b, 
                g(\ ρ\ \ cos\ θ) ≤ \ ρ\ sen\ θ ≤ \ h(ρ\ cos\ θ)\)}. Es importante describir bien el conjunto B porque para calcular la integral \(TraditionalForm\`∫\(∫\_B\ \(f(ρ\ cos\ θ, \ ρ\ sen\ θ)\) ρ\ \(d(ρ, θ)\)\)\) tienes que aplicar, naturalmente, el teorema de Fubini. Si, por ejemplo, el conjunto B es de tipo II,  B = {(ρ, θ) : 
                  c ≤ θ ≤ d,  
                g(θ) ≤ ρ ≤ h(θ)}, entonces
		\(TraditionalForm\`∫\(∫\_A\( f(x, y)\) \(d(x, 
                y)\)\) = \(∫\(∫\_B\ \(f(ρ\ cos\ θ, \ ρ\ sen\ θ)\) ρ\ \(d(ρ, θ)\)\) = ∫\_c\%d\([∫\_\(g(θ)\)\%\(h(θ)\)\(f(ρ\ cos\ θ, \ ρ\ sen\ θ)\) ρ\ dρ]\) dθ\)\)
Las coordenadas polares son especialmente útiles cuando el conjunto A es un círculo, o un sector circular o una corona circular, pues en estos casos el conjunto B es muy sencillo. Si, por ejemplo, A es el disco de centro el origen y radio R,  A  = D(0, 
            R) = {(x, y) : 
              (x^2) + (y^2) ≤ (R^( 2))} entonces \(TraditionalForm\`B = \({\((ρ, θ)\) : ρ ≤ \(R\^\ \)} = \([0, R]\)×\([\(-π\), π]\)\)\). Por tanto
	\(TraditionalForm\`∫\(∫\_\(D(0, R)\)\(f(x, y)\) \(d(
                x, y)\)\) = \(∫\_\(\(-π\)\(\ \)\)\%π\([∫\_0\%R\( 
                      f(ρ\ cos\ θ, \ ρ\ sen\ θ)\) ρ\ dρ]\) dθ = ∫\_0\%R\ \([∫\_\(-π\)\%π\( 
                      f(ρ\ cos\ θ, \ ρ\ sen\ θ)\) ρ\ dθ]\) dρ\)\)

El teorema del cambio de variable para integrales triples afirma que
			∫(∫_A) f(x, y, z) d(x, y, 
                z) = ∫(∫_B) f(g(u, v, w)) | 
          det (Subscript[J, g])(u, v, w) | d(u, v, w) 
donde se supone que la función g es una biyección de B sobre A de clase (C^1)(sus funciones componentes tienen derivadas parciales continuas) con determinante jacobiano distinto de cero, esto es, det (Subscript[J, g])(u, v, w) ≠  0 para todo (u, v, w) ∈ B.

Coordenadas esféricas. La función \(TraditionalForm\`g(ρ, θ, ϕ) = \((ρ\ senϕ\ cos\ θ, \ ρ\ senϕ\ sen\ θ, ρ\ cosϕ)\)\) es una biyección de Ω=\(TraditionalForm\`\(\(ℝ\^+\)\(×\)\)\)]-π, π[×]0, π[ sobre \(TraditionalForm\`ℝ\^3\{\((x, 0, z)\) : 
            x ≤ 0}\). Las componentes de g tienen derivadas parciales continuas y fácilmente se comprueba que \(TraditionalForm\`det\ \(\(J\_g\)(ρ, θ, ϕ)\) = \(-ρ\^2\) senϕ\). Los números ρ y θ y ϕ dados por x = ρ senϕ cos θ,  
      y = ρ senϕsen θ,  z = ρ cosϕ donde ρ > 0, -π<θ≤π, 0 ≤  ϕ ≤ π se llaman las coordenadas esféricas del punto (x, y, z).

esfericas;

Como aplicación del teorema del cambio de variables obtenemos la igualdad
		\(TraditionalForm
      \`∫\(∫\(∫\_A\ 
              \(f(x, y, z)\) \(d(x, y, z)\)\)\) = 
        ∫\(∫\(∫\_B\ 
              \(f(ρ\ senϕ\ cos\ θ, \ 
                  ρ\ senϕ\ sen\ θ, ρ\ cosϕ)\) 
                \(ρ\^2\) senϕ\ 
                \(d(ρ, θ, ϕ)\)\)\)\) 
donde B={(ρ, θ, ϕ):(ρ senϕ cos θ,  
                      ρ senϕ sen θ, ρ cosϕ)∈ A}. Las coordenadas esféricas son especialmente útiles cuando el conjunto A es una esfera, o un sector esférico o una corona esférica, pues en estos casos el conjunto B es muy sencillo. Si, por ejemplo, A es la esfera de centro el origen y radio R,  \(TraditionalForm
      \`A\  = \(B(0, R) = {
            \((x, y, z)\) : x\^2 + y\^2 + z\^2 ≤ R\^\(\ 2\)}
            \)\) entonces \(TraditionalForm
      \`B = \({\((ρ, θ, ϕ)\) : 
              ρ ≤ R\^\(\  \)} = 
          \([0, R]\)×\([\(-π\), π]\)×\([0, π]
              \)\)\). Por tanto
	\(TraditionalForm
      \`∫\(∫\(∫\_\(B(0, R)\)\(f(x, y, z)\) 
                \(d(x, y, z)\)\)\) = 
        ∫\_0\%R
            \([∫\_\(\(-π\)\ 
                  \)\%π\([
                    ∫\_0
                        \%π\( f(ρ\ senϕ\ cos\ θ, \ 
                          ρ\ senϕ\ sen\ θ, ρ\ cosϕ)
                          \) \(ρ\^2\) sen\ ϕ dϕ]
                    \) dθ]\) dρ\)
y la integral iterada puede hacerse en el orden que se quiera.

Ejercicio 4