Integración numérica con MATLAB
Integrales simples
La función 
>> f=@(x) -x.^2+14*x+21; >> q=integral(f,0,10) q = 576.6667
Integrales dobles
Queremos calcular una integral doble de la función f(x,y) en la región rectangular de abscisas a y b y de ordenadas c y d, tal como se muestra en la figura
MATLAB dispone de la función 
Creamos un script, en el que definimos la función anónima f(x,y) y llamamos a 
f=@(x,y) x.^2.*y;
res=integral2(f,1,2,0,3);
fprintf('El valor de la integral es: %2.3f\n',res)
En la ventana de comandos corremos el script
El valor de la integral es: 10.500
Como puede comprobarse fácilmente haciendo la integral doble.
Vamos a calcular la integral doble de la función f(x,y) en la región no rectangular comprendida entre las curvas y=c(x), y=d(x) y las rectas x=a y x=b
Al integrar la función f(x,y) respecto de la variable y entre los límites c(x) y d(x) nos queda la función g(x) que integramos entre los límites a y b.
Utilizamos la función MATLAB 
Definimos la función g(x) del siguiente modo
function z=g(x)
    n=length(x);
    z=zeros(size(x));
    for j=1:n
    %límites de y
        c=x(j)^3-x(j);
        d=x(j)^2+x(j);
    %integrando f(x,y)
        f=@(y) x(j)^3*y.^4+x(j)*y.^2;
        z(j)=integral(f,c,d);
    end
end
    En el script llamamos a la función 
res=integral(@g,0,2);
fprintf('El valor de la integral es: %3.3f\n',res)
    En la ventana de comandos corremos el script 
El valor de la integral es: 961.181
Ejemplos
Comprobar los resultados de las integrales siguientes:
Comprobamos la primera integral. Cargamos el fichero 
function z=g1(x)
    n=length(x);
    z=zeros(size(x));
    for j=1:n
    %límites de y
        c=x(j)^3;
        d=x(j)^2;
    %integrando f(x,y)
        f=@(y) exp(y./x(j));
        z(j)=integral(f,c,d);
    end
end
    Modificamos el script para llamar a la función 
res=integral(@g1,0.1,0.5);
fprintf('El valor de la integral es: %1.4f\n',res)
    En la ventana de comandos corremos el script 
El valor de la integral es: 0.0333
Comprobamos la segunda integral. Cargamos el fichero 
function z=g1(x)
    n=length(x);
    z=zeros(size(x));
    for j=1:n
    %límites de y
        c=x(j);
        d=2*x(j);
    %integrando f(x,y)
        f=@(y) x(j)^2+y.^3;
        z(j)=integral(f,c,d);
    end
end
    Modificamos el script para llamar a la función 
res=integral(@g2,0,1);
fprintf('El valor de la integral es: %1.4f\n',res)
    En la ventana de comandos corremos el script
El valor de la integral es: 1.0000
Ejemplos en el curso de Física
Flexión de una viga en voladizo (II)
Pandeo de una barra delgada empotrada en un extremo
Fuerza de rozamiento proporcional al cuadrado de la velocidad
Vaciado de un depósito cerrado
Presión producida por la curvatura de una superficie
Función de distribución de Weibull
 
