Gaussian Quadrature -Deriving a Formula-

Question

The following is an exercise in the problem section of the Gaussian Quadrature chapter.

The theorem:

Derive a formula of the form $$\int_{a}^{b} f(x)dx \approx w_0f(x_0) + w_1f(x_1) + w_2f'(x_2) + w_3f'(x_3)$$

I don't really know how to even start this. I have the Gaussian Quadrature theorem, but...I suppose I don't know how to use it. Can anyone show me how to do this?

We have four weights, $w_0, w_1,w_2,w_3$ and four nodes $x_0, x_1,x_2,x_3$. So our polynomial $q(x)$ will be of degree $4$. (Will the number of weights always determine the degree of the polynomial?). So from $$ \int_{a}^{b} x^kq(x)dx = 0$$ where $0 \leq k \leq n$ and the degree of $q(x)$ is $n+1$, we can write $$\int q(x)dx = \int xq(x)dx = \int x^2q(x) dx = 0$$ right?

From here I don't know what to do.

Answer 1

Let's look first at $\int_{-1}^{1} f(t) \ dt$.

We want the formula to evaluate $\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 $ and $\int_{-1}^{1} t^{7} \ dt$ exactly.

That leads to the following system of equations:

$$ 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}$$

Solve the system using a numerical solver.

Then use the the fact that $$\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

Old question and already has an accepted answer, but it's not clear exactly what answer the original questioner eventually arrived at. The problem is really difficult, the more so given that the book in question, Numerical Mathematics and Computing Sixth Edition, Ward Cheney and David Kincaid, Thomson Brooks/Cole, Redmont, CA 2008 has a rather cursory presentation on Gaussian quadrature. Indeed we find in Chapter 6, Problem 8, p. 240 that the reader is tasked with deriving a formula of the form $$\int_a^bf(x)dx\approx w_0f(a)+w_1f(b)+w_2f^{\prime}(a)+w_3f^{\prime}(b)$$ Note that in the actual problem statement the sample points are already specified and the problem is really asking about the Trapezoidal Rule with correction terms. It is known that the formula can't be exact for polynomials of degree $4$ because $$\int_a^b(x-a)^2(b-x)^2dx=\frac{(b-a)^5}{30}$$ whereas any formula will predict a value of $0$ for this integral. The known solution is $$\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]$$ Which you could get by applying the formula to $f(x)\in\{1,x,x^2,x^3\}$ resulting in the matrix equation $$\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}$$ But the actual question asked is much more interesting. Let's write it in Matlab style with $1$-based vectors: $$\int_{-1}^1f(x)dx\approx w_1f(x_1)+w_2f(x_2)+w_3f^{\prime}(x_3)+w_4f^{\prime}(x_4)$$ Also we take a fixed interval of integration because we could map any finite interval of integration onto it via a linear transformation. One approach might be to assume the solution has symmetry such that $x_1=-x_2$, $w_1=w_2$, $x_3=-x_4$, and $w_3=-w_4$. This automatically produces the right answer for odd functions so we only need apply the formula for $f(x)\in\{1,x^2,x^4,x^6\}$ to get a formula exact for $f(x)\in\mathcal{P}_7$. $8^{th}$ degree polynomials are out of the question as the formula would predict an integral of $0$ for the nonnegative trial function $f(x)=(x-x_1)^2(x-x_2)^2(x-x_3)^2(x-x_4)^2$. We arrive at the matrix equation $$\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}$$ Gaussian elimination from the bottom up yields $$\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}$$ And downward elimination produces $$\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}$$ From the first $3$ rows we get $$\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}$$ And substituting the latter value into the $4^{th}$ row we get $$x_1^8-\frac43x_1^6+\frac65x_1^4-\frac47x_1^2+\frac{37}{525}=0$$ If we let $$A = \sqrt{\frac{16\sqrt[3]{14}}{105}-\frac{16}{45}}$$ This factors to $$\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$$ The first factor has complex roots but the second has all real roots, so we get $2$ possible answers by the quadratic formula: $$\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}$$ And $$\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}$$ But the real problem is: are these all the solutions? Is it possible that there are asymmetric solutions? In that case we have to test all $7$ functions $f(x)\in\{1,x,x^2,x^3,x^4,x^5,x^6,x^7\}$ and our matrix equation is now $$\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}$$ As recommended by @Random Variable. Now this system was intractable for me so I wrote up a Matlab program:

% 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;

If we started off with x = [-0.75 0.75]';, y = [-0.25 0.25]'; we got the first solution above, while if we started with x = [-0.25 0.25]';, y = [-0.75 0.75]'; we reached the second, but searching for an asymmetric solution starting from x = [-0.75 -0.25]';, y = [0.25 0.75]'; turned up an asymmetric solution $$\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}$$ So I'm not sure that I've found all the real solutions and don't think I've found all the complex solutions to the problem. Thus the answer above must be considered only a partial answer that leaves open some more questions.