1
$\begingroup$

I want to numerically solve the advection-diffusion equation: \begin{equation} u_t(x,t) + cu_x(x,t) = v u_{xx}(x,t) \end{equation} for $x \in [0,1]$ and $t \geq 0$ subject to the boundary conditions $u(0,t) = u(1,t) = 0 $ and $u(x,0) = f(x)$. To compare the accuracy of the numerical solution, I will at first derive the analytical solution.

Analytical Solution The solution satisfies \begin{equation} u(x,t) = h(x,t) \exp \left( \alpha x + \beta t \right) \end{equation} If one sets $\alpha = c /2v$ and $\beta = - c^{2} / 4v$ the function $h(x,t)$ adheres to the conventional heat equation: $h_t(x,t) = vh_{xx}(x,t)$ suject to $h(0,t) = h(1,t) = 0$ and $h(x,0) = g(x)$. The solution is particularly inexpensive to calculate if one sets $g(x) = \sin ( 2\pi x)$. Then, the function $u$ is: \begin{equation} u(x,t) = \sin (2 \pi x) \exp \left( -4 \pi^2 v t - \dfrac{c^2 t}{4v} + \dfrac{cx}{2v} \right ) \end{equation}

Numerical Solution with BCTS and Crank-Nicolson To solve the advection-diffusion equation numerically, I use the BTCS and Crank-Nicolson algorithms. Define first \begin{align*} r &= \dfrac{v \Delta t}{\Delta x^{2}} \\ R &= \dfrac{c \Delta t}{\Delta x} \end{align*} The discretized version of the advection-diffusion equation adherring to the BTCS algorithm is: \begin{equation} u_{k}^{n}(1 + 2r) + u_{k-1}^{n} ( -R/2 - r) + u_{k+1}^{n} (R/2 -r) = u_{k}^{n-1} \end{equation} The differential equation for the Crank-Nicolson algorithm can be written as: \begin{align*} &u_{k}^{n} ( 1 + r) + u_{k+1}^{n} ( R/4 - r/2) + u_{k-1}^{n} ( -R/4 - r/2) = \\ & u_{k}^{n-1} ( 1 - r) + u_{k+1}^{n} ( -R/4 + r/2) + u_{k-1}^{n} ( R/4 + r/2) \end{align*}

Simulations When I solve the advection-diffusion algorithm with $c=1, v=2$ and $\Delta t = 0.001$ and $\Delta x \approx 0.0416 $ the solution looks as follows: enter image description here

The analytical solution is both positive and negative, the numerical solution however is an inverted parabola. With a reduced strength of the diffusion term $v=1/6$ the solution is very accurate:

enter image description here

How can I compute the numerical solution accurately over a large range of values for $v$?

I am particularly concerned about this question, as I solved the advection-diffusion equation to understand the algorithms and want to apply them to a non-linear PDE of which I don't know the analytical solution. The computer code used for the example is:

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt


class ConvectionDiffusion(object):
    """
    Class to construct solutions to the convection-diffusion equation
    """

    def __init__(self, DELTA_T, M, V, C):
        """
        Parameters:
        -----------
        DELTA_T: scalar(float):
            The time step Delta_t

        M: scalar(float):
            The number of grid-points for x

        V: scalar(float):
            The constant for the diffusion term

        C: scalar(float):
            The constant for the advection term
        """
        self.DELTA_T = DELTA_T
        self.M = M
        self.DELTA_X = 1 / (M - 1)
        self.xVec = np.linspace(0, 1., num=M)
        self.C = C
        self.V = V
        self.R = C * DELTA_T / self.DELTA_X
        self.r = V * DELTA_T / self.DELTA_X**2
        self.u0 = np.sin(2 * np.pi * self.xVec) * \
            np.exp(self.C * self.xVec / (2 * self.V))

    def Implicit(self, stencil, RHS, Sol, l_and_u, T):
        """
        Solving a PDE with implicit algorithms through (banded) matrix inversion

        Parameters:
        -----------
        stencil: fun
            Function defining for each value u^{n-1} the stencil such that
            the discretized PDE can be iterated forward

        RHS: fun
            Constructs the RHS of A x = b.

        l_and_u: tupel(int)
            The number of lower and upper diagonal elements used in stencil to
            use linalg.solve_banded

        T: scalar(float):
            The time at which the solution to the PDE is expressed.

        Sol: fun
            Function generating next periods u^{n} from the solution in
            the interior `x` and the boundary values

        Return:
        -------
        initial: array_like(float):
            The solution to the PDE at time T
        """
        DELTA_T, initial = self.DELTA_T, self.u0
        t = 0.0
        while t < T:
            A = stencil(initial)
            b = RHS(initial)
            x = linalg.solve_banded(l_and_u, A, b)
            initial = Sol(x)
            t = t + DELTA_T
        return initial

    def stencilCN(self, initial):
        """
        Stencil for the Crank-Nicolson algorithm; Stencil * u^n = b

        Parameters:
        ----------
        initial: array_like(float):
            The previous solution

        Returns:
        --------
        A: banded_matrix(float):
            The matrx in baned for to pass it to scipy.linalg.banded_solve
        """
        A = np.zeros((3, len(initial) - 2))
        R, r = self.R, self.r
        A[0, 1:] = R / 4 - r / 2
        A[1, :] = 1 + r
        A[2, :-1] = -R / 4 - r / 2
        return A

    def stencil(self, initial):
        """
        Stencil for the BCTS algorithm, ; Stencil * u^n = u^{n-1}

        Parameters:
        ----------
        initial: array_like(float):
            The previous solution

        Returns:
        --------
        A: banded_matrix(float):
            The matrx in baned for to pass it to scipy.linalg.banded_solve
        """
        A = np.zeros((3, len(initial) - 2))
        R, r = self.R, self.r
        A[0, 1:] = R / 2 - r  # uper diagonal
        A[1, :] = 1 + 2 * r
        A[2, :-1] = -R / 2 - r  # lower diagonal
        return A

    def RHS(self, initial):
        """
        Pepare the RHS for the BCTS algorithm.

        Parameters:
        -----------
        initial: array_like(float):
            The previous solution u^{n-1}

        Returns:
        --------
        b: array_like(float):
            Only the interior values of u^{n-1}
        """
        return initial[1: -1]

    def RHSCN(self, initial):
        """
        Pepare the RHS for the Crank-Nicolson algorithm.

        Parameters:
        -----------
        initial: array_like(float):
            The previous solution u^{n-1}

        Returns:
        --------
        b: array_like(float):
            A weighted sum of previous u^{n-1} values
        """
        R, r = self.R, self.r
        present = initial[1:-1] * (1 - r)
        past =  initial[:-2] * (R / 4 + r / 2)
        future = initial[2:] *(-R / 4 + r / 2)
        b =  present + past + future
        return b

    def Sol(self, x):
        """
        Add the boundary condition

        Parameters:
        x: array_like(float):
            The solution in the interior

        Returns:
        --------
        uNew: array_like(float):
            The entire solution
        """
        uNew = np.zeros(len(x) + 2)
        uNew[1:-1] = x
        return uNew

    def analyticalSol(self, x, T):
        """
        Analytical Solution of the advection-difussion equation

        Parameters:
        -----------
        x   array_like(float):
            The state space
        T:  scalar(float):
            The time at which the solution is evaluated

        Returns:
        v:  array_like(float):
            The solution for the PDE
        """
        C, V = self.C, self.V
        exponential = (- 4 * np.pi**2 * V * T
                       - C**2 * T / (4 * V) + C * x[1:-1] / (2 * V) )
        interior = np.sin(2 * np.pi * x[1:-1]) * np.exp(exponential)
        v = self.Sol(interior)
        return v

    def ComparisonSol(self, T):
        """
        Calculates the analytical solution and the its numercial approximation
        according to the BCTS, C-N and Explicit algorithm at time T

        Parameters:
        T:  scalar(float):
            The time value

        Returns:
        --------
        sol: array_like(float):
            The array with all solutions

        e:  array_like(float):
            The normalized second error norm of each approximation
        """
        xVec, DELTA_T, DELTA_X = self.xVec, self.DELTA_T, self.DELTA_X
        sol, E = np.empty((len(self.u0), 3)), np.zeros(2)

        ## == Anaytical == ##
        sol[:, 0] = self.analyticalSol(xVec, T)

        ## == BTCS == ##
        sol[:, 1] = self.Implicit(self.stencil, self.RHS, self.Sol, (1, 1), T)
        E[0] = np.linalg.norm(sol[:, 0] - sol[:, 1]) / np.linalg.norm(sol[:, 0])

        ## == CN == ##
        sol[:, 2] = self.Implicit(
            self.stencilCN, self.RHSCN, self.Sol, (1, 1), T)
        E[1] = np.linalg.norm(sol[:, 0] - sol[:, 2]) /np.linalg.norm(sol[:, 0])
        return sol, E

if __name__ == "__main__":
    Simulation = ConvectionDiffusion(DELTA_T=0.001, M=50, V=1/6, C=1.)
    y, e = Simulation.ComparisonSol(0.5)
    x = Simulation.xVec
    fig, ax1 = plt.subplots()
    ax1.plot(x, y[:, 0], 'b-', label='Anaytical Solution')
    ax1.set_xlabel('x')
    ax1.set_ylabel('Axis Anaytical Solution', color='b')
    ax1.tick_params('y', colors='b')
    ax1.legend(loc=2)

    ax2 = ax1.twinx()
    ax2.plot(x, y[:, 1], 'r--', label='BTCS ')
    ax2.plot(x, y[:, 2], 'r:', label='C-N ')
    ax2.set_ylabel('Axis Numerical Solutions', color='r')
    ax2.tick_params('y', colors='r')
    ax2.legend()
    fig.tight_layout()
    plt.show()
$\endgroup$
  • $\begingroup$ Have you tried rescaling your variables? $\endgroup$ – nicoguaro Jan 25 '17 at 16:37
  • $\begingroup$ @nicoguaro: Thank you for your comment! Can you explain briefly in which situations rescaling facilitates the computation of a numerical solution? $\endgroup$ – fabian Jan 26 '17 at 9:31
  • $\begingroup$ In general, nondimensionalization is a good idea in general. It might improve the condition number of your system. $\endgroup$ – nicoguaro Jan 31 '17 at 13:15
4
$\begingroup$

The solution you believe to be inaccurate is actually by far the more accurate one; you've simply plotted it in a very deceptive way. For $\nu=2$, the exact solution is actually no bigger than about $10^{-35}$ everywhere -- it's zero for all intents and purposes. Therefore the numerical solution is correct to 10 digits -- far better than the accuracy of your solution when $\nu=1/6$.

$\endgroup$
  • $\begingroup$ Thank you very much for your reply and putting the numerical solution into perspective. I was a bit worried that the solution doesn't fall below zero, which I interpreted as an error. I learn from your reply that this worry was ill-founded. $\endgroup$ – fabian Jan 26 '17 at 9:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.