Lo que he aprendido: integración numérica con Gauss

Esto no es técnicamente algo que haya aprendido ahora, llevo años integrando numéricamente como manda Gauss, pero creo que siempre lo había hecho o bien a mano o bien buscando los puntos de integración y los pesos en las tablas. La novedad aquí es que uso los polinomios de Legendre para calcular cuáles son los puntos de integración y los pesos correspondientes para cualquier número de puntos.

Recordemos que en el método de Gauss:

  • Los puntos de integración corresponden a las raíces del polinomio de Legendre de orden n
  • Los pesos se obtienen a partir del valor de la derivada en los puntos de integración
  • El método se define para el intervalo [-1,1], para cualquier otro intervalo tenemos que hacer un cambio de variable y ajustar el valor final de la integral

Con mis reducidas habilidades julianas, me queda una función con esta pinta (está in inglis porque la migré de Matlab de otra función que es parte del código de elementos finitos que estoy fabricando en el curro):

# Pkg.add("Polynomials")
using Polynomials

"""
legendre(x)

Create Legendre polynomial of degree `x`
"""
legendre(x) = 1//(2^x*factorial(x))*polyder(Poly([-1, 0, 1])^x,x)

"""
gaussGeneral(fun, points, interval)

Compute numerical integral of function `fun` using number
of `points` in `interval`
"""
function gaussGeneral(fun, points, interval)

a = interval[1]
b = interval[2]

# Compute roots of Legendre polynomial of order 'points'
p = roots(legendre(points))
d = polyder(legendre(points))

# Weights
w = 2*ones(length(p),1)./((ones(length(p),1)-p.^2).*d(p).^2)

# Variable change
p = (a+b)/2*ones(length(p),1) + (b-a)/2*p

# Sum of weighted values
return (b-a)/2*w'*fun(p)

end

Lo chachi de definir así la integral numérica es que me permite calcular la integral de la función con diferente número de puntos y comparar con el valor analítico. Yo soy ingeniera y no me creo las típicas frases la cuadratura de Gauss es un método para integrar de manera exacta polinimios hasta de orden 2n-1 donde n es el número de puntos utilizados si no lo pruebo.

Y eso era. Dejo una pregunta para los sabios: ¿me podéis recomendar una referencia en la que se explique de manera entendible por qué las raíces del polinomio de Legendre son los puntos óptimos de integración? Creo que nunca lo llegué a comprender.


Suena ahora:

4 pensamientos en “Lo que he aprendido: integración numérica con Gauss

  1. Pingback: Lo que he aprendido: leer matrices MMF en Julia | Onda Hostil

  2. Pingback: Compilación: los artículos del 2018 | Onda Hostil

¡Opina sin miedo! (Puedes usar Markdown)

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google photo

Estás comentando usando tu cuenta de Google. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s