Method of Lines Runge-Kutta nonlinear stability and behavior

Question

I have a system of 4 nonlinear 1st-order PDEs. I want to solve them numerically by method of lines, first discretizing space. This leads to the system of $N\times 4$ coupled ODEs.

$$ \mathbf{u}_{i} = (u^{1},u^{2},u^{3},u^{4})_{i},\quad i=1,\dots,N\\ \frac{\mathrm{d}u_{i}^{1}}{\mathrm{d}t} = q_{1}(\mathbf{u}_{i}) \\ \frac{\mathrm{d}u_{i}^{2}}{\mathrm{d}t} = q_{2}(\mathbf{u}_{i}) \\ \frac{\mathrm{d}u_{i}^{3}}{\mathrm{d}t} = F_{1}(\mathbf{u}_{i-k},\dots,\mathbf{u}_{i+l}) \\ \frac{\mathrm{d}u_{i}^{4}}{\mathrm{d}t} = F_{2}(\mathbf{u}_{i-k},\dots,\mathbf{u}_{i+l}) $$ where $F$ represents a choice of spatial discretization of the two conservation equations, which will depend on neighboring points. The solutions are traveling waves in all four fields. I plan to use a 4th-order Runge-Kutta for the time integration but with varying spatial schemes.

I have read that the stability for nonlinear FD schemes is governed by the local `amplification' matrix. I have no idea how to find the amplification matrix for this system, as the RK4 time step is rather complicated.

Also, can the behavior of the solution be expressed generally in terms of $F$? I expect that the amount of numerical dispersion and diffusion will be related to the accuracy of the spatial scheme, but $F$ by itself is just an expression over the spatial grid points and thus says nothing about how Fourier modes will be affected by the time-stepping scheme.

Answer 1

It is useful, as a first step, to analyze method-of-lines (MOL) treating the time integration as exact. For example, suppose we are solving the 1D advection equation

$ \frac{\partial{n}}{\partial{t}} =-c \frac{\partial{n}}{\partial{x}}, $

where $c$ is the advection speed.

Then, for example, for central difference, we'd have

$ \frac{d}{dt} n_i = - \frac{c}{2 h} [n_{i+1} - n_{i-1}], $

where $h$ is the spatial grid spacing, $i$ is the grid point index, and making the ansatz $n_i(t) = e^{I k h i} e^{-I \omega t}$, where $I$ is the imaginary unit, we'd obtain

$ \omega = c \frac{\sin(k h)}{h}, $

which is the dispersion relation for the spatially-discretized equation. From this analysis, one can see that there are no unstable modes here, and in the limit of small $kh$ it becomes the exact disperion relation, $\omega = kc$. On the other hand, for $k h = \pi$ the modes don't propagate, that's the infamous zigzag (Nyquist) mode.

For multiple nonlinear equations the analysis is conceptually similar but the procedure would be to linearize the system and spatially discretize it, and then find the eigenvalues; that would usually provide some useful insights.

Here is another illustration, for a system of linear PDEs describing the acoustic wave,

$ \frac{\partial{u}}{\partial{t}} =-\frac{c^2}{n_0} \frac{\partial{n}}{\partial{x}}, \\ % \frac{\partial{n}}{\partial{t}} =-n_0 \frac{\partial{u}}{\partial{x}}, $

Using MOL for the PDEs in the normalized form, and using the Fourier harmonics ansatz, the system discretized in space by central difference becomes

$ \frac{d}{dt} \begin{bmatrix} u_i \\ n_i \\ \end{bmatrix} = \begin{bmatrix} 0 & -I \frac{\sin(k h)}{k h} \\ -I \frac{\sin(k h)}{k h} & 0\\ \end{bmatrix} % \begin{bmatrix} u_i \\ n_i \\ \end{bmatrix} $

which leads to an eigenvalue problem for the frequency

Analyzing the scheme for a specific time integration algorithm can be carried out as well. For the acoustic wave example, for the first-order explicit Euler with step size $\tau$, we'd have

$ \frac{e^{-I \omega \tau} -1}{\tau} \begin{bmatrix} u_i \\ n_i \\ \end{bmatrix} = \begin{bmatrix} 0 & -I \frac{\sin(k h)}{k h} \\ -I \frac{\sin(k h)}{k h} & 0\\ \end{bmatrix} % \begin{bmatrix} u_i \\ n_i \\ \end{bmatrix} $

Solving for the eigenvalues $\lambda_k$ we'd arrive at the dispersion relation

$ \frac{e^{-I \omega \tau} -1}{\tau} = \lambda_k $

which is the relation between $\omega$ and $k$ describing time-evolution of individual Fourier modes, for the chosen spatial discretization and time-integration method.