# High order unconditionally stable discretization for a scalar hyperbolic PDE

In order to numerically solve the following differential equation: \begin{equation} \text{Fr}\{f\} := v(k)\dfrac{\partial f(z,k)}{\partial z} - F(z) \dfrac{\partial f(z,k)}{\partial k} = -\dfrac{f-f_0}{\tau} \end{equation} I have used the finite volume method, and have discretized the left hand side to: \begin{equation} \text{Fr}'\{f\} := v(k_j)\Delta k \Bigg[ f(z_{i+},k_j) - f(z_{i-},k_j) \Bigg] - F(z_i)\Delta z \Bigg[ f(z_{i},k_{j+}) - f(z_{i},k_{j-}) \Bigg] \end{equation} for every box.

The so-called flux averaging approximation, \begin{align} & f(z_{i+},k_j) = \dfrac{f(z_i,k_j)+f(z_{i+1},k_j)}{2} \\ & f(z_{i-},k_j) = \dfrac{f(z_i,k_j)+f(z_{i-1},k_j)}{2} \end{align} will lead to instability especially if the flux term is weak. This could be avoided by applying the upwind scheme. In this case: \begin{equation} f(z_{i+},k_j) = \begin{cases} f(z_{i},k_j) & v(k_j)>0 \\ f(z_{i+1},k_j) & v(k_j)<0 \end{cases} \end{equation} \begin{equation} f(z_{i-},k_j) = \begin{cases} f(z_{i-1},k_j) & v(k_j)>0 \\ f(z_{i},k_j) & v(k_j)<0 \end{cases} \end{equation} I have implemented the above upwinding method, and my results seem accurate for very small values of $F(z_i)$ throughout the system. However, when $F(z_i)$ gets large, the obtained results lose accuracy and deviate from the correct result.

The reason to this deviation, I guess, is that the analytical solution to $\text{Fr}\{f\}=0$ is an exponential function. Since the equations are linearly discretized, the discretization cannot follow the large exponential changes accurately.

Do you have any idea for increasing the accuracy of the above discretization while holding on to unconditional stability?

• Your equation looks similar to the one solved in the case of Lattice Boltzmann methods. Perhaps that is a way to go? Otherwise you may be interested in some schemes like QUICK? – nluigi Apr 29 '16 at 20:41
• Can you please give more information about boundary conditions, the number of grid points you used and exactly how inaccurate the solution is when $F(z_i)$ is increased? This may help with people being able to better answer your question. – Charles Apr 30 '16 at 16:55
• Also, it looks like you are missing $dA$ in your equations where $dA$ is the differential "area" from $dz$ and $dk$. – Charles Apr 30 '16 at 16:57

Let us think that $f=f(z,k,t)$, $f(z,k,0)=f_0(z,k)$ and we search for $f(z,k,\tau)$ that fulfills for $t \in (0,\tau)$ the advection equation (I use a shorter notation for partial derivatives): \begin{equation} \partial_t f + v(k) \partial_z f - F(z) \partial_k f = 0 \end{equation}
Then your upwind scheme is so called implicit full upwind (the first order accurate) method, in which you try to solve for $f(z,k,\tau)$ the advection equation in one time step. For small values of $v$ and $F$ it can work, i.e. the results can have satisfactory accuracy, but for large values the accuracy must be very poor.
Although the implicit upwind scheme has no restriction on the time step (say your $\tau$) due to stability reasons, it has a natural restriction on it due to accuracy. In fact, the standard CFL condition is pretty good guess for such accuracy restriction, it means $$\tau \le h \max \limits_{i,j} \left(|v(k_j)| + |F(z_i)| \right)^{-1}$$ where $h=\min \{\Delta k, \Delta z\}$.
If for any serious reasons your $\tau$ is given and fixed, you might look to your problem like the suggested advection equation and solve it with Lagrangian methods if you significantly violate the above condition for $\tau$. Such type of method can be implemented with arbitrary large $\tau$.
• Your first equation is one dimensional, the second one is two-dimensional, so you need a literature on two-dimensional advection equation. The much simpler one dimensional case is typically studied in the form $$\dfrac{\partial f}{\partial t} - c(x) \dfrac{\partial f}{\partial x}=0$$ and the general (linear) two-dimensional case $$\dfrac{\partial f}{\partial t} - a(x,y) \dfrac{\partial f}{\partial x} - b(x,y) \dfrac{\partial f}{\partial y}=0$$. – Peter Frolkovič May 3 '16 at 12:24