Cuadratura gaussiana -Derivación de una fórmula-

Question

El siguiente es un ejercicio de la sección de problemas del capítulo de cuadratura gaussiana.

El teorema:

Derivar una fórmula de la forma $$\int_{a}^{b} f(x)dx \approx w_0f(x_0) + w_1f(x_1) + w_2f'(x_2) + w_3f'(x_3)$$

La verdad es que no sé ni cómo empezar esto. Tengo el teorema de la cuadratura de Gauss, pero... supongo que no sé cómo usarlo. ¿Alguien puede mostrarme cómo hacerlo?

Tenemos cuatro pesos, $w_0, w_1,w_2,w_3$ y cuatro nodos $x_0, x_1,x_2,x_3$ . Así que nuestro polinomio $q(x)$ será de grado $4$ . (¿El número de pesos determinará siempre el grado del polinomio?). Entonces, a partir de $$ \int_{a}^{b} x^kq(x)dx = 0$$ donde $0 \leq k \leq n$ y el grado de $q(x)$ es $n+1$ podemos escribir $$\int q(x)dx = \int xq(x)dx = \int x^2q(x) dx = 0$$ ¿verdad?

A partir de aquí no sé qué hacer.

Answer 1

Veamos primero $\int_{-1}^{1} f(t) \ dt$ .

Queremos que la fórmula evalúe $\int_{-1}^{1} dt, \int_{-1}^{1} t \ dt, \int_{-1}^{1} t^{2} \ dt, \int_{-1}^{1} t^{3} \ dt, \int_{-1}^{1} t^{4} \ dt, \int_{-1}^{1} t^{5} \ dt, \int_{-1}^{1} t^{6} \ dt $ y $\int_{-1}^{1} t^{7} \ dt$ exactamente.

Esto nos lleva al siguiente sistema de ecuaciones:

$$ 2 = \omega_{0} + \omega_{1}$$

$$ 0 = \omega_{0} x_{0} + \omega_{1} x_{1}+ \omega_{2}+\omega_{3}$$

$$ \frac{2}{3} = \omega_{0} x_{0}^{2} + \omega_{1} x_{1}^{2} + 2 \omega_{2} x_{2} + 2 \omega_{3} x_{3}$$

$$ 0 = \omega_{0} x_{0}^{3} + \omega_{1} x_{1}^{3}+ 3 \omega_{2} x_{2}^{2} + 3 \omega_{3} x_{3}^{2} $$

$$ \frac{2}{5} = \omega_{0} x_{0}^{4} + \omega_{1} x_{1}^{4}+ 4 \omega_{2} x_{2}^{3} + 4 \omega_{3} x_{3}^{3} $$

$$0 = \omega_{0} x_{0}^{5} + \omega_{1} x_{1}^{5} + 5 \omega_{2} x_{2}^{4} + 5 \omega_{3} x_{3}^{4} $$

$$\frac{2}{7} = \omega_{0} x_{0}^{6} + \omega_{1} x_{1}^{6} + 6 \omega_{2} x_{2}^{5} + 6\omega_{3} x_{3}^{5} $$

$$ 0 = \omega_{0} x_{0}^{7} + \omega_{1} x_{1}^{7} + 7 \omega_{2} x_{2}^{6} + 7\omega_{3} x_{3}^{6}$$

Resuelve el sistema utilizando un solucionador numérico.

Entonces utiliza el hecho de que $$\int_{a}^{b} f(x) \ dx = \frac{b-a}{2} \int_{-1}^{1} f \Big(a+(1+t)\frac{b-a}{2} \Big) \ dt$$

Answer 2

Es una pregunta antigua y ya tiene una respuesta aceptada, pero no está claro a qué respuesta llegó finalmente el autor de la pregunta original. El problema es realmente difícil, tanto más cuanto que el libro en cuestión, Numerical Mathematics and Computing Sixth Edition, Ward Cheney y David Kincaid, Thomson Brooks/Cole, Redmont, CA 2008 tiene una presentación bastante somera sobre la cuadratura gaussiana. De hecho, en el capítulo 6, problema 8, p. 240, se pide al lector que derive una fórmula de la forma $$\int_a^bf(x)dx\approx w_0f(a)+w_1f(b)+w_2f^{\prime}(a)+w_3f^{\prime}(b)$$ Obsérvese que en el enunciado real del problema los puntos de muestra ya están especificados y el problema realmente pregunta por la Regla Trapezoidal con términos de corrección. Se sabe que la fórmula no puede ser exacta para polinomios de grado $4$ porque $$\int_a^b(x-a)^2(b-x)^2dx=\frac{(b-a)^5}{30}$$ mientras que cualquier fórmula predice un valor de $0$ para esta integral. La solución conocida es $$\int_a^bf(x)dx\approx\frac{b-a}2\left[f(a)+f(b)\right]+\frac{(b-a)^2}{12}\left[f^{\prime}(a)-f^{\prime}(b)\right]$$ Que se puede obtener aplicando la fórmula a $f(x)\in\{1,x,x^2,x^3\}$ resultando la ecuación matricial $$\begin{bmatrix}1&1&0&0\\ a&b&1&1\\ a^2&b^2&2a&2b\\ a^3&b^3&3a^2&3b^2\end{bmatrix} \begin{bmatrix}w_0\\w_1\\w_2\\w_3\end{bmatrix} =\begin{bmatrix}b-a\\\frac12(b^2-a^2)\\\frac13(b^3-a^3)\\\frac14(b^4-a^4)\end{bmatrix}$$ Pero la pregunta real que se hace es mucho más interesante. Escribámosla en estilo Matlab con $1$ -vectores: $$\int_{-1}^1f(x)dx\approx w_1f(x_1)+w_2f(x_2)+w_3f^{\prime}(x_3)+w_4f^{\prime}(x_4)$$ También tomamos un intervalo fijo de integración porque podríamos mapear cualquier intervalo finito de integración en él mediante una transformación lineal. Una aproximación podría ser asumir que la solución tiene simetría tal que $x_1=-x_2$ , $w_1=w_2$ , $x_3=-x_4$ y $w_3=-w_4$ . Esto produce automáticamente la respuesta correcta para las funciones Impares, por lo que sólo tenemos que aplicar la fórmula para $f(x)\in\{1,x^2,x^4,x^6\}$ para obtener una fórmula exacta para $f(x)\in\mathcal{P}_7$ . $8^{th}$ Los polinomios de grado están descartados ya que la fórmula predeciría una integral de $0$ para la función de prueba no negativa $f(x)=(x-x_1)^2(x-x_2)^2(x-x_3)^2(x-x_4)^2$ . Llegamos a la ecuación matricial $$\begin{bmatrix}2&0\\ 2x_1^2&4x_3\\ 2x_1^4&8x_3^3\\ 2x_1^6&12x_3^5\end{bmatrix} \begin{bmatrix}w_1\\w_3\end{bmatrix} =\begin{bmatrix}2\\2/3\\2/5\\2/7\end{bmatrix}$$ La eliminación gaussiana de abajo hacia arriba da como resultado $$\begin{bmatrix}2&0\\ 0&4x_3\\ 0&(8x_3^2-4x_1^2)x_3\\ 0&(12x_3^2-8x_1^2)x_3^3\end{bmatrix} \begin{bmatrix}w_1\\w_3\end{bmatrix} =\begin{bmatrix}2\\\frac23-2x_1^2\\\frac25-\frac23x_1^2\\\frac27-\frac25x_1^2\end{bmatrix}$$ Y la eliminación hacia abajo produce $$\begin{bmatrix}2&0\\0&4x_3\\0&0\\0&0\end{bmatrix}\begin{bmatrix}w_1\\w_3\end{bmatrix}=\begin{bmatrix}2\\\frac23-2x_1^2\\\frac25-2x_1^4-\left(\frac23-2x_1^2)\cdot2x_3^2\\\frac27-\frac25x_1^2+\left(\frac43x_1^2-4x_1^4\right)x_3^2-\left(2-6x_1^2\right)x_3^4\right)\end{bmatrix}$$ Desde la primera $3$ filas obtenemos $$\begin{align}w_1&=1\\w_3&=\frac{\frac23-2x_1^2}{4x_3}\\x_3^2&=\frac{\frac25-2x_1^4}{2\left(\frac23-2x_1^2\right)}\end{align}$$ Y sustituyendo este último valor en el $4^{th}$ fila obtenemos $$x_1^8-\frac43x_1^6+\frac65x_1^4-\frac47x_1^2+\frac{37}{525}=0$$ Si dejamos que $$A = \sqrt{\frac{16\sqrt[3]{14}}{105}-\frac{16}{45}}$$ Estos factores a $$\left[x_1^4+\left(A-\frac23\right)x_1^2+\frac12A^2-\frac13A+\frac{17}{45}+\frac{32}{945A}\right]\left[x_1^4-\left(A+\frac23\right)x_1^2+\frac12A^2+\frac13A+\frac{17}{45}-\frac{32}{945A}\right]=0$$ El primer factor tiene raíces complejas pero el segundo tiene todas las raíces reales, por lo que obtenemos $2$ posibles respuestas mediante la fórmula cuadrática: $$\begin{array}{rrl}x_1=&-x_2&=0.771864648647792\\ x_3=&-x_4&=0.543327162125273\\ w_1=&w_2&=1\\ w_3=&-w_4&=-0.241513512293663\end{array}$$ Y $$\begin{array}{rrl}x_1=&-x_2&=0.423175690707864\\ x_3=&-x_4&=0.737785505959038\\ w_1=&w_2&=1\\ w_3=&-w_4&=0.104539643894699\end{array}$$ Pero el verdadero problema es: ¿son estas todas las soluciones? ¿Es posible que haya soluciones asimétricas? En ese caso tenemos que probar todas $7$ funciones $f(x)\in\{1,x,x^2,x^3,x^4,x^5,x^6,x^7\}$ y nuestra ecuación matricial es ahora $$\begin{bmatrix}1&1&0&0\\ x_1&x_2&1&1\\ x_1^2&x_2^2&2x_3&2x_4\\ x_1^3&x_2^3&3x_3^2&3x_4^2\\ x_1^4&x_2^4&4x_3^3&4x_4^3\\ x_1^5&x_2^5&5x_3^4&5x_4^4\\ x_1^6&x_2^6&6x_3^5&6x_4^5\\ x_1^7&x_2^7&7x_3^6&7x_4^6\\ \end{bmatrix} \begin{bmatrix}w_1\\w_2\\w_3\\w_4\end{bmatrix} =\begin{bmatrix}2\\0\\2/3\\0\\2/5\\0\\2/7\\0\end{bmatrix}$$ Como recomienda @Random Variable. Ahora este sistema era intratable para mí, así que escribí un programa de Matlab:

% asym.m

clear all;
close all;
x = [-0.75 -0.25]';
y = [0.25 0.75]';
A = [ones(1,2) zeros(1,2);
    x' ones(1,2);
    x'.^2 2*y';
    x'.^3 3*y'.^2];
b = [2 0 2/3 0]';
w = A\b;
for k = 1:20,
    A = [ones(1,2) zeros(1,2);
        x' ones(1,2);
        x'.^2 2*y';
        x'.^3 3*y'.^2;
        x'.^4 4*y'.^3;
        x'.^5 5*y'.^4;
        x'.^6 6*y'.^5;
        x'.^7 7*y'.^6];
    b = [2 0 2/3 0 2/5 0 2/7 0]';
    res = A*w-b;
    F = [zeros(1,2) zeros(1,2);
        ones(1,2) zeros(1,2);
        2*x' 2*ones(1,2);
        3*x'.^2 6*y';
        4*x'.^3 12*y'.^2;
        5*x'.^4 20*y'.^3;
        6*x'.^5 30*y'.^4;
        7*x'.^6 42*y'.^5];
    err = [F*diag(w) A]\res;
    x = x-err(1:2);
    y = y-err(3:4);
    w = w-err(5:8);
    if norm(err) < 1.0e-12,
        break;
    end
end
err = norm(err)
format long;
x
y
w
format;

Si empezamos con x = [-0.75 0.75]'; , y = [-0.25 0.25]'; obtenemos la primera solución anterior, mientras que si empezamos con x = [-0.25 0.25]'; , y = [-0.75 0.75]'; llegamos a la segunda, pero buscando una solución asimétrica a partir de x = [-0.75 -0.25]'; , y = [0.25 0.75]'; resultó una solución asimétrica $$\begin{align}x_1&=-0.840793653320079\\ x_2&=-0.140661084689637\\ x_3&=0.194506019228943\\ x_4&=0.716583972644591\\ w_1&=0.402433842928781\\ w_2&=1.597566157071219\\ w_3&=0.437121762781648\\ w_4&=0.125957446751173 \end{align}$$ Así que no estoy seguro de haber encontrado todas las soluciones reales y no creo haber encontrado todas las soluciones complejas al problema. Por lo tanto, la respuesta anterior debe considerarse sólo una respuesta parcial que deja abiertas algunas preguntas más.