Integración numérica de $\int_0^1 \log(x) dx$

Question

¿Tienen alguna idea de cómo manejar algo como

$\int_0^1 log(x) dx$ numéricamente en Matlab (sólo me interesa la parte real, por cierto)?

He utilizado Quad, Quadl, "integral" etc todos ellos encuentran intfy en el origen.

He intentado cambiar el límite inferior a, por ejemplo, 0,00001, y entonces el resultado es finito, pero 0,00001 no es lo suficientemente bueno y si utilizo algo más pequeño que 0,00001 el resultado vuelve a ser infimo.

Debo añadir que el integrando no es un simple log(x) sino una expresión larga que va como/(expansión de Taylor) log(x) cerca de x=0. He definido este integrando en un archivo de función separado.

Ahora $\lim_{x->0} J(x) = -\infty$ donde $J$ es el integrando, PERO la integral debe ser finita en x=0.

¿Alguna idea de qué hacer?

Gracias.

Answer 1

Probablemente lo más fácil sea utilizar el resultado de la forma cerrada $\int_0^1 \log(x) \mathrm{d}x = -\int_0^\infty t \mathrm{e}^{-t} \mathrm{d}t = -1$ y aplicar las cuadraturas a $J(x)-\log(x)$ : $$ \int_0^1 J(x) \mathrm{d}x = -1 +\int_0^1 \left(J(x)-\log(x)\right)\mathrm{d}x $$

Answer 2

Pruebe los métodos abiertos de Newton-Cotes (el más sencillo, la regla del punto medio), o la cuadratura gaussiana. Estos métodos no necesitan ser evaluados en los límites del intervalo.

Answer 3

Para calcular numéricamente una integral de la forma

$$I = \int_0^1 \log(x) f(x) dx$$

donde $f(x)$ se supone que es regular en ambos extremos, se puede pasar por el formalismo de la cuadratura gaussiana como explicada aquí . Se pueden encontrar los polinomios ortogonales al elegir $\log(x)$ como función de peso utilizando la relación de recursión de 3 términos dada al final de esta sección . Las siguientes reglas de cuadratura pueden derivarse utilizando los métodos allí explicados. Si denotamos el $n$ aproximación de cuadratura gaussiana de orden a $I$ como $I_n$ entonces tenemos..:

$$I_n = \sum_{k=1}^{n}w_{n,k}f\left(x_{n,k}\right)$$

las abscisas vienen dadas por:

$$ \begin{split} x_{1,1} &= \frac{1}{4}\\ x_{2,1} &= \frac{1}{42} \left(15-\sqrt{106}\right)\\ x_{2,2} &= \frac{1}{42} \left(15+\sqrt{106}\right)\\ x_{3,1} &= 0.0638907930873254049961166031363\ldots\\ x_{3,2} &= 0.368997063715618765546197645857\ldots\\ x_{3,3} &= 0.766880303938941455423682659817\ldots\\ x_{4,1}&=0.0414484801993832208033213101564\ldots\\ x_{4,2}&=0.245274914320602251939675759523\ldots\\ x_{4,3}&=0.556165453560275837180184354376\ldots\\ x_{4,4}&=0.848982394532985174647849188085\ldots\\ x_{5,1}&=0.0291344721519720533037267621154\ldots\\ x_{5,2}&=0.173977213320897628701139710829\ldots\\ x_{5,3}&=0.411702520284902043174931924646\ldots\\ x_{5,4}&=0.677314174582820380701802667998\ldots\\ x_{5,5}&=0.894771361031008283638886204455\ldots \end{split} $$

Los pesos son:

$$ \begin{split} w_{1,1} &= -1\\ w_{2,1}&= -\frac{1}{2}-\frac{9}{4 \sqrt{106}}\\ w_{2,2}&= -\frac{1}{2}+\frac{9}{4 \sqrt{106}}\\ w_{3,1}&=-0.51340455223236332512930049757\ldots\\ w_{3,2}&=-0.39198004120148755480628718097\ldots\\ w_{3,3}&=-0.09461540656614912006441232147\ldots\\ w_{4,1}&=-0.38346406814513512485004652234\ldots\\ w_{4,2}&=-0.3868753177747626273360082346\ldots\\ w_{4,3}&=-0.1904351269501424153613600145\ldots\\ w_{4,4}&=-0.03922548712995983245258522856\ldots\\ w_{5,1}&=-0.29789347178289445727225787789\ldots\\ w_{5,2}&=-0.3497762265132241803750718703\ldots\\ w_{5,3}&=-0.234488290044052418886906858\ldots\\ w_{5,4}&=-0.098930459516633146976180711\ldots\\ w_{5,5}&=-0.0189115521431957964895826824\ldots \end{split} $$

Por ejemplo, la diferencia entre $I_5$ para $f(x) = \exp(x)$ y el valor exacto es aproximadamente $3\times 10^{-13}$ .

Answer 4

Sólo voy a lanzar esto a Wolfram por un segundo...

Ok, ahora aquí hay algunas pruebas..terminaras con errores de lo contrario..

min=1e-10;
max=1;
delta = 1e-2;
x = min:delta:max;
y = @(x) (log(x));
I = trapz(x,y(x))

I =

   -1.0829

min=1e-12;
max=1;
delta = 1e-2;
x = min:delta:max;
y = @(x) (log(x));
I = trapz(x,y(x))

I =

   -1.1059

hmmm..

min=1e-12;
max=1;
delta = 1e-6;
x = min:delta:max;
y = @(x) (log(x));
I = trapz(x,y(x))

arreglado...hecho..

I =

   -1.0000