I'm trying to numerically solve the following fairly simple non-linear PDE:
\begin{equation}
\frac{\partial \omega}{\partial t}=-V_g(\omega)\frac{\partial \omega}{\partial z}
\end{equation}
I've discretised this as follows:
\begin{align}
\frac{\omega(z,t+\Delta t) – \omega(z,t)}{\Delta t} &= -V_g(\omega(z,t))\frac{\omega(z+\Delta z,t) – \omega(z-\Delta z,t)}{2\Delta z}\\
\omega(z,t+\Delta t) &=\omega(z,t)-\Delta t V_g(\omega(z,t))\frac{\omega(z+\Delta z,t) – \omega(z-\Delta z,t)}{2\Delta z}
\end{align}
which is simple enough. However I've managed to derail myself when it comes to the boundary conditions. The initial conditions are essentially the delta function:
\begin{equation}
\omega(z,0) = \delta(z)
\end{equation}
This is where I come unstuck. With everything as written above, the system won't evolve, because $V_g(\omega)$ is, based on these boundary conditions, zero everywhere except $z=0$, and so no evolution occurs. I'm fairly sure the issue is with how I'm handling $V_g(\omega)$, but for the life of me I can't figure out what the solution is.
Best Answer
The equation $\partial_t\omega + V_g(\omega)\partial_z\omega = 0$ is a first-order quasilinear PDE. This nonlinear transport equation carries information at the characteristic speed $V_g(\omega)$. Let us consider a regular grid with uniform mesh size $\Delta z = z_{i+1} -z_{i}$, and variable time step $\Delta t = t_{n+1} -t_{n}$. We write $\omega_i^n \simeq \omega(z_i, t_n)$. The finite-difference scheme $$ \frac{\omega_i^{n+1}-\omega_i^n}{\Delta t} = V_g (\omega_i^n) \frac{\omega_{i+1}^n-\omega_{i-1}^n}{2\Delta x} $$ written in OP is unstable (see e.g. [1]). The small modification $$ \frac{\omega_i^{n+1}-\frac{1}{2}(\omega_{i-1}^n+\omega_{i+1}^n)}{\Delta t} = V_g (\omega_i^n) \frac{\omega_{i+1}^n-\omega_{i-1}^n}{2\Delta x} $$ makes it conditionally stable. This is the non-conservative Lax-Friedrichs method, stable under the Courant-Friedrichs-Lewy (CFL) condition $\max_i V_g (\omega_i^n)\frac {\Delta t}{\Delta x} < 1$. Viewed as a finite-volume method, one should initialize the data to the cell-averages of $\omega (z,0)=\delta (z)$, i.e. $\omega_i^{0} = \delta_i^{0}/\Delta z$, where $\delta_i^{0}$ is the Kronecker delta.
An alternative consists in writing the nonlinear transport equation above in conservative form. If we define the flux $F_g(\omega) = \int^\omega V_g(\nu)\, \text{d}\nu$ as an antiderivative of $V_g$, then the transport equation rewrites as a conservation law: $\partial_t\omega + \partial_z F_g(\omega) = 0$. Thus, one can implement straightforwardly finite-volume methods such as the conservative Lax-Friedrichs method $$ \frac{\omega_i^{n+1}-\frac{1}{2}(\omega_{i-1}^n+\omega_{i+1}^n)}{\Delta t} = \frac{F_g (\omega_{i+1}^n)-F_g (\omega_{i-1}^n)}{2\Delta x} . $$
[1] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, 1992.