Integración numérica con MATLAB

Integrales simples

La función integral(function, a, b); devuelve el resultado aproximado de la integral cuando se le pasa la función a integrar en el primer parámetro function, los límites a y b de la integral. Ejemplos:

0 10 ( t 2 +14t+21)·dt = t 3 3 +7 t 2 +21t | 0 10 = 1730 3 =576.6667m

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

a b c d f(x,y)·dy·dx

MATLAB dispone de la función integral2 que realiza esta tarea. Supongamos que queremos calcular la integral doble

1 2 0 3 x 2 y · d y · d x

Creamos un script, en el que definimos la función anónima f(x,y) y llamamos a integral2 pasándole la función, los límites de integración de la variable x, y los límites de integración de la variable y.

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

a b c(x) d(x) f(x,y)dy·dx

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.

a b ( c(x) d(x) f(x,y)·dy ) dx= a b g(x)·dx

Utilizamos la función MATLAB integral, para resolver ambas integrales, tal como vemos en el siguiente ejemplo.

0 2 x 3 x x 2 +x ( x 3 y 4 +x y 2 )·dy·dx = 0 2 ( x 3 x x 2 +x ( x 3 y 4 +x y 2 )·dy ) dx

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 integral para resolver la integral de g(x)

res=integral(@g,0,2);
fprintf('El valor de la integral es: %3.3f\n',res)

En la ventana de comandos corremos el script integral_doble1

   El valor de la integral es: 961.181

Ejemplos

Comprobar los resultados de las integrales siguientes:

0.1 0.5 x 3 x 2 e y/x dy·dx0.0333 0 1 x 2x ( x 2 + y 3 )dy·dx 1.0001

Comprobamos la primera integral. Cargamos el fichero g.m en el editor y lo guardamos (File/Save As..) con nombre g1.m. Cambiamos el nombre de la función g1 en vez de g. Modificamos los límites c y d de la integral así como la definición de la función f(x,y)

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 integral para calcular la integral de g1(x) entre los nuevos límites a y b

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 integral_doble2

El valor de la integral es: 0.0333

Comprobamos la segunda integral. Cargamos el fichero g1.m en el editor y lo guardamos (File/Save As..) con nombre g2.m. Cambiamos el nombre de la función g2 en vez de g1. Modificamos los límites c y d de la integral así como la definición de la función f(x,y)

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 integral para calcular la integral de g2(x) entre los nuevos límites de la integral a y b

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

El ascensor espacial

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

La radiación del cuerpo negro

Función de distribución de Weibull

Función de distribución de Rayleigh