# Solving Partial Differential Equation

## Poisson's equation

Poisson equation has the form of
\begin{equation}
\dfrac{\partial^2 u}{\partial x^2}+\dfrac{\partial^2 u}{\partial y^2} = -\dfrac{q}{K}
\tag{1}
\end{equation}
In this particular case $u$ is the steady-state temperature. Given the boundary condition we can apply finite difference into this PDE and solve it numerically. The finite difference has the following relation.
\begin{equation}
w_{ij} = \dfrac{1}{2+2\left( \dfrac{h}{j}\right)^2}
\left[ w_{i+1,j} + w_{i-1,j} + (\dfrac{h}{k})^2(w_{i,j+1}+w_{i,j-1})
-h^2f_{ij} \right]
\tag{2}
\end{equation}
Instead of solving this linear equation, we were using relaxation method to solve it which is let the solution iterate until it reach a desired tolerance. Follow figure show the temperature distribution in x-y plane with line represent the isotherm area. The area around x and y axis is what we expected that the center are the hottest and gradually cool down to another two boundary.

### Code

program poisson
    implicit none
    ! double precision, parameter :: h=.4, k=1./3.
    double precision, parameter :: a=0., b=6., c=0., d=5.
    double precision, parameter :: f=-1.5/1.04
    double precision, parameter :: tol=1.
    integer, parameter :: n=15, m=15
    double precision :: h, k
    double precision :: w((n+1), (n+1)), x((n+1)), y((n+1))
    double precision :: z
    double precision :: NORM, lamb, mu
    integer :: i, j, l, Nmax=100

    h = (b-a) / n
    k = (d-c) / m

    ! x and y
    do i = 1, n-1
        x(i) = a + dble(i) * h
    end do
    do j = 1, m-1
        y(j) = c + dble(j) * k
    end do
    ! init w
    w = 0.
    l = 1
    lamb = h * h / k / k
    mu = 2. * (1. + lamb)
    ! step 6
    do
        ! step 7
        z = (-h*h*f + y(m-1)*(6.-y(m-1)) + lamb * 0. + lamb * &
            w(1, m-2) + w(2, m-1)) / mu
        NORM = abs(z - w(1, m-1))
        w(1, m-1) = z
        ! step 8
        do i = 2, n-2
            z = (-h*h*f + lamb * 0. + w(i-1, m-1) + w(i+1, m-1)&
                + lamb * w(i, m-2)) / mu
            if (abs(w(i, m-1) - z).gt.NORM) then
                NORM = abs(w(i, m-1) - z)
            end if
            w(i, m-1) = z
        end do
        ! step 9
        z = (-h*h*f + 0. + lamb * 0. + w(n-2, m-1) + &
            lamb * w(n-1, m-2)) / mu
        if (abs(w(n-1, m-1) - z).gt.NORM) then
            NORM = abs(w(n-1, m-1) - z)
        end if
        w(n-1, m-1) = z
        ! step 10
        do j = m-2, 2, -1
            ! step 11
            z = (-h*h*f + y(j) * (5. - y(j)) + lamb * w(1, j+1)&
                + lamb * w(1, j-1) + w(2, j)) / mu
            if (abs(w(1, j) - z).gt.NORM) then
                NORM = abs(w(1, j) - z)
            end if
            w(1, j) = z
            ! step 12
            do i = 2, n-2
                z = (-h*h*f + w(i-1, j) + lamb * w(i, j+1) + &
                    w(i+1, j) + lamb * w(i, j-1)) / mu
                if (abs(w(i, j) - z).gt.NORM) then
                    NORM = abs(w(i, j) - z)
                end if
                w(i, j) = z
            end do
            ! step 13
            z = (-h*h*f + 0. + w(n-2, j) + lamb * w(n-1, j) &
                + lamb * w(n-1, j+1) + lamb * w(n-1, j-1)) / mu
            if (abs(w(n-1, j) - z).gt.NORM) then
                NORM = abs(w(n-1, j) - z)
            end if
            w(n-1, j) = z
        end do
        ! step 14
        z = (-h*h*f + y(1)*(5.-y(1)) + lamb * x(1)*(6.-x(1)) + &
            lamb * w(1, 2) + w(2, 1))/ mu
        if (abs(w(1, 1) - z).gt.NORM) then
            NORM = abs(w(1, 1) - z)
        end if
        w(1, 1) = z
        ! step 15
        do i = 2, n-2
            z = (-h*h*f + lamb * x(i) * (6. - x(i)) + w(i-1, 1) &
                + lamb * w(i, 2) + w(i+1, 1)) / mu
            if (abs(w(i, 1) - z).gt.NORM) then
                NORM = abs(w(i, 1) - z)
            end if
            w(i, 1) = z
        end do
        ! step 16
        z = (-h*h*f + 0. + lamb * x(n-1) * (6. - x(n-1)) + w(n-2, 1)&
            + lamb * w(n-1, 2)) / mu
        if (abs(w(n-1, 1) - z).gt.NORM) then
            NORM = abs(w(n-1, 1) - z)
        end if
        w(n-1, 1) = z
        if (NORM.le.tol) then
            open(66, file='data.dat', status='unknown')
            do i = 1, n+1
                do j = 1, m+1
                    write(66, *) w(i, j)
                end do
            end do
            close(66)
            print *, 'success'
            exit
        end if
        if (l.gt.Nmax) then
            exit
        end if
        l = l + 1
        print*, NORM
    end do
end program poisson

## Diffusion equation
Similarly, when solving diffusion equation which is another type of PDE. We use finite difference to approximate the derivative terms. The equation of stress-strain of a cylindrical material can be written as
\begin{equation}
\dfrac{\partial^2 T}{\partial r^2} + \dfrac{1}{r}\dfrac{\partial T}{\partial r} = \dfrac{1}{4K}\dfrac{\partial T}{\partial t}
\tag{3}
\end{equation}
In this case the finite difference equation can be written as
\begin{equation}
\left( 4kK/h^2 - 2kK/hr_i \right)w_{i-1j} + \left( -8kK/h^2 +1\right)w_{ij}+ \left( 4kK/h^2 + 2kK/hr_i\right)w_{i+1j} = w_{ij+1}
\tag{4}
\end{equation}
which can also be written as
\begin{equation}
\textbf{A}w^j = w^{j+1} \tag{5}
\end{equation}
However another approach is combining forward difference and backward difference when dealing with time derivative term. In this case the finite difference has the form of
\begin{equation}
\begin{split}
&\left( 2/h^2 - 1/4hr_i \right)w_{i-1j}+\left( -1/h^2 + 1/4kK \right)w_{ij}+\left( 1/2h^2 + 1/4hr_i \right)w_{i+1j}\\&=
\left( -2/h^2 + 1/4hr_i \right)w_{i-1j+1}+\left( 1/h^2 + 1/4kK \right)w_{ij+1}+\left( -1/2h^2 + 1/4hr_i \right)w_{i+1j+1}
\end{split}
\tag{6}
\end{equation}
which looks a lit bit intimidating, however we can transform this into matrix form of
\begin{equation}
\textbf{A}w^{j+1}=\textbf{B}w^{j} \tag{7}
\end{equation}
Consider a boundary condition is written as a array then this can be further write as
\begin{equation}
w^{j+1} = \textbf{A}^{-1}( \textbf{B} w^{j} + D^j - D^{j+1} ) \tag{8}
\end{equation}
Above two method has similar solution given by below figure.

Strain $I$ can be evaluated by
\begin{equation}
I = \int_{0.5}^{1}\alpha T(r,t)rdr
\tag{9}
\end{equation}
Two methods give following result.

| | Forward | Crank Nicolson |
| -------- | -------- | -------- |
| I | 1501.30 | 1368.58 |

### Code
### Forward finite-difference

program diffusion
    implicit none
    double precision, parameter :: diff=.1, k=.5, h=.1
    double precision, parameter :: rmin=.5, rmax=1., tmax=10.
    integer, parameter :: n=nint(tmax/k), m=nint((rmax-rmin)/h)
    double precision :: a, b, c
    double precision :: w((m-1), n), r((m-1)), D((m-1), n)
    integer :: i, j

    ! step 1: r
    do i = 1, m-1
        r(i) = rmin + h * dble(i)
    end do
    ! step 2: initial w(r, 0)
    do i = 1, m-1
        w(i, 1) = 200. * (r(i) - .5)
    end do
    do j = 1, n-1
        a = 4. * diff * k / h / h - 2. * diff * k / h / r(1)
        b = -8. * diff * k / h / h + 1.
        c = 4. * diff * k / h / h + 2. * diff * k / h / r(1)
        D(1, j+1) = a * (dble(j) * k) + b * w(1, j) + c * w(2, j)
        do i = 2, m-2
            a = 4. * diff * k / h / h - 2. * diff * k / h / r(i)
            b = -8. * diff * k / h / h + 1.
            c = 4. * diff * k / h / h + 2. * diff * k / h / r(i)
            D(i, j+1) = a * w(i-1, j) + b * w(i, j) + c * w(i+1, j)
        end do
        a = 4. * diff * k / h / h - 2. * diff * k / h / r(m-1)
        b = -8. * diff * k / h / h + 1.
        c = 4. * diff * k / h / h + 2. * diff * k / h / r(m-1)
        D(m-1, j+1) = a * w(m-2, j) + b * w(m-1, j) + c * (100. + 40. * (dble(j) * k))

        w(:, j+1) = D(:, j+1)
    end do
    ! step 6: write w(r, 10)
    open(66, file='out.dat', status='unknown')
    j = 20
    write(66, *) rmin, dble(j) * k
    do i = 1, m-1
        write(66, *) r(i), w(i, j) 
    end do
    write(66, *) rmax, 100. + 40. * (dble(j) * k)

program diffusion
    implicit none
    double precision, parameter :: diff=.1, k=.5, h=.1
    double precision, parameter :: rmin=.5, rmax=1., tmax=10.
    integer, parameter :: n=nint(tmax/k), m=nint((rmax-rmin)/h)
    double precision :: a, b, c
    double precision :: w((m-1), n), r((m-1)), D((m-1), n)
    integer :: i, j

    ! step 1: r
    do i = 1, m-1
        r(i) = rmin + h * dble(i)
    end do
    ! step 2: initial w(r, 0)
    do i = 1, m-1
        w(i, 1) = 200. * (r(i) - .5)
    end do
    do j = 1, n-1
        a = 4. * diff * k / h / h - 2. * diff * k / h / r(1)
        b = -8. * diff * k / h / h + 1.
        c = 4. * diff * k / h / h + 2. * diff * k / h / r(1)
        D(1, j+1) = a * (dble(j) * k) + b * w(1, j) + c * w(2, j)
        do i = 2, m-2
            a = 4. * diff * k / h / h - 2. * diff * k / h / r(i)
            b = -8. * diff * k / h / h + 1.
            c = 4. * diff * k / h / h + 2. * diff * k / h / r(i)
            D(i, j+1) = a * w(i-1, j) + b * w(i, j) + c * w(i+1, j)
        end do
        a = 4. * diff * k / h / h - 2. * diff * k / h / r(m-1)
        b = -8. * diff * k / h / h + 1.
        c = 4. * diff * k / h / h + 2. * diff * k / h / r(m-1)
        D(m-1, j+1) = a * w(m-2, j) + b * w(m-1, j) + c * (100. + 40. * (dble(j) * k))

        w(:, j+1) = D(:, j+1)
    end do
    ! step 6: write w(r, 10)
    open(66, file='out.dat', status='unknown')
    j = 20
    write(66, *) rmin, dble(j) * k
    do i = 1, m-1
        write(66, *) r(i), w(i, j) 
    end do
    write(66, *) rmax, 100. + 40. * (dble(j) * k)
end program diffusion

#### Crank-Nicolson

program diffusion
    implicit none
    double precision, parameter :: diff=.1, k=.5, h=.1
    double precision, parameter :: rmin=.5, rmax=1., tmax=10.
    integer, parameter :: n=nint(tmax/k), m=nint((rmax-rmin)/h)
    double precision :: w((m-1), n), r((m-1)), A((m-1), (m-1)), B((m-1), (m-1)), C((m-1))
    double precision :: work2((m-1)), pivot((m-1))
    integer :: i, j, rc

    ! step 1: r
    do i = 1, m-1
        r(i) = rmin + h * dble(i)
    end do
    ! step 2: initial w(r, 0)
    do i = 1, m-1
        w(i, 1) = 200. * (r(i) - .5)
    end do
    do j = 1, n-1
        A(1, 1) = 1. / h / h + 1. / 4. / diff / k
        A(1, 2) = - 1. / 2. / h / h - 1. / 4. / h / r(1)
        B(1, 1) = - 1. / h / h + 1. / 4. / diff / k
        B(1, 2) = 1. / 2. / h / h + 1. / 4. / h / r(1)
        do i = 2, m-2
            A(i, i-1) = - 1. / 2. / h / h + 1. / 4. / h / r(i)
            A(i, i) = 1. / h / h + 1. / 4. / diff / k
            A(i, i+1) = - 1. / 2. / h / h - 1. / 4. / h / r(i)
            B(i, i-1) = 1. / 2. / h / h - 1. / 4. / h / r(i)
            B(i, i) = - 1. / h / h + 1. / 4. / diff / k
            B(i, i+1) = 1. / 2. / h / h + 1. / 4. / h / r(i)
        end do
        A(m-1, m-2) = - 1. / 2. / h / h + 1. / 4. / h / r(m-1)
        A(m-1, m-1) = 1. / h / h + 1. / diff / k
        B(m-1, m-2) = 1. / 2. / h / h - 1. / 4. / h / r(m-1)
        B(m-1, m-1) = - 1. / h / h + 1. / 4. / diff / k
        call dgetrf(m-1, m-1, A, m-1, pivot, rc)
        call dgetri(m-1, A, m-1, pivot, work2, m-1, rc)
        C = matmul(B, w(:, j))
        C(1) = C(1) + (1. / 2. / h / h - 1. / 4. / h / r(1)) * dble(j) * k &
                    - (- 1. / 2. / h / h + 1. / 4. / h / r(1)) * dble(j+1) * k
        C(m-1) = C(m-1) + (1. / 2. / h / h + 1. / 4. / h / r(m-1)) * (100. + 40. * dble(j) * k) &
                        - (- 1. / 2. / h / h - 1. / 4. / h / r(m-1)) * (100. + 40. * dble(j+1) * k)
        w(:, j+1) = matmul(A, C)
    end do
    ! step 6: write w(r, 10)
    open(66, file='out2.dat', status='unknown')
    j = 20
    write(66, *) rmin, dble(j) * k
    do i = 1, m-1
        write(66, *) r(i), w(i, j) 
    end do
    write(66, *) rmax, 100. + 40. * (dble(j) * k)
end program diffusion

## Wave equation
Wave equation has similar approach again! Turning a PDE into finite difference equation and iterate it over and over. The finite difference equation is
\begin{equation}
w_{ij+1} = 2(1-\lambda^2)w_{ij} + \lambda^2(w_{i+1j}+w_{i-1j})-w_{ij-1}
\tag{10}
\end{equation}
where $\lambda \equiv \alpha k/h$.
Then it can be represent as matrix form of
\begin{equation}
w^{j+1}=\textbf{C}w^j-w^{j-1}
\tag{11}
\end{equation}
And since it required two step to evaluate the next step, we'll have difficulty staring, so we need another equation to start off.
\begin{equation}
w_{i1} = (1-\lambda^2)f(x_i) + dfrac{\lambda^2}{2}f(x_{i+1})+dfrac{\lambda^2}{2}f(x_{i-1})+kg(x_i)
\tag{12}
\end{equation}
The solution at $t=0.2$ and $t=0.5$ of $V$ and $i$ are given by below figure.

### Code

program wave
    implicit none
    double precision, external :: f, g
    double precision, parameter :: L=.3d0, C=.1d0, alpha=1.d0/(L*C)
    double precision, parameter :: lmax=200.d0, tmax=.2d0
    double precision, parameter :: h=10.d0, k=.1d0
    double precision, parameter :: lambda=k*alpha/h
    integer, parameter :: m=int(lmax/h), n=int(tmax/k)
    double precision :: w(0:m, 0:n)=0.d0
    integer :: i, j

    w(0, 0) = f(0.d0)
    w(m, 0) = f(lmax)

    do i = 1, m-1
        w(i, 0) = f(dble(i) * h)
        w(i, 1) = (1.d0 - lambda * lambda) * f(dble(i) * h) + &
                  + .5d0 * (lambda * lambda) * (f(dble(i+1) * h) &
                  + f(dble(i-1) * h)) + k * g(dble(i) * h)
    end do

    do j = 1, n-1
        do i = 1, m-1
            w(i, j+1) = 2.d0 * (1.d0 - lambda * lambda) * w(i, j) &
                        + lambda * lambda * (w(i+1, j) + w(i-1, j)) &
                        + w(i, j-1)
        end do
    end do

    ! output
    open(66, file='1_1.dat', status='unknown')
    j = n
    do i = 0, m
        write(66, *) dble(j-2) * k, dble(i-1) * h, w(i, j)
    end do
end program wave

function f(x)
    implicit none
    double precision, parameter :: pi=4.d0*atan(1.d0)
    double precision :: f, x
    f = 110.d0 * dsin(pi * x / 200.d0)
end function

function g(x)
    implicit none
    double precision :: g, x
    g = 0.d0 * x
end function