Convergence of Adaptive Filtered Schemes for First Order Evolutionary Hamilton-Jacobi Equations

TL;DR

This paper introduces an adaptive, automatic method for choosing the switching parameter in filtered schemes for first order Hamilton-Jacobi equations, improving accuracy and robustness through smoothness indicators and convergence analysis.

Contribution

It proposes a novel adaptive and automatic approach for selecting the filter parameter in filtered schemes, enhancing convergence and accuracy for Hamilton-Jacobi equations.

Findings

The adaptive scheme converges with proven error estimates.

Numerical tests show improved accuracy over non-adaptive methods.

The method effectively detects smoothness regions for optimal switching.

Abstract

We consider a class of "filtered" schemes for first order time dependent Hamilton-Jacobi equations and prove a general convergence result for this class of schemes. A typical filtered scheme is obtained mixing a high-order scheme and a monotone scheme according to a filter function $F$ which decides where the scheme has to switch from one scheme to the other. A crucial role for this switch is played by a parameter $\varepsilon=\varepsilon({\Delta t,\Delta x})>0$ which goes to 0 as the time and space steps $(\Delta t,\Delta x)$ are going to 0 and does not depend on the time $t_n$, for each iteration $n$. The tuning of this parameter in the code is rather delicate and has an influence on the global accuracy of the filtered scheme. Here we introduce an adaptive and automatic choice of $\varepsilon=\varepsilon ^n (\Delta t, \Delta x)$ at every iteration modifying the classical set up. The adaptivity is controlled by a smoothness indicator which selects the regions where we modify the regularity threshold $\varepsilon^n$. A convergence result and some error estimates for the new adaptive filtered scheme are proved, this analysis relies on the properties of the scheme and of the smoothness indicators. Finally, we present some numerical tests to compare the adaptive filtered scheme with other methods.

Tables

Table 1. Table 1 : (Example 1a.) Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		F-HC ($5\Delta x$)		AF-HC		AF-LW4ord		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$40$	$10$	$1.36$e-$02$		$1.70$e-$02$		$7.88$e-$03$		$8.02$e-$02$
$80$	$20$	$2.56$e-$03$	$2.41$	$2.56$e-$03$	$2.73$	$8.66$e-$06$	$9.83$	$2.62$e-$02$	$1.62$
$160$	$40$	$5.76$e-$04$	$2.15$	$5.76$e-$04$	$2.15$	$5.43$e-$07$	$4.00$	$4.50$e-$03$	$2.54$
$320$	$80$	$1.40$e-$04$	$2.04$	$1.40$e-$04$	$2.04$	$3.40$e-$08$	$4.00$	$1.95$e-$04$	$4.52$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$40$	$10$	$3.58$e-$02$		$3.29$e-$02$		$1.08$e-$02$		$2.07$e-$01$
$80$	$20$	$6.66$e-$03$	$2.43$	$6.66$e-$03$	$2.30$	$2.25$e-$05$	$8.90$	$4.14$e-$02$	$2.32$
$160$	$40$	$1.48$e-$03$	$2.17$	$1.48$e-$03$	$2.17$	$1.40$e-$06$	$4.01$	$5.09$e-$03$	$3.02$
$320$	$80$	$3.57$e-$04$	$2.05$	$3.57$e-$04$	$2.05$	$8.69$e-$08$	$4.01$	$3.08$e-$04$	$4.05$

Table 2. Table 2 : (Example 1b.) Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		F-HC ($10\Delta x$)		AF-HC		AF-LW4ord		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$50$	$50$	$3.46$e-$01$		$3.55$e-$01$		$3.31$e-$01$		$3.47$e-$01$
$100$	$100$	$1.41$e-$01$	$1.29$	$1.90$e-$01$	$0.90$	$1.72$e-$01$	$0.94$	$2.07$e-$01$	$0.75$
$200$	$200$	$9.69$e-$02$	$0.54$	$1.17$e-$01$	$0.70$	$9.72$e-$02$	$0.82$	$1.28$e-$01$	$0.70$
$400$	$400$	$7.29$e-$02$	$0.41$	$7.27$e-$02$	$0.69$	$5.47$e-$02$	$0.83$	$7.66$e-$02$	$0.74$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$50$	$50$	$4.34$e-$01$		$3.31$e-$01$		$2.68$e-$01$		$3.62$e-$01$
$100$	$100$	$1.41$e-$01$	$1.63$	$1.19$e-$01$	$1.47$	$9.27$e-$02$	$1.53$	$1.39$e-$01$	$1.39$
$200$	$200$	$4.24$e-$02$	$1.73$	$3.03$e-$02$	$1.98$	$1.30$e-$02$	$2.83$	$3.83$e-$02$	$1.86$
$400$	$400$	$1.38$e-$02$	$1.62$	$9.51$e-$03$	$1.67$	$3.07$e-$03$	$2.08$	$8.39$e-$03$	$2.19$

Table 3. Table 3 : (Example 1b.) CPU times in seconds.

$N_x$	$N_t$	F-HC	AF-HC	AF-LW4ord	WENO 2/3
$50$	$50$	$0.000s$	$0.001s$	$0.001s$	$0.002s$
$100$	$100$	$0.001s$	$0.004s$	$0.005s$	$0.006s$
$200$	$200$	$0.004s$	$0.016s$	$0.018s$	$0.026s$
$400$	$400$	$0.019s$	$0.061s$	$0.077s$	$0.095s$

Table 4. Table 4 : (Example 2a.) Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		F-LWR ($5\Delta x$)		AF-LWR		AF-LW4ord		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$40$	$8$	$1.96$e-$02$		$1.89$e-$02$		$1.95$e-$02$		$6.81$e-$02$
$80$	$16$	$4.48$e-$03$	$2.13$	$3.56$e-$03$	$2.41$	$1.04$e-$02$	$0.90$	$3.42$e-$02$	$1.00$
$160$	$32$	$1.06$e-$03$	$2.08$	$8.53$e-$04$	$2.06$	$1.45$e-$03$	$2.85$	$1.62$e-$02$	$1.08$
$320$	$64$	$2.56$e-$04$	$2.05$	$2.20$e-$04$	$1.96$	$2.31$e-$04$	$2.65$	$7.52$e-$03$	$1.11$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$40$	$8$	$1.52$e-$02$		$1.63$e-$02$		$1.28$e-$02$		$2.05$e-$02$
$80$	$16$	$3.78$e-$03$	$2.01$	$3.61$e-$03$	$2.17$	$1.11$e-$03$	$3.53$	$4.68$e-$03$	$2.13$
$160$	$32$	$8.94$e-$04$	$2.08$	$8.80$e-$04$	$2.04$	$7.48$e-$05$	$3.89$	$9.55$e-$04$	$2.29$
$320$	$64$	$2.09$e-$04$	$2.09$	$2.08$e-$04$	$2.08$	$7.14$e-$06$	$3.39$	$1.40$e-$04$	$2.78$

Table 5. Table 5 : (Example 2b.) Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		F-LWR ($5\Delta x$)		AF-LWR		AF-LW4ord		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$40$	$8$	$1.91$e-$02$		$2.35$e-$02$		$2.42$e-$02$		$2.33$e-$02$
$80$	$16$	$9.24$e-$03$	$1.04$	$3.37$e-$03$	$2.80$	$7.51$e-$03$	$1.69$	$1.02$e-$02$	$1.19$
$160$	$32$	$5.77$e-$03$	$0.68$	$1.58$e-$03$	$1.09$	$2.14$e-$03$	$1.81$	$4.10$e-$03$	$1.32$
$320$	$64$	$3.46$e-$03$	$0.74$	$7.09$e-$04$	$1.16$	$6.92$e-$04$	$1.63$	$1.22$e-$03$	$1.75$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$40$	$8$	$2.38$e-$02$		$2.24$e-$02$		$2.28$e-$02$		$2.96$e-$02$
$80$	$16$	$8.48$e-$03$	$1.49$	$5.70$e-$03$	$1.98$	$2.05$e-$03$	$3.48$	$7.04$e-$03$	$2.07$
$160$	$32$	$3.41$e-$03$	$1.32$	$1.82$e-$03$	$1.65$	$3.20$e-$04$	$2.68$	$1.43$e-$03$	$2.30$
$320$	$64$	$1.52$e-$03$	$1.17$	$5.84$e-$04$	$1.64$	$6.38$e-$05$	$2.33$	$2.82$e-$04$	$2.34$

Table 6. Table 6 : (Example 3.) T = 4 / ( 5 ​ π 2 ) 𝑇 4 5 superscript 𝜋 2 T=4/(5\pi^{2}) . Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		F-LWR ($10\Delta x$)		AF-LWR		AF-LW4ord		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$40$	$8$	$1.30$e-$02$		$9.61$e-$03$		$1.89$e-$03$		$1.04$e-$02$
$80$	$16$	$8.67$e-$03$	$0.59$	$2.77$e-$03$	$1.79$	$2.84$e-$04$	$2.73$	$2.12$e-$03$	$2.30$
$160$	$32$	$5.07$e-$03$	$0.77$	$7.24$e-$04$	$1.94$	$2.68$e-$05$	$3.41$	$1.82$e-$04$	$3.54$
$320$	$64$	$2.66$e-$03$	$0.93$	$1.83$e-$04$	$1.99$	$1.89$e-$06$	$3.83$	$2.05$e-$05$	$3.15$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$40$	$8$	$3.76$e-$03$		$3.13$e-$03$		$3.31$e-$04$		$3.67$e-$03$
$80$	$16$	$1.29$e-$03$	$1.54$	$8.20$e-$04$	$1.93$	$1.85$e-$05$	$4.16$	$6.57$e-$04$	$2.48$
$160$	$32$	$4.49$e-$04$	$1.52$	$2.04$e-$04$	$2.01$	$1.43$e-$06$	$3.70$	$5.43$e-$05$	$3.60$
$320$	$64$	$1.82$e-$04$	$1.30$	$5.09$e-$05$	$2.00$	$9.80$e-$08$	$3.86$	$2.98$e-$06$	$4.19$

Table 7. Table 7 : (Example 3.) T = 3 / ( 2 ​ π 2 ) 𝑇 3 2 superscript 𝜋 2 T=3/(2\pi^{2}) . Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		F-LWR ($10\Delta x$)		AF-LWR		AF-LW4ord		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$40$	$16$	$4.88$e-$02$		$5.53$e-$02$		$5.86$e-$02$		$3.89$e-$02$
$80$	$32$	$2.47$e-$02$	$0.98$	$2.50$e-$02$	$1.15$	$2.62$e-$02$	$1.16$	$1.61$e-$02$	$1.27$
$160$	$64$	$9.81$e-$03$	$1.33$	$9.99$e-$03$	$1.32$	$1.03$e-$02$	$1.34$	$5.12$e-$03$	$1.65$
$320$	$128$	$2.57$e-$03$	$1.93$	$2.59$e-$03$	$1.95$	$2.67$e-$03$	$1.95$	$8.40$e-$04$	$2.61$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$40$	$16$	$5.17$e-$03$		$5.38$e-$03$		$3.18$e-$03$		$3.69$e-$03$
$80$	$32$	$1.26$e-$03$	$2.03$	$1.28$e-$03$	$2.08$	$6.73$e-$04$	$2.24$	$6.94$e-$04$	$2.41$
$160$	$64$	$2.86$e-$04$	$2.14$	$2.87$e-$04$	$2.15$	$1.31$e-$04$	$2.36$	$8.67$e-$05$	$3.00$
$320$	$128$	$5.68$e-$05$	$2.33$	$5.68$e-$05$	$2.34$	$1.70$e-$05$	$2.95$	$6.40$e-$06$	$3.76$

Table 8. Table 8 : (Example 3.) T = 3 / ( 2 ​ π 2 ) 𝑇 3 2 superscript 𝜋 2 T=3/(2\pi^{2}) . CPU times in seconds.

$N_x$	$N_t$	F-LWR	AF-LWR	AF-LW4ord	WENO 2/3
$40$	$16$	$0.000s$	$0.000s$	$0.001s$	$0.000s$
$80$	$32$	$0.000s$	$0.001s$	$0.002s$	$0.001s$
$160$	$64$	$0.001s$	$0.004s$	$0.005s$	$0.004s$
$320$	$128$	$0.005s$	$0.016s$	$0.020s$	$0.016s$

Table 9. Table 9 : (Example 4.) Global Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		F-HC ($5\Delta x$)		AF-HC		AF-LW4ord		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$40$	$10$	$3.41$e-$02$		$1.87$e-$02$		$2.07$e-$02$		$1.40$e-$02$
$80$	$20$	$1.69$e-$02$	$1.01$	$8.08$e-$03$	$1.21$	$8.75$e-$02$	$1.24$	$4.88$e-$03$	$1.51$
$160$	$40$	$9.12$e-$03$	$0.89$	$3.07$e-$03$	$1.40$	$3.35$e-$03$	$1.38$	$1.32$e-$03$	$1.88$
$320$	$80$	$7.35$e-$03$	$0.31$	$2.89$e-$03$	$0.09$	$3.16$e-$03$	$0.09$	$2.08$e-$03$	$-0.65$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$40$	$10$	$6.06$e-$03$		$3.06$e-$03$		$2.89$e-$03$		$2.93$e-$03$
$80$	$20$	$1.97$e-$03$	$1.62$	$8.10$e-$04$	$1.92$	$6.85$e-$04$	$2.08$	$5.37$e-$04$	$2.45$
$160$	$40$	$6.46$e-$04$	$1.61$	$2.07$e-$04$	$1.97$	$1.68$e-$04$	$2.02$	$6.82$e-$05$	$2.98$
$320$	$80$	$2.21$e-$04$	$1.55$	$6.82$e-$05$	$1.60$	$5.07$e-$05$	$1.73$	$2.58$e-$05$	$1.40$

Table 10. Table 10 : (Example 4.) Local Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		F-HC ($5\Delta x$)		AF-HC		AF-LW4ord		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$40$	$10$	$5.52$e-$03$		$2.59$e-$03$		$4.04$e-$03$		$4.04$e-$03$
$80$	$20$	$1.58$e-$03$	$1.81$	$4.92$e-$04$	$2.40$	$1.37$e-$04$	$4.88$	$9.31$e-$04$	$2.12$
$160$	$40$	$2.97$e-$04$	$2.41$	$2.22$e-$04$	$1.15$	$1.46$e-$05$	$3.22$	$6.06$e-$05$	$3.94$
$320$	$80$	$8.02$e-$05$	$1.89$	$4.33$e-$05$	$2.36$	$1.32$e-$07$	$3.47$	$2.23$e-$06$	$4.77$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$40$	$10$	$2.09$e-$03$		$1.14$e-$03$		$6.23$e-$04$		$1.38$e-$03$
$80$	$20$	$4.47$e-$04$	$2.23$	$3.07$e-$04$	$1.91$	$2.08$e-$05$	$4.91$	$2.33e$-$04$	$2.57$
$160$	$40$	$9.08$e-$05$	$2.30$	$8.06$e-$05$	$1.93$	$1.16$e-$06$	$4.16$	$1.62$e-$05$	$3.84$
$320$	$80$	$2.13$e-$05$	$2.09$	$2.00$e-$05$	$2.01$	$5.13$e-$07$	$1.18$	$7.66$e-$07$	$4.40$

Table 11. Table 11 : (Example 5.) T = 4 / ( 5 ​ π 2 ) 𝑇 4 5 superscript 𝜋 2 T=4/(5\pi^{2}) . Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		LWR		F-LWR ($10\Delta x$)		AF-LWR		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$20$	$10$	$7.45$e-$02$		$7.75$e-$02$		$9.35$e-$02$		$8.66$e-$02$
$40$	$20$	$3.38$e-$02$	$1.14$	$5.12$e-$02$	$0.60$	$3.44$e-$02$	$1.44$	$3.59$e-$02$	$1.27$
$80$	$40$	$1.49$e-$02$	$1.18$	$3.25$e-$02$	$0.66$	$5.98$e-$03$	$2.52$	$1.30$e-$02$	$1.47$
$160$	$80$	$6.42$e-$03$	$1.22$	$1.94$e-$02$	$0.75$	$1.78$e-$03$	$1.75$	$4.87$e-$03$	$1.41$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$20$	$10$	$3.67$e-$02$		$4.42$e-$02$		$4.15$e-$02$		$3.71$e-$02$
$40$	$20$	$9.53$e-$03$	$1.94$	$1.21$e-$02$	$1.87$	$8.97$e-$03$	$2.21$	$1.00$e-$02$	$1.89$
$80$	$40$	$2.28$e-$03$	$2.06$	$4.29$e-$03$	$1.49$	$1.83$e-$03$	$2.29$	$1.95$e-$03$	$2.36$
$160$	$80$	$6.51$e-$04$	$1.81$	$1.70$e-$03$	$1.33$	$5.23$e-$04$	$1.81$	$4.50$e-$04$	$2.11$

Table 12. Table 12 : (Example 5.) T = 3 / ( 2 ​ π 2 ) 𝑇 3 2 superscript 𝜋 2 T=3/(2\pi^{2}) . Errors and orders in L ∞ superscript 𝐿 L^{\infty} and L 1 superscript 𝐿 1 L^{1} norms.

		LWR		F-LWR ($10\Delta x$)		AF-LWR		WENO 2/3
$N_x$	$N_t$	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord	$L^{\infty }$ Err	Ord
$20$	$20$	$1.69$e-$01$		$9.49$e-$02$		$1.12$e-$01$		$8.68$e-$02$
$40$	$40$	$6.39$e-$02$	$1.40$	$3.67$e-$02$	$1.37$	$3.66$e-$02$	$1.61$	$2.27$e-$02$	$1.93$
$80$	$80$	$3.23$e-$02$	$0.98$	$1.41$e-$02$	$1.38$	$1.46$e-$02$	$1.33$	$9.08$e-$03$	$1.32$
$160$	$160$	$2.64$e-$02$	$0.29$	$3.73$e-$03$	$1.91$	$3.86$e-$03$	$1.92$	$2.22$e-$03$	$2.03$
$N_x$	$N_t$	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord	$L^1$ Err	Ord
$20$	$20$	$6.18$e-$02$		$4.62$e-$02$		$5.40$e-$02$		$3.60$e-$02$
$40$	$40$	$1.74$e-$02$	$1.83$	$8.19$e-$03$	$2.50$	$8.79$e-$03$	$2.62$	$4.68$e-$03$	$2.94$
$80$	$80$	$4.54$e-$03$	$1.94$	$1.88$e-$03$	$2.13$	$1.82$e-$03$	$2.27$	$6.92$e-$04$	$2.76$
$160$	$160$	$1.13$e-$03$	$2.00$	$3.73$e-$04$	$2.33$	$3.83$e-$04$	$2.25$	$7.62$e-$05$	$3.18$

Table 13. Table 13 : (Example 5.) T = 3 / ( 2 ​ π 2 ) 𝑇 3 2 superscript 𝜋 2 T=3/(2\pi^{2}) . CPU times in seconds.

$N_x$	$N_t$	LWR	F-LWR	AF-LWR	WENO 2/3
$20$	$20$	$0.002s$	$0.002s$	$0.009s$	$0.010s$
$40$	$40$	$0.015s$	$0.016s$	$0.055s$	$0.054s$
$80$	$80$	$0.108s$	$0.127s$	$0.427s$	$0.324s$
$160$	$160$	$0.867s$	$0.902s$	$3.642s$	$2.569s$

Equations

$$\left\{\begin{array}[]{ll}v_{t}+H(v_{x})=0,&(t,x)\in[0,T]\times\mathbb{R},\\ v(0,x)=v_{0}(x),&x\in\mathbb{R},\end{array}\right.$$

$$u^{n+1}_{j}\equiv S^{F}(u^{n})_{j}:=S^{M}(u^{n})_{j}+\varepsilon\Delta tF\bigg{(}\frac{S^{A}(u^{n})_{j}-S^{M}(u^{n})_{j}}{\varepsilon\Delta t}\bigg{)},\quad j\in\mathbb{Z},$$

$$u_j^{n+1}=S^{AF}(u^n{)}_j:=S^M(u^n{)}_j+{\varphi }_j^n{\epsilon }^n\Delta tF\left(\frac{S^A(u^n{)}_j-S^M(u^n{)}_j}{{\epsilon }^n\Delta t}\right),$$

$$u_j^{n+1}\equiv S^M(u_j^n):=u_j^n-\Delta t\text{leavevmode}{h}^M(D^{-}{u}_j^n,D^{+}{u}_j^n)$$

$${\mathcal{E}}_M(v)(t,x):=\left|\frac{v(t+\Delta t,x)-S^M(v(t,\cdot ))(x)}{\Delta t}\right|\le C_M\left(\Delta t||v_{tt}|{|}_{\infty }+\Delta x||v_{xx}|{|}_{\infty }\right),$$

$$\frac{\partial h^M}{\partial p^{-}}(p^{-},p^{+})\ge 0,\frac{\partial h^M}{\partial p^{+}}(p^{-},p^{+})\le 0,$$

$$\frac{\Delta t}{\Delta x}\left(\frac{\partial h^M}{\partial p^{-}}(p^{-},p^{+})-\frac{\partial h^M}{\partial p^{+}}(p^{-},p^{+})\right)\le 1.$$

$$||S^M(u)-S^M(v)|{|}_{\infty }\le ||u-v|{|}_{\infty }.$$

$$v_t+|v_x|=0,$$

$$h^M(p^{-},p^{+}):=\max \{p^{-},-p^{+}\}.$$

$$h^M(p^{-},p^{+}):=\frac{1}{a^{+}-a^{-}}\left[a^{-}H(p^{+})-a^{+}H(p^{-})-a^{+}{a}^{-}(p^{+}-p^{-})\right],$$

$$h^M(p^{-},p^{+}):=H\left(\frac{p^{-}+p^{-}}{2}\right)-\frac{\theta }{2}(p^{+}-p^{-})$$

$$u_j^{n+1}=S^A(u^n{)}_j:=u_j^n-\Delta t{h}^A(D^{k,-}{u}_j,\dots ,D^{-}{u}_j^n,D^{+}{u}_j^n,\dots ,D^{k,+}{u}_j^n),$$

$${\mathcal{E}}_A(v)(t,x):=$$

$$\le C_{A,l}\left(\Delta t^l||{\partial }_t^{l+1}v|{|}_{\infty }+\Delta x^l||{\partial }_x^{l+1}v|{|}_{\infty }\right).$$

$$\frac{{\partial }^{k}v(t,x)}{\partial t^k}$$

$$=(-1{)}^k\frac{{\partial }^{k-2}}{\partial x^{k-2}}\left(H_p^{k-1}(v_x(t,x))\frac{\partial }{\partial x}H(v_x(t,x))\right),$$

$$v_{tt}=\frac{\partial }{\partial t}(-H(v_x))=-H_p(v_x)v_{xt}=-H_p(v_x)\frac{\partial }{\partial x}(-H(v_x))=H_p^2(v_x)v_{xx},$$

$$\frac{{\partial }^{k+1}v}{\partial t^{k+1}}$$

$$=\frac{\partial }{\partial t}\left((-1{)}^k\frac{{\partial }^{k-2}}{\partial x^{k-2}}\left(H_p^k(v_x))v_{xx}\right)\right)\text{ by inductive hypothesis }$$

$$=(-1{)}^k\frac{{\partial }^{k-2}}{\partial x^{k-2}}\left(\frac{\partial }{\partial t}\left(H_p^k(v_x))v_{xx}\right)\right)$$

$$=(-1{)}^k\frac{{\partial }^{k-2}}{\partial x^{k-2}}\left(\frac{\partial }{\partial p}\left(H_p^k(v_x)\right)v_{xt}{v}_{xx}+H_p^k(v_x)v_{xxt}\right)$$

$$=(-1{)}^k\frac{{\partial }^{k-2}}{\partial x^{k-2}}\left(\frac{\partial }{\partial x}(H_p^k(v_x))v_{xt}+H_p^k(v_x)\frac{\partial }{\partial x}(v_{xt})\right)$$

$$=(-1{)}^k\frac{{\partial }^{k-1}}{\partial x^{k-1}}\left(H_p^k(v_x)v_{tx}\right)$$

$$=(-1{)}^{k+1}\frac{{\partial }^{k-1}}{\partial x^{k-1}}\left(H_p^{k+1}(v_x)v_{xx}\right),$$

$$v(t+\Delta t,x)$$

$$=v(t,x)-\Delta tH(v_x(t,x))+$$

$$\sum _{k=2}^r\frac{(-\Delta t{)}^k}{k!}\frac{{\partial }^{k-2}}{\partial x^{k-2}}\left(H_p^k(v_x(t,x))v_{xx}(t,x)\right)+O(\Delta t^{r+1}),$$

$$v(t+\Delta t,x)=v(t,x)-\Delta tH(v_x(t,x))+\frac{\Delta t^2}{2}{H}_p^2(v_x(t,x))v_{xx}(t,x)+O(\Delta t^3).$$

$${\mathcal{E}}_A(v)(t,x):=$$

$$\le C_{A,l}\left(\Delta t^l||{\partial }_t^{l+1}v|{|}_{\infty }+\Delta x^l||{\partial }_x^{l+1}v|{|}_{\infty }\right).$$

$$u_t=h_{\ast }^A(D^{-}u(t,x),D^{+}u(t,x))),$$

$$h_{\ast }^A(D^{-}{v}_j^n,D^{+}{v}_j^n)=H(v_x(t^n,x_j))+O(\Delta x^2),$$

$$h_{\ast }^A(D^{-}{u}_j^n,D^{+}{u}_j^n)=H\left(\frac{D^{-}{u}_j^n+D^{+}{u}_j^n}{2}\right),$$

$$\left\{\begin{array}[]{l}u^{*}_{j}=u^{n}_{j}-\Delta th_{*}^{A}(D^{-}u^{n}_{j},D^{+}u^{n}_{j})\\ u^{n+1}_{j}=\frac{1}{2}u^{n}_{j}+\frac{1}{2}u^{*}_{j}-\frac{\Delta t}{2}h_{*}^{A}(D^{-}u^{*}_{j},D^{+}u^{*}_{j}).\end{array}\right.$$

$$h^A(D^{-}{u}_j^n,D^{+}{u}_j^n)=\frac{1}{2}\left[h_{\ast }^A(D^{-}{u}_j^n,D^{+}{u}_j^n)+h_{\ast }^A(D^{-}{u}_j^{\ast },D^{+}{u}_j^{\ast })\right].$$

Full text

**Keywords: ** High-order Filtered schemes, Hamilton-Jacobi equations, Convergence, Smoothness indicators

1 Introduction

Here we propose and analyze a new adaptive filter scheme and prove its convergence to the viscosity solution of the scalar evolutionary Hamilton-Jacobi equation

[TABLE]

where Hamiltonian $H$ and the initial data $v_{0}$ are Lipschitz continuous functions. A precise result of existence and uniqueness in the framework of weak viscosity solutions can be found in [4] and the precise setting of assumptions will be given in Sect. 2.

The accurate numerical solution of Hamilton-Jacobi (HJ) equations is a challenging topic of growing importance in many fields of application, e.g. control theory, KAM theory, image processing and material science. Due to the lack of regularity of viscosity solutions, this issue is delicate and the construction of high-order methods can be rather complicated and the proof of convergence is challenging. It is well known that simple monotone schemes are at most first order accurate as shown in [9] so monotonicity should be abandoned to get high-order convergence. Our goal is to present a rather simple way to construct convergent schemes to the viscosity solution $v$ of (1) with the property to be of high-order in the region of regularity.

In recent years a general approach to the construction of high-order methods using filters has been proposed by Lions and Souganidis in [21] and further developed by Oberman and Salvador [22]. It is also interesting to mention that filtered schemes were also used for second order problems by Froese and Oberman [15] for the Monge-Ampere equation and, more recently, by Bokanowski, Picarelli and Reisinger [7] who studied second order time dependent HJB equations. Let us remind that a typical feature of a filtered scheme $S^{F}$ is that at the node $x_{j}$ the scheme is combination of a high-order scheme $S^{A}$ and a monotone scheme $S^{M}$ according to a filter function $F$. The scheme is written as

[TABLE]

where $u^{n+1}_{j}:=u(t_{n+1},x_{j})$ is the numerical approximation at time $t_{n+1}$ and node $x_{j}$, $\varepsilon=\varepsilon(\Delta t,\Delta x)>0$ is a fixed parameter going to 0 as $(\Delta t,\Delta x)$ is going to 0 and does not depend on $n$. Filtered schemes are high-order accurate where the solution is smooth, monotone otherwise, and this feature is crucial to prove a convergence result as in [6]. Note that the choice of the parameter $\varepsilon$ is delicate because it plays a crucial role in the switching so its tuning is rather important (see [6] for a detailed discussion of this point). Then it seems natural to adapt its choice to the regularity of the solution in the cell via a smoothness indicator. Here we improve the filtered scheme (2) introducing an adaptive and automatic choice of the parameter $\varepsilon=\varepsilon^{n}$ at every iteration. To this end, we introduce a smoothness indicator to select the regions where we have to update the regularity threshold $\varepsilon^{n}$, this indicator is chosen according to the analysis proposed in [19] although other proposals with similar properties can be applied.

To set this paper into perspective let us remind that the construction of high order methods for hyperbolic equations has been a very active research area started by the seminal paper [18]. Several techniques have been proposed to improve the accuracy leading to essentially non oscillatory schemes ENO and weighted ENO (so called WENO) for conservation laws as in [17, 2, 16, 1], for a survey on these high-order techniques we refer to [24, 25]. More recently a centered and more efficient version (called CWENO) has been proposed in [10]. We should also mention that high-order methods have been proposed for Hamilton-Jacobi either extending the ENO approach as in [19, 20, 8] or by semi-lagrangian techniques as extensively discussed in [12]. For a recent survey on the numerical approximation of Hamilton-Jacobi equations we refer the interested reader to [11].

The paper is organized as follows:

In Sect. 2 we construct the new adaptive filtered scheme and present in detail all its building blocks, the main assumptions are given there. Sect. 3 is focused on the analysis of the smoothness indicators in one dimension. In Sect. 4 we state and prove the main convergence result (that was announced in [13]), some technical lemmas are proved in Appendix A at the end of this paper. Finally in Sect. 5 we present several tests to show the effectiveness of the adaptive scheme with respect to the basic filtered scheme and to other state-of-the-art methods. Sect. 6 contains the conclusions with final comments.

2 A new Adaptive Filtered scheme

Consider the first order evolutionary Hamilton-Jacobi equation (1) where the hamiltonian $H$ and the initial data $v_{0}$ are Lipschitz continuous functions. It is well known that with these assumptions we have the existence and uniqueness of the viscosity solution. Notice that to keep the ideas clear we are considering the most simple scalar case with the hamiltonian depending only on the derivative of the solution, with more general situations following directly.

Our aim is to present a rather simple way to construct convergent schemes to the viscosity solution $v$ of (1) with the property to be of high-order whenever some regularity is detected. Starting from the ideas of [6] on filtered schemes, we proceed in this study introducing a procedure to compute the regularity threshold $\varepsilon$ in an automatic way, in order to exploit the local regularity of the solution.

Let us begin defining a uniform grid in space $x_{j}=j\Delta x$, $j\in\mathbb{Z}$, and in time $t_{n}=t_{0}+n\Delta t$, $n\in[0,N]\cap\mathbb{N}$, with $(N-1)\Delta t<T\leq N\Delta t$. Then, we compute the numerical approximation $u^{n}_{j}=u(t_{n},x_{j})$ with the simple formula

[TABLE]

where $u^{n+1}_{j}:=u(t_{n+1},x_{j})$, $S^{M}$ and $S^{A}$ are respectively the monotone and the high-order scheme, $F$ is the filter function needed to switch between the two schemes, $\varepsilon^{n}$ is the switching parameter at time $t_{n}$ and $\phi_{j}^{n}$ is the smoothness indicator function at the node $x_{j}$ and time $t_{n}$. More details on the components of the schemes will be given in the following sections.

Notice that if $\varepsilon^{n}\equiv\varepsilon\Delta x$, with $\varepsilon>0$ and $\phi_{j}^{n}\equiv 1$, we get the Basic Filtered Scheme (2).

2.1 Assumptions on the schemes

In this section we present the basic components of our scheme, which are a monotone finite difference scheme $S^{M}$ and a high-order, possibly unstable, scheme $S^{A}$. Let us begin by giving the assumptions on the monotone scheme.

Assumptions on $S^{M}$.

(M1)

The scheme can be written in differenced form

[TABLE]

for a function $h^{M}(p^{-},p^{+})$, with $D^{\pm}u^{n}_{j}:=\pm\frac{u^{n}_{j\pm 1}-u^{n}_{j}}{\Delta x}$;

(M2)

$h^{M}$ is a Lipschitz continuous function;

(M3)

(Consistency) $\forall v$, $h^{M}(v,v)=H(v)$;

(M4)

(Monotonicity) for any functions $u,v$ $\quad u\leq v\quad\Rightarrow\quad S^{M}(u)\leq S^{M}(v)$.

Under assumption (M2), the consistency property (M3) is equivalent to say that for all functions $v\in C^{2}([0,T]\times\mathbb{R})$, there exists a constant $C_{M}\geq 0$ independent on $\Delta=(\Delta t,\Delta x)$ such that

[TABLE]

where $\mathcal{E}_{M}$ is the consistency error. The last relation clearly shows the bound on the accuracy of the monotone schemes, which are at most first order accurate even for regular solutions.

Remark 2.1**.**

As pointed out in [6], under the Lipschitz assumption (M2) the monotonicity property (M4) can be restated in terms of some quantities that can be easily computed. In fact, it is enough to require, for a.e. $(p^{-},p^{+})\in\mathbb{R}^{2}$,

[TABLE]

and the CFL condition

[TABLE]

We call the CFL number, dependent on the hamiltonian of the considered problem, the constant ratio $\lambda:=\frac{\Delta t}{\Delta x}$ such that (6) is satisfied. Notice that working with explicit finite difference schemes this number can always be computed.

An important consequence of property (M4) is the nonexpansivity in $L^{\infty}$ of the mapping $S^{M}$ (see **[9]**, page 8), that is, for any functions $u,v$,

[TABLE]

Example 2.2**.**

We give some examples of monotone schemes in differenced form which satisfy (M1)-(M4). Other examples may be found in the pioneering work [9] or in [24].

• •

For the eikonal equation,

[TABLE]

we can use the simple numerical hamiltonian

[TABLE]

• •

For general equations, instead, we recall the Central Upwind scheme of **[20]**

[TABLE]

with $a^{+}=\max\{H_{p}(p^{-}),H_{p}(p^{+}),0\}$ and $a^{-}=\min\{H_{p}(p^{-}),H_{p}(p^{+}),0\}$, using the usual notation $H_{p}$ for the derivative of $H$ with respect to $v_{x}$.

• •

Another numerical hamiltonian we could use is the Lax-Friedrichs hamiltonian

[TABLE]

where $\theta>0$ is a constant. The scheme is monotone under the restrictions $\max_{p}|H_{p}(p)|<\theta$ and $\theta\lambda\leq 1$.

Next, we define the requirements on the high-order scheme.

Assumptions on $S^{A}$.

(A1)

The scheme can be written in differenced form

[TABLE]

for some function $h^{A}(p^{-},p^{+})$ (in short), with $D^{k,\pm}u^{n}_{j}:=\pm\frac{u^{n}_{j\pm k}-u^{n}_{j}}{k\Delta x}$;

(A2)

$h^{A}$ is a Lipschitz continuous function.

(A3)

(High-order consistency) Fix $k\geq 2$ order of the scheme, then for all $l=1,\dots,k$ and for all functions $v\in C^{l+1}$, there exists a constant $C_{A,l}\geq 0$ such that

[TABLE]

It is interesting to notice that we are not making any assumption on the stability of the high-order scheme, that is because filtered schemes are able to stabilize a possibly unstable scheme.

Before giving some examples of high-order schemes satisfying (A1)-(A3), let us state an interesting property of the solution $v$ of (1) in case of enough regularity. Notice that we are considering the simplest case of $H$ dependent only on the derivative of $v$.

Lemma 2.1**.**

Let $v$ be the solution of (1). Then, if $v\in C^{r}\left(\Omega_{(t,x)}\right)$, $r\geq 2$, where $\Omega_{(t,x)}$ is a neighborhood of a point $(t,x)\in\Omega:=[0,T]\times\mathbb{R}$, it holds

[TABLE]

for $k=2,\dots,r$.

Proof.

Let us proceed by induction on $2\leq k\leq r$, omitting the dependence on $(t,x)$ to simplify the notation. For $k=2$, we have

[TABLE]

and the statement holds in this case. Suppose now that (11) holds for $2<k<r-1$, then we can compute

[TABLE]

as we wanted. ∎

Let us now consider the value of the solution at $v(t+\Delta t,x)$, with $\Delta t>0$ and its Taylor expansion of order $r\geq 2$ around the point $(t,x)$. Using Lemma 2.1, we can rewrite

[TABLE]

which for $r=2$ simply reads

[TABLE]

Remark 2.3**.**

Using this last relation we could show that, assuming (A1)-(A2), the consistency property is equivalent to require that for $l=2,\dots,k$, and for all $v\in C^{l+1},$

[TABLE]

Now, let us give some examples of high-order schemes satisfying (A1)-(A3) with $l=2$.

Example 2.4**.**

As a first example let us consider the class of schemes obtained combining a high-order in space numerical hamiltonian $h^{A}_{*}$ and the second order Runge-Kutta SSP (or Heun scheme). To explain the simple procedure, let us consider the semidiscrete problem

[TABLE]

where $h_{*}^{A}$, is a high-order in space numerical hamiltonian of second order,

[TABLE]

such as the simple second order central approximation

[TABLE]

then to obtain the same accuracy in time we discretize using the second order SSP Runge-Kutta scheme,

[TABLE]

The scheme can be written in differenced form in the sense of (A1)-(14) defining

[TABLE]

Notice that through this procedure the stencil of the scheme (16) becomes doubled for $h^{A}$. Notice also that this procedure can be easily extended to the case of hamiltonian dependent on the space variable $x$.

Example 2.5**.**

Then we propose a couple of numerical hamiltonians $h^{A}$ obtained discretizing directly the formula (13) or, equivalently, obtained from the same Lax-Wendroff schemes for conservation laws by the substitution $u^{n}_{j}=\frac{v^{n}_{j+1}-v^{n}_{j}}{\Delta x}$. The first is the original Lax-Wendroff scheme

[TABLE]

and the second is its variation proposed by Richtmyer,

[TABLE]

Example 2.6**.**

Following the approach of the Lax-Wendroff schemes and making use of the expansion (2.1), we can easily write higher order schemes, in both space and time, using very compact stencils. The idea is simply to discretize directly the above expansion using finite difference approximations of the right order. For example, if we want to write a fourth order Lax-Wendroff scheme using only five points, one of the possibilities is to define

[TABLE]

and then compute

[TABLE]

It is straightforward to verify that, if the solution $v$ is regular enough, using Taylor expansion we have

• •

$H_{1}=H(v_{x})+O(\Delta x^{4})$,

• •

$H_{2}=H^{2}_{p}(v_{x})v_{xx}+O(\Delta x^{4})$,

• •

$H_{3}=\frac{\partial}{\partial x}\left(H_{p}^{3}(v_{x})v_{xx}\right)+O(\Delta x^{2})$,

• •

$H_{4}=\frac{\partial^{2}}{\partial x^{2}}\left(H_{p}^{4}(v_{x})v_{xx}\right)+O(\Delta x^{2})$,

and that the resulting scheme satisfies (A1)-(A3) with $l=4$. Notice that to obtain fourth order it would have been enough to have approximations of one order lower for $H_{2}$ and $H_{4}$, but thanks to the symmetry of the discretizations we can get higher orders without increasing the number of points in the stencil.

2.2 Filter function

In order to couple the schemes and their properties, we need to define a function $F$, called filter function F, such that

(F1)

$F(x)\approx x$ for $|x|\leq 1$,

(F2)

$F(x)=0$ for $|x|>1$,

which implies that

• •

If $|S^{A}-S^{M}|\leq\Delta t\varepsilon^{n}$ and $\phi_{j}^{n}=1\Rightarrow S^{AF}\approx S^{A}$

• •

If $|S^{A}-S^{M}|>\Delta t\varepsilon^{n}$ or $\phi^{n}_{j}=0\Rightarrow S^{AF}=S^{M}$.

It is clear that, with just these two requirements, several filter functions can be considered, which differ for regularity properties. Four examples are reported in Fig. 1. The first filter, $F_{1}$, which we use in our numerical tests, has been defined in [6] as

[TABLE]

which is clearly discontinuous at $x=-1,1$.

As a second possibility, we propose the family of regular filter functions given by the formula

[TABLE]

for appropriate choices of the parameters $a$, $b$ and $c$. In Fig. 1, the filter $F_{2}$ belongs to that family with $a=0.25,\ b=20,\ c=4$. Functions of that kind are very regular ($F\in C^{\infty}$) and developing with Taylor we can see that they satisfy the properties (F1)-(F2).

Another example of filter functions satisfying (F1)-(F2), continuous but not necessarily derivable, is the following family of functions

[TABLE]

varying the parameters $a$ and $b$. In Fig. 1 we show $F_{3}(x)$ with $a=0.001$ and $b=1.05$, for which the function approach better the value $1$ at $x=-1,1$.

Finally, we recall also the filter defined in [15] as

[TABLE]

which satisfies (F1)-(F2) with a rather wide transition phase.

After extensive computations, we noticed that the results obtained with our adaptive filtered (AF) scheme are not sensitive with respect to changes in regularity of the filter function, even with very large transition phases. That is probably because, as we will see in the next section, the parameter $\varepsilon^{n}$ is designed to obtain the property (F1) whenever possible, then in regions of regularity of the solution the argument of $F$ lies most probably in $[-1,1]$, where all the filter functions are practically the same. Some major differences, instead, can be seen in the results obtained with the basic filtered scheme, for which the threshold $\varepsilon$ is fixed at the beginning, as it is highlighted in the introduction of [6].

2.3 Tuning of the parameter $\varepsilon^{n}$

The last step is to show how to compute the switching parameter $\varepsilon^{n}$, which is the real core of the adaptivity of our scheme. Then, if we want the scheme (3) to switch to the high-order scheme when some regularity is detected, we have to choose $\varepsilon^{n}$ such that

[TABLE]

in the region of regularity at time $t_{n}$, that is

[TABLE]

For the moment, to simplify the presentation we assume the existence of a function $\phi$ such that

[TABLE]

where $I_{j}=(x_{j-1},x_{j+1})$, referring to the next section for some examples of practical computation of the function $\phi$.

Assuming $v$ sufficiently regular and proceeding by Taylor expansions as in [6], we have for the monotone scheme

[TABLE]

where we used the relation

[TABLE]

whereas for the high-order scheme, by the consistency property,

[TABLE]

Whence, from (25) we obtain

[TABLE]

Finally, we use a numerical approximation of the lower bound on the right hand side of the previous inequality to obtain the following formula for $\varepsilon^{n}$,

[TABLE]

with $K>\frac{1}{2}$, $\lambda:=\frac{\Delta t}{\Delta x}$ and $D\,u^{n}_{j}:=\frac{u^{n}_{j+1}-u^{n}_{j-1}}{2\Delta x}$. Notice that if we assume enough regularity on the solution $v$, then (2.3) gives a second order approximation of the right hand side of (28) multiplied by $2K$.

3 Smoothness indicator function

In the previous section we assumed the existence of a smoothness indicator function $\phi$, in the sense that

[TABLE]

where $I_{j}=(x_{j-1},x_{j+1})$ and $\omega^{n}_{j}$ is the smoothness indicator at the node $x_{j}$ depending on the values of the approximate solution $u^{n}$. The aim of this section is precisely to show a simple construction of a function satisfying (27) which makes use of smoothness indicators widely known in literature. Moreover, in the process we review the theory of the smoothness indicators of [19], defined for the construction of the WENO schemes for (1),

[TABLE]

for $k=0,\dots,r-1$, where $P_{k}$ is the Lagrange polynomial of degree $r$ interpolating the values of $u^{n}$ on the stencil $\mathcal{S}_{j+k}=\{x_{j+k-r},\dots,x_{j+k}\}$.

Then, before proceeding with the construction of $\phi$, let us state a fundamental result on the behavior of the indicators (31).

Proposition 3.1**.**

Assume $f\in C^{r+1}\left(\Omega\setminus\{x_{s}\}\right)$, with $\Omega$ a neighborhood of $x_{s}$, and $f^{\prime}(x_{s}^{-})\not=f^{\prime}(x_{s}^{+})$. Moreover, just for simplicity, let $f^{\prime\prime}(x)\neq 0$, $\forall x\in\Omega\setminus{x_{s}}$. Then, for $k=0,\dots,r-1$ and $j\in\mathbb{Z}$, the followings are true:

1. i)

If $x_{s}\in\Omega\ \setminus\stackrel{{\scriptstyle\circ}}{{\mathcal{S}}}_{j+k}\quad\Rightarrow\beta_{k}(f)=O(\Delta x^{2})$, 2. ii)

If $x_{s}\in\ \stackrel{{\scriptstyle\circ}}{{\mathcal{S}}}_{j+k}\quad\Rightarrow\beta_{k}(f)=O(1)$,

where $\stackrel{{\scriptstyle\circ}}{{\mathcal{S}}}_{j+k}=(x_{j-r+k},x_{j+k})$.

We skip the proof, which is rather technical, but the interested reader can find it in Appendix A.

Remark 3.2**.**

Notice that we could avoid the restrictions on $f^{\prime\prime}$ at points of regularity by adding a small quantity $\sigma_{\Delta x}:=\sigma\Delta x^{2}$, for some constant $\sigma>0$, to the indicators $\beta_{k}$ and consider instead

[TABLE]

This is necessary in order to avoid a reduction of accuracy at points such that $f^{\prime\prime}=0$, as it has been thoroughly discussed in [3] in the case of discontinuous functions. We will use this assumption in our numerical tests, choosing $\sigma=1$.

Our aim is to identify the points (or the intervals) in which a function $f$ presents a singularity in the first derivative using only its nodal values $f_{j}$, $j\in\mathbb{Z}$. Let us focus the attention on a point $x_{j}$ of the grid and consider the simplest case of $r=2$, which is enough for our purpose. Let us consider separately the intervals $(x_{j-1},x_{j}]$ and $[x_{j},x_{j+1})$ defining

[TABLE]

for $k=0,1$, where $P_{0}$, $P_{1}$ are the polynomials interpolating the function, respectively, on the stencils $\{x_{j-2},x_{j-1},x_{j}\}$ and $\{x_{j-1},x_{j},x_{j+1}\}$, and symmetrically

[TABLE]

for $k=0,1$, where now $P_{0}$, $P_{1}$ are the interpolating polynomials on the stencils $\{x_{j-1},x_{j},x_{j+1}\}$ and $\{x_{j},x_{j+1},x_{j+2}\}$. From the definition it is clear that $(\beta^{+})_{j}=(\beta^{-})_{j+1}$ so we have to compute the quantities just once (note that these quantities are always nonnegative). Then, we define as in [19]

[TABLE]

with $\sigma_{\Delta x}:=\sigma\Delta x^{2}$ the parameter we introduced in Remark 3.2, and focus on the information given by the interpolating polynomial on $\{x_{j-1},x_{j},x_{j+1}\}$ defining

[TABLE]

to inspect the regularity, respectively, on $[x_{j},x_{j+1})$ and for $(x_{j-1},x_{j}]$.

By Prop. 3.1 and Remark 3.2 we know that $\widetilde{\beta}_{k}=O(\Delta x^{2})$ if there is no singularity in the stencil, and $\widetilde{\beta}_{k}=O(1)$ otherwise, so in presence of a singularity we can only fall in one of the following cases:

• •

If $x_{j-2}<x_{s}\leq x_{j-1}$, then ${\widetilde{\beta}}^{-}_{0}=O(1)$, ${{\widetilde{\beta}}}^{-}_{1}={\widetilde{\beta}}^{+}_{0}=O(\Delta x^{2})$, ${\widetilde{\beta}}^{+}_{1}=O(\Delta x^{2})$,

• •

If $x_{j-1}<x_{s}<x_{j}$, then ${\widetilde{\beta}}^{-}_{0}=O(1)$, ${\widetilde{\beta}}^{-}_{1}={\widetilde{\beta}}^{+}_{0}=O(1)$, ${\widetilde{\beta}}^{+}_{1}=O(\Delta x^{2})$,

• •

If $x_{s}=x_{j}$, then ${\widetilde{\beta}}^{-}_{0}=O(\Delta x^{2})$, ${\widetilde{\beta}}^{-}_{1}={\widetilde{\beta}}^{+}_{0}=O(1)$, ${\widetilde{\beta}}^{+}_{1}=O(\Delta x^{2})$,

• •

If $x_{j}<x_{s}<x_{j+1}$, then ${\widetilde{\beta}}^{-}_{0}=O(\Delta x^{2})$, ${\widetilde{\beta}}^{-}_{1}={\widetilde{\beta}}^{+}_{0}=O(1)$, ${\widetilde{\beta}}^{+}_{1}=O(1)$,

• •

If $x_{j+1}\leq x_{s}<x_{j+2}$, then ${\widetilde{\beta}}^{-}_{0}=O(\Delta x^{2})$, ${\widetilde{\beta}}^{-}_{1}={\widetilde{\beta}}^{+}_{0}=O(\Delta x^{2})$, ${\widetilde{\beta}}^{+}_{1}=O(1)$,

with $x_{s}$ point of singularity. Now, we can compute

[TABLE]

which, noticing that, if the function is smooth in both stencils of $\beta_{0}^{\pm}$ and $\beta_{1}^{\pm}$, we have

[TABLE]

leads to

[TABLE]

Whence we can deduce that if the solution is regular enough in both stencils

[TABLE]

On the other hand, if there is a singularity in at least one of the stencils, by Prop. 3.1 and the definition (35) we have that

[TABLE]

then it is easy to verify that the behavior of our $\omega_{\pm}$ falls in the following cases:

• •

If $x_{j-2}<x_{s}\leq x_{j-1}$, then $\omega_{-}=1+O(\Delta x^{4})$, $\omega_{+}=1/2+O(\Delta x)$

• •

If $x_{j-1}<x_{s}<x_{j}$, then $\omega_{-}=O(1)$, $\omega_{+}=O(\Delta x^{4})$

• •

If $x_{s}=x_{j}$, then $\omega_{-}=O(\Delta x^{4})$, $\omega_{+}=O(\Delta x^{4})$

• •

If $x_{j}<x_{s}<x_{j+1}$, then $\omega_{-}=O(\Delta x^{4})$, $\omega_{+}=O(1)$

• •

If $x_{j+1}\leq x_{s}<x_{j+2}$, then $\omega_{-}=1/2+O(\Delta x)$, $\omega_{+}=1+O(\Delta x^{4})$,

where with $\omega_{\pm}=O(1)$ we mean a number dependent on the jump of the derivative. Now, defining $\omega_{j}:=\min\{\omega_{-},\omega_{+}\}$ we can rewrite

[TABLE]

Finally, what is left is to define the function $\phi$ such that $\phi=1$ if $\omega$ is close to $\frac{1}{2}$ and $\phi=0$, otherwise. Notice that in the latter are included both cases in which the function has a singularity in the first derivative ($\omega=O(\Delta x^{4})$) and when the second derivative is discontinuous ($\omega=O(1)$). The simplest choice is to take

[TABLE]

with $M<\frac{1}{2}$, a number possibly dependent on $\Delta x$.

Remark 3.3**.**

Notice that to construct the function $\phi$ using the indicators (31) with $r=2$ we need only five points to inspect the regularity in $I_{j}$.

Next, we show that if we make a particular choice for $M$ we are able to prove the following result, which can be seen as an “inverse” of Prop. 3.1 for numerical solutions and, if we use the previous simple construction for $\omega_{j}$, gives a useful tool for the analysis of the next section.

Before proceeding, let us remind that we are working with structured grids, then if we consider a one-parameter family of grid values $\{f_{j}(\Delta x)\}_{j\in J(\Delta x)}$, as $\Delta x$ goes to 0, the indexed family of sets of indices $J(\Delta x)$ is expanding, in the sense that if $\Delta x_{2}<\Delta x_{1}$, then $J(\Delta x_{1})\subset J(\Delta x_{2})$, where $J(\Delta x)\subseteq\mathbb{Z}$, for all $\Delta x>0$. Moreover, we define $I_{s}(\Delta x)$ as the set of indices $j$ such that $\phi_{j}=0$ and assume, for simplicity, $|I_{s}(\Delta x)|<\overline{I}_{s}$, where $\overline{I}_{s}$ is a positive constant.

Lemma 3.1**.**

Let $\omega$ be computed using (33)-(34) and $\phi$ be defined by (42) with $M(\Delta x)=\frac{1}{2}-C\Delta x$, for some constant $C$ such that $0<M(\Delta x)<\frac{1}{2}$. Consider a one-parameter family of sequences $\{f_{j}(\Delta x)\}_{j\in J(\Delta x)}$ with compact support in the interval $[-b,b]$, and a partition $\{R_{i}\}_{i=0,\dots,|I_{s}|}$ of the regularity set $\mathcal{R}=\{j\in\mathbb{Z}:\phi_{j}=1\}=\bigcup_{i}R_{i}$, and $\mathcal{R}=\mathbb{Z}$ if $I_{s}=\emptyset$. Then, if for all $i=0,\dots,|I_{s}|$, there exists $j_{i}\in R_{i}$, such that $|D^{2}f_{j_{i}}(\Delta x)|<\infty$, we have that

[TABLE]

for a constant $B$ independent of $\Delta x$.

Proof.

Since $\{f_{j}\}$ has compact support we have $|I_{s}|<\infty$ and it will be enough to prove the assertion just for one $i\in I_{s}$. More simply, in the regular case we have $\mathcal{R}=\mathbb{Z}$ and we want to show that the statement is true if the discrete derivative is bounded at some point, so if there exist an index $\widehat{j}\in\mathcal{R}$ and a positive constant $\delta$ (independent of $\Delta x$) such that $|D^{2}f_{\widehat{j}}(\Delta x)|\leq\delta$ (note that for the nodes outside the support of $f_{j}$ we can even set $\delta=0$). In the following, we simplify the notation dropping the dependence of $f_{j}$ on $\Delta x$. By definition of $\phi$ and $\omega$, if $\phi_{j}=1$ then both $\omega_{\pm}>M$. Moreover, (33)-(34) imply that the coefficients $\beta^{\pm}$ are always nonnegative as well as $\omega_{\pm}$.

Let us consider the case $j<\widehat{j}$. Then, by definition,

[TABLE]

which leads by simple computations to

[TABLE]

then, dividing by $\Delta x^{2}$ and recalling that $\sigma_{\Delta x}=\sigma\Delta x^{2}$ , we get

[TABLE]

Now let us iterate (44) on $j$ till $\widehat{j}$ and define $L_{j}\equiv\widehat{j}-j$, we have

[TABLE]

For $j>\widehat{j}$, we can use the relation $\omega_{-}>M$ and iterate back to $\widehat{j}$ redefining $L_{j}\equiv j-\widehat{j}$, the calculations are similar also for this case.

Since $L_{j}\Delta x$ is bounded by $b$, we have $L_{j}\leq\frac{b}{\Delta x}$, $\forall j\in\mathcal{R}$. Recalling that $M=\frac{1}{2}-C\Delta x$, we can use the previous bound on $|D^{2}f_{j}|^{2}$ to proceed

[TABLE]

We get the final bound passing to the limit for $\Delta x$ going to 0 in (45), in conclusion we get

[TABLE]

and the statement follows simply taking $B:=\sqrt{e^{3Cb}(\sigma+\delta^{2})-\sigma}$. ∎

Unfortunately, we noticed through numerical tests that the $O(\Delta x)$ term in regular regions may produce heavy oscillations around the optimal value $\overline{\omega}=1/2$. To increase the accuracy, we can use higher order smoothness indicator ($r>2$), but we would need a bigger reconstruction stencil. Otherwise, if we want to keep the compactness of the stencil, we can use the mappings defined in [16],

[TABLE]

which have the properties that $g(0)=0$, $g(1)=1$, $g(\overline{\omega})=\overline{\omega}$, $g^{\prime}(\overline{\omega})=0$ and $g^{\prime\prime}(\overline{\omega})=0$. Then, we define

[TABLE]

Notice that with respect to the definition in [16] we avoided the second weighting which seems unnecessary in our case. More explicitly, the mapping we use is

[TABLE]

It is important to remind that, at the moment, Lemma 3.1 is valid only for indicators $\omega$ using the standard construction for $r=2$, without the possibility to introduce any modification, or higher order indicators. Moreover, as it will be briefly discussed in Remark 4.4, it introduces some limitations in the applicability even when using the standard indicators, testifying the necessity of some improvements in the argument used. Notice that the previous lemma strongly relies on the fact that $\omega$ is computed using (33)-(34) without introducing the mappings (48). In fact, if we were to use (48), we could develop the algebra until the inequality

[TABLE]

but, by definition, $g^{-1}$ cannot be expanded in Taylor series around the point $\frac{1}{2}$, whence we could not use the notable limit to conclude.

Therefore, we are forced to add a “technical” assumption in order to justify the proof of Prop. 4.1. More precisely, when using the alternative constructions for $\omega$ (using the mapping (48)), we define the region of regularity $\mathcal{R}$ detected by the function $\widetilde{\phi}$ as the set

[TABLE]

for some constant $B\gg 0$. Notice that with this definition, which, we recall, is needed only for theoretical reasons, it is not necessary to require $M(\Delta x)\to 0$, then we can simply choose a constant $M>0$ small enough (e.g. $M=0.1$), as we will do in the numerical tests of Sect. 5.

4 Convergence result

We are now able to present our main result, but before doing so let us state a useful proposition about the numerical solution and the parameter $\varepsilon^{n}$.

Proposition 4.1**.**

Let $u^{n}$ be the solution obtained by the scheme (3)-(2.3) and assume that $v_{0}$ and $H$ are Lipschitz continuous functions. Assume also that $\mathcal{R}^{n}$ is defined by (26) or (49), with $\phi$ given by (42), and that $\lambda=\Delta t/\Delta x$ is a constant such that (6) is satisfied. Then, $\varepsilon^{n}$ is well defined and $u^{n}$ satisfies, for any $i$ and $j$, the discrete Lipschitz estimate

[TABLE]

for some constant $L>0$, for $0\leq n\leq T/\Delta t$. Moreover, there exists a constant $C>0$ such that

[TABLE]

Proof.

Before proceeding with the proof let us notice that, if $u^{n}$ satisfies (50) for a constant $L_{n}>0$, calling for brevity

[TABLE]

we have that

[TABLE]

whence we can conclude

[TABLE]

where $L_{h^{M}}$ is the Lipschitz constant of $h^{M}$, whereas $L_{H}$ and $L_{H2}$ are the local Lipschitz constant of $H$ on $[-L_{n},L_{n}]$ and $[-2L_{n}-\Delta tL_{H}B,2L_{n}+\Delta tL_{H}B]$, respectively. Notice that, if the smoothness indicators are computed using the definitions (33)-(34), we have

[TABLE]

Then, in such case by Lemma 3.1, $x_{j}\in\mathcal{R}^{n}\Rightarrow D^{2}u_{j}^{n}<B$, for some constant $B>0$ independent on $n$. Otherwise, we can obtain the same estimate by the definition (49) of $\mathcal{R}^{n}$.

Notice also that if the function $H$ is globally Lipschitz continuous we have the same estimate with $L_{H2}=L_{H}$, where now $L_{H}$ is the global Lipschitz constant of $H$. Consequently, the last statement would follow with $C=KB(\lambda L^{2}_{H}+2L_{h^{M}})$.

Let us now prove the main statement proceeding, as usual, by induction on $n\geq 0$ and noticing that it is sufficient to prove (50) for $i$ and $j$ such that $i=j\pm 1$.

For $n=0$, as we take $u^{0}_{j}=v_{0}(x_{j})$ for $j\in\mathbb{Z}$, we have that (50) is satisfied by the Lipschitz continuity assumption on $v_{0}$, with constant $L_{0}$.

Now, assuming that (50) is satisfied for $n-1>0$ so that $\varepsilon^{k}$ for $k=0,\dots,n-1$ are bounded by (4), we can compute

[TABLE]

then, iterating back and using the same arguments,

[TABLE]

where $C$ is well defined by (4). Notice that we have used the nonexpansivity in $L^{\infty}$ of $S^{M}$ and the fact that $|F|\leq 1$, $|\phi|\leq 1$. ∎

Therefore, it is clear that by construction our scheme is $\varepsilon$-monotone, in the sense of the following

Definition 4.2** ($\varepsilon$-monotonicity).**

A numerical scheme $S$ is $\varepsilon$-monotone if for any functions $u,v$,

[TABLE]

where $C$ is constant and $\varepsilon\rightarrow 0$ as $\Delta=(\Delta t,\Delta x)\rightarrow 0$.

Thanks to that property, by applying the Barles-Souganidis result [5], the convergence follows directly. We conclude this section with the following theorem, which gives us the order of convergence for the Adaptive Filtered Schemes.

Theorem 4.3**.**

Let the assumptions on $S^{M}$ and $S^{A}$ be satisfied. Assume that $v_{0}$ and $H$ are Lipschitz continuous functions, $u_{j}^{n+1}$ is computed by (3)-(2.3), with $K>1/2$ and $\lambda=\frac{\Delta t}{\Delta x}$, a constant such that (6) is satisfied. Assume also that $\mathcal{R}^{n}$ is defined by (26) or (49), with $\phi$ given by (42). Let us denote by $v_{j}^{n}:=v(t^{n},x_{j})$ the values of the viscosity solution on the nodes of the grid. Then,

1. i)

the AF scheme (3) satisfies Crandall-Lions estimate **[9]**

[TABLE]

for some constant $C_{1}>0$ independent of $\Delta x$. 2. ii)

(First order convergence for regular solutions) Moreover, if $v\in C^{2}([0,T]\times\mathbb{R})$, then

[TABLE]

for some constant $C_{2}>0$ independent of $\Delta x$. 3. iii)

(High-order local consistency) Let $k\geq 2$ be the order of the scheme $S^{A}$. If $v\in C^{l+1}$ in some neighborhood of a point $(t,x)\in[0,T]\times\mathbb{R}$, then for $1\leq l\leq k$,

[TABLE]

for $t^{n}-t$, $x_{j}-x$, $\Delta t$, $\Delta x$ sufficiently small.

Proof.

i) Let us proceed as has been done in [6] defining $w^{n+1}_{j}=S^{M}(w^{n})_{j}$, the solution computed with the monotone scheme alone with $w^{0}_{j}=v_{0}(x_{j})$. Then by definition,

[TABLE]

whence, exploiting the nonexpansivity in $L^{\infty}$ of $S^{M}$, the definition of $\varepsilon^{n}$ and that $|F|\leq 1$,

[TABLE]

Then, proceeding recursively on $n\leq N$ and recalling that by Prop. 4.1 there exists a constant $C>0$ such that $\varepsilon^{n}\leq C\Delta x:=\varepsilon$ for each $n$,

[TABLE]

At this point, by the triangular inequality

[TABLE]

whence we have that

[TABLE]

with $C_{CL}>0$ given by the Crandall-Lions estimate for $S^{M}$.

ii) Let us recall that by (4), in the case of $v\in C^{2}$ the consistency error for the monotone scheme is such that $\mathcal{E}_{M}(v^{n})_{j}\leq C_{M}(\Delta t+\Delta x)$. Then we can compute

[TABLE]

whence, by recursion on $n\leq N$ and recalling what we have done in the previous point,

[TABLE]

To finish this proof what is left is to use the estimate on $\mathcal{E}_{M}$ and Prop. 4.1 .

iii) In order to show that $S^{AF}(v^{n})_{j}=S^{A}(v^{n})_{j}$ for $\Delta t$ e $\Delta x$ small enough it is sufficient to prove that

[TABLE]

which follows directly from the computation we have done in Sect. 2.3 for the tuning of the parameter $\varepsilon^{n}$. In fact, if we plug (2.3) inside the previous inequality, we can deduce that

[TABLE]

which, using that $K>1/2$ by assumption, leads to the thesis as $(\Delta t,\Delta x)\to 0$. Notice that we have used the property $\varepsilon^{n}=O(\Delta x)$ and exploited the CFL condition. ∎

Remark 4.4**.**

Notice that the assumption $M(\Delta x)=\frac{1}{2}-C\Delta x$, for some constant $C>0$ such that $M(\Delta x)>0$, needed to apply Lemma 3.1, may give some problems in the proof of third assertion of the previous theorem. In fact, applying the standard definition (39) to the viscosity solution $v$ at a point $x_{j}$ and recalling the computations that led to (37), we get that

[TABLE]

Consequently, in order to be sure that if $v\in C^{3}$, then $j\in\mathcal{R}$, we have to choose the constant $C$ such that

[TABLE]

or require additional smoothness assumptions on $v$, for example $v^{\prime\prime\prime}_{j}\ll v^{\prime\prime}_{j}$. This in fact poses a strong limitation on the applicability of Lemma 3.1, at least in the present formulation.

5 Numerical Tests

In this section we will present some one-dimensional examples designed to show the properties of our scheme, stated by Theorem 4.3. Our goal is also to compare the performances of our Adaptive Filtered Schemes $S^{AF}$ with those of the Filtered Scheme $S^{F}$ introduced in [6] and of the WENO scheme of second/third order of [19]. Regarding the basic filtered scheme, we decided to avoid the introduction of the limiter used in [6] in all the numerical tests here presented for a more direct comparison. For all our numerical examples, we will use the function $\phi$ defined in (42), with $\beta_{k}$ given by (33)-(34), the mapping (48), and $M=0.1$, the parameter $\varepsilon^{n}$ defined in (2.3), and we will compute the errors and orders in $L^{\infty}$ and $L^{1}$ norm. For each test, we will specify the monotone and high-order schemes composing the filtered scheme. As already stated in Sect. 2.2, in all our numerical simulations we will use the discontinuous filter function defined in (22). This choice is justified by comparison reasons, since in [6] this is the only filter function used and we suppose the authors in [6] used it since it gives the best performances for their $S^{F}$ scheme. Since our scheme is not sensitive to the choice of the filter function, we use the same as in [6] for best comparisons. At the end of the section, we will also show briefly how to use these schemes in order to approximate simple two dimensional problems. To be precise, in the following examples we will refer to the standard CFL condition

[TABLE]

to define $\lambda$ , which is alternative to (6) and more easily computed.

All the numerical tests have been implemented in language C++, with plots generated by using MATLAB. The computer used for the simulations is a Notebook Asus F556U Intel Core i7-6500U with speed of 2.59 GHz and 12 GB of RAM.

Example 1: Transport equation. In order to test the capability of our scheme to handle both regular and singular regions, let us begin with a simple linear example and consider the problem

[TABLE]

with periodic boundary conditions, in two different situations. At first, aiming to test the full accuracy of the schemes, we consider the regular initial data (Case a),

[TABLE]

with $\Omega=[-2,2]$ and $T=0.9$. Then, as a second test, we take the mixed initial datum (Case b),

[TABLE]

with $\Omega=[-1.5,3.5]$ and $T=2$. The latter problem models the transport of a function composed by two peaks, the first with one point of singularity whereas the second is in $C^{2}$. For these tests we use the Central Upwind scheme (9) as monotone scheme and the simple Heun-Centered (HC) scheme (16)-(17) as high-order scheme, with $\lambda=0.9$ for Case a and $\lambda=0.4$ for Case b. We also compare the results obtained using $S^{AF}$ with the $4$th order Lax-Wendroff scheme (21) as high-order scheme. We recall that the latter high-order scheme has a very compact $5$-points stencil, whereas the WENO scheme of second/third order (coupled with the third order Runge Kutta scheme) has a stencil of nine points.

In the first case (Case a) of this test, all the schemes are very accurate and achieve optimal order in both norms, as shown in Tab. 1. In this case, both filtered schemes have the same numerical results, except for a slight difference with the coarstest grid, and coincide with the simple HC high-order scheme, as expected (we avoided to add another column in the table to report also the results for the HC high-order scheme since they are the same). Moreover, we can see that our fourth order scheme is much more accurate even than the WENO scheme, despite the smaller stencil required. In Fig. 2 we reported only the AF-HC scheme and the WENO scheme, avoiding to show all the schemes since no differences are visible for that case.

For the second case (Case b), looking at Fig. 3 we can observe that the adaptive tuning of $\varepsilon^{n}$ is able to contain the oscillations behind the peaks produced by the unstable HC scheme, which are clearly visible instead in the case of $S^{F}$ with $\varepsilon=10\Delta x$. We can also see that our scheme coupled with the fourth order scheme produces again almost always the best results in terms of errors and orders in both norms (see Tab. 2) and gives the best resolution of the peaks, preserving better the kink of the singularity and the feet of the regular part, without introducing any oscillation. In Tab. 3 we reported the CPU times for this Case b, in which the evolution lasts longer. All the schemes are very fast and complete the computations in less than $0.1$ s for all the refinements. Note that our two adaptive filtered schemes perform faster than the WENO scheme, even in the case of the fourth-order scheme. On the other hand, as could be expected, the adaptive procedure increases the cost of the filtering process three/four times depending on the refinement with respect to the basic filtered scheme.

Example 2: Eikonal equation. As a first nonlinear problem let us consider the eikonal equation

[TABLE]

where $v_{0}$ is a Lipschitz continuous initial datum with high regularity (Case a). Then, we repeat the simulation with the “reversed” initial datum (Case b)

[TABLE]

which presents also a major problem in the origin because of the saddle point in the hamiltonian, where two directions of propagation occur. Here the aim is mainly to compare the results obtained by the unfiltered high-order schemes with their filtered versions, in order to show the stabilization property of the filtering process. For the monotone scheme we use the numerical hamiltonian (8), whereas to achieve high-order we use the Lax-Wendroff-Richtmyer (LWR) scheme (20). Moreover, as in the previous example, we present also the results obtained with the AF scheme coupled with the fourth order LW scheme. The CFL number is set to $0.375$ for both simulations.

Let us first point out that, as Figs. 4 - 5 clearly show, the LWR scheme is unstable in the origin in both situations, whereas the AF scheme (and the simple filtered scheme) is stable. Then, for the first case, looking at Tab. 4 we can see that the filtered-LWR schemes give almost the same results, are of high-order in both norms and get lower errors with respect to the WENO scheme in almost all simulations. Moreover, we can recognize the typical improvements and drawbacks of the fourth order LW scheme, which has a slightly wider stencil. In fact, as will be shown also in the following examples, the scheme has bigger errors in the $L^{\infty}$ norm with respect to the second order AF scheme whereas has way better errors and orders in the $L^{1}$ norm, achieving almost optimal order, which testifies the overall improvement.

For Case b, looking at Tab. 5 we can repeat almost the same considerations made for Case a, but this time the improvements given by the adaptive filtering are more evident. The AF-LWR scheme is again of high-order especially in the $L^{1}$ norm, without the need to introduce any limiter as has been done in [6], and the numerical results are always comparable to those obtained by the WENO scheme of second/third order, whereas the AF-LW4ord scheme produces again bigger errors in $L^{\infty}$ with respect to the second-order AF-LWR scheme and better orders in the $L^{1}$ norm.

Example 3: Burgers’ equation. Let us consider now the Burgers’ equation for HJ with a regular initial datum

[TABLE]

which is a test case widely used in literature. In order to test the full accuracy of the schemes even in the nonlinear case, we first run the simulation for $T=\frac{4}{5\pi^{2}}$, when the solution is still regular, with $\lambda=\frac{2}{\pi^{2}}\approx 0.2<{\max|H_{p}|}^{-1}=0.5$. Then, we consider the final time $T=\frac{3}{2\pi^{2}}$ when a moving (to the right) singularity appears, taking $\lambda=\frac{15}{8\pi^{2}}\approx 0.19$. For both simulations we use the the Central Upwind monotone scheme and the LWR scheme for both the filtered schemes and compare the results as before with the WENO scheme and the fourth order AF scheme. In Fig. 6 we report the intial datum of the problem and the solution produced by the AF-LWR scheme at the two different times in order to show the different behavior.

This example summarizes all the behaviors already seen in the previous cases. In fact, as displayed by Tabs. 6-7, if the solution is still regular the fourth order AF scheme gives the best results and achieves the optimal order in both norms, whereas when the singularity appears, it gets bigger errors in $L^{\infty}$ norm but lower errors and better orders in the $L^{1}$ norm with respect to the second order filtered schemes. Here we have to notice that the WENO scheme has better errors and orders in the second simulation with respect to all the filtered schemes. Moreover, we can clearly see that the simple filtered scheme depends heavily on the choice of $\varepsilon$, in fact after extensive computations we noticed that choosing for example $\varepsilon=5\Delta x$ we get worse results in both cases, whereas if we increase the constant we get better results in the regular case and worse in the latter. In the tables we presented the results for the choice that gives the best results in the singular case, whereas it has clearly problems in the first situation. This is the main advantage of the adaptive $\varepsilon^{n}$ which is able to tune itself in the right way depending on the local (in time) regularity of the solution.

In order to give a visual evidence of that latter property, in Fig. 7 we reported the regions of activity of the schemes composing the two second order filtered schemes. There we can clearly see that our procedure is able to better localize the presence of the singularity, whereas when the solution is still regular the high-order scheme is always active. On the other hand, if we look at the computational times in Tab. 8, we can see that the basic filtered scheme is of course the fastest scheme, whereas the other three schemes have very similar CPU times, with the fourth order scheme performing slightly slower.

Example 4: Nonconvex Hamiltonian. In this example we consider a well known test case for nonconvex Hamiltonians (see e.g. [19]), that is

[TABLE]

with periodic boundary conditions and final time $T=3/(2\pi^{2})$, when two singularities appear in the solution, as can be seen in Fig. 8.

In order to define the monotone scheme for this test, we use the Lax-Friedrichs hamiltonian (10) with $\theta=1$ as in [6], whereas the CFL number is set to $0.31$. A reference solution is computed by using the AF-LW4ord scheme with $10240$ points.

We reported the global errors in Tab. 9 and the errors far away from singular points (regularity region) in Tab. 10. More precisely, for the second table we consider the set of grid points $x$ such that $|x-x_{i}|\geq 0.05$, for $i=1,2$, where $x_{1}=-0.895$ and $x_{2}=0.245$ are approximately the position of the singularities.

Looking at Tab. 9 we can note that all the tested schemes suffer a sort of “saturating” effect in $L^{\infty}$ norm, showing some difficulties in dropping the error in the last refinement, especially in the case of the WENO scheme. On the other hand, the high-order convergence rate in $L^{1}$ norm testifies the reliability of the schemes also in this situation. The best results in Tab. 9 are clearly given by the WENO scheme, nevertheless, the AF-HC scheme performs better in terms of error and orders in both norms with respect to its basic version with $\varepsilon=5\Delta x$, whereas the fourth order scheme presents the usual behavior in $L^{\infty}$ norm with respect to the second order AF-HC scheme, however maintaining the same order of errors. If instead we look at the errors in regions of regularity reported in Tab. 10, we can acknowledge that all the schemes achieve optimal order in both norms, with best results now given by the AF-LW4ord scheme, also with respect to the WENO scheme.

As already seen in the previous Example 3, in Fig. 9 we can recognize the ability of the AF scheme to better localize the regions of singularity with respect to the basic procedure which uses the monotone scheme way more than necessary.

Example 5: Evolution in 2D by dimensional splitting. We conclude this section on numerical simulations showing a convenient procedure to solve simple two-dimensional problems by making use of the one-dimensional schemes defined in the previous sections. Let us consider a classical problem similar to the Burgers’ equation, which is strictly connected to (65),

[TABLE]

with $\Omega=[0,2]^{2}$ and periodic boundary conditions. As done for problem (65), we consider the final time $T=\frac{4}{5\pi^{2}}$, when the solution is still smooth, and then $T=\frac{3}{2\pi^{2}}$, time at which an interesting set of singularities develops. The exact solution is computed by the Hopf-Lax formula,

[TABLE]

with $A=[-5,5]$.

In this situation, since the hamiltonian can be expressed as a sum of one-dimensional hamiltonian, depending on the evolution along the $x$ and $y$ direction, respectively, we can use a dimensional splitting to solve the problem. More precisely, if we write $H(v_{x},v_{y})=H_{1}(v_{x})+H_{2}(v_{y})$, we can approximate the solution by solving sequentially the problems in one space dimension

[TABLE]

keeping each time the other space variable constant. Since the hamiltonians trivially commute, we can use the simple Lie-Trotter splitting

[TABLE]

where $S^{\Delta t}_{x}$ and $S^{\Delta t}_{y}$ are numerical schemes of time step $\Delta t$ for the problems in the $x$ and $y$ direction, respectively, without introducing errors in the time evolution. For more details about dimensional splitting techniques we refer the reader to [23] and the references therein.

We use the same schemes as in Example 3 and a slightly more restrictive CFL number with respect to problem (65) in order to use coarser grids, which is set to $\lambda=\frac{4}{5\pi^{2}}\approx 0.08$ for the first test, and $\lambda=\frac{3}{4\pi^{2}}\approx 0.076$ for the latter.

As we could expect, in this example we have analogous result with respect to Example 3, with the AF scheme performing well in both situations and better than the F scheme in the regular case (see Tabs. 11-12). Here again the simple filtered scheme has slightly better results after the singularities develop, due to the action of the $\phi$ function in the regions of singularity, but the loss of accuracy is in fact minimal. Moreover, our scheme performs as good as the WENO scheme when the solution is still regular, whereas the latter performs better in the second case. Concerning the CPU times, looking at Tab. 13 we can see that, differently from the one-dimensional case, the WENO scheme is faster than the AF scheme in the last two refinements, whereas the F-LWR scheme is comparable to the simple LWR high-order scheme and both are three/four times faster than the other two considered schemes.

6 Conclusions

We have presented a rather simple way to construct convergent schemes coupling a monotone and a high-order scheme via a filter function. A typical feature of filtered schemes is their high-order accuracy in the regions of regularity for the solution. In fact, the filter function can stabilize an otherwise unstable (high-order) scheme, still preserving its accuracy. The main novelty here is the adaptive and automatic choice of the parameter $\varepsilon^{n}$ which improves the scheme in [6]. The computation of the switching parameter $\varepsilon^{n}$, although more expensive, is still affordable in low dimension. The adaptive scheme is able to reduce the oscillations which may appear choosing a constant $\varepsilon$ and, as shown by the numerical tests, gives always better results. Finally, we note that the accuracy of adaptive filtered schemes is close to WENO schemes of the same order but filtered schemes are easier to implement, give a rather flexible way to couple different schemes and, as we proved, converge to the viscosity solution. A fully 2D scheme has been used in our recent paper [14], in which new formulas for the 2D smoothness indicators are proposed, instead of applying a splitting techniques, for the resolution of the image segmentation problem. Thanks to recent computations, we are able to prove a result similar to Prop. 3.1 for our new 2D smoothness indicators. Thanks to such a result, we believe that the convergence and high-order consistency could be proven also in fully two-dimensional problems, but this is a future work.

Acknowledgements

We would like to thank the national group INdAM-GNCS for the financial support given to this research.

Appendix A APPENDIX: Technical results

For completeness and reader’s convenience we give the proofs of Prop. 3.1 and of the properties of the undivided differences and the binomial coefficients involved. This analysis follows the ideas in [1] where a similar analysis is developed for conservation laws.

of Proposition 3.1.

Let us take $r>1$ and without loss of generality, let $x_{s}=0$ (to simplify the notation). Moreover, we will use the convention $h:=\Delta x$. Let us start by reminding that, using the Newton form of the interpolating polynomial, for $k=0,\dots,r-1$ and $j\in\mathbb{Z}$, we get

[TABLE]

where $\omega_{i}(x)=(x-x_{j-r+k})\cdots(x-x_{j-r+k+i})$ and $f[\cdot]$ denotes the divided difference of $f$.

We proceed with the proof of i). In this case it is sufficient to observe that, since the function $f$ is regular in $\Omega\setminus\{x_{s}\}$, the properties of the interpolating polynomial directly give

[TABLE]

Moreover, expanding with Taylor, it holds

[TABLE]

where $m=\max\{s+1,l\}$ and $s=\max\{k:f^{(i)}(x_{j})=0,\ \forall i\leq k\}$ ($s\leq r$). Then, integrating (remembering that by hypothesis $s=0\Rightarrow m=l$), we get

[TABLE]

as we wanted.

Let us continue with the proof of ii). In this case the proof is a little more complicated and it is better to treat separately the following two cases:

1. (a)

[math] is a point of the grid $\{x_{i}\}$, $i\in\mathbb{Z}$; 2. (b)

$0\in I_{i}=(x_{i-1},x_{i})$ for some $i\in\mathbb{Z}$.

Case a. By hypothesis $0\in\mathcal{S}_{j+k}$ for at least one $k=0,\dots,r-1$, then, for each fixed $k$, there exists an integer $j_{s}\in\{k-r+1,\dots,k-1\}$ such that $x_{j}=-j_{s}h$ (for $j_{s}=k-r$ and $j_{s}=k$ we fall in the case treated previously). Substituting in (68),

[TABLE]

with $\omega_{i}(x)=(x+(j_{s}-k+r)h)\cdots(x+(j_{s}-k+r-i)h)$. Moreover, if we define the function $f_{h}(y):=f(x_{j}+hy)=f(h(y-j_{s}))$, we can write

[TABLE]

where $f_{h}[\cdot]$ denotes the undivided difference of $f_{h}$. Now, defining the polynomial

[TABLE]

where $q_{i}(y)=(y-(k-r))\cdots(y-(k-r-i))$, $f_{h}(y)=f(x_{j}+hy)$, we can rewrite

[TABLE]

Then, applying the change of variable $y=(x-x_{j})/h$ in the integral in (31), we have

[TABLE]

where

[TABLE]

At this point, it is useful to notice that (75) for $l=1$ reads, for $i=1,\dots,r$,

[TABLE]

In order to simplify the notation let us call $i_{s}:=j_{s}-k+r$, that is to say the index $i_{s}\in\{1,\dots,r-1\}$ such that $x_{j}+(k-r+i_{s})h=0$. Then, by hypothesis, we can write for all $i>t:=\max\{i_{s},l-1\}$,

[TABLE]

and, noticing that for $h\to 0$

[TABLE]

we can conclude that

[TABLE]

having exploited the relations $\sum_{j=0}^{i}\binom{i}{j}(-1)^{i-j}=0$ and $\sum_{j=0}^{n}\binom{i}{j}(-1)^{i-j}=\binom{i-1}{n}(-1)^{i-n}$, for $0\leq n<i$ by Lemma A.1. Furthermore for $l\leq i\leq i_{s}$, using the relation (76), we can conclude

[TABLE]

From what we have done so far we can deduce, recalling that $t:=\max\{i_{s},l-1\}$,

[TABLE]

where $C_{rk}=\int_{-1}^{0}\left(\sum_{i=t+1}^{r}\binom{i-2}{i_{s}-1}(-1)^{i-i_{s}+1}\frac{q_{i-1}^{(l)}}{i!}\right)^{2}$, which is the thesis for Case a.

Case b. By hypothesis there exists an integer $j_{s}\in\{k-r+1,\dots,k\}$ and a number $0<a_{s}<1$ such that $x_{j}=(-j_{s}+a_{s})h$. It is clear now that we can repeat the same constructions of the previous case defining the function $f_{a_{s},h}(y):=f(h(y-j_{s}+a_{s}))$ and using it in place of $f_{h}$; so, to obtain (72) it will suffice to apply the change of variables $y=\frac{x}{h}+j_{s}-a_{s}$. Then, naming $i_{s}=j_{s}-k+r$, for $i\geq t:=\max\{i_{s},l\}$,

[TABLE]

whence, noticing that

[TABLE]

and that

[TABLE]

we can infer that if $i=i_{s}$ (in this case in (74) on the right side of the equation we have only the second term), then $f_{a_{s},h}[k-r,\dots,(k-r+i)]\approx a_{s}h\left[f^{\prime}(0^{+})-f^{\prime}(0^{-})\right]\not=0$, whereas if $i>i_{s}$,

[TABLE]

The last quantity, as it can be easily shown, it is null if and only if $a_{s}=\frac{i-i_{s}}{i-1}$; more precisely, for $k$ fixed there exists an integer $i\geq t$ such that $f_{a_{s},h}[k-r,\dots,(k-r+i)]\approx Ch\left[f^{\prime}(0^{+})-f^{\prime}(0^{-})\right]$ with $C\not=0$, whence the thesis even in the last case. ∎

Lemma A.1**.**

Let us assume $i\geq 1$ and write $f[\cdot]$ for the undivided difference of a function $f$. Then, it holds

[TABLE]

Moreover, we have that

[TABLE]

Proof.

Let us start from the proof of (75) and let us proceed by induction on $i$.

Firstly, let us notice that for $l=i$ the identity is trivially satisfied, whence the case $i=1$ follows directly. Then, for any $l=0,\dots,i-1$, suppose that the statement holds for $i-1$ and for $i>0$ let us compute,

[TABLE]

as we wanted.

Remark A.1**.**

To simplify the notation we have stated the result for $f[0,\dots,i]$ but the proof clearly holds for $f[k,\dots,k+i]$, $\forall k$. In the second identity of the previous chain we have assumed this fact applying the inductive hypothesis on both terms.

Let us focus now on the second relation of the lemma (76) and proceed again by induction, but this time on $n:0\leq n<i$. For $n=0$ we have $(-1)^{i}=(-1)^{i}$, then the identity holds. Suppose that (76) holds for $n-1<i-1$ and compute

[TABLE]

For $n=i$ instead, from what we have just seen we can easily compute

[TABLE]

∎