Solving Partial Differential Equation

Poisson's equation

Poisson equation has the form of

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}=-{\displaystyle \frac{q}{K}}\end{array}$$
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{array}{}\text{(2)}& w_{ij}={\displaystyle \frac{1}{2+2{\left({\displaystyle \frac{h}{j}}\right)}^2}}[w_{i+1,j}+w_{i-1,j}+({\displaystyle \frac{h}{k}}{)}^2(w_{i,j+1}+w_{i,j-1})-h^2{f}_{ij}]\end{array}$$
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{array}{}\text{(3)}& {\displaystyle \frac{{\partial }^{2}T}{\partial r^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial T}{\partial r}}={\displaystyle \frac{1}{4K}}{\displaystyle \frac{\partial T}{\partial t}}\end{array}$$
In this case the finite difference equation can be written as

$$\begin{array}{}\text{(4)}& (4kK/h^2-2kK/h{r}_i)w_{i-1j}+(-8kK/h^2+1)w_{ij}+(4kK/h^2+2kK/h{r}_i)w_{i+1j}=w_{ij+1}\end{array}$$
which can also be written as

$$\begin{array}{}\text{(5)}& \mathbf{\text{A}}{w}^j=w^{j+1}\end{array}$$
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{array}{}\text{(6)}& \begin{array}{rl}& (2/h^2-1/4h{r}_i)w_{i-1j}+(-1/h^2+1/4kK)w_{ij}+(1/2h^2+1/4h{r}_i)w_{i+1j}\\ & =(-2/h^2+1/4h{r}_i)w_{i-1j+1}+(1/h^2+1/4kK)w_{ij+1}+(-1/2h^2+1/4h{r}_i)w_{i+1j+1}\end{array}\end{array}$$
which looks a lit bit intimidating, however we can transform this into matrix form of

$$\begin{array}{}\text{(7)}& \mathbf{\text{A}}{w}^{j+1}=\mathbf{\text{B}}{w}^j\end{array}$$
Consider a boundary condition is written as a array then this can be further write as

$$\begin{array}{}\text{(8)}& w^{j+1}={\mathbf{\text{A}}}^{-1}(\mathbf{\text{B}}{w}^j+D^j-D^{j+1})\end{array}$$
Above two method has similar solution given by below figure.

Strain

$I$ can be evaluated by

$$\begin{array}{}\text{(9)}& I={\int }_{0.5}^1\alpha T(r,t)rdr\end{array}$$
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{array}{}\text{(10)}& w_{ij+1}=2(1-{\lambda }^2)w_{ij}+{\lambda }^2(w_{i+1j}+w_{i-1j})-w_{ij-1}\end{array}$$
where

$\lambda \equiv \alpha k/h$.
Then it can be represent as matrix form of

$$\begin{array}{}\text{(11)}& w^{j+1}=\mathbf{\text{C}}{w}^j-w^{j-1}\end{array}$$
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{array}{}\text{(12)}& w_{i1}=(1-{\lambda }^2)f(x_i)+dfrac{\lambda }^{2}2f(x_{i+1})+dfrac{\lambda }^{2}2f(x_{i-1})+kg(x_i)\end{array}$$
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