Implementing Dirichlet BC for the Advection-Diffusion equation using a second-order Upwind Scheme finite difference discretization

Question

i am implementing a Matlab code to solve the following equation numerically :

$$ (\frac{\partial c}{\partial t} =-D_{e} \frac{\partial^2 c}{\partial z^2} +U_{z}\frac{\partial c}{\partial z}) $$ with the folllowing boundary conditions:
$$ at \ \ z = 0 \ ; \ c= c_{i0} \\ \\ at \ \ z = L \ ; \ \frac{\partial c_{i}}{\partial z }=0 $$ what i am doing right now is applying following discretization scheme:

central difference approximation for the second derivative

$$ \displaystyle \frac{\partial^2 c}{\partial z^2 } \ \approx \ \frac{c^{i-1}-2c^{i}+c^{i+1}}{h^2} - O (h^2) \\ \\ $$

and forward difference for the first derivative

$$ \displaystyle \frac{\partial c}{\partial h } \ \approx \ \frac{c^{i} - c^{i-1}} { h} - O(h) $$

with all BC's implemented i get the following system of equations:

$$ \displaystyle \frac{\partial c}{\partial t } = \frac{-D_{e}}{h^2} \begin{bmatrix} 0 & 0 & \dots & 0 & 0 \\ 1 & -2 & 1 & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 1 & -2 & 1 \\ 0 & \dots & 0 & -2 & 2 \\ \end{bmatrix} +\frac{U_{z}}{h} \begin{bmatrix} 0 & 0 & \dots & 0 & 0 \\ -1 & 1 & 0 & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & -1 & 1 & 0 \\ 0 & \dots & 0 & -1 & 1 \\ \end{bmatrix} \\ $$ I plugged this system o equqtion into matlab's ODE 15s and it worked fine.

Now my problem is the following: because of the inconsistency of the error in $O(h^{2})$ and $O(h)$ i wanted to try out a higher order discretization scheme for the first derivative

$$ \displaystyle \frac{\partial c}{\partial z } \ \approx \ \frac{3c^{i} - 4c^{i-1} + c^{i-2}} { 2h} - O(h^2) $$

So my questions are the following:

1. As far as i understood, with this type of matrix organisation i can set the first row of the matrix all zero because $\frac{\partial c}{\partial t } = 0 $ for my boundary nodes(z=0).

2. How to handle this equation for the second gridpoint, especially the $c^{i-2}$ in this equation is their an approach to incorporate this higher order discretization scheme properly $$ \displaystyle \frac{\partial c}{\partial z } \ \approx \ \frac{3c^{2} - 4c^{1} + *c^{i-2}*} { 2h} - O(h^2) $$

3. Is there a beter and/or more accurate way to implement the ditrichlet bc

Thank You!

Answer 1

So if I get you right, what you are recommending is to reorganise the system of equations like this:

$$ \displaystyle \frac{\partial c}{\partial t } = \frac{-D_{e}}{h^2} \begin{bmatrix} 0 & 0 & \dots & 0 & 0 \\ 1 & -2 & 1 & \dots & 0 \\ 0 & 1 & -2 & 1 & \dots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 1 & -2 & 1 \\ 0 & \dots & 0 & 0 & 0 \\ \end{bmatrix} \begin{bmatrix} c_{1} \\ c_{2} \\ 0 \\ 0 \\c_{N-1} \\ c_{N} \end{bmatrix} + \begin{bmatrix} c_{10} \\ 0 \\ 0 \\ 0 \\0 \\ X_{1} \end{bmatrix} +\frac{U_{z}}{2h} \begin{bmatrix} 0 & 0 & \dots & 0 & 0 \\ -4 & 3 & 0 & \dots & 0 \\ 1 & -4 & 3 & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 1 & -4 & 3 & 0 \\ 0 & \dots & 0 & 0 & 0 \\ \end{bmatrix} \begin{bmatrix} c_{1} \\ c_{2} \\ 0 \\ 0 \\c_{N-1} \\ c_{N} \end{bmatrix} + \begin{bmatrix} c_{10} \\ X_2 \\ 0 \\ 0 \\ 0\\ X_3 \end{bmatrix} $$ As far as I understand with this notation all the boundary nodes would not be incorporated into the state vector but into the inhomogeneous part of the system of equations and this is what you are recomending?

If this is what you meant, I don't know how to derive the equations for the values for $X_1$, $X_2$ and $X_3$