Stability Analysis

Question

The partial differential equation, \begin{align} \dfrac{\partial f}{\partial x} + a(x)\dfrac{\partial f}{\partial y} = 0 \qquad & f(0,y) = f(L_1,y) = c_0e^{-y} \\ & f(x,0) = c_0 \;,\; f(x,L_2) = c_0e^{-L_2} \end{align} for $f=f(x,y)$ has the nontrivial solution, \begin{equation} f(x,y) = Ae^{(-y-h(x))} \;,\quad a(x)=-\dfrac{dh(x)}{dx} \;. \end{equation} If I discretize the equation by box integration over $[y_{i},y_{i+1}]$, I get: \begin{equation} \dfrac{\Delta y}{2} \left[\dfrac{\partial f(x,y_{i})}{\partial x} + \dfrac{\partial f(x,y_{i+1})}{\partial x}\right] + a(x)\big[f(x,y_{i+1}) - f(x,y_i)\big] = 0 \end{equation} which converges to the nontrivial solution as $\Delta y\rightarrow 0$. However, I have noticed such discretizations are numerically unstable for $a(x)\neq 0$ and result in oscillations when solved on a computer. How can I investigate the above equation in terms of numerical stability? Does this equation belong to a specific class of PDEs, and are there techniques to improve its stability?

Answer 1

This is a first-order PDE but you are trying to impose Dirichlet BC on the whole boundary. This is not well-posed in general but maybe your data is such that you expect it to have a unique solution. You can write this PDE as $$ v \cdot \nabla f = \nabla\cdot (v f) = 0 $$ where $$ v(x,y) = (1, a(x)) $$ This is a hyperbolic PDE. If $n$ is outward normal to your domain boundary, then you can specify Dirichlet BC only on those portions of the boundary where $v \cdot n < 0$. So you can specify BC on the left boundary $(x=0)$ but not the right boundary. On the bottom and top boundaries, you have to check the sign of $a(x)$ to decide this.

Draw the characteristic curves $$ dy/dx = a(x) $$ to see how BC will flow into your domain. This can identify a time-like direction along which you can march to find the solution.

For example, if $a(x)=x$ and your domain is $[0,L_1]\times[0,L_2]$ so that $a(x) \ge 0$ in your domain, then you can specify BC on the left and bottom sides. You can march along $x$. As Wolfgang wrote, for $\partial f/\partial y$ you can use backward-difference type scheme. An FD scheme may look like this $$ \frac{f_{i+1,j}-f_{i,j}}{\Delta x} + a_i \frac{f_{i,j} - f_{i,j-1}}{\Delta y} = 0, \qquad a_i \ge 0 $$ Note that $f_{0,j}$ and $f_{i,0}$ are the given boundary values and you determine the remaining values using the scheme above.

This is only first-order accurate but a good starting point.