# Numerical Integration

This problem try to solve following integral and it's analytic solution also given,
\begin{equation}
\int_0^2 e^{2x}\sin{3x}dx = \dfrac{1}{13}\left( -3e^4\cos(6) + 2e^4\sin(6) + 3 \right) \tag{1}
\end{equation}

## Trapeziodal method
Trapezoidal method divide the interval into $n$ portion, and averaging the left and the right Riemann sums.
\begin{equation}
\int_a^b f(x)dx \approx \sum_{k=1}^{n} \dfrac{f(x_{k-1}) + f(x_k)}{2}h
\tag{2}
\end{equation}
This method give the error of
\begin{equation}
error \ = \ -\dfrac{(b-a)^3}{12n^2}f''(\zeta)
\tag{3}
\end{equation}

## 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{equation}
\int_a^b f(x)dx \approx \dfrac{h}{3}
\sum_{j=1}^{n/2}\left[ f(x_{2j-2})+4f(x_{2j-1}) +
f(x_{2j}) \right]
\tag{4}
\end{equation}
This method give the error of

\begin{equation}
error \ = \ -\dfrac{h^4}{180}(b-a)f^{(4)}(\zeta)
\tag{5}
\end{equation}

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{equation}
\varepsilon_{trap} = -\dfrac{(b-a)h^2}{12}f''(\zeta)
\tag{6}
\end{equation}
\begin{equation}
\varepsilon_{simp} =
-\dfrac{(b-a)h^4}{180}f^{(4)}(\zeta)
\tag{7}
\end{equation}
where $h = \dfrac{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{equation*}
\left| \varepsilon_{trap} \right| =
\dfrac{(b-a)h^2}{12}f''(\zeta)
= \dfrac{(2-0)^3}{12n^2}f''(\zeta)
< 10^{-6}
\tag{8}
\end{equation*}
\begin{equation*}
n > 21 \times 10^3
\tag{9}
\end{equation*}
So if $n$ is greater than $2.1\times10^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{equation*}
\left| \varepsilon_{simp} \right| =
\dfrac{(b-a)h^4}{180}f^{(4)}(\zeta)
= \dfrac{(2-0)^5}{180n^4}f^{(4)}(\zeta)
< 10^{-6}
\tag{10}
\end{equation*}
\begin{equation*}
n > 170 \tag{11}
\end{equation*}
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{equation}
\int_0^{48}\sqrt{1+(\cos x)^2} dx
\tag{12}
\end{equation}

\begin{equation}
\int_{-1}^{-1} f(x)dx \approx \sum_{i=1}^n w_i f(x_i)
\tag{13}
\end{equation}
where $x_i$ is the $i$th zeros of $P_n(x)$, and $w_i=\dfrac{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{equation}
\int_1^{\infty} \dfrac{\sin {1/x}}{x^{3/2}} dx \tag{14}
\end{equation}
But in order to let the computer calculate this integral, we can not simply let the upper range go to $\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}$.

\begin{equation*}
dt = -x^{-2}dx
\end{equation*}
\begin{equation*}
dx = -x^{2}dt = -t^{-2}dt
\end{equation*}
so
\begin{equation}
\int_1^{\infty} \dfrac{\sin {1/x}}{x^{3/2}} dx
= \int_0^1 \dfrac{\sin t}{\sqrt{t}}dt
\tag{15}
\end{equation}
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{equation}
\begin{split}
I &= \int_0^{\infty} \dfrac{1}{1+x^4}dx \\
&= \int_0^{1} \dfrac{1}{1+x^4}dx + \int_1^{\infty} \dfrac{1}{1+x^4}dx \\
&= \int_0^{1} \dfrac{1}{1+x^4}dx + \int_0^{1} \dfrac{x^2}{x^4+1}dx
\end{split}
\end{equation}
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{equation}
\int_a^b \int_{c(x)}^{d(x)} f(x, y)dy dx
\tag{16}
\end{equation}
For the step size of $h = (b-a) / 2$ with the step size of y varies with x is
\begin{equation}
k(x) = \dfrac{d(x)-c(x)}{2}
\tag{17}
\end{equation}
Simpson's rule can be write as
\begin{equation}
\begin{split}
\int_a^b \int_{c(x)}^{d(x)} f(x, y)dy dx
\approx
&\int_a^b \dfrac{k(x)}{3}[f(x, c(x)) + 4f(x, c(x)+
k(x)) + f(x, d(x))]dx \\
\approx
&\dfrac{h}{3}\{
\dfrac{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))]\\
&+\dfrac{k(b)}{3}[f(b,c(b))+4f(b,c(b)+k(b))+
f(b,d(b))]
\}
\end{split}
\tag{18}
\end{equation}
Equation (18) can be easily expand to (n, m) point formula.
On the other hand, Gaussian quadrature is given by following formula.
\begin{equation}
\begin{split}
\int_a^b \int_{c(x)}^{d(x)} & f(x, y)dy dx \\
& \approx \int_a^b \dfrac{d(x)-c(x)}{2}\sum_{j-1}^{n}
c_{n,j}f\left(x, \dfrac{(d(x)-c(x))r_{n,j}+d(x)+
c(x)}{2}\right)dx
\end{split}
\tag{19}
\end{equation}
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}\%$). As you can see the Gaussian quadrature cost roughly
10 time less than the Simpson's rule did.
\begin{equation}
\int_0^{\pi /4}\int_{\sin x}^{\cos x}
(2y\sin x + \cos x)dydx
= \dfrac{1}{6}(5\sqrt{2}-4)
\end{equation}

## 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