Implementation of Composite Simpson Rule and Gaussian Quadrature in Octave

John Bryan

• Two integration approximation techniques
  1. Composite Simpson Rule
  2. Gaussian Quadrature
  are implemented in Octave.
  
  
• 1. 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) \end{equation}
  2. Gaussian Quadrature.
     • Legendre-Gauss Integration implemented.
     • a=lower limit of integration
     • b=upper limit of integration
     • n=number of points
     • $x_{k,n}$=abscissa=kth root location of the Legendre Polynomial $P_n$ or order n.
     • $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) \end{equation} 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}) \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})} \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} \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} \label{eq:aw} \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} \end{equation} In \eqref{eq:aw}, \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)!} \end{equation}
     
     
• Octave implementation.
  
  
• Results.
  
  

  $$\int_0^1 sin(x) dx$$
  Exact	0.459697694131860
  Six-Point LG (N=6)	0.459697694131861
  Composite Simpson N=600	0.459697694131861
  Six-Point LG (N=6) difference from exact	9.436895709313831e-16
  Composite Simpson N=600 difference from exact	9.992007221626409e-16

  
  
  
  
  
  

  $$\int_1^3 x^5 dx$$
  Exact	121.3333333333333
  Six-Point LG (N=6)	121.3333333333338
  Composite Simpson N=600	121.3333333333539
  Six-Point LG (N=6) difference from exact	4.405364961712621e-13
  Composite Simpson N=600 difference from exact	2.061995019175811e-11

  
  
  
  
  
  

  Composite Quadrature for $\int_0^{50} sin(x) dx$
  Exact	3.503397150788667e-02
  Composite Six-Point LG (N-6), 50 segments	3.503397150788601e-02
  Composite Simpson N=12500	3.503397150788788e-02
  Adaptive Simpson	3.503397150788623e-02
  Composite Six-Point LG (N=6), 50 segments, difference from exact	6.661338147750939e-16
  Composite Simpson N=12500, difference from exact	1.207367539279858e-15
  Adaptive Simpson difference from exact	4.440892098500626e-16