Quadrature rules for the weight function $w(x) = \frac{1}{1+x^2}$

Question

I would like to ask a question regarding the integral $\int_a^b w(x)f(x)dx$ where $w(x)$ is a weight function and $f(x)$ is well-approximated by a polynomial of degree $2n − 1$ or less on [a, b]. In particular, I am looking for a quadrature rule for the weight function $w(x) = \frac{1}{1+x^2}$. I had a look at https://en.wikipedia.org/wiki/Gaussian_quadrature. However, I could not find the weight function that I need. Therefore, could anyone please tell me whether there exists a quadrature rule for $w(x) = \frac{1}{1+x^2}$? Thank you very much in advance!

Answer 1

The general procedure is to use the Gram-Schmidt process to create orthogonal polynomials.

Let $$ w(x)=\frac{1}{1+x^2}, $$ define $$ \langle a, b\rangle \equiv \int_a^b a(x)b(x)w(x)dx. $$ Goal is to create a set of orthogonal polynomials $\left\{p_n(x)|n\in \mathbb{Z}\right\}$ starting with $p_0(x)=1$ and creating the subsequent by ensuring the next are orthogonal to all earlier polynomials (you can visualize this in a large dimensional vector space, subtracting off the components of the vectors).

To do this, the general formula is $$ p_n(x) = x^n - \sum_{i=0}^{n-1} \frac{\langle x^n,p_i\rangle}{\langle p_i, p_i\rangle} $$ which you can do up to whatever order of quadrature you want, say $N^\textrm{th}$ order quadrature.

Quadrature approximates the integral $$ I[f] = \int_a^b f(x)w(x)dx $$ at a set of abscissa with carefully chosen weights. Define the error in quadrature as $$ E[f] = \sum_{i=0}^{N-1} w_i f(x_i) - I[f] $$

Once you have the polynomials up to $N-1$ you are really free to choose the abscissa $\left\{x_n|n\in [0,N-1] \right\}$wherever you would like, and then solve the $N$ tower equations for the weights $w_i$ $$ \sum_{i=0}^{N-1}w_i p_n(x_i) = \delta_{n,0}\langle p_0,p_0\rangle $$ which ensures that the polynomials are integrated exactly up to degree $N-1$, and as they form a basis for all polynomials, this also ensures that any polynomial up to degree $N-1$ is also exactly integrated.

However, if you choose the abscissa to be the zeros of $p_N(x)$ then you get a quadrature method that is exact up until order $2N-1$. To see this, suppose you have a polynomial $d(x)$ of order less between $N$ and $2N-1$, then it can be polynomial divided by $p_N(x)$ to produce $$ d(x)=p_N(x)q(x) + r(x) $$ where the remainder polynomial $r(x)$ is of order $N-1$ or less. Applying the quadrature operator on $d(x)$ $$ I[d] = \sum_{i=0}^{N-1} w_ip_i(x_i)p_N(x_i)q(x_i) + \sum_{i=0}^{N-1} w_ip_i(x_i)r(x_i) $$ The first term has $p_N(x)$ evaluated exactly at its zeros, hence it vanishes. The second term is exact because of the low order of the remainder.

It is much easier if the interval is symmetric, in which case the polynomials all have alternating parity and the roots come in pairs, however this is not a requirement.

Choosing the interval to be symmetric $[-a,a]$, the first four polynomials are $$ p_0(x)=1\\ p_1(x)=x\\ p_2(x)= x^{2} - \frac{a - \arctan\left(a\right)}{\arctan\left(a\right)}\\ p_3(x)=x^3 - \frac{{\left(a^{3} - 3 \, a + 3 \, \arctan\left(a\right)\right)}}{3 \, {\left(a - \arctan\left(a\right)\right)}}x $$ from which we can read off the abscissa for third order quadrature $$ x_i\in \left\{-\sqrt{\frac{{\left(a^{3} - 3 \, a + 3 \, \arctan\left(a\right)\right)}}{3 \, {\left(a - \arctan\left(a\right)\right)}}},0,\sqrt{\frac{{\left(a^{3} - 3 \, a + 3 \, \arctan\left(a\right)\right)}}{3 \, {\left(a - \arctan\left(a\right)\right)}}} \right\} $$ Define $c=\sqrt{\frac{{\left(a^{3} - 3 \, a + 3 \, \arctan\left(a\right)\right)}}{3 \, {\left(a - \arctan\left(a\right)\right)}}}$ to ease notation. The resulting weights are $$ w_0 = \frac{a-\arctan(a)}{c^2} \\ w_1 = 2\arctan(a)-\frac{2(a-\arctan(a))}{c^2}\\ w_2 = \frac{a-\arctan(a)}{c^2}. $$

I tested this out in sage math (iPython) with the following code:

def intw(f,a):
    c = sqrt((a^3-3*a+3*arctan(a))/(3*(a-arctan(a))))
    wi = vector([(a-arctan(a))/c^2, 2*arctan(a)-2*(a-arctan(a))/c^2, (a-arctan(a))/c^2])
    fx = vector([f(-c), f(0), f(c)])
    return expand(wi.dot_product(fx))