Integración numérica

Antes de introducir al lector en al cálculo del área comprendida entre una función y=f(x) y el eje X en un intervalo dado [a,b]. Vamos a ver como se calcula el área de un polígono como suma de las áreas de trapecios.

Areas

El área del trapecio de la figura de la izquierda es la suma de dos áreas: un rectángulo y un triángulo

(x2-x1y1+(x2-x1)·(y2-y1)/2=(x2-x1)·(y2+y1)/2

El área de la figura de la derecha formada por los puntos (x1,y1), (x2,y2) y (x3,y3) es

(x2-x1)·(y2+y1)/2+(x3-x2)·(y3+y2)/2

Para n puntos, el área es

A= 1 2 i=1 n1 ( x i+1 x i ) ( y i+1 + y i )

Supongamos que en una experiencia hemos medido la velocidad de un móvil en función del tiempo,

t(s) 1 2 3 4 5 6 7 8 9 10
v(m/s) 5.0 6.0 5.5 7.0 8.3 7.6 6.2 6.1 7.0 5.7

calculamos el desplazamiento aproximado del móvil sumando áreas y utilizando la función trapz.

x=1:10;
y=[5.0,6.0,5.5,7.0,8.3,7.6,6.2,6.1,7.0,5.7];
area=0;
for i=1:length(x)-1
     area=area+(x(i+1)-x(i))*(y(i)+y(i+1))/2;
end
disp(['Area: ', num2str(area)]);
area=trapz(x,y);
disp(['Area: ', num2str(area)]);
Area 1: 59.05
Area 2: 59.05

Area de un polígono

El área del triángulo formado por tres puntos (x1,y1), (x2,y2) y (x3,y3) es

(x2-x1)·(y2+y1)/2+(x3-x2)·(y3+y2)/2-(x3-x1)·(y3+y1)/2=(x2-x1)·(y2+y1)/2+(x3-x2)·(y3+y2)/2+(x1-x3)·(y1+y3)/2

Para un polígono de n lados, formado por los puntos (x1,y1), (x2,y2) ... (xn,yn) el área se calcula mediante la fórmula

A= 1 2 ( ( x 1 x n )( y 1 + y n )+ 1 n1 ( x i+1 x i )( y i+1 + y i ) )

Para cerrar la figura poligonal, añadimos el vértice n+1, que coincide con el primer vértice (x1,y1), (x2,y2) ... (xn,yn), (x1,y1). El cálculo del área se simplifica

A= 1 2 ( 1 n ( x i+1 x i )( y i+1 + y i ) )

Consideremos primero un polígono regular, un pentágono. Calculamos el área mediante la función MATLAB polyarea y también, mediante la fórmula que hemos obtenido anteriormente

n=5; %número de lados
radio=3; %radio de la circunferencia 
x=radio*cos((0:n)*2*pi/n);
y=radio*sin((0:n)*2*pi/n);
%las coordenadas del último vértice son las del primer vértice
hold on
fill(x,y,'y')
plot(x,y,'-ro', 'markersize',4,'markeredgecolor','r','markerfacecolor','r')
hold off
xlabel('x')
ylabel('y')
axis equal
area=0;
for i=1:n
    area=area+(x(i+1)-x(i))*(y(i+1)+y(i))/2;
end
area=abs(area)
title (['Area del polígono regular: ' num2str(polyarea(x,y))])

area =   21.3988

El área de un polígono regular de n inscrito una circunferencia de radio r es la suma del área de n triángulos isósceles de vértice θ=2π/n.

A=n·(a·h)/2

Como vemos en la figura

a=2·r·sin( θ 2 )h=r·cos( θ 2 ) A=n ah 2 = 1 2 n r 2 sinθ= 1 2 n r 2 sin( 2π n )

>> area=n*radio^2*sin(2*pi/n)/2
area =   21.3988

Consideramos ahora, un polígono cualesquiera

x = [15,13,5,4,0,0,1,4,5,7,9,12,14,13,15,15];
y = [1,0,1,2,4,4,7,8,10,11,10,9,8,4,4,1];
%las coordenadas del último vértice son las del primer vértice
hold on
fill(x,y,'y')
plot(x,y,'-ro', 'markersize',4,'markeredgecolor','r','markerfacecolor','r')
hold off
xlabel('x')
ylabel('y')
axis equal
area=0;
for i=1:length(x)-1
    area=area+(x(i+1)-x(i))*(y(i+1)+y(i))/2;
end
title (['Area del polígono: ' num2str(area])

La función MATLAB polyarea solamente precisa los n vértices del polígono, dos vectores x e y de dimensión n

>> x = [15,13,5,4,0,0,1,4,5,7,9,12,14,13,15];
>> y = [1,0,1,2,4,4,7,8,10,11,10,9,8,4,4];
>> polyarea(x,y)
ans =   108

Método del punto medio

El método del punto medio es muy fácil de entender, ya que aproxima la integral definida

a b f(x)dx

a la suma del área de pequeños rectángulos. Como se ve en la figura el área de cada rectángulo es f(xm) ·Δxi, donde xm=(xi+1+xi)/2 es la posición intermedia entre xi y xi+1.

velocidad

Definimos la función integral_2 para que realice las siguientes tareas:

  1. Calcula la posición intermedia xm=(xi+1+xi)/2 entre xi y xi+1
  2. Calcula el área del rectángulo f(xm) ·Δxi, donde Δxi=xi+1-xi
  3. Calcula el área total aproximada en el intervalo a, b.
function suma=integral_2(f,x)
    suma=0;
    for i=1:length(x)-1
        xm=(x(i+1)+x(i))/2;
        suma=suma+f(xm)*(x(i+1)-x(i));
    end
end

En el script llamamos a integral_2

t0=input('tiempo inicial, t0: ');
tf=input('tiempo final, tf: ');
n=input('número intervalos, n: ');

vel=@(t) -t.^2+14*t+21;  %función velocidad
t1=linspace(t0,tf,n+1);
res=integral_2(vel,t1); %calcula la integral
fprintf('El desplazamiento es: %3.2f\n',res)

En la ventana de comandos corremos el script

tiempo inicial, t0: 0
tiempo final, tf: 10
número intervalos, n: 10
El desplazamiento es: 577.50

Fórmula del trapecio

Para calcular la integral definida de la función f(x) en el intervalo comprendio entre x0 y x, dividimos este intervalo en pequeños intervalos de longitud h=xi+1-xi. Sustituimos la función por la recta que une los puntos (xi, yi) y (xi+1, yi+1). El área sombreada en la figura es la suma del área de un rectángulo más el área de un triángulo, vale

trapecio

h· y i + 1 2 h·( y i+1 y i )= 1 2 h·( y i+1 + y i )

El área total aproximada bajo la curva es

h 2 ( y 1 + y 2 )+ h 2 ( y 2 + y 3 )+ h 2 ( y 3 + y 4 )+...+ h 2 ( y n2 + y n1 )+ h 2 ( y n1 + y n )= h 2 ( y 1 +2 y 2 +2 y 3 +2 y 4 +...+2 y n2 +2 y n1 + y n )

Definimos en la función trapecio el procedimiento. Donde f es la función integrando y x el vector de n+1 datos comprendidos entre la abscisa inicial x0 y la final xf ambas incluidas, si el número de intervalos es n.

function sum=trapecio(f,x0, xf, nInterv)
    %el número de puntos es n intervalos más uno
    x=linspace(x0,xf,nInterv+1);
    n=length(x);    
    sum=(f(x(1))+f(x(n)))/2;
    for i=2:n-1
        sum=sum+f(x(i));
    end
    sum=sum*(x(2)-x(1));
end

Dada la función integrando f(x), calculamos la aproximación a la integral definida en el intervalo comprendido entre x0 y x.

x0=input('abscisa inicial, x0: ');
xf=input('abscisa final, xf: ');
n=input('número intervalos, n: ');

deriv=@(x) -x^2+14*x+21;  %definición del integrando 
res=trapecio(deriv,x0,xf,n); %calcula la integral
fprintf('El valor aproximado de la integral: %3.2f\n',res)

En la ventana de comandos corremos el script

abscisa inicial, x0: 0
abscisa final, xf: 10
número intervalos, n: 10
El valor aproximado de la integral: 575.00

El método de Simpson (1/3)

En este procedimiento, se toma el intervalo de anchura 2h, comprendido entre xi y xi+2, y se sustituye la función f(x) por la parábola que pasa por tres puntos (xi, yi), (xi+1, yi+1), y (xi+2, yi+2).

simpson

Calculamos la contribución a la integral del primer intervalo (x0, x0+2h) y generalizaremos para el resto de los intervalos.

La ecuación de la parábola y=ax2+bx+c que pasa por los puntos (x0, y0), (x0+h, y1), (x0+2h, y2) es

y 0 =a x 0 2 +b x 0 +c y 1 =a ( x 0 +h ) 2 +b( x 0 +h)+c y 2 =a ( x 0 +2h ) 2 +b( x 0 +2h)+c

Este sistema de tres ecuaciones con tres incógnitas, se reduce a

y 1 = y 0 +2a x 0 h+a h 2 +bh y 2 = y 0 +4a x 0 h+4a h 2 +2bh

Despejamos el coeficiente a, y 2ax0+b

a= y 2 2 y 1 + y 0 2 h 2 2a x 0 +b= 4 y 1 3 y 0 y 2 2h

Sustituimos la curva por la porción de parábola en el intervalo (x0, x0+2h). La integral vale.

I= x 0 x 0 +2h (a x 2 +bx+c)dx= a 3 ( 6 x 0 2 +12 x 0 h 2 +8 h 3 )+ b 2 ( 4 x 0 h+4 h 2 )+c( 2h )= 2h( a x 0 2 +b x 0 +c )+2( 2a x 0 +b ) h 2 + 8 3 a h 3 = h 3 ( y 0 +4 y 1 + y 2 )

En general, el valor del área aproximada, en el intervalo (xi, xi+2h) sombreada en la figura, es

h 3 ( y i +4 y i+1 + y i+2 )

El área aproximada en el intervalo (a, b) es

a b f(x)·dx h 3 ( y 0 +4 y 1 + y 2 )+ h 3 ( y 2 +4 y 3 + y 4 )+....+ h 3 ( y n2 +4 y n1 + y n )

o bien, agrupando términos

a b f(x)·dx h 3 ( ( y 0 + y n )+4( y 1 + y 3 +.... y n1 )+2( y 2 + y 4 +.... y n2 ) )

El primer paréntesis, contiene la suma de los extremos, el segundo, la suma de los términos de índice impar, y el tercero la suma de los términos de índice par. En el método de Simpson, el número de divisiones n debe de ser par.

Definimos en la función simpson el procedimiento. Donde f es la función integrando y x el vector de n+1 datos comprendidos entre la abscisa inicial x0 y la final xf ambas incluidas, si el número de intervalos es n que tiene que ser un número PAR.

function suma=simpson(f,x0,xf,n)
%n número par de intervalos, n+1 número de puntos en el vector    
    x=linspace(x0,xf,n+1);   
    h=x(2)-x(1);
    suma=f(x(1))+f(x(n+1));
    for i=2:2:n
            suma=suma+4*f(x(i));
        end
        for i=3:2:n-1
            suma=suma+2*f(x(i));
        end
        suma=suma*h/3;
end

Dada la función integrando f(x), calculamos la aproximación a la integral definida en el intervalo comprendido entre x0 y x

x0=input('abscisa inicial, x0: ');
xf=input('abscisa final, xf: ');
n=input('número intervalos (par), n: ');
if rem(n,2)==1
  disp('El número intervalos tiene que ser par ');
  break
end

deriv=@(x) -x^2+14*x+21;  %definición del integrando 
res=simpson(deriv,x0,xf,n); %calcula la integral
fprintf('El valor aproximado de la integral: %3.2f\n',res)

En la ventana de comandos corremos el script

abscisa inicial, x0: 0
abscisa final, xf: 10
número intervalos, n: 10
El valor aproximado de la integral: 576.67

Apreciamos la mejora en el resultado de la integral con el método de Simpson.

Existe una versión del método de Simpson denominada 3/8 en el que se emplea una ecuación cúbica para concectar cuatro puntos.