# Stability Optimization of Explicit Runge-Kutta Methods with Higher-Order Derivatives
## 5 Numerical Examples
### The Problem for the Lorenz System with Diffusion  
The following system of nonlinear parabolic PDEs is considered [45]:  
\[\frac{\partial u}{\partial t}=D_{1}\frac{\partial^{2}u}{\partial x^{2}}-\xi u+ \xi v,\frac{\partial v}{\partial t}=D_{2}\frac{\partial^{2}v}{\partial x^{2}}+ \eta u-v-uv,\frac{\partial w}{\partial t}=D_{3}\frac{\partial^{2}w}{\partial x ^{2}}-\zeta w+uv, \tag{20}\]  
where \(t>0\), \(x\in(a,b)\), with the following boundary conditions:  
\[\frac{\partial u}{\partial x}(t,a)=\frac{\partial u}{\partial x}(t,b)=\frac{ \partial v}{\partial x}(t,a)=\frac{\partial v}{\partial x}(t,b)=\frac{ \partial w}{\partial x}(t,a)=\frac{\partial w}{\partial x}(t,b)=0,\]  
and initial conditions:  
\[u(0,x)=v(0,x)=w(0,x)=\varphi_{0}.\]  
The following parameter values are used [45]: \(a=0\), \(b=1\), \(\xi=10\), \(\eta=24.8\), \(\zeta=8/3\), \(\varphi_{0}=8.1\). Computations are performed for the following values of diffusion coefficients: \(D_{1}=D_{2}=D_{3}=10^{-2}\). The time interval \([0,150]\) is considered. A space grid with \(N=200\) nodes is considered. For the construction of finite-difference schemes for (20), we use the approximations of second derivatives, as for Equation (19). Derivatives in boundary conditions are approximated by the first-order finite differences, so, in boundary nodes, the extrapolated boundary conditions are used. It is easy to demonstrate that Jacobian and Hessian matrices have a sparse structure, so only vectors, which represent products \(f^{\prime}f\), \((f^{{}^{\prime}})^{2}f\), and \(f^{{}^{\prime\prime}}(f,f)\), are calculated. In Table 8, the values of parabolic CFL number \(\tilde{\gamma}\), obtained in computational experiments, are presented. The largest values of \(\tilde{\gamma}\) are observed for the methods corresponding to optimal parameters, which maximize length \(\mathcal{L}\). This is explained by the behavior of the spectra of Jacobian matrices, presented for the final time moment in Figure 8. As can be seen, despite the complex structure, the eigenvalues are situated along the real negative axis. Figure 7: Spectrum of matrix \(\mathbf{B}\), corresponding to sixth-order approximations of the second derivative. In order to analyze the effect of the obtained parameters on the considered methods, computations on the practical estimation of local error are performed. For the experimental estimation of the error, [PERSON]'s extrapolation, which proposes the estimation of the  
Figure 8: Spectrum of Jacobian matrix, corresponding to high-order approximations of the second derivative: (**a**) fourth-order approximation; (**b**) fifth-order approximation; (**c**) sixth-order approximation. leading term of error, is used [27]. For this purpose, in the computation with constant step size in every node, we obtain a solution with step size \(\Delta t\) and a solution with two step sizes \(\Delta t/2\). For these solutions, notations \(y_{\Delta t}\) and \(y_{\Delta t/2}\) are used. The norm of the principal term of the local error of the \(p\)-th order method is estimated as [27]  
\[R=\frac{||y_{\Delta t}-y_{\Delta t/2}||}{2^{p}-1}.\]  
Computations are performed on the time interval \([0,10]\) on the same spatial grid. The values of \(\gamma\) are chosen in order to guarantee the stability for every parameter set. For GJRKMs, it is equal to 0.2. For TDRKMs, it is equal to 0.6. For ThDRKMs with \(s=2\) and \(s=2\), 0.7 and 0.8 are used, respectively. In Figure 9, plots of the logs of the absolute values of differences in experimental local error estimations are presented. As can be seen, for all methods, the values of errors corresponding to optimized and non-optimized methods are very close to each other (maximums correspond to order \(10^{-7}\)), so stability optimization does not realize any serious effect on the method accuracy.