# Integración numérica con MAXIMA

## Integración numérica

Podemos obtener un  valor aproximado de una integral definida mediante métodos numéricos. 

Estos métodos se emplean, o bien cuando no conocemos la primitiva de la función o bien porque 
sólo conocemos la función en un conjunto de valores discretos.

Recordemos que la integral definida de una función se obtiene como un proceso de aproximación. Construíamos una familia de rectángulos, calculábamos el área de cada rectángulo y sumábamos para obtener la aproximación. Para mejorar esta aproximación se utiliza
un mayor número de rectángulos.

En general estos rectángulos pueden ser de cualquier tamaño, pero por simplicidad, supondremos que el intervalo [a,b] está dividido en  n
partes iguales:

    a= x_0   <  x_1  <  x_2  <  ....  <  x_n=b

donde cada x_k está definido como

    x_k  = a + k*(b-a) / n

### Sumas de Riemann

Para la definición de sumas inferiores y superiores, necesitaríamos el mínimo y máximo de la función f(x) en cada uno de los subintervalo de la partición, no obstante, podemos utilizar la suma de Riemann que nos permite elegir un punto cualquiera dentro de dicho intervalo.

Supongamos que queremos calcular  la siguiente integral, usando sumas de Riemann

'integrate((x+1)*x*(x+3),x,1,3);

integrate((x+1)*x*(x+3),x,1,3);

Definimos los extremos del intervalo, la función y el número de intervalos de la partición

a : 1.0;
b : 3.0;
n : 10;
f(x) := (x+1)*x*(x+3);

Definimos los puntos de la partición

xl : makelist( a + k*(b-a)/n,k,0,n);

xl[11];

Podemos elegir un punto al azar dentro de cada intervalo

s:make_random_state(445)$                                           /* Esta función y la siguiente inician el sistema para la generación de números aleatorios */
set_random_state(s)$
xa : makelist( xl[k] + random(1.0)*(b-a)/n,k,1,n);             /* Con estas instrucciones generamos un punto aleatorio en cada subintervalo */
fp:makelist(f(xa[k]),k,1,n);                                                   /* Aquí se evalúa la función en cada uno de los puntos anteriores */

Y calcular su suma

sum(fp[k]*(b-a)/n,k,1,n);

O bien, podemos elegir el extremo inferior del intervalo

inferior : makelist(f(xl[k]),k,1,n);
sum(inferior[k]*(b-a)/n,k,1,n);

O el extremo superior

superior : makelist(f(xl[k]),k,2,n+1);
sum(superior[k]*(b-a)/n,k,1,n);

wxplot2d((x+1)*x*(x+3),[x,0,3]);

Si aumentamos el número de subintervalo

a : 1.0;
b : 3.0;
n : 20000;
f(x) := (x+1)*x*(x+3);
xl : makelist( a + k*(b-a)/n,k,0,n)$
s:make_random_state(445)$                                        
set_random_state(s)$
xa : makelist( xl[k] + random(1.0)*(b-a)/n,k,1,n)$          
fp:makelist(f(xa[k]),k,1,n)$     
sum(fp[k]*(b-a)/n,k,1,n);
inferior : makelist(f(xl[k]),k,1,n)$
sum(inferior[k]*(b-a)/n,k,1,n);
superior : makelist(f(xl[k]),k,2,n+1)$
sum(superior[k]*(b-a)/n,k,1,n);

200/3.0;

### Métodos rectángulo, trapecio y Simpson.

Regla del rectángulo

Es simplemente considerar las sumas de Riemann asociadas a la partición P, tomando 
como punto el extremo inferior, el superior o el intermedio

Supongamos que queremos calcular la siguiente integral

y supongamos que conocemos los valores de f(x) en un conjunto de puntos

f(x_0), f(x_1), ...., f(x_n)

kill(all)$
'integrate(exp(x)/x,x,1,3);

Creamos la partición en el intervalo [1,3]

para ello hacemos uso del comando makelist

makelist(expr; var; a; b);

donde 

    - expr es la fórmula a emplear
    - var es la variable usada
    - a: extremo inferior
    - b: extremo superior

en este caso los puntos del intervalo son

    x_k  = a + k*(b-a)/n = 1 + k*(3-1)/n = 1 + k*2/n

siendo n el número de puntos que vamos a utilizar, por ejemplo si n=5 tendremos

x_0, x_1, x_2, x_3, x_4, x_5

a:1$
b:3$
d:(b-a)$
n : 5$
p:makelist(a+k*d/n,k,0,n);

Supongamos que usamos el método del Rectángulo, eligiendo el extremo inferior de cada intervalo 

f(x0)*(x1-x0) + f(x1)*(x2-x1) + ... + f(x_(n-1))*(x_n-x_n-1)

Y si son equidistantes cada intervalo tiene longitud (b-a)/n

que podemos sacar factor común para obtener

(b-a)*(f(x0)+f(x1)+....f(x_n-1))/n

f(x):=exp(x)/x;
n:6$
a:1$
b:3$
fp : makelist(f(a+k*(b-a)/n),k,0,n)$
sum(fp[k],k,1,n)*(b-a)/n,numer;

El cálculo de las sumas superiores

Si queremos utilizar el extremo superior del intervalo

(b-a)*(f(x1)+...f(xn))/n

n:6$
fp : makelist(f(a+k*(b-a)/n),k,0,n)$
sum(fp[k],k,2,n+1)*(b-a)/n,numer;

Y si elegimos el punto medio del intervalo, entonces

f((x0+x1)/2) + f( (x1+x2)/2) + ......

n:6$;
fp:makelist(f(a+(b-a)*(2*k-1)/2/n),k,1,n)$
sum(fp[k],k,1,n)*(b-a)/n,numer;

El valor exacto de la integral es:

integrate(exp(x)/x,x,1,3),numer;

Método Trapecio

En lugar de elegir el rectángulo de base [xk-1,xk] se elige como valor aproximado el área del trapecio
que tiene por vértices 

[xk-1,xk,f(xk-1),f(xk)]

que está dado por

AT =( f(xk-1)+f(x))*(xk-x_k-1)/2

de forma que la suma de las áreas de todos los trapecios sería

A:( f(x0)+f(x1))*(x1-x0)/2 + ( f(x2)+f(x1))*(x2-x1)/2 + .... +( f(xn-1)+f(xn))*(xn-x_n-1)/2

y suponiendo equidistantes

A: (b-a)*(f(x0) + 2 f(x1) + .....2*f(xn-1) + f(xn))/n

n:6$;
fp:makelist(f(a+(b-a)*k/n),k,0,n)$
suma1:fp[1]+fp[n+1]$
suma2: sum(fp[k],k,2,n)$
suma2:2*suma2$
suma : suma1 + suma2$
suma : suma*(b-a)/(2*n)$
suma,numer
;

Regla de Simpson

En esta regla se utilizan 3 puntos y dos intervalos. Necesitamos que n sea par y en cada intervalo

x0,x1,x2         ==> f(x0) + 4f(x1) + f(x2)
x2,x3,x4         ==> f(x2) + 4f(x3) + f(x4)
x4,x5,x6
..
xn-2, xn-1,xn

La suma total sería

f(x0) + 4I + 2P + f(xn)

donde 

I: suma de los valores de la función en los nodos impares
P: suma de los valores de la función en los nodos pares

kill(all)$
n : 6$
a : 1$ 
b : 3$ 
h : float((b-a)/n)$
p : makelist(a + k *h, k, 0, n)$
f(x):=exp(x)/x;
q : f(p)$
I : sum(q[2*k],k,1,n/2)$
P: sum(q[2*k+1],k,1,n/2-1)$
suma: float((q[1]+4*I+2*P+q[n+1])*h/3);
integrate(f(x),x,1,3),numer;

### Métodos de MAXIMA

En MAXIMA hay implementados 3 comandos:

quad_quags(f(x),x,a,b)
quad_qagi(f(x),x,a,b)
romberg(f(x),x,a,b)

que se utilizan para la integración numérica

quad_qags

devuelve 4 valores

 - El valor aproximado de la integral
 - Una cota del error
- El número de veces que se ha evaluado la función del integrando
- El código de error que puede tomar estos valores
    ‘0’ Si no ha habido problemas;
    ‘1’ Si se utilizaron demasiados intervalos;
    ‘2’ Si se encontró un número excesivo de errores de redondeo;
    ‘3’ Si el integrando ha tenido un comportamiento extraño frente a la integración;
    ‘4’ Fallo de convergencia;
    ‘5’ La integral es probablemente divergente o de convergencia lenta;
    ‘6’ Si los argumentos de entrada no son válidos.

f(x):=exp(x)/x;
'integrate(f(x),x,1,3);
quad_qags(f(x),x,1,3);

La función romberg sólo devuelve el valor aproximado de la integral

romberg(f(x),x,1,3);

INTEGRAGIÓN IMPROPIA:

Las dos funciones anteriores se utilizan para intervalos finitos. 

Para integrales impropias de primera especie utilizamos quad_qagi

f(x):=cos(x)/x^2;
'integrate(f(x),x,1,inf);

Función quad_qags

quad_qags(f(x),x,1,inf);

Función romberg

romberg(f(x),x,1,inf);

Función quad_qagi

quad_qagi(f(x),x,1,inf);

INTEGRACIÓN DE SEGUNDA ESPECIE

g(x):=sin(1/x)/x;

'integrate(g(x),x,0,1);

Función quad_qags

quad_qags(g(x),x,0,1);

Función romberg

romberg(g(x),x,0,1);

trigreduce(%);