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

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!

• Thank's ! This was of course a false notation!
– Ivan
Jan 20, 2021 at 10:45
• Why the backwards difference and not the central difference quotient? The backwards and forwards formulas are usually only used to implement boundary conditions. // You should separate the computation for the inner nodes from the one for the boundary nodes. Better even, take the boundary nodes out of the state vector and use the BC to compute their values for the difference quotients. Jan 20, 2021 at 11:52
• 1. I am using the backwards difference because the problem I am investigating is driven by the advection term, I found sources that recommend using the upwind scheme for advection driven problems. en.wikipedia.org/wiki/Upwind_scheme (Of course I found other Sources as well, this is just the fastest I could find).
– Ivan
Jan 20, 2021 at 12:48
• 2. how do I separate the computation of the inner cell from that of the boundary node? Could you provide some literature recommendation of how to realise that for a non-steady-state problem? If I take the boundary nodes out of the state vector and compute them separately how would the matrix look like? thank you for your thoughts!
– Ivan
Jan 20, 2021 at 12:56
• You would get an inhomogeneous term in $\dot u_{1..N}=Lu_{0..N+1}=\bar Lu_{1..N}+b(u_0,u_{N+1})$. As the BC here are homogeneous linear, $u_0=0$, $u_{N+1}=(4u_N-u_{N-1})/3$, one can fold the second term back into the first to get $\dot u_{1..N}=\tilde Lu_{1..N}$. // If they are finite in number, the global order is not influenced if the local truncation order drops by $1$. With that you should be able to solve the $i=1$ problem, just leave it at the order 1 difference quotient there. Or use the central formula in this node only. Jan 20, 2021 at 13:28

$$\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$$