Numerical Integration

This problem try to solve following integral and it's analytic solution also given,

$$\begin{array}{}\text{(1)}& {\int }_0^2{e}^{2x}\sin 3xdx={\displaystyle \frac{1}{13}}(-3e^4\cos (6)+2e^4\sin (6)+3)\end{array}$$

Trapeziodal method

Trapezoidal method divide the interval into

$n$ portion, and averaging the left and the right Riemann sums.

$$\begin{array}{}\text{(2)}& {\int }_a^{b}f(x)dx\approx \sum _{k=1}^n{\displaystyle \frac{f(x_{k-1})+f(x_k)}{2}}h\end{array}$$
This method give the error of

$$\begin{array}{}\text{(3)}& error=-{\displaystyle \frac{(b-a{)}^3}{12n^2}}{f}^{″}(\zeta )\end{array}$$

Composite Simpson's rule

Another method is using Composite Simpson's rule which split up the interval into n subinterval. Then, the formula is given by

$$\begin{array}{}\text{(4)}& {\int }_a^{b}f(x)dx\approx {\displaystyle \frac{h}{3}}\sum _{j=1}^{n/2}[f(x_{2j-2})+4f(x_{2j-1})+f(x_{2j})]\end{array}$$
This method give the error of

$$\begin{array}{}\text{(5)}& error=-{\displaystyle \frac{h^4}{180}}(b-a)f^{(4)}(\zeta )\end{array}$$

I wrote the following code to compare the result of two subroutine. Follow by that, there are two subroutine which use trapezoidal rule and Simpson's rule respectly. With

$n=100$, Simpson's rule got a upper hand both in accuracy and time consuming.

Code

Main

PROGRAM MAIN
!========================================================
!   Purpose: Compare different integration method
!
!   Methods: trapezoidal + simpson's
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/20/20       MorrisH        Original Code
!========================================================
    IMPLICIT NONE
    DOUBLE PRECISION, EXTERNAL :: f
    DOUBLE PRECISION :: start, finish, ans, exact

    ! exact value
    exact = (-3.d0 * dexp(4.d0) * dcos(6.d0) + &
              2.d0 * dexp(4.d0) * dsin(6.d0) + 3.d0) / 13.d0
    print *, "Exact value of this integral is: ", exact
    write(*, '(a)') repeat('=', 60)

    ! evaluate trapezoidal rule
    CALL CPU_TIME(start)
    CALL TRAPEZOIDAL(f, 0.d0, 2.d0, 100, ans)
    print *, "Trapezoidal method gives: ", ans
    print *, "Error = ", abs(ans - exact) / exact
    write(*, '(a)') repeat('=', 60)
    ! print *, "Cost", finish-start, 'sec'
    CALL CPU_TIME(finish)

    ! evaluate simpson's rule
    CALL CPU_TIME(start)
    CALL SIMPSON(f, 0.d0, 2.d0, 100, ans)
    print *, "Simpson's method gives: ", ans
    print *, "Error = ", abs(ans - exact) / exact
    write(*, '(a)') repeat('=', 60)
    ! print *, "Cost", finish-start, 'sec'
    CALL CPU_TIME(finish)
END PROGRAM MAIN

Trapezoidal

SUBROUTINE TRAPEZOIDAL(func, a, b, n, ans)
!========================================================
!   Purpose: Evaluate integral
!
!   Methods: Trapezoidal rule
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/20/20       MorrisH        Original Code
!========================================================
    IMPLICIT NONE
    DOUBLE PRECISION, EXTERNAL :: func
    DOUBLE PRECISION :: a, b, ans
    DOUBLE PRECISION :: h, u, d
    INTEGER :: n, i

    h = (b - a) / DBLE(n)
    ! divide [a, b] into n portion
    ans = 0.d0
    do i = 1, n
        d = func( (DBLE(i) - 1.d0) * h + a )
        u = func( DBLE(i) * h + a )
        ans = ans + ( u + d )
    end do
    ans = ans * h / 2.d0
END SUBROUTINE TRAPEZOIDAL

Simpson

SUBROUTINE SIMPSON(func, a, b, n, ans)
!========================================================
!   Purpose: Evaluate integral
!
!   Methods: Simpson's rule
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/20/20       MorrisH        Original Code
!========================================================
    IMPLICIT NONE
    DOUBLE PRECISION, EXTERNAL :: func
    DOUBLE PRECISION :: a, b, ans
    DOUBLE PRECISION :: h, u, m, d
    INTEGER :: n, i

    ans = 0.d0
    h = ( b - a ) / n
    do i = 1, n/2
        d = func( a + ( 2.d0 * DBLE(i) - 2.d0) * h )
        m = func( a + ( 2.d0 * DBLE(i) - 1.d0) * h )
        u = func( a + 2.d0 * DBLE(i) * h)
        ans = ans + ( d + 4.d0 * m + u )
    end do
    ans = h / 3.d0 * ans
END SUBROUTINE SIMPSON

Error analysis

The error formula of these two methods are the following,

$$\begin{array}{}\text{(6)}& {\epsilon }_{trap}=-{\displaystyle \frac{(b-a)h^2}{12}}{f}^{″}(\zeta )\end{array}$$

$$\begin{array}{}\text{(7)}& {\epsilon }_{simp}=-{\displaystyle \frac{(b-a)h^4}{180}}{f}^{(4)}(\zeta )\end{array}$$
where

$h={\displaystyle \frac{b-a}{n}}$ and the maximum value of second degree and fourth degree derivatives both happened at

$x=2$.
In order to reduce the error down to

${10}^{-6}$ for trapzoidal, it required following relation.

$$\begin{array}{}\text{(8)}& \left|{\epsilon }_{trap}\right|={\displaystyle \frac{(b-a)h^2}{12}}{f}^{″}(\zeta )={\displaystyle \frac{(2-0{)}^3}{12n^2}}{f}^{″}(\zeta )<{10}^{-6}\end{array}$$

$$\begin{array}{}\text{(9)}& n>21\times {10}^3\end{array}$$
So if

$n$ is greater than

$2.1\times {10}^4$ the accracy must less than

${10}^{-6}$ (neglecting machine error).
Similarly, in order to reduce the same error for Simpson's method, it must satified following condition

$$\begin{array}{}\text{(10)}& \left|{\epsilon }_{simp}\right|={\displaystyle \frac{(b-a)h^4}{180}}{f}^{(4)}(\zeta )={\displaystyle \frac{(2-0{)}^5}{180n^4}}{f}^{(4)}(\zeta )<{10}^{-6}\end{array}$$

$$\begin{array}{}\text{(11)}& n>170\end{array}$$
Following result verified above approximation. Also, the cpu_time for both subroutine to achieve same accracy is quite different.

Gaussian Quadratures

In this problem I try to appoximate following difinite integral using Gauss-Lengendre quadratures which given by (13).

$$\begin{array}{}\text{(12)}& {\int }_0^{48}\sqrt{1+(\cos x{)}^2}dx\end{array}$$

$$\begin{array}{}\text{(13)}& {\int }_{-1}^{-1}f(x)dx\approx \sum _{i=1}^n{w}_{i}f(x_i)\end{array}$$
where

$x_i$ is the

$i$th zeros of

$P_n(x)$, and

$w_i={\displaystyle \frac{2}{[p_n^1(x_i){]}^2}}$.

As you can see I evaluate with

$N=4,8,16,48,100$, and the integral converges very fast compare to Simpson's rule. Moreover, it was much more easier to program!

Code

PROGRAM QAULEG
!========================================================
!   Purpose: Evaluate integral
!
!   Methods: Gaussian quadrature
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/21/20       MorrisH        Original Code
!========================================================
    IMPLICIT NONE
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE ::  x, w
    DOUBLE PRECISION, EXTERNAL :: f
    DOUBLE PRECISION :: ans=0.d0
    INTEGER :: s(5)=(/4, 8, 16, 48, 100 /)
    INTEGER :: i, j, n
    CHARACTER(len=10) :: txt

    ! write result to a file
    open(99, file='result.dat', status='unknown')

    do j = 1, 5  

        ! number of node
        n = s(j)

        ! open file
        write(txt, '(I3.3)') n
        open(66, file=TRIM(txt)//'.wx', status='old', action='read')

        ! allocate memory
        allocate(x(n))
        allocate(w(n))

        ! read file into x, w array
        do i = 1, n
            read(66, *) x(i), w(i)
        end do

        ! evaluate integral
        ans = 0.d0
        do i = 1, n
            ans = ans + w(i) * f(x(i))
        end do

        ! deallocate the memory
        deallocate(x)
        deallocate(w)

        ! show the result
        print *, "For n = ", n, ":", ans
        write(*, '(A)') repeat("=", 50)

        ! write result
        write(99, *) n, ans

        close(66)
    end do
END PROGRAM QAULEG


FUNCTION f(x)
    IMPLICIT NONE
    DOUBLE PRECISION :: f, x
    f = DSQRT(1 + DCOS(24.d0 * x + 24.d0)**2) * 24.d0
END FUNCTION

Improper Integral

We try to evalute following imporper integral.

$$\begin{array}{}\text{(14)}& {\int }_1^{\mathrm{\infty }}{\displaystyle \frac{\sin 1/x}{x^{3/2}}}dx\end{array}$$
But in order to let the computer calculate this integral, we can not simply let the upper range go to

$\mathrm{\infty }$. Instead, we use change of variable and letting the singularity occurs in the lower of integration.
Since there is

$x$ at the deniminator, we let

$t=x^{-1}$.

$$dt=-x^{-2}dx$$

$$dx=-x^{2}dt=-t^{-2}dt$$
so

$$\begin{array}{}\text{(15)}& {\int }_1^{\mathrm{\infty }}{\displaystyle \frac{\sin 1/x}{x^{3/2}}}dx={\int }_0^1{\displaystyle \frac{\sin t}{\sqrt{t}}}dt\end{array}$$
The relation between two integral is given by below figure.

Now we can use the Simpson's rule to calculate this integral with

$N=10,50,100$. And the result is given below.

If I plot the integral value as the n increase, we can see that it's value oscillatint bwtween two converging value. As the n increase, it doesn't settle to a single value.

When evaluate an improper integral, like following two example,
we can make a twist that transform the sigularity which
happened at the upper bound of the integral into the lower bound.

$$\begin{array}{rl}I& ={\int }_0^{\mathrm{\infty }}{\displaystyle \frac{1}{1+x^4}}dx\\ & ={\int }_0^1{\displaystyle \frac{1}{1+x^4}}dx+{\int }_1^{\mathrm{\infty }}{\displaystyle \frac{1}{1+x^4}}dx\\ & ={\int }_0^1{\displaystyle \frac{1}{1+x^4}}dx+{\int }_0^1{\displaystyle \frac{x^2}{x^4+1}}dx\end{array}$$
Then we can apply Simpson's rule to evaluate the integral, the result
show in below figures.
As you can see, both integral oscillate and gradually converge to
a single value as the

$n$ increase.

Double Integral

Double integral has the form of

$$\begin{array}{}\text{(16)}& {\int }_a^b{\int }_{c(x)}^{d(x)}f(x,y)dydx\end{array}$$
For the step size of

$h=(b-a)/2$ with the step size of y varies with x is

$$\begin{array}{}\text{(17)}& k(x)={\displaystyle \frac{d(x)-c(x)}{2}}\end{array}$$
Simpson's rule can be write as

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\int }_a^b{\int }_{c(x)}^{d(x)}f(x,y)dydx\approx & {\int }_a^b{\displaystyle \frac{k(x)}{3}}[f(x,c(x))+4f(x,c(x)+k(x))+f(x,d(x))]dx\\ \approx & {\displaystyle \frac{h}{3}}\{{\displaystyle \frac{k(a)}{3}}[f(a,c(a))+4f(a,c(a)+k(a))+f(a,d(a))]\\ & +\frac{4k(a+h)}{3}[f(a+h,c(a+h))+4f(a+h,c(a+h)\\ & +k(a+h))+f(a+h,d(a+h))]\\ & +{\displaystyle \frac{k(b)}{3}}[f(b,c(b))+4f(b,c(b)+k(b))+f(b,d(b))]\}\end{array}\end{array}$$
Equation (18) can be easily expand to (n, m) point formula.
On the other hand, Gaussian quadrature is given by following formula.

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\int }_a^b{\int }_{c(x)}^{d(x)}& f(x,y)dydx\\ & \approx {\int }_a^b{\displaystyle \frac{d(x)-c(x)}{2}}\sum _{j-1}^n{c}_{n,j}f(x,{\displaystyle \frac{(d(x)-c(x))r_{n,j}+d(x)+c(x)}{2}})dx\end{array}\end{array}$$
where

$r_{n,j}$ is the roots of Legendre polynomial and

$c_{n,j}$ is the weight to each roots.

Following two subroutine compute the integrl of eq (19)
and compare the cpu_time for the same accuracy(

${10}^{-7}\mathrm{\%}$). As you can see the Gaussian quadrature cost roughly
10 time less than the Simpson's rule did.

$${\int }_0^{\pi /4}{\int }_{\sin x}^{\cos x}(2y\sin x+\cos x)dydx={\displaystyle \frac{1}{6}}(5\sqrt{2}-4)$$

Code

DSIMP

SUBROUTINE DSIMP(func, a, b, c, d, n, m, ans)
!========================================================
!   Purpose: Evaluate double integral
!
!   Methods: Simpson's rule
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/24/20       MorrisH        Original Code
!========================================================
    IMPLICIT NONE
    DOUBLE PRECISION, EXTERNAL :: func, c, d
    DOUBLE PRECISION :: a, b, ans
    DOUBLE PRECISION :: h
    INTEGER :: n, m, i, j
    DOUBLE PRECISION :: j1, j2, j3
    DOUBLE PRECISION :: k1, k2, k3
    DOUBLE PRECISION :: x, y, hx, Q, L

    ! initialize
    ans = 0.d0
    h = ( b - a ) / DBLE(n)
    j1 = 0.d0
    j2 = 0.d0
    j3 = 0.d0

    do i = 0, n
        x = a + DBLE(i) * h
        hx = (d(x) - c(x)) / DBLE(m)
        k1 = func(x, c(x)) + func(x, d(X))
        k2 = 0.d0
        k3 = 0.d0
        do j = 1, m-1
            y = c(x) + DBLE(j) * hx
            Q = func(x, y)
            if (mod(j, 2) == 0) then
                k2 = k2 + Q
            else
                k3 = k3 + Q
            end if
        end do
        L = (k1 + 2.d0 * k2 + 4.d0 * k3) * hx / 3.d0
        if ((i == 0) .or. (i == n)) then
            j1 = j1 + L
        else if (mod(i, 2) == 0) then
            j2 = j2 + L
        else
            j3 = j3 + L
        end if
    end do
    ans = h * (j1 + 2.d0 * j2 + 4.d0 * j3) / 3.d0
END SUBROUTINE DSIMP

DGauLeg

SUBROUTINE dGAULEG(func, a, b, c, d, n, m, ans)
!========================================================
!   Purpose: Evaluate integral
!
!   Methods: Gaussian quadrature
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/21/20       MorrisH        Original Code
!========================================================
    IMPLICIT NONE
    DOUBLE PRECISION, EXTERNAL :: func, c, d
    DOUBLE PRECISION :: a, b, ans
    INTEGER :: n, m
    DOUBLE PRECISION :: h1, h2, jx
    DOUBLE PRECISION :: x, y, d1, c1, k1, k2, Q
    INTEGER :: i, j
    CHARACTER(len=10) :: txt
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: coef, r

    write(txt, '(I3.3)') max(n, m)
    open(66, file=TRIM(txt)//'.wx', status='old', &
         action='read')

     allocate(coef(max(n, m)))
     allocate(r(max(n, m)))

    ! read weight and node
    do i = 1, n
        read(66, *) coef(i), r(i)
    end do

    h1 = (b - a) / 2.d0
    h2 = (b + a) / 2.d0
    ans = 0.d0

    do i  = 1, m
        jx = 0
        x = h1 * r(i) + h2 
        d1 = d(x)
        c1 = c(x)
        k1 = (d1 - c1) / 2.d0
        k2 = (d1 + c1) / 2.d0
        do j = 1, n
            y = k1 * r(j) + k2
            Q = func(x, y)
            jx = jx + coef(j) * Q
        end do
        ans = ans + coef(i) * k1 * jx
    end do
END SUBROUTINE dGAULEG