Implementation of Composite Simpson Rule and Gaussian Quadrature in Octave
John Bryan
-
Two integration approximation techniques
- Composite Simpson Rule
- Gaussian Quadrature
are implemented in Octave.
-
- Composite Simpson Rule.
- a=lower limit of integration
- b=upper limit of integration
- N=number of 1-4-1 calculations
- 2N+1=number of abscissae
- h=(b-a)/2N=distance between adjacent abscissae
\begin{equation}
\int_a^b f(x) dx \approx \frac{h}{3}\sum_{k=1}^N f(a+[2k-2]h)+4f(a+[2k-1]h)+f(a+2kh)
\tag{eq. 1}
\end{equation}
- Gaussian Quadrature. Legendre-Gauss Integration is implemented.
- a=lower limit of integration
- b=upper limit of integration
- n=number of points
- $x_{k,n}$=abscissa
- $w_{k,n}$=weight
- c=(b-a)/2
- d=(b+a)/2
\begin{equation}
\int_a^b f(x) dx \approx c\sum_{k=1}^n w_{k,n}f(cx_{k,n}+d)
\tag{eq. 2}
\end{equation}
$x_{k,n}$ is the kth root location of the Legendre Polynomial $P_n$ of order n. Tricomi initial approximation of the roots in descending order is given by
\begin{equation}
x_{k,n} \approx (1-\frac{1}{8n^2}+\frac{1}{8n^3})cos(\pi\frac{4k-1}{4n+2}) \tag{eq. 3}
\end{equation}
Refinement of $x_{k,n}$ is obtained using Newton-Raphson iteration. With $\nu$=iteration,
\begin{align}
x_{{k,n}_{\nu+1}} & = x_{{k,n}_\nu} - \frac{P(x_{{k,n}_\nu})}{P'(x_{{k,n}_\nu})} \\[10pt]
& = x_{{k,n}_\nu} - \frac{1-x^2_{{k,n}_\nu}}{n} \frac{P_n(x_{{k,n}_\nu})}{P_{n-1}(x_{{k,n}_\nu})-x_{{k,n}_\nu}P_n(x_{{k,n}_\nu})} \tag{eq. 4}
\end{align}
using $$ P'(x_{{k,n}_\nu})=\frac{n}{1-x^2_{{k,n}_\nu}}(P_{n-1}(x_{{k,n}_\nu})-x_{{k,n}_\nu}P_n(x_{{k,n}_\nu})). $$
The weight,
\begin{equation}
w_{k,n} = \frac{2(1-x^2_{k,n})}{(n+1)^2[P_{n+1}(x_{k,n})]^2} \tag{eq. 5}
\end{equation}
where
\begin{equation}
P_{n+1}(x)= 2^{n+1}\sum_{k=0}^{n+1} x^k \binom{n+1}{k} \binom{\frac{n+k}{2}}{n+1} \tag{eq. 6}
\end{equation}
using
\begin{equation}
P_n(x)= 2^n\sum_{k=0}^{n} x^k \binom{n}{k} \binom{\frac{n+k-1}{2}}{n} \tag{eq. 7}
\end{equation}
In eq. 6,
\begin{equation}
\binom{\frac{n+k}{2}}{n+1} = \frac{(\frac{n+k}{2})(\frac{n+k}{2}-1)(\frac{n+k}{2}-2) \cdots (\frac{n+k}{2}-k+1)}{(n+1)!} \tag{eq. 8}
\end{equation}
-
Octave implementation.
-
Results.
$$ \int_1^4 e^x cos(x) dx $$
$$
\begin{array}{c|lcr}
\text{Exact} & -40.382 \\
\text{Composite Simpson (N=5)} & -40.375 \\
\text{Composite Simpson (N=10)} & -40.381 \\
\text{Gaussian Quadrature (N=5)} & -40.382
\end{array}
$$
$$ \int_3^5 (x^8-x^7) dx $$
$$
\begin{array}{c|lcr}
\text{Exact} & 166818.889 \\
\text{Composite Simpson (N=6)} & 166822.541 \\
\text{Composite Simpson (N=12)} & 166819.118 \\
\text{Gaussian Quadrature (N=6)} & 166818.889
\end{array}
$$