# Stability Optimization of Explicit Runge-Kutta Methods with Higher-Order Derivatives
## 5 Numerical Examples
### The Problem for the Discrete Kinetic Equations  
The following nonlinear hyperbolic system, used in the lattice Boltzmann method for the simulation of compressible gas flows, is considered [46]:  
\[\frac{\partial f_{i}}{\partial t}+V_{i}\frac{\partial f_{i}}{\partial x}=- \frac{1}{\lambda}\Big{(}f_{i}-f_{i}^{(eq)}\Big{)},\quad i=1,\ldots,n, \tag{21}\]  
where \(f_{i}\) are distribution functions of particles with velocities \(V_{i}\), \(\lambda\) is the relaxation time, and \(f_{i}^{(eq)}\) are equilibrium distribution functions. System (21) describes the 1D dynamics of an ideal gas with the equation of state, written as \(T=(\kappa-1)e\), where \(T\) is the temperature, \(e\) is the internal energy, and \(\kappa\) is a specific-heat ratio. The macrovariables, such as density \(\rho\), velocity \(u\), and full energy, are computed as  
\[\rho(t,x)=\sum_{i=1}^{n}f_{i}(t,x),\quad\rho(t,x)u(t,x)=\sum_{i=1}^{n}f_{i}(t,x)V_{i}, \tag{22}\]  
\[\rho(t,x)\bigg{(}e(t,x)+\frac{u^{2}(t,x)}{2}\bigg{)}=\frac{1}{2}\sum_{i=1}^{ n}f_{i}(t,x)(V_{i}^{2}+\eta_{i}^{2}), \tag{23}\]  
Figure 10: Work-precision plots for fourth-order GJRKM, applied to the problem for the Lorenz system. where \(\eta_{i}\) are the parameters for control of the value of \(\kappa\left(\frac{\eta_{i}^{2}}{2}\right)\) play the roles of energy levels for extra degrees of freedom). For the simulation of inviscid flow, the model with \(n=5\) can be used [46]. For the numerical simulations, the so-called D2Q5 lattice is used, which is defined by the following discrete velocities:  
\[V_{i}= \begin{cases}0,\quad i=1,\\ V_{1}\cos(\pi i),\quad i=2,3,\\ V_{2}\cos(\pi i),\quad i=4,5,\end{cases} \]  
and the following values of \(\eta_{i}\) are used:  
\[\eta_{i}= \begin{cases}\eta_{0},\quad i=1,\\ 0,\quad i=2,\ldots,5.\end{cases} \]  
The equilibrium distribution functions are presented as [46]  
\[f_{i}^{(eq)}=\rho(A_{i}+B_{i}uV_{i}),\quad i=1,\ldots,5, \tag{24}\]  
where \(A_{i}\) and \(B_{i}\) are computed as  
\[A_{i}= \begin{cases}\frac{b-1}{\eta_{0}^{2}}T,i=1,\\ \frac{1}{2(V_{1}^{2}-V_{2}^{2})}\left(-V_{2}^{2}+\left((b-1)\frac{V_{2}^{2}}{ \eta_{0}^{2}}+1\right)T+u^{2}\right),\quad i=2,3,\\ \frac{1}{2(V_{2}^{2}-V_{1}^{2})}\left(-V_{1}^{2}+\left((b-1)\frac{V_{1}^{2}}{ \eta_{0}^{2}}+1\right)T+u^{2}\right),\quad i=4,5.\end{cases} \]  
\[B_{i}= \begin{cases}\frac{-V_{2}^{2}+(b+2)T+u^{2}}{2V_{1}^{2}V_{2}^{2}-V_{2}^{2 }},\quad i=2,3,\\ \frac{-V_{1}^{2}+(b+2)T+u^{2}}{2V_{2}^{2}(V_{2}^{2}-V_{1}^{2})},\quad i=4,5, \end{cases} \]  
where \(b=\frac{2}{\kappa-1}\). For computational experiments, the test problem from [47] is chosen. The interval \(x\in(0,1)\) is considered with periodic boundary conditions for \(f_{i}\). Initial conditions for macrovariables are stated as  
\[\rho(0,x)=1+0.2\sin(2\pi x),\quad u(0,x)=0.2,\quad p(0,x)=1,\]  
where \(p\) is a pressure. Initial values for \(f_{i}\) are computed from the values of macrovariables as \(f_{i}^{(eq)}\). The value of \(\kappa\) is stated as 1.4. For computations, the time interval \([0,1]\) is considered. The following values of model parameters are used: \(\eta_{0}=2\), \(V_{1}=V\), \(V_{2}=2V\), \(V=1\), \(\lambda=10^{-3}\). It should be noted that, according to the nonlinear dependence of \(f_{i}^{(eq)}\) on all functions \(f_{i}\) (see (22)-(24)), analytic formulas for the partial derivatives of \(f_{i}^{(eq)}\) on \(f_{i}\) are quite cumbersome, and their computations are expensive. So, for this problem, we do not consider methods based on the third derivative. In the group of methods, based on the second derivative, only GJRKMs are considered, according to only one computation of \(f^{{}^{\prime}}\).