How do we find finite difference schemes to numerically solve this B.V.P?

Question

\begin{equation} \frac{\partial P}{\partial t}=\frac{\partial}{\partial x}\left[D(x)\frac{\partial P}{\partial x}\right]-v\frac{\partial P}{\partial x}-\tau P+S \end{equation} and \begin{equation} \frac{\partial Q}{\partial t}=\frac{\partial}{\partial x}\left[D(x)\frac{\partial Q}{\partial x}\right]-v\frac{\partial Q}{\partial x}-\tau P, \end{equation} where $t$ is time, $x$ is the spatial variable along the column, $D(x)$ is the diffusion rate, $v$ is the advection rate, $\tau$ is the decay rate, and $S(x)$ is the source rate. We will be looking for the late-time, steady state solution to these equations with $\partial P/\partial t=\partial Q/\partial t = 0$. The diffusion rate and source rate may be assumed to be smooth and non-negative. We will also assume $v$ is non-negative.

The $x$ domain is $x\in[0,L]$. Boundary conditions for $P$ and $Q$ are mixed at $x=0$, representing zero influx \begin{equation} \left.\left(D\frac{\partial P}{\partial x}-vP\right)\right|_{x=0}=\left.\left(D\frac{\partial Q}{\partial x}-vQ\right)\right|_{x=0}=0. \end{equation} At $x=L$ the boundary condition will be $P=Q=c$, where $c$ is a constant.

The question is we wish to find a finite difference scheme (break up the domain into $N$ gridpoints ranging from $0$ to $N-1$ and form index equations for each point) in order to create matrices such as $A$ where $AP = S$ then numerically solve this equation. I'm wondering which schemes would be most useful. I have thought of using centred/forward/backward finite difference scheme but I'm very confused about how to input the boundary conditions and which to use?

We require a tridiagonal matrix A. e.g as produced in this example. values of D(x) and S(x) are provided at all points on the index grid.

https://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/node9.html

Answer 1

You can translate that almost literally. For the second derivative apply twice central difference quotients with step size $\Delta x$, that is, half-steps from the center. This will lead to $$ \left.\frac{\partial}{\partial x}\left(D(x)\frac{\partial P(x,t)}{\partial x}\right)\right|_{(x,t)=(x_i,t_j)} \approx \frac{D_{i+\frac12}(P^j_{i+1}-P^j_{i})-D_{i-\frac12}(P^j_{i}-P^j_{i-1})}{Δx^2} $$ For the first derivative take the central difference quotient with step size $2Δx$, $$ \left.\frac{\partial P(x,t)}{\partial x}\right|_{(x,t)=(x_i,t_j)} \approx \frac{P^j_{i+1}-P^j_{i-1}}{2Δx} $$ For the right boundary conditions you set that component as constant and remove it from the active computation. For the left condition it is probably easiest to introduce ghost cells at index $i=-1$ so that the PDE can be applied there and also the boundary condition can be constructed with central difference quotients.

Answer 2

1) Premise

Consider the function $P(x,t)$ to correspond to the matrix $$ P(x,t) \Rightarrow {\bf P} = \left( {\matrix{ {P(0,0)} & {P(0,\Delta t)} & {P(0,2\Delta t)} & \cdots \cr {P(\Delta x,0)} & {P(\Delta x,\Delta t)} & {P(\Delta x,2\Delta t)} & \cdots \cr {P(2\Delta x,0)} & {P(2\Delta x,\Delta t)} & {P(2\Delta x,2\Delta t)} & \cdots \cr \vdots & \vdots & \vdots & \ddots \cr } } \right) $$ where the starting point has been taken to be $(0,0)$, but could of course be shifted.

Then consider ${\bf E}, \overline {\bf E}$ to be the shift matrices $$ {\bf E} = \left( {\matrix{ 0 & 0 & 0 & \cdots \cr 1 & 0 & 0 & \cdots \cr 0 & 1 & 0 & \cdots \cr \vdots & \vdots & \ddots & \ddots \cr } } \right)\quad \overline {\bf E} = \left( {\matrix{ 0 & 1 & 0 & \cdots \cr 0 & 0 & 1 & \cdots \cr 0 & 0 & 0 & \ddots \cr \vdots & \vdots & \vdots & \ddots \cr } } \right) $$ (we denote the transpose with an overline)

So we will have $$ \left\{ \matrix{ {\partial \over {\partial x}}P(x,t) \Rightarrow {{P_{\,n,m} - P_{\,n - 1,m} \,} \over {\Delta x}} = {1 \over {\Delta x}}\nabla _x P(x,t) = {1 \over {\Delta x}}\left( {{\bf I} - {\bf E}} \right){\bf P} \hfill \cr {\partial \over {\partial x}}P(x,t) \Rightarrow {{P_{\,n + 1,m} - P_{\,n,m} \,} \over {\Delta x}} = {1 \over {\Delta x}}\Delta _x P(x,t) = {1 \over {\Delta x}}\left( {\overline {\bf E} - {\bf I}} \right){\bf P} \hfill \cr {\partial \over {\partial t}}P(x,t) \Rightarrow {1 \over {\Delta t}}\nabla _t P(x,t) = {1 \over {\Delta t}}{\bf P}\left( {{\bf I} - \overline {\bf E} } \right) \hfill \cr {\partial \over {\partial t}}P(x,t) \Rightarrow {1 \over {\Delta t}}\Delta _t P(x,t) = {1 \over {\Delta t}}{\bf P}\left( {{\bf E} - {\bf I}} \right) \hfill \cr} \right. $$

It is to be underlined that, with matrices of finite dimension ($[0,h]^2$) we get $$ \left( {{\bf I} - {\bf E}} \right){\bf P} = \left( {\matrix{ {P_{\,0,\,0} - 0} & {P_{\,0,\,1} - 0} & \cdots & {P_{\,0,\,h} - 0} \cr {P_{\,1,\,0} - P_{\,0,\,0} } & {P_{\,1,\,1} - P_{\,0,\,1} } & \cdots & {P_{\,1,\,h} - P_{\,0,\,h} } \cr \vdots & \vdots & \ddots & \vdots \cr {P_{\,h,\,0} - P_{\,h - 1,\,0} } & {P_{\,h,\,1} - P_{\,h - 1,\,1} } & \cdots & {P_{\,h,\,h} - P_{\,h - 1,\,h} } \cr } } \right) \Rightarrow {\partial \over {\partial x}} H(x)P(x,t) $$ meaning that it actually reproduces the Backward difference of $P(x,t)$ with $P(x,t)=0 \; |x<0$, i.e. for $P(x,t)$ multiplied by the Heaviside step function.
While the Forward difference is suited to $P(x,t)$ cropped for $h < x$ $$ \left( {\overline {\bf E} - {\bf I}} \right){\bf P} = \left( {\matrix{ {P_{\,1,\,0} - P_{\,0,\,0} } & {P_{\,1,\,1} - P_{\,0,\,1} } & \cdots & {P_{\,1,\,h} - P_{\,0,\,h} } \cr {P_{\,2,\,0} - P_{\,1,\,0} } & {P_{\,2,\,1} - P_{\,1,\,1} } & \cdots & {P_{\,2,h} - P_{\,1,\,h} } \cr \vdots & \vdots & \ddots & \vdots \cr {0 - P_{\,h,\,0} } & {0 - P_{\,h,\,1} } & \cdots & {0 - P_{\,h,\,h} } \cr } } \right) \Rightarrow{\partial \over {\partial x}} \left( {1 - H(x - h - 1)} \right)P(x,t) $$ and analogously for the differences in $t$.

The Central difference expressed as $$ \eqalign{ & \mathop \diamondsuit \nolimits_x P(x,t) = {{P_{\,n + 1,m} - P_{\,n - 1,m} \,} \over {2\Delta x}} = {1 \over 2}\left( {\nabla _x + \Delta _x } \right)P(x,t)\; \Rightarrow \cr & \Rightarrow \;{1 \over {2\Delta x}}\left( {\overline {\bf E} - {\bf E}} \right){\bf P} = {1 \over {2\Delta x}}\left( {\matrix{ {P_{\,1,\,0} - 0} & {P_{\,1,\,1} - 0} & \cdots & {P_{\,1,\,h} - 0} \cr {P_{\,2,\,0} - P_{\,0,\,0} } & {P_{\,2,\,1} - P_{\,0,\,1} } & \cdots & {P_{\,2,h} - P_{\,0,\,h} } \cr \vdots & \vdots & \ddots & \vdots \cr {0 - P_{\,h - 1,\,0} } & {0 - P_{\,h - 1,\,1} } & \cdots & {0 - P_{\,h - 1,\,h} } \cr } } \right)\; \cr} $$ uses the matrix $\overline {\bf E} - {\bf E}$, which however is not invertible, and is suitable for a "boxed" $P(x,t)$ .

I am not going to list and discuss here the other formulations ("stencils") that may be adopted, and the important aspects of stability and convergence which are well summarized in these three Wikipedia articles:Finite diff., FD method, FD coefficients.
I will focus on the matrix formulation.

2) Example

Before coming to your specific equation let's analyze the simpler case of the [Heat equation] $$ \left\{ \matrix{ {\partial \over {\partial t}}T(x,t) = {{\partial ^{\,2} } \over {\partial x^{\,2} }}T(x,t) \hfill \cr T(0,t) = T(h\Delta x,t) = 0 \hfill \cr T(x,0) = T_{\,0} (x) \hfill \cr} \right. $$

Solution by the Explicit Method

Using the so called FTCS scheme, we put $$ \left\{ \matrix{ {\partial \over {\partial t}}T(x,t) \Rightarrow {1 \over {\Delta t}}\Delta _t T(x,t) = {1 \over {\Delta t}}{\bf T}\left( {{\bf E} - {\bf I}} \right) \hfill \cr {{\partial ^{\,2} } \over {\partial x^{\,2} }}T(x,t) \Rightarrow {{T_{\,n + 1,m} - 2T_{\,n,m} + T_{\,n - 1,m} \,} \over {\Delta x^{\,2} }} = \hfill \cr = {{\Delta _x T(x,t) - \nabla _x T(x,t)} \over {\Delta x^{\,2} }} = {1 \over {\Delta x^{\,2} }}\left( {\overline {\bf E} - 2{\bf I} + {\bf E}} \right){\bf T} \hfill \cr} \right. $$ thereby the PDE translates to $$ {\bf T}\left( {{\bf E} - {\bf I}} \right) = {{\Delta t} \over {\Delta x^{\,2} }}\left( {\overline {\bf E} - 2{\bf I} + {\bf E}} \right){\bf T} = r\left( {\overline {\bf E} - 2{\bf I} + {\bf E}} \right){\bf T} $$ which, due to the presence of both left and right multiplications on the variable $\bf T$, is solved by progressively determining the column vectors of $\bf T$ as $$ {\bf TE} = r\left( {\overline {\bf E} - 2{\bf I} + {\bf E}} \right){\bf T} + {\bf T} = \left( {1 - 2r} \right)\left( {{\bf I} + {r \over {\left( {1 - 2r} \right)}} \left( {\overline {\bf E} + {\bf E}} \right)\,} \right){\bf T} $$ that is $$ {\bf T}_{\,m + 1} \left( x \right) = \left( {1 - 2r} \right) \left( {{\bf I} + {r \over {\left( {1 - 2r} \right)}}\left( {\overline {\bf E} + {\bf E}} \right)\,} \right) {\bf T}_{\,m} \left( x \right)\quad \left| {\,0 \le m \le h - 1} \right. $$ where the two ends of the vector $$ {\bf T}_{\,m + 1} \left( 0 \right),\;{\bf T}_{\,m + 1} \left( {h\Delta x} \right) $$ are not actually computed, yet set to the imposed boundary values (here taken to be null).

Solution by the Implicit Method

This method , that takes for ${\partial / {\partial t}}$ the backward difference, would lead to $$ {\bf T}_{\,m + 1} \left( x \right) = {1 \over {1 + 2r}}\left( {{\bf I} - {r \over {\left( {1 + 2r} \right)}} \left( {\overline {\bf E} + {\bf E}} \right)} \right)^{\, - \,{\bf 1}} {\bf T}_{\,m} \left( x \right) $$ 3) Conclusion

Coming to the specific equations you are interested in, we put $D(x)$ as a diagonal matrix $$ D(x) \Rightarrow \left( {\matrix{ {D(0)} & 0 & 0 & \cdots \cr 0 & {D(1\,\Delta x)} & 0 & \cdots \cr 0 & 0 & {D(2\,\Delta x)} & \cdots \cr \vdots & \vdots & \vdots & \ddots \cr } } \right) = \left( {D(n\,\Delta x) \circ {\bf I}} \right) $$ then $$ D(x){\partial \over {\partial x}}P(x,t) \Rightarrow {1 \over {\Delta x}}\left( {D(n\,\Delta x) \circ {\bf I}} \right) \left( {\overline {\bf E} - {\bf I}} \right){\bf P} $$ if the forward difference is taken.

For the term $$ {\partial \over {\partial x}}\left( {D(x){\partial \over {\partial x}}P(x,t)} \right) $$ we have a even wider selection of "stencils" that can be applied, noting that there are slight differences between them, because for instance we have $$ \eqalign{ & {{\partial ^{\,2} } \over {\partial x^{\,2} }} \Rightarrow {1 \over {\Delta x^{\,2} }}\left( {\overline {\bf E} - 2{\bf I} + {\bf E}} \right) = {1 \over {\Delta x}}\left( {{1 \over {\Delta x}}\left( {\overline {\bf E} - {\bf I}} \right) - {1 \over {\Delta x}} \left( {{\bf I} - {\bf E}} \right)} \right) \approx \cr & \approx {1 \over {\Delta x}}\left( {{\bf I} - {\bf E}} \right){1 \over {\Delta x}} \left( {\overline {\bf E} - {\bf I}} \right) = {1 \over {\Delta x^{\,2} }} \left( {\overline {\bf E} - {\bf I} - {\bf E}\overline {\bf E} + {\bf E}} \right) \approx \cr & \approx {1 \over {\Delta x}}\left( {\overline {\bf E} - {\bf I}} \right){1 \over {\Delta x}}\left( {{\bf I} - {\bf E}} \right) = {1 \over {\Delta x^{\,2} }}\left( {\overline {\bf E} - {\bf I} - \overline {\bf E} {\bf E} + {\bf E}} \right) \cr} $$ since ${\bf E}\overline {\bf E}$ and $\overline {\bf E} {\bf E}$ differ from $\bf I$ by having one diagonal term null (the first at bottom, the second on top).
In practice one adopts the stencil whose (a-)symmetry matches that of $D(x)$ and of (the expected) $P(x)$.

Also, if $D(x)$ is known as a continuous function, the term can be expanded into $$ {\partial \over {\partial x}}\left( {D(x){\partial \over {\partial x}}P(x,t)} \right) = \left( {{\partial \over {\partial x}}D(x)} \right){\partial \over {\partial x}}P(x,t) + D(x){{\partial ^{\,2} } \over {\partial x^{\,2} }}P(x,t) $$

Assuming we can do that, and taking the forward difference for the first derivative in $t$, while for that in $x$ we take the central one supposing that we expect an almost central-symmetric $P$ vs. $x$, the equation $$ {\partial \over {\partial t}}P(x,t) = \left( {D'(x) - v} \right){\partial \over {\partial x}}P(x,t) + D(x){{\partial ^{\,2} } \over {\partial x^{\,2} }}P(x,t) - \tau P(x,t) + S\left( x \right) $$ translates into $$ \eqalign{ & {\bf P}_{m + 1} (x) = \left( {{{\Delta t} \over {2\Delta x}}\left( {\left( {D'(n\Delta x) - v} \right) \circ {\bf I}} \right) \left( {\overline {\bf E} - {\bf E}} \right) + {{\Delta t} \over {\Delta x^{\,2} }} \left( {D(n\Delta x) \circ {\bf I}} \right)\left( {\overline {\bf E} - 2{\bf I} + {\bf E}} \right) + \left( {1 - \tau \,\Delta t} \right)\,{\bf I}} \right){\bf P}_m (x) + {\bf S}\left( x \right) = \cr & = \left( {\left( {a(n) + b(n) \circ {\bf I}} \right)\overline {\bf E} + \left( {c(n) \circ {\bf I}} \right) + \left( { - a(n) + b(n) \circ {\bf I}} \right){\bf E}} \right){\bf P}_m (x) + {\bf S}\left( x \right) \cr} $$ where, as in the example , the edges ${\bf P}_{m + 1} (0),\;{\bf P}_{m + 1} (h\Delta x)$ shall be set to the prescribed boundary values.