Implicit-Explicit multistep methods for hyperbolic systems with multiscale relaxation

TL;DR

This paper develops high-order Implicit-Explicit multistep numerical methods for hyperbolic systems with relaxation, capable of accurately capturing multiscale behavior and asymptotic limits without restrictive time step constraints.

Contribution

It introduces novel high-order space-time discretizations based on IMEX multistep methods that effectively handle multiple scales and asymptotic regimes in hyperbolic systems with relaxation.

Findings

Methods successfully capture hyperbolic and parabolic limits

Numerical examples confirm theoretical accuracy and stability

Approach avoids time step restrictions in multiscale regimes

Abstract

We consider the development of high order space and time numerical methods based on Implicit-Explicit (IMEX) multistep time integrators for hyperbolic systems with relaxation. More specifically, we consider hyperbolic balance laws in which the convection and the source term may have very different time and space scales. As a consequence the nature of the asymptotic limit changes completely, passing from a hyperbolic to a parabolic system. From the computational point of view, standard numerical methods designed for the fluid-dynamic scaling of hyperbolic systems with relaxation present several drawbacks and typically lose efficiency in describing the parabolic limit regime. In this work, in the context of Implicit-Explicit linear multistep methods we construct high order space-time discretizations which are able to handle all the different scales and to capture the correct asymptotic behavior, independently from its nature, without time step restrictions imposed by the fast scales. Several numerical examples confirm the theoretical analysis.

Figures

Tables

Table 1. Table 1. Coefficients c r ​ j subscript 𝑐 𝑟 𝑗 c_{rj} for the 5 5 5 -th order WENO reconstruction on equispaced grid points.

r	$j=0$	$j=1$	$j=2$
-1	11/6	-7/6	1/3
0	1/3	5/6	-1/6
1	-1/6	5/6	1/3
2	1/3	-7/6	11/6

Table 2. Table 2. L 1 superscript 𝐿 1 L^{1} error and estimated convergence rates for u 𝑢 u in the AP-explicit case.

		$\epsilon =1$		$\epsilon =0.1$		$\epsilon =0.01$		$\epsilon =0.001$
IMEX	$N$	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate
IMEX-SG(3,2)	128	9.9402e-05	–	8.2324e-05	–	0.00019196	–	1.0662e-07	–
	256	3.0565e-05	1.7014	2.5215e-05	1.707	5.7846e-05	1.7305	4.262e-08	1.3229
	512	8.3499e-06	1.8721	6.8755e-06	1.8748	1.5671e-05	1.8841	1.1871e-08	1.8441
IMEX-BDF2	128	4.4737e-05	–	3.6789e-05	–	8.9909e-05	–	9.2432e-08	–
	256	1.4423e-05	1.6331	1.1759e-05	1.6455	2.8204e-05	1.6726	2.4676e-08	1.9053
	512	4.0033e-06	1.8491	3.2494e-06	1.8555	7.7212e-06	1.869	6.0694e-09	2.0235
IMEX-TVB(3,3)	128	7.3451e-08	–	1.2994e-07	–	5.0813e-06	–	3.1308e-08	–
	256	1.2953e-08	2.5035	1.8201e-08	2.8358	7.5974e-07	2.7416	3.8147e-09	3.0369
	512	1.8478e-09	2.8094	2.4055e-09	2.9197	1.0251e-07	2.8898	4.4859e-10	3.0881
IMEX-BDF3	128	8.7902e-08	–	1.6714e-07	–	6.6227e-06	–	5.8099e-08	–
	256	1.5388e-08	2.514	2.4439e-08	2.7738	9.8974e-07	2.7423	6.0486e-09	3.2638
	512	2.1902e-09	2.8127	3.2815e-09	2.8968	1.3293e-07	2.8964	6.5821e-10	3.2
IMEX-TVB(4,4)	128	2.2234e-08	–	6.9585e-09	–	1.2029e-06	–	2.0263e-07	–
	256	8.4006e-10	4.7262	1.7448e-10	5.3177	9.1633e-08	3.7145	5.8212e-09	5.1214
	512	4.7195e-11	4.1538	3.7747e-12	5.5305	6.0883e-09	3.9118	2.2316e-10	4.7052
IMEX-BDF4	128	2.3209e-08	–	7.4749e-09	–	4.5554e-07	–	2.449e-08	–
	256	6.8151e-10	5.0898	2.1477e-10	5.1212	3.3431e-08	3.7683	1.1948e-09	4.3573
	512	2.9394e-11	4.5352	5.499e-12	5.2875	2.2255e-09	3.909	6.3684e-11	4.2297
IMEX-TVB(5,5)	128	1.043e-08	–	4.1161e-08	–	9.4337e-09	–	1.9208e-07	–
	256	3.2924e-10	4.9854	1.2742e-09	5.0137	3.3621e-10	4.8104	2.4637e-09	6.2848
	512	1.4718e-11	4.4835	4.7629e-11	4.7416	1.0529e-11	4.9969	5.4809e-11	5.4903
IMEX-BDF5	128	4.9951e-08	–	3.49e-08	–	4.906e-09	–	1.7525e-08	–
	256	1.5618e-09	4.9992	1.0894e-09	5.0016	1.5768e-10	4.9594	4.2623e-10	5.3616
	512	5.6577e-11	4.7869	3.2187e-11	5.0809	4.6748e-12	5.076	1.1341e-11	5.232

Table 3. Table 3. L 1 superscript 𝐿 1 L^{1} error and estimated convergence rates for v 𝑣 v in the AP-explicit case.

		$\epsilon =1$		$\epsilon =0.1$		$\epsilon =0.01$		$\epsilon =0.001$
IMEX	$N$	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate
IMEX-SG(3,2)	128	4.1297e-05	–	0.0015868	–	0.30746	–	0.017621	–
	256	1.3046e-05	1.6624	0.00049005	1.6951	0.092824	1.7278	0.006956	1.341
	512	3.6092e-06	1.8538	0.00013416	1.869	0.025173	1.8826	0.0019759	1.8158
IMEX-BDF2	128	5.3468e-05	–	0.0010541	–	0.14541	–	0.015341	–
	256	1.7643e-05	1.5996	0.00034308	1.6194	0.045795	1.6669	0.0040821	1.91
	512	4.9529e-06	1.8328	9.5659e-05	1.8426	0.012564	1.8658	0.001048	1.9616
IMEX-TVB(3,3)	128	7.2288e-07	–	8.2789e-06	–	0.0083555	–	0.0044838	–
	256	1.1051e-07	2.7096	1.2742e-06	2.6999	0.0012579	2.7317	0.00064242	2.8031
	512	1.5103e-08	2.8712	1.7412e-07	2.8714	0.0001703	2.8849	8.423e-05	2.9311
IMEX-BDF3	128	6.3076e-07	–	8.0827e-06	–	0.010784	–	0.0087392	–
	256	9.8098e-08	2.6848	1.2676e-06	2.6727	0.0016207	2.7342	0.00099636	3.1328
	512	1.3503e-08	2.861	1.7436e-07	2.862	0.0002183	2.8923	0.00011713	3.0885
IMEX-TVB(4,4)	128	8.1836e-08	–	8.5455e-08	–	0.0019902	–	0.033128	–
	256	3.443e-09	4.571	7.8768e-09	3.4395	0.00015204	3.7104	0.00099626	5.0554
	512	1.5925e-10	4.4343	6.285e-10	3.6476	1.0128e-05	3.9079	4.1573e-05	4.5828
IMEX-BDF4	128	6.8621e-08	–	2.6004e-08	–	0.00076076	–	0.0045931	–
	256	2.4028e-09	4.8359	2.607e-09	3.3183	5.5621e-05	3.7737	0.00022733	4.3366
	512	5.3325e-11	5.4938	2.1977e-10	3.5683	3.7004e-06	3.9099	1.2649e-05	4.1676
IMEX-TVB(5,5)	128	1.7977e-08	–	2.1366e-07	–	9.0073e-06	–	0.030281	–
	256	4.6141e-10	5.2839	6.7708e-09	4.9798	3.4995e-07	4.6858	0.00039906	6.2457
	512	6.7507e-11	2.773	1.182e-10	5.84	1.1336e-08	4.9481	9.4581e-06	5.3989
IMEX-BDF5	128	1.2276e-07	–	1.533e-07	–	1.6471e-06	–	0.0022391	–
	256	3.8413e-09	4.9981	4.7927e-09	4.9993	5.7491e-08	4.8405	5.2309e-05	5.4197
	512	2.5154e-10	3.9327	1.9215e-10	4.6405	1.4584e-09	5.3009	1.3736e-06	5.251

Table 4. Table 4. L 1 superscript 𝐿 1 L^{1} error and estimated convergence rates for u 𝑢 u in the AP-implicit case.

		$\epsilon =1$		$\epsilon =0.1$		$\epsilon =0.01$		$\epsilon =0.001$
IMEX	$N$	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate
IMEX-SG(3,2)	128	0.0001009	–	8.4701e-05	–	0.00019279	–	0.015388	–
	256	3.1168e-05	1.6948	2.6145e-05	1.6958	5.8207e-05	1.7278	0.0047058	1.7093
	512	8.5342e-06	1.8687	7.1564e-06	1.8692	1.5785e-05	1.8826	0.0012828	1.8751
IMEX-BDF2	128	4.5983e-05	–	3.8769e-05	–	9.0615e-05	–	0.0074794	–
	256	1.4986e-05	1.6175	1.2626e-05	1.6185	2.8544e-05	1.6666	0.0023976	1.6413
	512	4.1817e-06	1.8414	3.5217e-06	1.8421	7.8321e-06	1.8657	0.00066348	1.8535
IMEX-TVB(3,3)	128	5.6622e-08	–	1.5426e-07	–	5.1657e-06	–	0.00024864	–
	256	9.4982e-09	2.5756	2.2969e-08	2.7476	7.7718e-07	2.7326	3.7993e-05	2.7103
	512	1.3103e-09	2.8577	3.121e-09	2.8796	1.0522e-07	2.8849	5.1683e-06	2.878
IMEX-BDF3	128	7.3382e-08	–	1.8956e-07	–	6.7012e-06	–	0.0002871	–
	256	1.2326e-08	2.5737	2.8971e-08	2.7099	1.0064e-06	2.7352	4.443e-05	2.6919
	512	1.7158e-09	2.8448	3.9775e-09	2.8647	1.3555e-07	2.8923	6.0575e-06	2.8747
IMEX-TVB(4,4)	128	2.9453e-08	–	2.5367e-08	–	1.2204e-06	–	2.9521e-05	–
	256	1.2448e-09	4.5645	7.248e-10	5.1292	9.3626e-08	3.7044	2.2475e-06	3.7153
	512	5.0075e-11	4.6356	1.6962e-11	5.4172	6.2417e-09	3.9069	1.5677e-07	3.8416
IMEX-BDF4	128	2.582e-08	–	7.3488e-09	–	4.6271e-07	–	1.2261e-05	–
	256	8.5367e-10	4.9187	2.0595e-10	5.1572	3.4155e-08	3.7599	8.3808e-07	3.8709
	512	1.6737e-11	5.6726	1.69e-12	6.9291	2.2839e-09	3.9025	5.5639e-08	3.9129
IMEX-TVB(5,5)	128	3.6171e-08	–	3.5169e-08	–	9.4689e-09	–	1.1944e-06	–
	256	1.1333e-09	4.9962	1.1385e-09	4.9491	3.3274e-10	4.8308	7.5898e-09	7.298
	512	5.1332e-11	4.4646	6.4841e-11	4.1341	1.2796e-11	4.7006	3.7672e-10	4.3325
IMEX-BDF5	128	2.0507e-08	–	2.7439e-08	–	4.9118e-09	–	1.3221e-06	–
	256	6.2502e-10	5.0361	8.5434e-10	5.0053	1.5808e-10	4.9575	2.1912e-08	5.915
	512	2.5331e-11	4.625	3.7038e-11	4.5277	1.5295e-12	6.6914	1.0639e-09	4.3643

Table 5. Table 5. L 1 superscript 𝐿 1 L^{1} error and estimated convergence rates for v 𝑣 v in the AP-implicit case.

		$\epsilon =1$		$\epsilon =0.1$		$\epsilon =0.01$		$\epsilon =0.001$
IMEX	$N$	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate	${\Vert E_{\Delta x,\Delta t}^k\Vert }_1$	Rate
IMEX-SG(3,2)	128	4.4648e-05	–	0.00089276	–	0.30567	–	0.11134	–
	256	1.3846e-05	1.6891	0.00027568	1.6953	0.092283	1.7278	0.032643	1.7702
	512	3.7984e-06	1.866	7.5475e-05	1.8689	0.025026	1.8826	0.0087961	1.8918
IMEX-BDF2	128	2.3413e-05	–	0.00042114	–	0.14368	–	0.060252	–
	256	7.6904e-06	1.6062	0.00013723	1.6177	0.045257	1.6667	0.0174	1.7919
	512	2.1539e-06	1.8361	3.8288e-05	1.8416	0.012418	1.8657	0.0046647	1.8992
IMEX-TVB(3,3)	128	2.0289e-07	–	1.2756e-06	–	0.0081773	–	0.001989	–
	256	3.279e-08	2.6294	1.9536e-07	2.707	0.0012317	2.731	0.00026696	2.8974
	512	4.5301e-09	2.8557	2.6669e-08	2.8729	0.00016679	2.8845	3.4663e-05	2.9451
IMEX-BDF3	128	2.3141e-07	–	1.3649e-06	–	0.010612	–	0.0023402	–
	256	3.8064e-08	2.6039	2.1348e-07	2.6766	0.0015951	2.7339	0.00032067	2.8674
	512	5.2868e-09	2.848	2.9369e-08	2.8617	0.00021487	2.8921	4.2074e-05	2.9301
IMEX-TVB(4,4)	128	7.1151e-08	–	9.5652e-08	–	0.0019469	–	0.00028374	–
	256	2.6963e-09	4.7218	2.8851e-09	5.0511	0.0001488	3.7098	1.6152e-05	4.1348
	512	1.0634e-10	4.6642	3.3063e-10	3.1253	9.9148e-06	3.9076	1.0431e-06	3.9527
IMEX-BDF4	128	6.5438e-08	–	1.7232e-08	–	0.00074557	–	9.082e-05	–
	256	2.1846e-09	4.9047	4.6192e-10	5.2213	5.4526e-05	3.7733	5.8265e-06	3.9623
	512	3.0347e-11	6.1697	3.0551e-11	3.9183	3.6281e-06	3.9097	3.7459e-07	3.9592
IMEX-TVB(5,5)	128	8.0907e-08	–	1.4931e-07	–	8.7961e-06	–	3.1526e-05	–
	256	2.3712e-09	5.0926	4.7065e-09	4.9875	3.4175e-07	4.6859	2.7811e-07	6.8247
	512	1.4198e-10	4.0619	5.5157e-11	6.415	1.1008e-08	4.9563	9.1116e-09	4.9318
IMEX-BDF5	128	4.6717e-08	–	1.0372e-07	–	1.6077e-06	–	9.3079e-06	–
	256	1.4517e-09	5.0081	3.1885e-09	5.0237	5.6036e-08	4.8425	1.6743e-07	5.7968
	512	8.9081e-12	7.3484	1.0165e-10	4.9711	1.4002e-09	5.3226	7.4625e-09	4.4878

Table 6. Table 6. Examples of IMEX-LM methods. The numbers between brackets denote respectively: s 𝑠 s number of steps, p 𝑝 p order of the method. For BDF methods the number of steps is equivalent to the order of the method.

IMEX	Explicit		Implicit
SG(3,2)	$a^T$	$(-\frac{3}{4},0,-\frac{1}{4},0,0)$	$c^T$	$(0,0,\frac{1}{2},0,0)$
	$b^T$	$(\frac{3}{2},0,0,0,0)$	$c_{-1}$	$1$
BDF2	$a^T$	$(-\frac{4}{3},\frac{1}{3},0,0,0)$	$c^T$	$(0,0,0,0,0)$
	$b^T$	$(\frac{4}{3},-\frac{2}{3},0,0,0)$	$c_{-1}$	$\frac{2}{3}$
TVB(3,3)	$a^T$	$(-\frac{3909}{2048},\frac{1367}{1024},-\frac{873}{2048},0,0)$	$c^T$	$(-\frac{1139}{12288},-\frac{367}{6144},\frac{1699}{12288},0,0)$
	$b^T$	$(\frac{18463}{12288},-\frac{1271}{768},\frac{8233}{12288},0,0)$	$c_{-1}$	$1089/2048$
BDF3	$a^T$	$(-\frac{18}{11},\frac{9}{11},-\frac{2}{11},0,0)$	$c^T$	$(0,0,0,0,0)$
	$b^T$	$(\frac{18}{11},-\frac{18}{11},\frac{6}{11},0,0)$	$c_{-1}$	$\frac{6}{11}$
TVB(4,4)	$a^T$	$(-\frac{21531}{8192},\frac{22753}{8192},-\frac{12245}{8192},\frac{2831}{8192},0)$	$c^T$	$(-\frac{3567}{8192},-\frac{697}{24576},\frac{4315}{24576},-\frac{41}{384},0)$
	$b^T$	$(\frac{13261}{8192},-\frac{75029}{24576},\frac{54799}{24576},-\frac{15245}{24576},0)$	$c_{-1}$	$\frac{4207}{8192}$
BDF4	$a^T$	$(-\frac{48}{25},\frac{36}{25},-\frac{16}{25},\frac{3}{25},0)$	$c^T$	$(0,0,0,0,0)$
	$b^T$	$(\frac{48}{25},-\frac{72}{25},\frac{48}{25},-\frac{12}{25},0)$	$c_{-1}$	$\frac{12}{25}$
TVB(5,5)	$a^T$	$(-\frac{13553}{4096},\frac{38121}{8192},-\frac{7315}{2048},\frac{6161}{4096},-\frac{2269}{8192})$	$c^T$	$(-\frac{4118249}{5898240},\frac{768703}{2949120},\frac{47849}{245760},-\frac{725087}{2949120},\frac{502321}{5898240})$
	$b^T$	$(\frac{10306951}{5898240},-\frac{13656497}{2949120},\frac{1249949}{245760},-\frac{7937687}{2949120},\frac{3387361}{5898240})$	$c_{-1}$	$\frac{4007}{8192}$
BDF5	$a^T$	$(-\frac{300}{137},\frac{300}{137},-\frac{200}{137},\frac{75}{137},-\frac{12}{137})$	$c^T$	$(0,0,0,0,0)$
	$b^T$	$(\frac{300}{137},-\frac{600}{137},\frac{600}{137},-\frac{300}{137},\frac{60}{137})$	$c_{-1}$	$\frac{60}{137}$

Equations

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}u+\partial_{x}v=0,\\ \displaystyle\partial_{t}v+\frac{1}{{\varepsilon}^{2\alpha}}\partial_{x}p(u)=-\frac{1}{{\varepsilon}^{1+\alpha}}\left(v-f(u)\right),\quad\alpha\in[0,1]\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}u+\partial_{x}v=0,\\ \displaystyle\partial_{t}v+\partial_{x}p(u)=-\frac{1}{{\varepsilon}}\left(v-f(u)\right),\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}u+\partial_{x}v=0,\\ \displaystyle\partial_{t}v+\frac{1}{{\varepsilon}^{2}}\partial_{x}p(u)=-\frac{1}{{\varepsilon}^{2}}\left(v-f(u)\right).\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle v=f(u)-{\varepsilon}^{1-\alpha}\partial_{x}p(u)+{\varepsilon}^{1+\alpha}f^{\prime}(u)^{2}\partial_{x}u+{\mathcal{O}}({\varepsilon}^{2}),\\ \displaystyle\partial_{t}u+\partial_{x}f(u)={\varepsilon}^{1+\alpha}\partial_{x}\Bigg{[}\left(\frac{p^{\prime}(u)}{{\varepsilon}^{2\alpha}}-f^{\prime}(u)^{2}\right)\partial_{x}u\Bigg{]}+{\mathcal{O}}({\varepsilon}^{2}).\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle v=f(u),\\ \displaystyle\partial_{t}u+\partial_{x}f(u)=0,\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle v=f(u)-\partial_{x}p(u),\\ \displaystyle\partial_{t}u+\partial_{x}f(u)=\partial_{xx}p(u).\\ \end{array}\right.$$

$$f^{\prime }(u{)}^2<\frac{p^{\prime }(u)}{{\epsilon }^{2\alpha }},$$

$$\frac{u^{n+1}-u^n}{\Delta t}$$

$${\epsilon }^{1+\alpha }\frac{v^{n+1}-v^n}{\Delta t}$$

$$v^{n+1}=\frac{{\epsilon }^{1+\alpha }}{{\epsilon }^{1+\alpha }+\Delta t}{v}^n-\frac{\Delta t}{{\epsilon }^{1+\alpha }+\Delta t}\left({\epsilon }^{1-\alpha }{\partial }_{x}p(u^n)-f(u^n)\right).$$

$$\frac{u^{n+1}-u^n}{\Delta t}+\frac{{\epsilon }^{1+\alpha }}{{\epsilon }^{1+\alpha }+\Delta t}{v}_x^n+\frac{\Delta t}{{\epsilon }^{1+\alpha }+\Delta t}{\partial }_{x}f(u^n)$$

$$\frac{v^{n+1}-v^n}{\Delta t}+\frac{{\epsilon }^{1-\alpha }}{{\epsilon }^{1+\alpha }+\Delta t}{\partial }_{x}p(u^n)$$

$$\frac{u^{n+1}-u^n}{\Delta t}+{\partial }_{x}f(u^n)$$

$${\partial }_{t}u+\frac{{\epsilon }^{1+\alpha }}{{\epsilon }^{1+\alpha }+\Delta t}{\partial }_{x}v+\frac{\Delta t}{{\epsilon }^{1+\alpha }+\Delta t}{\partial }_{x}f(u)$$

$${\partial }_{t}v+\frac{{\epsilon }^{1-\alpha }}{{\epsilon }^{1+\alpha }+\Delta t}{\partial }_{x}p(u)$$

$$\displaystyle\frac{u^{n+1}-u^{n}}{\Delta t}=\partial_{t}u\big{|}_{t=t^{n}}+\mathcal{O}(\Delta t),\quad\frac{v^{n+1}-v^{n}}{\Delta t}=\partial_{t}v\big{|}_{t=t^{n}}+\mathcal{O}(\Delta t).$$

$${\lambda }_{\pm }^{\alpha }(\Delta t,\epsilon )=\frac{1}{2}\left(\gamma (1-{\theta }_{\alpha })\pm \sqrt{{\gamma }^2(1-{\theta }_{\alpha }{)}^2+4{\epsilon }^{-2\alpha }{\theta }_{\alpha }^2}\right),$$

$${\theta }_{\alpha }(\Delta t,\epsilon ):=\frac{{\epsilon }^{1+\alpha }}{{\epsilon }^{1+\alpha }+\Delta t},$$

$${\partial }_{x}f(u)=f^{\prime }(u){\partial }_{x}u=\gamma {\partial }_{x}u,{\partial }_{x}p(u)=p^{\prime }(u){\partial }_{x}u={\partial }_{x}u.$$

$${\lambda }_{\pm }^{\alpha }(0,\epsilon )=\pm \frac{1}{{\epsilon }^{\alpha }},$$

$${\lambda }_{\pm }^{\alpha }(\Delta t,0)=\frac{1}{2}\left(\gamma \pm |\gamma |\right).$$

$$u^{n+1}$$

$$v^{n+1}$$

$$Explicit$$

$$b^T=(b_0,b_1,\dots ,b_{s-1})$$

$$Implicit$$

$$u^{n+1}$$

$$v^{n+1}$$

$$=-a^T\cdot V+\left(a^T-\frac{{\epsilon }^2{a}^T}{{\epsilon }^2+\Delta t{c}_{-1}}-\frac{\Delta t{c}^T}{{\epsilon }^2+\Delta t{c}_{-1}}\right)\cdot V-\frac{\Delta t{b}^T}{{\epsilon }^2+\Delta t{c}_{-1}}\cdot \left({\partial }_{x}p(U)-f(U)\right),$$

$$v^{n+1}=-a^T\cdot V-\frac{\Delta t}{{\epsilon }^2+\Delta t{c}_{-1}}\left(c^T-c_{-1}{a}^T\right)\cdot V-\frac{\Delta t}{{\epsilon }^2+\Delta t{c}_{-1}}{b}^T\cdot \left({\partial }_{x}p(U)-f(U)\right).$$

$$u^{n+1}$$

$$+\frac{\Delta t^2{c}_{-1}}{{\epsilon }^2+\Delta t{c}_{-1}}\left(c^T-c_{-1}{a}^T\right)\cdot {\partial }_{x}V+\frac{\Delta t^2{c}_{-1}}{{\epsilon }^2+\Delta t{c}_{-1}}{b}^T\cdot {\partial }_x\left({\partial }_{x}p(U)-f(U)\right)$$

Full text

The work of G. Albi has been partially supported by the project Computer Science for Industry 4.0, ’MIUR Departments of Excellence 2018-2022’. The authors acknowledge the partial support of GNCS-INdAM 2019 project Numerical approximation of hyperbolic problems and applications.

Computer Science Department, University of Verona, (giacomo.albi@univr.it).

Mathematics Department, University of Ferrara, Italy (giacomo.dimarco@unife.it).

Mathematics Department, University of Ferrara, Italy (lorenzo.pareschi@unife.it).

Key words. Implicit-explicit methods, linear multistep methods, hyperbolic balance laws, fluid-dynamic limit, diffusion limit, asymptotic-preserving schemes.

AMS subject classification. 35L65, 65L06, 65M06, 76D05, 82C40.

1. 1 Introduction

2. 2 First order IMEX discretization

3. 3 AP-explicit and AP-implicit IMEX-LM methods

4. 3.1 AP-explicit methods in the diffusive case: $\alpha=1$

5. 3.1.1 Asymptotic preserving property for $\alpha=1$.

6. 3.2 AP-explicit methods in the general case: $\alpha\in[0,1)$.

7. 3.2.1 Asymptotic preserving property for $\alpha\in[0,1)$.

8. 3.3 Removing the parabolic stiffness: AP-implicit methods

9. 3.3.1 Asymptotic preserving property for the AP-implicit methods.

10. 4 Linear stability analysis

11. 4.1 AP-explicit methods

12. 4.2 AP-implicit methods

13. 5 Space discretization

14. 6 Numerical validation & applications

15. 6.1 Test 1. Numerical convergence study for a linear problem.

16. 6.2 Test 2: Riemann problem for the linear model.

17. 6.3 Test 3: Barenblatt solution for the porous media equation.

18. 6.4 Test 4: Applications to the Ruijgrook–Wu model.

19. 7 Conclusions

20. A Order conditions for IMEX-LM methods and examples

1. Introduction

The goal of the present work is to develop high order numerical methods based on IMEX linear multistep (IMEX-LM) schemes for hyperbolic systems with relaxation [14, 38, 42]. These systems often contain multiple space-time scales which may differ by several orders of magnitude. In fact, the various parameters characterizing the models permit to describe different physical situations, like flows which pass from compressible to incompressible regimes or flows which range from rarefied to dense states. This is the case, for example, of kinetic equations close to the hydrodynamic limits [5, 15, 13, 49]. In such regimes these systems can be more conveniently described in terms of macroscopic equations since these reduced systems permit to describe all the features related to the space-time scale under consideration [5, 49]. However, such macroscopic models can not handle all the possible regimes one is frequently interested in. For such reason one has to resort to the full kinetic models. They permit to characterize a richer physics but on the other hand they are computationally more expensive and limited by the stiffness induced by the scaling under consideration [16].

The prototype system we will use in the rest of the paper is the following [8, 41]

[TABLE]

where $\varepsilon$ is the scaling factor and $\alpha$ characterizes the different type of asymptotic limit that can be obtained. The condition $p^{\prime}(u)>0$ should be satisfied for hyperbolicity to hold true since the eigenvalues of (1.1) are given by $\pm\sqrt{p^{\prime}(u)}/{\varepsilon}^{\alpha}$. Note that, except for the case $\alpha=0$, the eigenvalues are unbounded for small values of ${\varepsilon}$.

System (1.1) is obtained from a classical $(2\times 2)$ $p$-system with relaxation under the space-time scaling $t\to t/{\varepsilon}^{1+\alpha}$, $x\to x/{\varepsilon}$ and by the change of variables $v=\tilde{v}/\varepsilon^{\alpha}$, $f(u)=\tilde{f}(u)/\varepsilon^{\alpha}$, where $\tilde{v}$ is the original unknown and $\tilde{f}(u)$ the original flux associated to the variable $u$ in the non rescaled $p$-system. For $\alpha=0$, system (1.1) reduces to the usual hyperbolic scaling

[TABLE]

whereas for $\alpha=1$ yields the so-called diffusive scaling

[TABLE]

More in general, thanks to the Chapman-Enskog expansion [15], for small values of ${\varepsilon}$ we get from (1.1) the following nonlinear convection-diffusion equation

[TABLE]

In the limit $\varepsilon\rightarrow 0$, for $\alpha\in[0,1)$, we are led to the conservation law

[TABLE]

while, when $\alpha=1$, in the asymptotic limit we obtain the following advection-diffusion equation

[TABLE]

Note that, the main stability condition for system (1.4) corresponds to

[TABLE]

and it is always satisfied in the limit ${\varepsilon}\rightarrow 0$ when $\alpha>0$, whereas for $\alpha=0$ the function $p(u)$ and $f(u)$ must satisfy the classical sub-characteristic condition [14, 38].

The space-time scaling just discussed, in classical kinetic theory, is related to the hydrodynamical limits of the Boltzmann equation. In particular, for $\alpha=0$ it corresponds to the compressible Euler scaling, whereas for $\alpha\in(0,1)$ to the incompressible Euler limit. In the case $\alpha=1$ dissipative effects become non-negligible and we get the incompressible Navier-Stokes scaling. We refer to [13, Chapter 11] and [49] for further details and the mathematical theory behind the hydrodynamical limits of the Boltzmann equation. Moreover, we refer to [36, 39] for theoretical results on the diffusion limit of a system like (1.3).

The development of numerical methods to solve hyperbolic systems with stiff source terms has attracted many researches in the recent past [11, 24, 24, 22, 44, 43, 10, 28, 31, 33, 34, 17, 37]. The main computational challenge is related to the presence of the different scales that require a special care to avoid loss of stability and spurious numerical solutions. In particular, in diffusive regimes the schemes should be capable to deal with the very large characteristic speeds of the system by avoiding a CFL condition of the type $\Delta t={\mathcal{O}}(\varepsilon^{\alpha})$. A particular successful class of schemes is represented by the so-called asymptotic-preserving (AP) schemes which aims at preserving the correct asymptotic behavior of the system without any loss of efficiency due to time step restrictions related to the small scales [18, 30, 32, 20].

In the large majority of these works, the authors focused on the specific case $\alpha=0$, where a hyperbolic to hyperbolic scaling is studied, or to the case $\alpha=1$, where a hyperbolic to parabolic scaling is analysed. Very few papers have addressed the challenging multiscale general problem for the various possible values of $\alpha\in[0,1]$ (see [40, 32, 8, 22]) and all of them refer to one-step IMEX methods in a Runge-Kutta setting. We refer to [1, 2, 4, 19, 21, 26, 45, 47] for various IMEX-LM methods developed in the literature and we mention that comparison between IMEX Runge-Kutta methods and IMEX-LM methods have been presented in [26].

In the present work, following the approach recently introduced by Boscarino, Russo and Pareschi in [8], we analyze the construction of IMEX-LM for such problems that work uniformly independently of the choices of $\alpha$ and ${\varepsilon}$. By this, we mean that the schemes are designed in such a way as to be stable for all different ranges of the scaling parameters independently of the time step. At the same time, they should ensure high order in space and time and should be able to accurately describe the various asymptotic limits. Moreover, whenever possible, the above described properties must be achieved without the need of an iterative solver for non linear equations.

In [8], using a convenient partitioning of the original problem, the authors developed IMEX Runge-Kutta schemes for a system like (1.1) which work uniformly with respect to the scaling parameters. Here, we extend these results to IMEX-LM, previous results for IMEX-LM methods refer to the case $\alpha=0$ (see [19, 26]). Among others, there are two main reasons to consider the development of such schemes. First, in contrast to the IMEX Runge-Kutta case, for IMEX-LM it is relatively easy to construct schemes up to fifth order in time and typically they show a more uniform behavior of the error with respect to the scaling parameters. Second, thanks to the use of BDF methods, it is possible consider only one evaluation of the source term per time step independently of the scheme order. The latter feature is particularly significant in term of computational efficiency for kinetic equations where often the source term represents the most expensive part of the computation [18, 19].

The rest of the paper is organized as follows. In Section 2, we discuss the discretization of these multiscale problems and motivate our partitioning choice of the system by analyzing a simple first order IMEX scheme. Next, in Section 3, we introduce the general IMEX-LM methods and discuss the asymptotic-preserving properties of our approach. Two classes of schemes are considered, AP-explicit and AP-implicit accordingly to the way the diffusion term in the limit equation is treated. In Section 4, we perform a linear stability analysis for $2\times 2$ linear systems in the case of IMEX-LM methods based on a backward-differentiation formula (BDF). In Section 5, the space discretization is briefly discussed. Several numerical examples are reported in Section 6 which confirm the theoretical findings. Final considerations and future developments are discussed at the end of the article.

2. First order IMEX discretization

In this part, we discuss a first order IMEX time-discretization of the relaxation system (1.1) and we analyze its relationship with a reformulated system in which the eigenvalues are bounded for any value of the scaling parameter $\varepsilon$. To this aim, following [8], we consider the following implicit-explicit first order partitioning of system (1.1)

[TABLE]

One can notice that in system (2.1) besides its implicit form, the second equation can be solved explicitly by inversion of the linear term $v^{n+1}$. This gives

[TABLE]

Then, making use of the above relation and inserting it in the first equation, one gets

[TABLE]

while a simple rewriting of the second equation gives

[TABLE]

Therefore, the IMEX scheme can be recast in an equivalent fully explicit form. Similarly to the continuous case, depending on the choice of $\alpha$, as ${\varepsilon}\to 0$, we have different limit behaviors. For $\alpha\in[0,1)$ we obtain

[TABLE]

whereas in the case $\alpha=1$ we get

[TABLE]

For small values of $\Delta t$, the scheme (2.3)-(2.4) corresponds up to first order in time to the system

[TABLE]

where the following Taylor expansion has been employed at $t=t^{n}=n\Delta t$

[TABLE]

The main feature of system (2.7) is that its left-hand side has bounded characteristic speeds. These are given by

[TABLE]

with

[TABLE]

and where, for simplicity, we considered $f^{\prime}(u)=\gamma$, $\gamma\in\mathbb{R}$ and $p^{\prime}(u)=1$ so that

[TABLE]

If we fix ${\varepsilon}$ and send $\Delta t\to 0$ we obtain the usual characteristic speeds of the original hyperbolic system, i.e.

[TABLE]

while for a fixed $\Delta t$, the characteristic speeds $\lambda^{\alpha}_{+}$ and $\lambda^{\alpha}_{-}$ are respectively decreasing and increasing functions of $\varepsilon$ and, as ${\varepsilon}\to 0$, they converge to

[TABLE]

In Figure 1, we show the shape of the eigenvalues (2.8) for $\gamma=1$ and different values of the scaling parameter $\alpha$ and the time step $\Delta t$. We observe that when $\varepsilon$ grows the absolute value of the eigenvalues diminish accordingly and when $\varepsilon$ diminishes the eigenvalues grow but remain bounded by the finite time step.

Thus, for a given $\Delta t$, if we denote by $\Delta x$ the space discretization parameter, from the left hand side of (2.7) we expect the hyperbolic CFL condition $\Delta t\leq\Delta x/|\gamma|$ in the limit ${\varepsilon}\to 0$. On the other hand, the stability restriction coming from the parabolic term requires $\Delta t={\mathcal{O}}(\Delta x^{2})$ when $\alpha=1$.

In the next Section, we will show how to generalize the above arguments to the case of high order IMEX multistep methods.

3. AP-explicit and AP-implicit IMEX-LM methods

In this part, we focus our attention on order $s$, $s$-step IMEX-LM methods with $s\geq 2$ (see Appendix A.1 for derivation and order conditions). First we discuss methods that yield a fully explicit discretization in the limit $\varepsilon\to 0$. For clarity of presentation, we separate the discussion of the diffusive case $\alpha=1$ from the general case $\alpha\in[0,1)$. In the second part, we discuss IMEX-LM discretizations which deal with the stiffness caused by the parabolic term in the asymptotic limit.

3.1. AP-explicit methods in the diffusive case: $\alpha=1$

In this case, we can write the $s-$step IMEX-LM for the original hyperbolic system (1.1) as follows

[TABLE]

where we introduced the following coefficients

[TABLE]

Methods for which $c_{j}=0$, $j=0,\ldots,s-1$ are referred to as implicit-explicit backward differentiation formula, IMEX-BDF in short. Another important class of LM is represented by implicit-explicit Adams methods, for which $a_{0}=-1$, $a_{j}=0$, $j=1,\ldots,s-1$. We refer to Appendix A for a brief survey of some IMEX multistep methods, and to [1, 2, 4, 19, 26, 21, 45, 47] for further details and additional schemes.

In what follows, we rely on the equivalent vector-matrix notation,

[TABLE]

where $U=(u^{n},\ldots,u^{n-s+1})^{T}$, $V=(v^{n},\ldots,v^{n-s+1})^{T}$, $\partial_{x}p(U)=(\partial_{x}p(u^{n}),\ldots,\partial_{x}p(u^{n-s+1}))^{T}$ and $f(U)=(f(u^{n}),\ldots,f(u^{n-s+1}))^{T}$ are $s$-dimensional vectors.

Similarly to the one step scheme (2.1), we proceed by rewriting the multistep methods in the fully explicit vector form. To that aim, we observe that the second equation in (3.3) can be explicitly solved in terms of $v$ due to the linearity in the relaxation part. This gives

[TABLE]

and thus

[TABLE]

Substituting equation (3.5) into the first one of (3.3) leads to

[TABLE]

Hence, we obtain the following system

[TABLE]

which is the generalization of system (2.4) to an $s$-step IMEX scheme.

Our aim now is to show that, similarly to the simple first order method analyzed in Section 2, the above discretization for a small value of $\Delta t$ corresponds, up to first order in time, to a modified hyperbolic problem where the characteristic speeds are bounded even in the limit $\varepsilon\to 0$. More precisely, see also [3], for a smooth solution in time, by Taylor series expansion about $t=t^{n}$, we have

[TABLE]

where $e$ is a vector of ones in $\mathbb{R}^{s}$, $J=(0,\ldots,s-1)^{T}$, $\partial^{q}_{t}$, $q=1,\ldots,s$ denotes the $q$-derivative and the vector powers must be understood component-wise. Similarly, we can expand $u^{n+1}$ and $v^{n+1}$ about $t=t^{n}$.

Therefore, we obtain

[TABLE]

From the order conditions we have

[TABLE]

and subsequentely

[TABLE]

Thus, scheme (3.7) for small values of $\Delta t$ can be considered as a first order approximation of the following modified system

[TABLE]

where the factor $\beta$ simplifies in all the terms.

Note that, the above system has exactly the same structure as system (2.7) (with $c_{-1}=1$ and $\alpha=1$) which was derived from to the first order time discretization. As a consequence, under the same simplification assumptions (2.9) on the fluxes $f(u)$ and $p(u)$, the eigenvalues of the hyperbolic part correspond to

[TABLE]

where

[TABLE]

Thus, the bounds for the characteristic velocities are the same as for the first order scheme and we get the limit cases

[TABLE]

3.1.1. Asymptotic preserving property for $\alpha=1$.

Now, we study the capability of the schemes (3.3) to become a consistent discretization of the limit system (1.6). To that aim, letting $\varepsilon\rightarrow 0$, in the reformulated scheme (3.7), we get from the first equation

[TABLE]

which corresponds to the explicit multistep scheme applied to the limiting convection-diffusion equation (1.6). For this reason, from now on, we refer to this class of IMEX-LM schemes as AP-explicit methods. On the other hand, we have for the second equation

[TABLE]

or equivalently

[TABLE]

Let us observe that in order to have at time $t^{n+1}$ a consistent projection over the asymptotic limit a condition over the states $U=(u^{n},..,u^{n-s+1})^{T}$ and $V=(v^{n},..,v^{n-s+1})^{T}$ should be imposed. This can be resumed by saying that the vector of the initial data should be well prepared to the asymptotic state. This means that for the first variable

[TABLE]

where $\bar{u}^{n-j}$ is a consistent solution of the limit system (1.6) while $\tilde{u}^{n-j}_{\varepsilon}$ is a perturbation that disappears in the limit. An analogous relation should hold true for the second variable $v$, i.e.

[TABLE]

where $\bar{v}^{n-j}=f(u^{n-j})-\partial_{x}p(u^{n-j})$, $j=0,\ldots,s-1$ is a consistent projection of the asymptotic limit, while $\tilde{v}^{n-j}_{\varepsilon}$ is a perturbation that disappears in the limit. Under this assumption, as a consequence of the order conditions, relation (3.13) is an $s$-order approximation of the asymptotic limit $v=f(u)-\partial_{x}p(u)$ and therefore, at subsequent time steps, the numerical solution is guaranteed to satisfy (3.14)-(3.15). If such conditions are not imposed on the initial values, then the numerical solution may present a spurious initial layer and deterioration of accuracy is observed. In particular, for IMEX-BDF methods, expression (3.13) simplifies to

[TABLE]

This shows that, even for non well prepared initial data in $v$ but only in $u$ we obtain an $s$-order approximation of the equilibrium state and subsequently of the numerical solution for all times. This stronger AP property of IMEX-BDF methods is satisfied also in the asymptotic limits analyzed for the various schemes in the sequel of the manuscript. We refer to [19, 18, 43] for more detailed discussions.

3.2. AP-explicit methods in the general case: $\alpha\in[0,1)$.

The IMEX-LM scheme in vector form reads

[TABLE]

Again, we rewrite the second equation by solving it in terms of $v^{n+1}$ as follows

[TABLE]

and using this solution in the first equation, we obtain the explicit scheme

[TABLE]

Considering, the same Taylor approximations (LABEL:eq:IMEXest2) the above schemes corresponds up to first order in time to the modified system

[TABLE]

Clearly, system (3.19) has again the same structure as (2.7). As a consequence, under the same simplification assumptions (2.9), the eigenvalues of the hyperbolic part correspond to

[TABLE]

where $\theta_{\alpha}$ is defined as follows

[TABLE]

and the bounds for the characteristic velocities for ${\varepsilon}=0$ and $\Delta t=0$ are

[TABLE]

3.2.1. Asymptotic preserving property for $\alpha\in[0,1)$.

We consider now the analogous asymptotic preserving property proved for the schemes (3.3) in the case of the schemes (3.17). Namely, we want to show that (3.17) becomes a consistent discretization of the limit system (1.5) when $\varepsilon\rightarrow 0$. Taking scheme (3.18), we get from the first equation

[TABLE]

which is a standard explicit multistep discretization of the asymptotic hyperbolic limit, i.e. of equation (1.5). On the other hand, we have for the second equation

[TABLE]

or equivalentely

[TABLE]

As a consequence of the order conditions, equation (3.23) defines an $s$-order consistent approximation of the asymptotic limit $v=f(u)$ provided that the vector of the initial data is well prepared. These conditions are the analogous of (3.14) and (3.15) except that now $\bar{v}^{n-j}=f(u^{n-j})$, $j=0,\ldots,s-1$.

For IMEX-BDF methods (3.23) reduces to

[TABLE]

so that, even for non well prepared initial data in $v$ but only in $u$ we obtain an $s$-order approximation of the equilibrium state and subsequently of the numerical solution for all times.

3.3. Removing the parabolic stiffness: AP-implicit methods

Although, the schemes developed in the previous Section overcome the stiffness related to the scaling factor $\varepsilon$, there is another stiffness that may appear in the equations close to the asymptotic limit. In fact, as shown previously, all the schemes originate a fully explicit scheme in the limit.

In diffusive regimes, this typically leads to the time step restriction $\Delta t={\mathcal{O}}(\Delta x^{2})$ when

[TABLE]

(see the diffusion coefficient in equation (3.19)). Therefore, for small $\varepsilon$ and in the case of $\alpha\simeq 1$, the main stability restriction is due to the second order term of the Chapmann–Enskog expansion (1.4). Note that, beside the case $\alpha=1$ and ${\varepsilon}\to 0$ where we obtain a parabolic problem in the limit, the above time step limitation may occur also in transient regimes for $\alpha\neq 1$ as soon as (3.25) holds true.

For this reason, we modify the partitioning of the system taking also $\partial_{x}p(u)$ implicit in the second equation as follows

[TABLE]

We can still solve the second equation in $v$, to get

[TABLE]

which, inserted into the first equation of (LABEL:eq:SPs_vec1) yields the IMEX formulation

[TABLE]

Note that, except for the case in which $p(u)$ is linear, in general, the first equation in (3.27) requires the adoption of a suitable solver for nonlinear problems to compute $u^{n+1}$. By the same arguments of in the previous Sections, for small values of $\Delta t$, the scheme (3.27) corresponds up to first order to the modified system (3.19). Thus, under the same simplification assumptions (2.9), the eigenvalues of the hyperbolic part are given by (3.20).

3.3.1. Asymptotic preserving property for the AP-implicit methods.

Finally, we conclude our analysis by studying the asymptotic preserving property. As we will see the main difference is that in the asymptotic limit the diffusive terms are integrated implicitly. For this reason we refer to this class of IMEX-LM schemes as AP-implicit methods. We first consider the case $\alpha=1$. Taking the reformulated IMEX scheme (3.27) and letting $\varepsilon\rightarrow 0$ with $\alpha=1$, gives

[TABLE]

which correspond to the IMEX multistep scheme applied to the limiting convection diffusion problem where the diffusion term is treated implicitly [1, 2]. Note that, for the second equation we have

[TABLE]

or equivalently

[TABLE]

which, under the assumption of well prepared initial values, as a consequence of the order conditions, corresponds to a $s$-order approximation of the equilibrium projection $v=f(u)-\partial_{x}p(u)$.

On the other hand, in the case $\alpha\in[0,1)$, we get the same asymptotic limit (3.22) of the AP-explicit method (see Section 3.2.1)

[TABLE]

We will not discuss further this limit system, but we emphasize that (3.30) is obtained as the limit of the implicit-explicit scheme (3.27) whereas (3.22) is obtained as the limit of the explicit scheme (3.18).

4. Linear stability analysis

Monotonicity properties for IMEX-LM have been previously studied in [4, 19, 21, 26, 27]. Due to the well-known difficulties in extending the usual stability analysis for linear systems to the implicit-explicit setting most results are limited to the single scalar equation. In our case, however, the schemes are specifically designed to deal with systems in the form (1.1), and we need therefore to tackle the stability properties in such a case. Here we show that, in the case of IMEX-BDF methods, we can generalize some of the stability results for the single scalar equations to linear multiscale systems of the form

[TABLE]

where $\varepsilon>0$, $\gamma>0$ and $\alpha\in[0,1]$. Note that, for small values of $\varepsilon$, the above system when $\alpha=1$ reduces to the convection-diffusion equation $\partial_{t}u+\gamma\partial_{x}u=\partial_{xx}u$, whereas when $\alpha=0$, if $\gamma<1$, yields the simple advection equation $\partial_{t}u+\gamma\partial_{x}u=0$. To simplify notations, in the sequel we will assume $\gamma=1$ and $\alpha>0$. The case $\alpha=0$ is rather classical and follows straightforwardly from our analysis. Under these assumptions, the Chapman-Enskog expansion for small values of $\varepsilon$ gives the limiting convection-diffusion equation

[TABLE]

In Fourier variables we get

[TABLE]

where, for example, $\xi=\sin(2k)/\Delta x$ if we use central differences and $k$ is the frequency of the corresponding Fourier mode.

The change of variables $y={\hat{u}}$, $z=\varepsilon^{\alpha}{\hat{v}}$, $\lambda_{I}=i\xi/\varepsilon^{\alpha}$, $\lambda_{R}=1/\varepsilon^{1+\alpha}$, $\lambda=\lambda_{I}+\lambda_{R}\in\mathbb{C}$ transforms the system into the problem

[TABLE]

Let us note that the above problem is equivalent to the second order differential equation

[TABLE]

4.1. AP-explicit methods

We then apply an IMEX-LM method to system (4.4) as follows

[TABLE]

In the case of IMEX-BDF methods the first equation (4.6) permits to write

[TABLE]

and more in general for $j=0,\ldots,s-1$

[TABLE]

Thus, by direct substitution into the second equation (4.7) we obtain a discretization to (4.5) in the form

[TABLE]

Finally, we can rewrite the resulting scheme as

[TABLE]

Note that, the above IMEX-LM in the limit $\varepsilon\to 0$ for $\alpha=1$ leads to the reduced scheme

[TABLE]

which corresponds to an explicit LMM for the convection-diffusion equation.

The characteristic equation for scheme (LABEL:eq:scstab) reads

[TABLE]

with $z_{R}=\lambda_{R}\Delta t$, $z_{I}=\lambda_{I}\Delta t$ and

[TABLE]

Stability corresponds to the requirement that all roots of (4.11) have modulus less or equal one and that all multiple roots have modulus less than one.

In Figure 2 we plot the stability regions of AP-explicit IMEX-BDF schemes with respect to the variable $z_{R}$ and $z_{I}$. The contour lines represent different values of the scaling parameter ${\varepsilon}^{\alpha}$. Note that, since we are assuming to use central differences, the stability regions are inversely proportional to the value of ${\varepsilon}^{\alpha}$, since ${\varepsilon}^{1-\alpha}$ measures the strength of the diffusive term in agreement with (4.2). As expected, as the order of the methods increase the corresponding stability requirements become stronger and the various stability regions are reduced.

4.2. AP-implicit methods

Next we apply an IMEX-LM method to (4.4) in the AP-implicit form

[TABLE]

Thus, restricting to IMEX-BDF methods, by direct substitution of (4.8) into the second equation (4.13) we get

[TABLE]

or equivalently

[TABLE]

Now, the above LMM method in the limit $\varepsilon\to 0$ for $\alpha=1$ leads to the scheme

[TABLE]

which corresponds to an implicit-explicit IMEX-BDF scheme for the convection-diffusion equation.

The characteristic equation associated to scheme (LABEL:eq:scstab2) takes the form

[TABLE]

We report in Figure 3 the stability regions of the AP implicit IMEX-BDF methods with respect to the variable $z_{R}$ and $z_{I}$. The contour lines, as for the AP explicit case, represent different values of the scaling parameter ${\varepsilon}^{\alpha}$. The second order method is uniformly stable when ${\varepsilon}^{\alpha}<0.5$, all other methods show better stability properties compared to the AP-explicit case, in particular for large values of $z_{R}$. Again the stability regions diminish for increasing values of ${\varepsilon}^{\alpha}$, in agreement with the limit problem (4.2), and as the order of the methods increases.

5. Space discretization

In this section we briefly discuss the space discretization adopted in the numerical examples. For the hyperbolic fluxes, we consider a WENO method of order five [48] combined with a Rusanov flux. We stress that the space discretization is not constructed over the original discretized systems, namely (3.1), (3.17) and (LABEL:eq:SPs_vec1). Instead, we introduce the space discretization on the reformulated systems (3.7) for the AP-explicit with $\alpha=1$, on (3.18) for the AP-explicit with $\alpha\in[0,1)$ and on (3.27) for the AP-implicit case. In fact, the adopted IMEX partitioning of the system, which guarantees boundedness of the eigenvalues, is of paramount importance to avoid instabilities of the fluxes and excessive numerical dissipation typical of diffusive scaling limits [31, 40, 41]. As a consequence, the numerical diffusion is chosen accordingly to (3.20) in the numerical fluxes reported below.

Given a generic flux function $F(Q)$ of $Q\in\mathbb{R}^{n}$, we first reconstruct the unknown values at the interfaces $Q_{-},Q_{+}$ and successively we employ the numerical Rusanov flux defined as follows

[TABLE]

where $\max_{q\in[Q^{-},Q^{+}]}\{|\lambda(F^{\prime}(Q))|\}$ represents the maximum modulus of the eigenvalues of the Jacobian matrix $F^{\prime}(Q)$ and $S\in\mathbb{R}^{n\times n}$ a transformation matrix. Hence, for systems (3.7), (3.18), (3.27) and according to (3.10), (3.20), this value will depend on the scaling parameter $\varepsilon$ and the choice of the discretization steps. In particular for the general hyperbolic system (1.1), independently on the scaling factor $\alpha$, we have two unknowns $Q=(u,v)^{T}$ and three fluxes

[TABLE]

and

[TABLE]

where to comply with (5.1) we consider

[TABLE]

The generic variables $w$ reconstructed at the grid interfaces $(i+\frac{1}{2})$ and $(i-\frac{1}{2})$ respectively on the right $w^{-}_{i+\frac{1}{2}}$ and on the left side $w^{+}_{i-\frac{1}{2}}$ are given by

[TABLE]

with weights

[TABLE]

and with standard smoothness indicators

[TABLE]

where ${\varepsilon}=10^{-8}$, $d_{0}=3/10=\tilde{d}_{2}$, $d_{1}=3/5=\tilde{d}_{1}$ and $d_{2}=1/10=\tilde{d}_{0}$. Finally, the values $w^{(r)}_{i+\frac{1}{2}}$ and $w^{(r)}_{i+\frac{1}{2}}$ represent the third order reconstructions of the pointwise values $\bar{w}_{i}$. These are obtained through the formulas

[TABLE]

where $\bar{w}_{i-r+j}$ are the pointwise values of the unknown evaluated at the points $S_{r}(i)=\{x_{i-r},..,x_{i-r+2}\}$ $\ r=0,1,2$. Since, we use equispaced grid points, the coefficients $c_{rj}$ can be precomputed and their values are reported in Table 1.

In addition, we have a second order term $\partial_{xx}p(u)$ which, for the AP-explicit case, may be discretized by two consecutive application of the Rusanov flux with WENO reconstruction of the state variables and with numerical diffusion $\Theta(u,v)$ fixed equal to zero, or by a specific space discretization which is consistent with the limit problem. For example, by a sixth order finite difference formula

[TABLE]

with $a=1/90,b=-3/20,c=3/2,d=-49/18$. This latter approach has been adopted in the case of AP-implicit schemes, since the term $\partial_{xx}p(u)$ is implicit and we want to avoid nonlinearities induced by the WENO reconstructions.

6. Numerical validation & applications

In this Section, we present different numerical tests to validate the analysis performed in the previous Sections. In particular, we report results for the IMEX linear multistep from order two up to order five both for the AP-explicit and for the AP-implicit formulations. For the details about the IMEX-LM methods used we refer to Appendix A. We remark that other IMEX-LM methods can be included as well in the present formulation, see for example [45, 47]. In all test cases, the initial data has been chosen well prepared and the IMEX-LM methods have been initialized with a third order IMEX Runge-Kutta scheme (see [7]) with a time step which satisfies the accuracy constraints.

6.1. Test 1. Numerical convergence study for a linear problem.

We consider the following linear hyperbolic model with diffusive scaling for $\alpha=1$

[TABLE]

where $\gamma>0$. In the diffusive limit $\varepsilon\to 0$ the second equation relaxes to the local equilibrium

[TABLE]

substituting into the first equation this gives the limiting advection-diffusion equation

[TABLE]

In particular, we consider the model (6.1) solved on the domain $x\in[0,1]$, with $\gamma=1$ and periodic boundary conditions and with smooth initial data given by

[TABLE]

Note that, the initial data is well prepared, in the sense that $v(x,0)=u(x,0)-\partial_{x}u(x,0)$. For this specific problem, we numerically estimate the order of convergence of the schemes for various values of $\varepsilon=1,0.1,0.01,0.001$ by measuring the space and time $L_{1}$-error of the numerical solutions computed by using as reference solution the thinner grids. Namely, given a coarser grid $\Delta x_{1}=1/N$ with $N=128$ we consider

[TABLE]

The time step for AP-explicit methods is chosen as $\Delta t=\lambda\Delta x\max\left\{\varepsilon,\Delta x\right\}$, with $\lambda=0.25$. Namely, the largest between the CFL condition imposed by the hyperbolic part and by the limiting parabolic part of the equations. Instead, for implicit schemes we choose $\Delta t=\lambda\Delta x\max\left\{\varepsilon,1\right\}$ since the diffusion term is integrated implicitly in the limit.

The local truncation error is measured for the two components $u$ and $v$ as follows

[TABLE]

and the order of convergence is estimated by computing the rate between two $L_{1}$ errors of distinct numerical solutions

[TABLE]

The analysis is performed for several different IMEX linear multi-step schemes from second to fifth order (see Appendix A). In particular, we focus on the BDF and the TVB classes of IMEX multistep methods thanks to their favorable stability properties (see [26, 27] for details and derivation). We report in Table 2 the space-time $L_{1}$ errors and the relative rates of convergence for increasing size of the meshes considering $N=2^{k}$ points with $k=8,..,11$ for the $u$ variable, while the space-time $L_{1}$ errors and the relative rates of convergence for the $v$ variable are shown in Table 3 for the AP-explicit schemes. In Table 4 and 5, we report the corresponding $L_{1}$ errors and rates of convergence for the AP-implicit schemes. In all cases we can conclude that the expected orders of convergence are achieved by the schemes for the different values of the asymptotic parameter ${\varepsilon}$ and that the behavior of the schemes outperforms the corresponding IMEX Runge-Kutta methods for the non conserved quantity $v$ (see Table 2 in [8] for example). In particular, we observe the tendency of $4$-th order methods to achieve higher then expected convergence rates on the conserved quantity $u$ for moderate values of the stiffness parameter. The same tendency, on both variables $u$ and $v$, is observed for $5$-th order methods close to the diffusion limit particularly in the AP-implicit case. On the contrary, the $4$-th order schemes in the AP-explicit implementation produce a slight deterioration on the non conserved quantity $v$ which is not observed in the AP-implicit setting. The SG($3$,$2$) scheme also suffers of a slight deterioration of accuracy close to the diffusion limit in the AP-explicit form. It should be noted that scheme in AP-explicit form close to the fluid limit have a smaller truncation error with respect to time thanks to the CFL condition $\Delta t=O(\Delta x^{2})$. This can be observed by comparing the $L_{1}$-errors of the AP-explicit and AP-explicit formulations for $\varepsilon=0.001$. Finally, in Figure 4 and 5, we summarized in a plot the order of convergence for the mesh corresponding to 256 nodes as a function of $\varepsilon$ for the AP-explicit and AP-implicit methods. These plots emphasize that the convergence rate for all schemes is almost uniform.

6.2. Test 2: Riemann problem for the linear model.

Next, we consider a Riemann problem defined on the space interval $[0,4]$ with discontinuous initial data as follows

[TABLE]

For the above initial data, the linear hyperbolic system in the form (1.1) with zero-flux boundary conditions is solved and comparisons with different values of the relaxation parameter $\varepsilon$ are shown. The same problem has been studied in [8] using IMEX Runge-Kutta schemes.

In the limit $\varepsilon\to 0$, the exact solution of the corresponding advection-diffusion equation is known and it reads

[TABLE]

with $\textrm{erf}(x)$ the error function. We report in Figure 6 the numerical solution for $u$ at final time $T=0.25$ computed using two different time integration schemes, namely the IMEX-BDF2, and the IMEX-TVB(4,4). We choose $N=80$ points in space, and compare the AP-explicit approach (3.18), with the AP-implicit one (3.27). The reference solution is computed with the IMEX-BDF5 scheme using $\Delta t_{\textrm{ref}}=\Delta t/10$ and $N=200$ space points. In order to preserve the CFL conditions the time steps for the different regimes of $\varepsilon$ are selected as follows: in the AP-explicit case $\Delta t=\lambda\Delta x\max\{\varepsilon,\Delta x\}$ and in the AP-implicit case $\Delta t=\lambda\Delta x\max\{\varepsilon,1\}$ with $\lambda=0.4$. In Figure 6, for the diffusive limit, we observe that the different schemes agree well with the exact solution (6.7). In the hyperbolic regime, $\varepsilon=0.5$, the shock is correctly captured compared to the reference solution by both schemes with a slightly better resolution in the case of IMEX-TVB(4,4) for the AP-implicit approach. Note that, no spurious oscillations are observed for IMEX-BDF2 scheme, even if it does not satisfy any specific TVB stability property in the hyperbolic regime.

6.3. Test 3: Barenblatt solution for the porous media equation.

We then consider a nonlinear diffusion limit, by considering the following hyperbolic system with diffusive relaxation

[TABLE]

This problem has been previously studied in [33, 41]. Note that, when $k(u)=1$, the system (6.8) is a model of relaxing heat flow and, as ${\varepsilon}\to 0$, it relaxes towards the heat equation

[TABLE]

On the other hand, by choosing $k(u)=(2u)^{-1}$, the limiting equation for this model results in the porous media equation

[TABLE]

For this problem we consider the IMEX-BDF2 and the IMEX-BDF5 schemes for $\alpha=0$ by computing the numerical solution with $N=80$ mesh points, using the AP-explicit approach and selecting the time step as $\Delta t=\lambda\Delta x\max\{\varepsilon,\Delta x\}$ with $\lambda=0.4$.

The numerical solution in the limit $\varepsilon\to 0$ is compared with the analytical Barenblatt solution for the porous media equation [6]

[TABLE]

where $r(t)=[12(t+1)]^{1/3},~{}~{}t\geq 0$ and $x\in[-10,10]$. Note that the above solution defines also the initial state of the system where in addition we considered $v(x,t=0)=0$. In Figure 7, this comparison is shown using ${\varepsilon}=10^{-6}$ at time $T=3$. The sharp front of the Barenblatt solution is very well captured by both schemes without observing any relevant different by the increase order of accuracy of the IMEX-BDF5 method. This is somewhat expected, since the CFL condition $\Delta t=O(\Delta x^{2})$, in practice, gives to the IMEX-BDF2 scheme the same accuracy of a fourth order method. Therefore, the accuracy barrier here is represented by the $5$-th order space discretization method.

6.4. Test 4: Applications to the Ruijgrook–Wu model.

We consider, in this last test case, an application of the schemes to the so-called Ruijgrook–Wu model of rarefied gas dynamic [46, 23, 22, 41, 8, 33]

[TABLE]

where $f^{+}$ and $f^{-}$ denote the particle density distribution at time $t$, position $x$ and with velocity $+1$ and $-1$ respectively. Here $Kn$ is the Knudsen number, $M$ is the Mach number of the system and $a$,$b$ and $c$ are positive constants which characterize the microscopic interactions . The local (Maxwellian) equilibrium $f_{\infty}^{\pm}$ is characterized by

[TABLE]

The macroscopic variables for the model are the density $u$ and momentum $v$ defined by

[TABLE]

The nondimensional multiscale problem is obtained taking $M={\varepsilon}^{\alpha}$ and $Kn={\varepsilon}$, the Reynolds number of the system is then defined as usual according to $Re=M/Kn=1/{\varepsilon}^{1-\alpha}$. In macroscopic variables taking $a=b=1/2$, $c=M={\varepsilon}^{\alpha}$ this can be written as [22]

[TABLE]

The model has the nice feature to provide nontrivial limit behaviors for several values of $\alpha$ including the corresponding compressible Euler ($\alpha=0$) limit and the incompressible Euler ($\alpha\in(0,1)$) and Navier-Stokes ($\alpha=1$) limits (see [23, 22]).

In [41] the above model has been used under a similar but different scaling. If we denote with the pair ${\tilde{\alpha}},{\tilde{\varepsilon}}$ the scaling parameters in [41], these correspond to take $M={\tilde{\varepsilon}}$, $Kn={\tilde{\varepsilon}}^{1+{\tilde{\alpha}}}$. The two scaling can be made equivalent for $\alpha\neq 0$, if we map the pair $\alpha,{\varepsilon}$ into ${\tilde{\alpha}},{\tilde{\varepsilon}}$ as follows

[TABLE]

Note, that the non-linearity on the source term depends both on $u$ and $v$, so that we have $f=f(u,v)$ in (1.1). Nevertheless, following the same strategy described in the previous Sections, our methods can be applied in a straightforward way also in this situation.

For the Ruijgrook–Wu model (6.15) it can be shown, via Chapman-Enskog expansion, that for $\alpha\in(1/3,1]$ so that $2\alpha>1-\alpha$ and small values of ${\varepsilon}$ we have

[TABLE]

Then, the solution behavior is characterized by the viscous Burgers equation

[TABLE]

In the sequel, for Test 4a we apply the AP-explicit approach with $\Delta t=\lambda\Delta x\max\{{\varepsilon},\min\{1,\Delta x/{\varepsilon}^{1-\alpha}\}\}$, $\lambda=0.4$ which again corresponds to the maximum between the CFL of the hyperbolic part and the one originated by the limiting viscous Burger equation, whereas for Test 4b and Test 4c we use the AP-implicit approach with $\Delta t=\lambda\Delta x\max\{{\varepsilon},1\}$ and $\lambda=0.1$. Therefore, for ${\varepsilon}\leq 1$ the AP-implicit approach uses the same CFL condition in all test cases independently of ${\varepsilon}$ and $\alpha$. The reference solution is computed using the IMEX-BDF5 scheme with time step $\Delta t_{\textrm{ref}}=\Delta t/10$ and $N_{\textrm{ref}}=400$ space points. The initial data for the non conserved quantity $v$ has been taken well prepared accordingly to (6.17). In the sequel we will restrict to ${\tilde{\alpha}}\in[0,1]$ as in [41], therefore, from (6.16), we have $\alpha\in[0.5,1]$. Note, however, that the schemes can be applied for any value of $\alpha>0$. We refer to [22] for results in the fluid limit $\alpha=0$.

Test 4a. Riemann problem

We select the Riemann problem characterized by (6.6) as initial data. The numerical parameters are also the same of the previous test case, therefore we fix the final time $T=0.25$, and we compute the solution $u(x)$ using two different time AP-explicit integration schemes, BDF2 and TVB(4,4) with $N=100$ points in space.

In Figure 8 we report in the left column the behavior in the rarefied (non stiff) regime $\varepsilon=0.5$ and in the right-column the limit behavior $\varepsilon=10^{-6}$, whereas top row depicts the behavior for $\alpha=1$ and bottom row $\alpha=2/3$. Even if the rarefied solutions for the different values of $\alpha$ are similar, we remark the difference of the solution profiles in the limit ${\varepsilon}\to 0$. For $\alpha=1$ we obtain the classical viscous Burger equation where the discontinuity is smeared out by the diffusion term, while for $\alpha=2/3$ we have the sharp shock formation of the inviscid Burger equation. In both cases the methods yield an accurate description of the dynamic without spurious oscillations or excessive numerical dissipation.

Test 4b. Propagation of a square wave

Next, we consider the Ruijgrook–Wu model, in the space interval $[-0.5,0.5]$ with initial data defined as follows

[TABLE]

where we account for reflecting boundary conditions, i.e. $v=0,\partial_{x}u=0$ on the boundaries $x=\pm 0.5$.

This test problem has been previously studied in [8, 41] with IMEX Runge-Kutta methods. We study the solution to (6.15) for three different regimes of the parameters $\alpha$ and ${\varepsilon}$ and we solve the model with $N=100$ space grid points using the AP-implicit IMEX-BDF4 method. We report in Figure 9 the evolution of the solution for the density $u$ and the momentum $v$, respectively on the top and bottom rows. The first column represents the rarefied regime, for ${\varepsilon}=0.7,\alpha=1$. In this regime, the transport part dominates and the initial data propagates in the directions of the particles. This behavior is well described by the method without spurious numerical oscillations. In the second column we have the hyperbolic limit for ${\varepsilon}=10^{-12},\alpha=2/3$, corresponding to the inviscid Burger equation with a shock propagating in the right direction. Even in this case the shock profile is well captured. Finally, in the last column we report the parabolic limit for ${\varepsilon}=10^{-10},\alpha=4/5$, corresponding to a viscous Burger equation. As expected, the shock profile is regularized by the presence of the diffusive term.

Test 4c. Anisotropy of the multiscale parameter $\alpha$.

In the last test case, we solve the model (6.15) considering a multiscale parameter $\alpha$ to be a function of the space $x$, whereas the relaxation parameter is fixed to ${\varepsilon}=10^{-8}$. This test aims at reproducing a realistic situation where the scaling terms may depend on the physical quantities and vary in the different regions of the computational domain. We report the results obtained using the AP-implicit scheme with IMEX-BDF3 and the same discretization parameters of the previous test case. The initial date is defined in the space interval $[-0.5,0.5]$ as follows

[TABLE]

In the left column of Figure 10, we report the value of function $\alpha(x)$ as a function of the space, varying between $0.5$ (hyperbolic regime) and $1$ (parabolic regime). The central and right columns depict respectively the evolution of $u(x,t)$ and $v(x,t)$ showing the initial and final profiles, a similar test case has been presented in [8] for IMEX Runge-Kutta methods.

In the first row, we account a single variation of the regime from the hyperbolic to the parabolic

[TABLE]

At final time $T=0.05$ we observe a rarefaction wave moving to the left (in the hyperbolic regime) and a smooth profile on the right (in the parabolic regime). Note that, the method captures well the complicated shock structure even at the interface between the two regions.

Second row considers two variations of the the regime from hyperbolic to parabolic, choosing $\alpha(x)$ as

[TABLE]

at final time is $T=0.1$ we observe that the discontinuous initial data $u_{0}(x)$ is blunted within the parabolic regime, whereas a shock and rarefaction waves emerge in the hyperbolic one, both well described by the numerical method.

7. Conclusions

In this work we have developed a unified IMEX multistep approach for hyperbolic balance laws under different scalings. These problems, inspired by the classical hydrodynamical limits of kinetic theory [13], are challenging for numerical methods since the nature of the asymptotic behavior is not known a-priori and depends on the scaling parameters. Therefore, the schemes should be able to capture correctly asymptotic limits characterized by hyperbolic conservation laws as well as diffusive parabolic equations. A major difficulty in the schemes construction is represented by the unbounded growth of the characteristic speeds of the system in diffusive regimes. For these problems, we developed two different kind of approaches, originating a problem reformulation with bounded characteristic speeds, and which lead respectively to explicit (AP-explicit) or explicit-implict (AP-implicit) time discretizations of the asymptotic limit. Several numerical results for linear and non linear hyperbolic relaxation systems have confirmed that the IMEX multistep methods are capable to describe correctly the solution for a wide range of relaxation parameters and for different values of the scaling coefficient $\alpha$. Compared to the IMEX Runge-Kutta approach developed in [8] the IMEX multistep schemes here constructed present several advantages. In particular, it is possible to achieve easily high order accuracy and in general a more uniform behavior of the error with respect to the scaling parameters is observed. In addition, when dealing with computationally challenging problems such as the case of kinetic equations with stiff collision terms its is possible to strongly reduce the number of evaluations of the most expensive part of the computation represented by the source term. Future research will be in the direction of extending the present results to the more difficult case of diffusion limits for non linear kinetic equations and, more in general, to the case of low Mach number limits and all Mach number flows [7, 28, 29, 37, 25].

Appendix A Order conditions for IMEX-LM methods and examples

In this appendix we give the details of the particular IMEX-LM methods used in the manuscript. Let us recall that an order $p$, $s$-step, IMEX-LM scheme is obtained provided that

[TABLE]

Moreover the following theorem holds true [4]

Theorem 1**.**

For an $s$-step IMEX scheme we have

1. (1)

If $p\leq s$ the $2p+1$ constraints of (A.1) are linearly independent, therefore there exist $s$-step IMEX-LM schemes of order $s$. 2. (2)

An $s$-step IMEX-LM scheme has accuracy at most $s$. 3. (3)

The family of $s$-step IMEX-LM schemes of order $s$ has $s$ parameters.

Listed below the IMEX-LM methods analyzed along the paper, for further details and additional methods we refer to [1, 2, 19, 26, 45, 47].

Acknowledgments

The authors are grateful to the unknown referees for having reported a series of inaccuracies in the first version of the manuscript and for the indications that led to this improved revision.