Analysis of the spectral symbol function for spectral approximation of a differential operator

TL;DR

This paper investigates the spectral symbol function for discretized differential operators, establishing conditions under which it accurately approximates the spectrum and challenging the effectiveness of uniform sampling methods.

Contribution

It proves that the spectral symbol provides a necessary condition for spectral approximation and shows that uniform sampling of the symbol is generally insufficient for accurate spectrum approximation.

Findings

Spectral symbol is necessary for spectral approximation.

Uniform sampling of spectral symbol often fails to approximate spectra accurately.

Proper discretization and grid refinement improve spectral approximation.

Abstract

Given a differential operator $\mathcal{L}$ along with its own eigenvalue problem $\mathcal{L}u = \lambda u$ and an associated algebraic equation $\mathcal{L}^{(n)} \mathbf{u}_n = \lambda\mathbf{u}_n$ obtained by means of a discretization scheme (like Finite Differences, Finite Elements, Galerkin Isogeometric Analysis, etc.), the theory of Generalized Locally Toeplitz (GLT) sequences serves the purpose to compute the spectral symbol function $\omega$ associated to the discrete operator $\mathcal{L}^{(n)}$ We prove that the spectral symbol $\omega$ provides a necessary condition for a discretization scheme in order to uniformly approximate the spectrum of the original differential operator $\mathcal{L}$. The condition measures how far the method is from a uniform relative approximation of the spectrum of $\mathcal{L}$. Moreover, the condition seems to become sufficient if the discretization method is paired with a suitable (non-uniform) grid and an increasing refinement of the order of approximation of the method. On the other hand, despite the numerical experiments in many recent literature, we disprove that in general a uniform sampling of the spectral symbol $\omega$ can provide an accurate relative approximation of the spectrum, neither of $\mathcal{L}$ nor of the discrete operator $\mathcal{L}^{(n)}$.

Tables

Table 1. Table 1: for α → 0 → 𝛼 0 \alpha\to 0 , we confront the sup-norm of the difference between the monotone rearrangment ω ~ α subscript ~ 𝜔 𝛼 \prescript{}{\alpha}{\tilde{\omega}} and the spectral symbol ω ​ ( θ ) = 2 − 2 ​ cos ⁡ ( θ ​ π ) 𝜔 𝜃 2 2 𝜃 𝜋 \omega(\theta)=2-2\cos(\theta\pi) , which characterizes the 3-point FD scheme discretization of the 1 ​ d 1 𝑑 1d Dirichlet Laplacian over [ 0 , 1 ] 0 1 [0,1] . The pointwise evaluation of ω ~ α subscript ~ 𝜔 𝛼 \prescript{}{\alpha}{\tilde{\omega}} and ω 𝜔 \omega is made over the uniform grid { k n + 1 } k = 1 n superscript subscript 𝑘 𝑛 1 𝑘 1 𝑛 \left\{\frac{k}{n+1}\right\}_{k=1}^{n} with n = 10 3 𝑛 superscript 10 3 n=10^{3} , whereas the approximation of ω ~ α subscript ~ 𝜔 𝛼 \prescript{}{\alpha}{\tilde{\omega}} is obtained evaluating ϕ − 1 α subscript superscript italic-ϕ 1 𝛼 \prescript{}{\alpha}{\phi^{-1}} from ( 4.5 ) by means of the fzero() function from Matlab r2018b.

	$\alpha =1$	$\alpha =10e-02$	$\alpha =10e-05$	$\alpha =10e-10$
${\Vert {}_{\alpha }{}^{}\tilde{\omega }-\omega \Vert }_{\infty }$	5.8049	0.3912	0.0103	4.2392e-06

Table 2. Table 2: For every fixed k 𝑘 k and α 𝛼 \alpha , the analytic relative error err ~ k ( n ) α subscript superscript subscript ~ err 𝑘 𝑛 𝛼 \prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{k}^{(n)}} converges to the lower bound c α , k subscript 𝑐 𝛼 𝑘 c_{\alpha,k} as n 𝑛 n increases, where c α , k subscript 𝑐 𝛼 𝑘 c_{\alpha,k} is given in ( 4.7 ). Observe that err ~ k ( n ) α subscript superscript subscript ~ err 𝑘 𝑛 𝛼 \prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{k}^{(n)}} seems to be monotone decreasing of order O ​ ( n − 2 ) 𝑂 superscript 𝑛 2 O(n^{-2}) . The approximation of ω ~ α subscript ~ 𝜔 𝛼 \prescript{}{\alpha}{\tilde{\omega}} is obtained evaluating ϕ − 1 α subscript superscript italic-ϕ 1 𝛼 \prescript{}{\alpha}{\phi^{-1}} from ( 4.5 ) by means of the fzero() function from Matlab r2018b.

		$\left|{\widetilde{\text{err}}}_k^{\left(n\right)}/c_{\alpha ,k}-1\right|$
		$n={10}^2$	$n={10}^3$	$n={10}^4$
$\alpha =0.1$	$c_{\alpha ,1}=0.0025$	0.0326	3.3223e-04	3.3283e-06
	$c_{\alpha ,5}=1.0131e-04$	20.3811	0.2076	0.0021
	$c_{\alpha ,10}=2.5330e-05$	325.3811	3.3222	0.0333
$\alpha =1$	$c_{\alpha ,1}=0.0247$	0.0041	4.1363e-05	4.1438e-07
	$c_{\alpha ,5}=0.0010$	2.5395	0.0259	2.5899e-04
	$c_{\alpha ,10}=2.5324e-04$	40.6422	0.4136	0.0041
$\alpha =2$	$c_{\alpha ,1}=0.0482$	0.0026	2.6120e-05	2.6167e-07
	$c_{\alpha ,5}=0.0020$	1.6056	0.0163	1.6354e-04
	$c_{\alpha ,10}=5.0635e-04$	25.7979	0.2612	0.0026
$\alpha =5$	$c_{\alpha ,1}=0.1124$	0.0020	2.0008e-05	2.0044e-07
	$c_{\alpha ,5}=0.0050$	1.2389	0.0125	1.2528e-04
	$c_{\alpha ,10}=0.0013$	20.4017	0.2002	1.2528e-04

Table 3. Table 3: In this table it is possible to check that for every fixed α 𝛼 \alpha , the n 𝑛 n -th term of the analytic relative error err ~ r , n ( n ) α subscript superscript subscript ~ err 𝑟 𝑛 𝑛 𝛼 \prescript{}{\alpha}{\tilde{\textnormal{err}}_{r,n}^{(n)}} converges to zero as n 𝑛 n increases, and that the quotient between ( n + 1 ) 2 ​ λ n ​ ( ℒ ^ dir , α ​ x 2 ( n ) 1 e α ) superscript 𝑛 1 2 subscript 𝜆 𝑛 superscript subscript superscript subscript ^ ℒ dir 𝛼 superscript 𝑥 2 𝑛 1 superscript e 𝛼 (n+1)^{2}\lambda_{n}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\hat{\mathcal{L}}_{\textnormal{dir},\alpha x^{2}}^{(n)}}\right) and λ n ​ ( ℒ dir , α ​ x 2 1 e α ) subscript 𝜆 𝑛 superscript subscript subscript ℒ dir 𝛼 superscript 𝑥 2 1 superscript e 𝛼 \lambda_{n}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right) converges to max ( x , θ ) ∈ [ 1 , e α ] × [ 0 , π ] ⁡ ω α ​ ( x , θ ) / π 2 = ω ~ α ​ ( 1 ) / π 2 subscript 𝑥 𝜃 1 superscript e 𝛼 0 𝜋 subscript 𝜔 𝛼 𝑥 𝜃 superscript 𝜋 2 subscript ~ 𝜔 𝛼 1 superscript 𝜋 2 \max_{(x,\theta)\in[1,\textrm{e}^{\sqrt{\alpha}}]\times[0,\pi]}\prescript{}{\alpha}{\omega(x,\theta)}/\pi^{2}=\prescript{}{\alpha}{\tilde{\omega}(1)}/\pi^{2} . In this case we choose r = n 𝑟 𝑛 r=n to keep computational costs at minimum.

		$n={10}^2$	$n={10}^3$	$n=3\cdot {10}^3$
$\alpha =0.5$	${}_{\alpha }{}^{}\widetilde{\text{err}}_{r,n}^{\left(n\right)}$	0.0210	0.0036	0.0016
	$\left|\frac{{\left(n+1\right)}^2{\lambda }_n^{\left(n\right)}/{\lambda }_n}{{}_{\alpha }{}^{}\tilde{\omega }\left(1\right)/{\pi }^2}-1\right|$	0.0482	0.0128	0.0065
$\alpha =1.2$	${}_{\alpha }{}^{}\widetilde{\text{err}}_{r,n}^{\left(n\right)}$	0.0249	0.0044	0.0020
	$\left|\frac{{\left(n+1\right)}^2{\lambda }_n^{\left(n\right)}/{\lambda }_n}{{}_{\alpha }{}^{}\tilde{\omega }\left(1\right)/{\pi }^2}-1\right|$	0.0615	0.0157	0.0079
$\alpha =4$	${}_{\alpha }{}^{}\widetilde{\text{err}}_{r,n}^{\left(n\right)}$	0.0297	0.0052	0.0024
	$\left|\frac{{\left(n+1\right)}^2{\lambda }_n^{\left(n\right)}/{\lambda }_n}{{}_{\alpha }{}^{}\tilde{\omega }\left(1\right)/{\pi }^2}-1\right|$	0.0765	0.0191	0.0095

Table 4. Table 4: In this table we check numerically the validity of Theorem 5.1.1 for different values of α 𝛼 \alpha and n 𝑛 n . It can be seen that for every α 𝛼 \alpha , as n 𝑛 n increases then the relative error between max k = 1 , … , n ⁡ | λ k ​ ( ℒ dir , α ​ x 2 ( n ) 1 e α ) / λ k ​ ( ℒ dir , α ​ x 2 1 e α ) | subscript 𝑘 1 … 𝑛 subscript 𝜆 𝑘 superscript subscript superscript subscript ℒ dir 𝛼 superscript 𝑥 2 𝑛 1 superscript e 𝛼 subscript 𝜆 𝑘 superscript subscript subscript ℒ dir 𝛼 superscript 𝑥 2 1 superscript e 𝛼 \max_{k=1,\ldots,n}\left|\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n)}}\right)/\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right)\right| and max x ∈ [ 0 , 1 ] ⁡ | ω ~ α ​ ( x ) / x 2 ​ π 2 | subscript 𝑥 0 1 subscript ~ 𝜔 𝛼 𝑥 superscript 𝑥 2 superscript 𝜋 2 \max_{x\in[0,1]}\left|\prescript{}{\alpha}{\tilde{\omega}}(x)/x^{2}\pi^{2}\right| decreases, confirming ( 4.13 ). In the table is reported as well the quotient k ¯ / n ¯ 𝑘 𝑛 \bar{k}/n , where k ¯ ¯ 𝑘 \bar{k} is the k 𝑘 k -th eigenvalue which achieves the maximum relative error between λ k ​ ( ℒ dir , α ​ x 2 ( n ) 1 e α ) subscript 𝜆 𝑘 superscript subscript superscript subscript ℒ dir 𝛼 superscript 𝑥 2 𝑛 1 superscript e 𝛼 \lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n)}}\right) and λ k ​ ( ℒ dir , α ​ x 2 1 e α ) subscript 𝜆 𝑘 superscript subscript subscript ℒ dir 𝛼 superscript 𝑥 2 1 superscript e 𝛼 \lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right) . We can notice that k ¯ / n ¯ 𝑘 𝑛 \bar{k}/n is always bounded and it tends to a finite value in ( 0 , 1 ] 0 1 (0,1] as n 𝑛 n increases. The approximation of ω ~ α subscript ~ 𝜔 𝛼 \prescript{}{\alpha}{\tilde{\omega}} is obtained evaluating ϕ − 1 α subscript superscript italic-ϕ 1 𝛼 \prescript{}{\alpha}{\phi^{-1}} from ( 4.5 ) by means of the fzero() function from Matlab r2018b.

		$n={10}^2$	$n={10}^3$	$n=5\cdot {10}^3$
$\alpha =0.5$	$\left|\frac{\max \left|{\lambda }_k^{\left(n\right)}/{\lambda }_k\right|}{\max \left|{}_{\alpha }{}^{}\tilde{\omega }\left(x\right)/x^2{\pi }^2\right|}-1\right|$	0.0104	0.0010	2.0853e-04
	$\overline{k}/n$	0.7900	0.7880	0.7878
$\alpha =1$	$\left|\frac{\max \left|{\lambda }_k^{\left(n\right)}/{\lambda }_k\right|}{\max \left|{}_{\alpha }{}^{}\tilde{\omega }\left(x\right)/x^2{\pi }^2\right|}-1\right|$	0.0158	0.0016	3.1754e-04
	$\overline{k}/n$	0.6700	0.6680	0.6676
$\alpha =1.2$	$\left|\frac{\max \left|{\lambda }_k^{\left(n\right)}/{\lambda }_k\right|}{\max \left|{}_{\alpha }{}^{}\tilde{\omega }\left(x\right)/x^2{\pi }^2\right|}-1\right|$	0.0180	0.0018	3.6226e-04
	$\overline{k}/n$	0.64	0.6310	0.6302
$\alpha =3$	$\left|\frac{\max \left|{\lambda }_k^{\left(n\right)}/{\lambda }_k\right|}{\max \left|{}_{\alpha }{}^{}\tilde{\omega }\left(x\right)/x^2{\pi }^2\right|}-1\right|$	0.0518	0.0097	0.0032
	$\overline{k}/n$	1	1	1

Table 5. Table 5: Comparison between the maximum of the eigenvalues relative errors of the discrete differential operators ℒ dir , α ​ τ ​ ( x ) 2 ( n , η ) 1 e α superscript subscript superscript subscript ℒ dir 𝛼 𝜏 superscript 𝑥 2 𝑛 𝜂 1 superscript e 𝛼 \prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha\tau(x)^{2}}^{(n,\eta)}} obtained by means of ( 2 ​ η + 1 ) 2 𝜂 1 (2\eta+1) -points central FD discretization on uniform and non-uniform grids, for different values of η 𝜂 \eta . The parameter α 𝛼 \alpha is fixed, with α = 1 𝛼 1 \alpha=1 . We observe that in the uniform grid case, as η 𝜂 \eta increases the maximum increases as well. On the contrary, in the non-uniform grid case given by ( 4.14 ), the maximum decreases significantly as η 𝜂 \eta increases. See Figure 8 for a general overview of the error distribution.

	${\max }_{k=1,\dots ,n}\left|\frac{{\lambda }_k^{(n,\eta )}}{{\lambda }_k}-1\right|$
	$\eta =1$	$\eta =10$	$\eta =15$
	$n={10}^3$	$n={10}^3$	$n={10}^3$
uniform grid	0.3201	0.9057	1.0000
non-uniform grid	0.5939	0.2210	0.1814

Table 6. Table 6: In this table we check numerically the validity of Theorem 5.1.1 for different values of α 𝛼 \alpha and n 𝑛 n . The discretization has been made by means of central FD of order η = 5 𝜂 5 \eta=5 . It can be seen that for every α 𝛼 \alpha , as n 𝑛 n increases then the relative error between max k = 1 , … , n ⁡ | λ k ​ ( ℒ dir , α ​ τ ​ ( x ) 2 ( n , η ) 1 e α ) / λ k ​ ( ℒ dir , α ​ x 2 1 e α ) | subscript 𝑘 1 … 𝑛 subscript 𝜆 𝑘 superscript subscript superscript subscript ℒ dir 𝛼 𝜏 superscript 𝑥 2 𝑛 𝜂 1 superscript e 𝛼 subscript 𝜆 𝑘 superscript subscript subscript ℒ dir 𝛼 superscript 𝑥 2 1 superscript e 𝛼 \max_{k=1,\ldots,n}\left|\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha\tau(x)^{2}}^{(n,\eta)}}\right)/\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right)\right| and max x ∈ [ 0 , 1 ] ⁡ | ω ~ η , r α ​ ( x ) / x 2 ​ π 2 | subscript 𝑥 0 1 subscript subscript ~ 𝜔 𝜂 𝑟 𝛼 𝑥 superscript 𝑥 2 superscript 𝜋 2 \max_{x\in[0,1]}\left|\prescript{}{\alpha}{\tilde{\omega}_{\eta,r}}(x)/x^{2}\pi^{2}\right| decreases, confirming ( 4.13 ). In the table is reported as well the quotient k ¯ / n ¯ 𝑘 𝑛 \bar{k}/n , where k ¯ ¯ 𝑘 \bar{k} is the k 𝑘 k -th eigenvalue which achieves the maximum relative error between λ k ​ ( ℒ dir , α ​ τ ​ ( x ) 2 ( n , η ) 1 e α ) subscript 𝜆 𝑘 superscript subscript superscript subscript ℒ dir 𝛼 𝜏 superscript 𝑥 2 𝑛 𝜂 1 superscript e 𝛼 \lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha\tau(x)^{2}}^{(n,\eta)}}\right) and λ k ​ ( ℒ dir , α ​ x 2 1 e α ) subscript 𝜆 𝑘 superscript subscript subscript ℒ dir 𝛼 superscript 𝑥 2 1 superscript e 𝛼 \lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right) . We can notice that k ¯ / n ¯ 𝑘 𝑛 \bar{k}/n is always bounded and it tends to a finite value in ( 0 , 1 ] 0 1 (0,1] as n 𝑛 n increases. The approximation of ω ~ η , r α subscript subscript ~ 𝜔 𝜂 𝑟 𝛼 \prescript{}{\alpha}{\tilde{\omega}_{\eta,r}} is obtained by means of Algorithm I with r 𝑟 r fixed, r = n 𝑟 𝑛 r=n .

		Uniform grid $\eta =5$		Nonuniform grid $\eta =5$
		$n={10}^2$	$n={10}^3$	$n={10}^2$	$n={10}^3$
$\alpha =0.1$	$\left|\frac{\max \left|{\lambda }_k^{(n,\eta )}/{\lambda }_k\right|}{\max \left|{}_{\alpha }{}^{}\tilde{\omega }_{\eta ,r}^{}\left(x\right)/x^2{\pi }^2\right|}-1\right|$	0.0369	0.0018	2.7821e-04	2.7510e-06
	$\overline{k}/n$	0.9100	0.9050	1	1
$\alpha =3$	$\left|\frac{\max \left|{\lambda }_k^{(n,\eta )}/{\lambda }_k\right|}{\max \left|{}_{\alpha }{}^{}\tilde{\omega }_{\eta ,r}^{}\left(x\right)/x^2{\pi }^2\right|}-1\right|$	0.0628	0.0096	0.0086	8.2742e-05
	$\overline{k}/n$	1	1	1	1

Table 7. Table 7: Comparison between the maximum of the eigenvalues relative errors of the discrete differential operators ℒ dir , α ​ τ ​ ( x ) 2 ( n , η ) 1 e α superscript subscript superscript subscript ℒ dir 𝛼 𝜏 superscript 𝑥 2 𝑛 𝜂 1 superscript e 𝛼 \prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha\tau(x)^{2}}^{(n,\eta)}} obtained by means of IgA discretization of order η 𝜂 \eta on uniform and non-uniform grids, for different values of η 𝜂 \eta . The parameter α 𝛼 \alpha is fixed, with α = 1 𝛼 1 \alpha=1 . We observe that in the non-uniform grid case given by ( 4.14 ), the maximum decreases significantly as η 𝜂 \eta increases. See Figure 11 for a general overview of the error distribution. Let us notice that we did not take in consideration the outliers, see Definition 3.1.5 .

	${\max }_{k=1,\dots ,n}\left|\frac{{\lambda }_k^{(n,\eta )}}{{\lambda }_k}-1\right|$
	$\eta =1$	$\eta =5$	$\eta =10$
	$n={10}^2$	$n={10}^2$	$n={10}^2$
uniform grid	1.7653	1.0688	1.1465
non-uniform grid	0.4433	0.0483	0.0265

Table 8. Table 8: In this table we check numerically the validity of Theorem 5.1.1 for different values of α 𝛼 \alpha and n 𝑛 n . The discretization has been made by means of IgA of order η = 4 𝜂 4 \eta=4 . It can be seen that for every α 𝛼 \alpha , as n 𝑛 n increases then the relative error between max k = 1 , … , n ⁡ | λ k ​ ( ℒ dir , α ​ x 2 ( n + η − 1 , η ) 1 e α ) / λ k ​ ( ℒ dir , α ​ x 2 1 e α ) | subscript 𝑘 1 … 𝑛 subscript 𝜆 𝑘 superscript subscript superscript subscript ℒ dir 𝛼 superscript 𝑥 2 𝑛 𝜂 1 𝜂 1 superscript e 𝛼 subscript 𝜆 𝑘 superscript subscript subscript ℒ dir 𝛼 superscript 𝑥 2 1 superscript e 𝛼 \max_{k=1,\ldots,n}\left|\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n+\eta-1,\eta)}}\right)/\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right)\right| and max x ∈ [ 0 , 1 ] ⁡ | ω ~ η , r α ​ ( x ) / x 2 ​ π 2 | subscript 𝑥 0 1 subscript subscript ~ 𝜔 𝜂 𝑟 𝛼 𝑥 superscript 𝑥 2 superscript 𝜋 2 \max_{x\in[0,1]}\left|\prescript{}{\alpha}{\tilde{\omega}_{\eta,r}}(x)/x^{2}\pi^{2}\right| decreases, confirming ( 4.13 ). In the table is reported as well the quotient k ¯ / n ¯ 𝑘 𝑛 \bar{k}/n , where k ¯ ¯ 𝑘 \bar{k} is the k 𝑘 k -th eigenvalue which achieves the maximum relative error between λ k ​ ( ℒ dir , α ​ x 2 ( n + η − 1 , η ) 1 e α ) subscript 𝜆 𝑘 superscript subscript superscript subscript ℒ dir 𝛼 superscript 𝑥 2 𝑛 𝜂 1 𝜂 1 superscript e 𝛼 \lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n+\eta-1,\eta)}}\right) and λ k ​ ( ℒ dir , α ​ x 2 1 e α ) subscript 𝜆 𝑘 superscript subscript subscript ℒ dir 𝛼 superscript 𝑥 2 1 superscript e 𝛼 \lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right) . We can notice that k ¯ / n ¯ 𝑘 𝑛 \bar{k}/n is always bounded and it tends to a finite value in ( 0 , 1 ] 0 1 (0,1] as n 𝑛 n increases. The approximation of ω ~ η , r α subscript subscript ~ 𝜔 𝜂 𝑟 𝛼 \prescript{}{\alpha}{\tilde{\omega}_{\eta,r}} is obtained by means of Algorithm I with r 𝑟 r fixed, r = n 𝑟 𝑛 r=n . Let us notice that we did not take in consideration the outliers, see Definition 3.1.5 .

		Uniform grid $\eta =4$		Nonuniform grid $\eta =4$
		$n={10}^2$	$n={10}^3$	$n={10}^2$	$n={10}^3$
$\alpha =0.1$	$\left|\frac{\max \left|{\lambda }_k^{(n,\eta )}/{\lambda }_k\right|}{\max \left|{}_{\alpha }{}^{}\tilde{\omega }_{\eta ,r}^{}\left(x\right)/x^2{\pi }^2\right|}-1\right|$	0.3397	0.0881	0.0016	1.7241e-04
	$\overline{k}/n$	1	1	0.9400	0.9300
$\alpha =3$	$\left|\frac{\max \left|{\lambda }_k^{(n,\eta )}/{\lambda }_k\right|}{\max \left|{}_{\alpha }{}^{}\tilde{\omega }_{\eta ,r}^{}\left(x\right)/x^2{\pi }^2\right|}-1\right|$	0.1735	0.0425	0.0016	1.7243e-04
	$\overline{k}/n$	1	1	0.9400	0.9300

Table 9. Table 9: Comparison between an equispaced sampling of the (approximated) monotone rearrangement ω ~ r subscript ~ 𝜔 𝑟 \tilde{\omega}_{r} and the eigenvalues λ k subscript 𝜆 𝑘 \lambda_{k} of the discrete differential operator of Problem ( 4.17 ). ω ~ r subscript ~ 𝜔 𝑟 \tilde{\omega}_{r} is computed according to Algorithm I with r = n 𝑟 𝑛 r=n . In the first row it is calculated the absolute error: it increases as n 𝑛 n increases as stated in ( 4.21 ). In the second row it is calculated the maximum of the analytic relative error: it saturates at a lower bound c > 0 𝑐 0 c>0 . In the third and fourth row are validated the estimates ( 4.18 ) and ( 4.19 ), respectively. In the sixth row both the left-hand side and the right-hand side of Equation ( 4.20 ) are compared.

	$r=n={10}^2$	$r=n={10}^3$	$r=n=2\cdot {10}^3$
${\Vert {\tilde{\omega }}_{r,k}^{\left(n\right)}-{\lambda }_k\Vert }_{\infty }$	11.7523	37.0283	52.3532
${\max }_{k=1,\dots ,n}\left\{{\widetilde{\text{err}}}_{r,k}^{\left(n\right)}\right\}$	0.4133	0.4136	0.4136
${\tilde{\omega }}_{r,n}^{\left(n\right)}/\sqrt{n+1}$	3.9990	4	4
${\lambda }_n/\sqrt{n+1}$	2.8296	2.8296	2.8296
$n^{-1}{\sum }_{k=1}^n{\lambda }_k$	3.7663	3.9200	3.9428
\hdashline[1pt/1pt] ${\int }_0^1\tilde{\omega }\left(x\right)𝑑{\mu }_1\left(x\right)$	4	4	4

Equations

$$\{\begin{array}{ll}-{\partial }_x\left(p(x){\partial }_{x}u(x)\right)+q(x)u(x)=f(x)& x\in (a,b),\\ \text{boundary conditions},& \end{array}$$

$$\{\begin{array}{ll}-{\partial }_x\left(p(x){\partial }_{x}u(x)\right)+q(x)u(x)=\lambda w(x)u(x)& x\in (a,b),\\ \text{boundary conditions}.& \end{array}$$

$${\mathcal{L}}^{(n)}{\mathbf{u}}_n={\lambda }^{(n)}{\mathbf{u}}_n,$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{k=1}^{n}F\left({\lambda }_k\left({\hat{\mathcal{L}}}^{(n)}\right)\right)=\frac{1}{{\mu }_m(D)}{\int }_{D}F(\omega (\mathbf{x}))d{\mu }_m(\mathbf{x}),$$

$$\text{prescript}ba{\hat{\mathcal{L}}}_{\text{dir},p(x),q(x),w(x)}^{(n,\eta )}.$$

$$\{\begin{array}{ll}-{\partial }_x\left(p(x){\partial }_{x}u(x)\right)+q(x)u(x)=0& x\in (a,b)\subset \mathbb{R},\\ {\sigma }_{1}u(a)-{\sigma }_{2}p(a){\partial }_{x}u(x{)}_{|x=a}=0& {\sigma }_1^2+{\sigma }_2^2>0,\\ {\zeta }_{1}u(b)+{\zeta }_{2}p(b){\partial }_{x}u(x{)}_{|x=b}=0& {\zeta }_1^2+{\zeta }_2^2>0,\end{array}$$

$$\{\begin{array}{ll}-{\partial }_x\left(p(x){\partial }_{x}u(x)\right)+q(x)u(x)=\lambda w(x)u(x)& x\in (a,b)\subset \mathbb{R},\\ {\sigma }_{1}u(a)-{\sigma }_{2}p(a){\partial }_{x}u(x{)}_{|x=a}=0& {\sigma }_1^2+{\sigma }_2^2>0,\\ {\zeta }_{1}u(b)+{\zeta }_{2}p(b){\partial }_{x}u(x{)}_{|x=b}=0& {\zeta }_1^2+{\zeta }_2^2>0.\end{array}$$

$$\mathcal{L}(\cdot ):=-w(x{)}^{-1}\left[{\partial }_x\left(p(x){\partial }_x(\cdot )\right)+q(x)(\cdot )\right],$$

$$\underset{n\to \infty }{\lim }\frac{{\lambda }_n}{n^2}={\pi }^2{\left({\int }_a^b\sqrt{\frac{w(x)}{p(x)}}dx\right)}^{-2}.$$

$$y:={\int }_a^x\sqrt{\frac{w(t)}{p(t)}}dt,$$

$$v(y)=u(x(y))(w(x(y))p(x(y)){)}^{1/4},$$

$$\{\begin{array}{ll}-\Delta v(y)+V(y)v(y)=\lambda v(y)& y\in (0,B)\subset \mathbb{R},\\ {\Sigma }_{1}v(0)-{\Sigma }_2{\partial }_{y}v(y{)}_{|y=0}=0& {\Sigma }_1^2+{\Sigma }_2^2>0,\\ Z_{1}v(B)+Z_2{\partial }_{y}v(y{)}_{|y=B}=0& Z_1^2+Z_2^2>0,\end{array}$$

$$B={\int }_a^b\sqrt{\frac{w(t)}{p(t)}}dt,$$

$$V(y)=\frac{g^{\prime \prime }(y)}{g(y)}+\frac{q(x(y))}{w(x(y))}\in C([0,B]),$$

$${\Sigma }_1=\frac{{\sigma }_1}{(w(a)p(a){)}^{1/4}}+{\sigma }_2\left(\frac{(wp{)}^{\prime }(a)}{4(p(a){)}^{1/4}(w(a){)}^{5/4}}\right),$$

$$Z_1=\frac{{\zeta }_1}{(w(b)p(b){)}^{1/4}}-{\zeta }_2\left(\frac{(wp{)}^{\prime }(b)}{4(p(b){)}^{1/4}(w(b){)}^{5/4}}\right),$$

$$v(0,\lambda )={\Sigma }_2,v^{\prime }(0,\lambda )={\Sigma }_1,$$

$$u(a,\lambda )={\sigma }_2,p(a)u^{\prime }(a,\lambda )={\sigma }_1.$$

$$\Psi (\lambda )={\zeta }_{1}u(b,\lambda )+{\zeta }_{2}p(b)u^{\prime }(b,\lambda )=Z_{1}v(B,\lambda )+Z_2{v}^{\prime }(B,\lambda ).$$

$$\mathbf{x}:=(x^1,\cdots ,x^{\nu },{\theta }^{\nu +1},\cdots ,{\theta }^m),\nu \in \{0,1\cdots ,m\}.$$

$$\{T^{(n)}{\}}_n{\sim }_{\lambda }\omega ,$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{k=1}^{n}F({\lambda }_k(T^{(n)}))=\frac{1}{{\mu }_m(D)}{\int }_{D}F(\omega (\mathbf{x}))d{\mu }_m(\mathbf{x}).$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{k=1}^n{\lambda }_k(T^{(n)})=\frac{1}{{\mu }_m(D)}{\int }_D\omega (\mathbf{x})d{\mu }_m(\mathbf{x}).$$

$${\lambda }_k\left(T^{(n)}(\omega )\right)=\omega \left(\frac{k\pi }{n+1}\right)+O(n^{-1})\text{mbox}foreveryk=1,\cdots ,n,$$

$${\Omega }_{ϵ}:=\bigcup _{z\in \Omega }{B}_{ϵ}^{\mathbb{C}}(z).$$

$$\underset{n\to \infty }{\lim }\frac{\left|\left\{k=1,\dots ,n:{\lambda }_k\left(T^{(n)}\right)\notin {\Omega }_{ϵ}\right\}\right|}{n}=0.$$

$$\left|{\lambda }_1\left(T^{(n)}\right)-z\right|\le \left|{\lambda }_2\left(T^{(n)}\right)-z\right|\le \dots \le \left|{\lambda }_n\left(T^{(n)}\right)-z\right|.$$

$$\underset{n\to \infty }{\lim }\left|{\lambda }_k\left(T^{(n)}\right)-z\right|=0\text{mbox}foreveryfixedk.$$

$$t\in R_{\omega }\iff {\mu }_m\left(\left\{\mathbf{x}\in D:\left|\omega (\mathbf{x})-t\right|<ϵ\right\}\right)>0\forall ϵ>0.$$

$${\left\{T^{(n)}\right\}}_{n\in \mathbb{N}}{\sim }_{\lambda }\omega (\mathbf{x}),\mathbf{x}\in D\subset {\mathbb{R}}^m.$$

$$\{B^{(n,m)}{\}}_n\to \{A^{(n)}{\}}_n\text{mbox}a.c.s.,$$

$$A^{(n)}=B^{(n,m)}+R^{(n,m)}+N^{(n,m)},\mathrm{rank}\left(R^{(n,m)}\right)\le c_1(m)n,\Vert N^{(n,m)}\Vert \le c_2(m),$$

Full text

We prove that the spectral symbol $\omega$ provides a necessary condition for a discretization scheme in order to uniformly approximate the spectrum of the original differential operator $\mathcal{L}$. The condition measures how far the method is from a uniform relative approximation of the spectrum of $\mathcal{L}$. Moreover, the condition seems to become sufficient if the discretization method is paired with a suitable (non-uniform) grid and an increasing refinement of the order of approximation of the method.

On the other hand, despite the numerical experiments in many recent literature, we disprove that in general a uniform sampling of the spectral symbol $\omega$ can provide an accurate relative approximation of the spectrum, neither of $\mathcal{L}$ nor of the discrete operator $\mathcal{L}^{(n)}$.

Keywords Spectral approximation differential operators $\cdot$ Spectral symbol $\cdot$ Eigenvalue problem $\cdot$ Sturm-Liouville $\cdot$ Matrix methods $\cdot$ Generalized locally Toeplitz

Mathematics Subject Classification 65N25 $\cdot$ 65N06 $\cdot$ 65N30 $\cdot$ 65N35 $\cdot$ 65L15

1 Introduction

For simplicity, throughout this paper our main model problem will be the one dimensional Sturm-Liouville equation

[TABLE]

and its associated (regular) Sturm-Liouville eigenvalues problem (SLP)

[TABLE]

As it will be clear in Subsection 5.2, a generalization to other kind of differential operators in higher dimensions can be handled by the same techniques we are going to discuss for this model case.

There are several many approaches to compute the eigenvalues $\lambda_{k}$ related to Problem (1.2), such as matrix methods, Prüfer’s methods, sampling methods, etc. We refer to [40, 47, 55] and [13, 2] for a general basic overview and we cite some recent interesting papers about transmutation operator methods [37, 36].

Matrix methods consist in discretizing Problem (1.2) in the form

[TABLE]

where $\mathcal{L}^{(n)}$ is a matrix of order $n\times n$, with $n$ a mesh fineness parameter, and then solving the associated algebraic eigenvalue problem. They generally suffer of a (fast) growth of the absolute (and relative) error, that is the quality of the numerical approximation they provide for the $k$-th eigenvalue deteriorates significantly as the index $k$ increases. Even if there have been developed suitable correction techniques for improving the accuracy of the computed eigenvalues (see as an example for the FD case the first paper appeared in this sense [39], and a more recent development [1] for the Numerov’s method), if the interest relies only on the computation of the eigenvalues/eigenfunctions, then other methods typically can be a better choice, such as shooting-type methods for example, which can provide more uniformly accurate estimates of the eigenvalues with respect to the index $k$, at the cost of a fair more effort in the implementation, see [40, Chapter 5].

Nevertheless, many other problems demand to discretize equation (1.1) in such a way that the spectrum of the discrete operator $\mathcal{L}^{(n)}$ preserves the spectrum of the continuous operator $\mathcal{L}$ uniformly with respect to $k=1,\cdots,n$. Within this regard, see for example the spectral gap problem of the $1d$ wave equation for the uniform observability of the control waves, [30] and [21, 7], or structural engineering problems, see [29, Section 3.1], or finally the discretization of the Laplace-Beltrami operator by means of the graph-Laplacian, see [42, 10, 33].

In this setting, the theory of Generalized Locally Toeplitz (GLT) sequences provides the necessary tools to understand whether the matrix methods used to discretize the operator $\mathcal{L}$ are efficient or not to spectrally approximate it. The GLT theory originated from the seminal work of Tilli on Locally Toeplitz (LT) sequences [50] and then afterwards developed by Serra-Capizzano in [45, 46]. It was devised to compute and analyze the spectral distribution of matrices arising from the numerical discretization of integral equations and differential equations. One of the purposes of this spectral analysis is the design of efficient numerical methods for computing the related numerical solutions by fast iterative solvers (especially, multigrid and preconditioned Krylov methods), in this sense see for example the related works [17, 18, 19].

It often happens that the discrete operator $\mathcal{L}^{(n)}$ from (1.3), properly weighted by a power of $n$ (depending on the dimension of underlying space and on the maximum order of derivatives involved) and which we will denote as $\hat{\mathcal{L}}^{(n)}$, enjoys an asymptotic spectral distribution as $n\to\infty$, i.e., as the mesh is progressively refined. More precisely, for any test functions $F(t)\in C_{c}(\mathbb{C})$,

[TABLE]

where $\lambda_{k}\left(\hat{\mathcal{L}}^{(n)}\right)$ is the $k$-th eigenvalue of the weighted operator $\hat{\mathcal{L}}^{(n)}$ and $\omega:D\subset\mathbb{R}^{m}\to\mathbb{C}$ is referred to as the spectral symbol of the sequence $\left\{\hat{\mathcal{L}}^{(n)}\right\}_{n\in\mathbb{N}}$, see [9, Definition 2.6] in relation with Toeplitz operators.

The GLT theory serves the purpose to compute the spectral symbol $\omega$ related to matrix methods employed to get the matrix equation (1.3), especially if the numerical method belongs to the family of the so-called local methods, such as FD methods, Finite Element (FE) methods and collocation methods with locally supported basis functions.

In several recent papers it was advanced the suggestion to exploit the sampling of the spectral symbol for approximating the spectrum of the discrete operator $\hat{\mathcal{L}}^{(n)}$, see as a main reference the paper [26] with all the references therein, which is a full review up to the state-of-the-art of the symbol-based analysis for the eigenvalues distribution carried on in the framework of the isogeometric Galerkin approximation (IgA). Unfortunately, this does not apply in general, as we are going to prove.

Nevertheless, the spectral symbol provides to a matrix method a necessary condition for the uniform spectral approximation of the continuous operator $\mathcal{L}$, in the sense of the relative error. The condition is easily computable and it seems to be sufficient if the discretization method is paired with a suitable (non-uniform) grid and an increasing refinement of the order of approximation of the method.

The paper is organized as follows.

• •

Section 2 is a summary about the main properties concerning the SLP model equations (1.1) and (1.2).

• •

In Section 3, the theory of GLT sequences is briefly introduced. In Subsections 3.2 and 3.3 we present the spectral symbols concerning the discretization of (2.2) by means of central FD and IgA schemes of varying order, respectively. In Subsection 3.4 it is introduced the monotone rearrangement of the spectral symbol and it is explained how to compute an approximation of it, whenever it is not feasible to get an exact analytical expression. We highlight the connection with the theory of uniformly distributed sequences, see Theorem 3.4.1.

• •

Section 4 is mainly dedicated to numerical experiments: we analyze the SLP in the Euler-Cauchy equation case. Since the eigenvalues for this regular SLP are computable in closed form, it is a suitable toy-model example with variable coefficients. In Subsection 4.1, by easy calculations we show that a uniform sampling of the monotone rearrangement of the spectral symbol produces a positive lower bound in the relative error approximation, independent by the mesh finesse parameter $n$, which prevents to obtain an accurate approximation of the eigenvalues of the Euler-Cauchy differential operator. In Subsection 4.4, furthermore we present a negative result in the case of spectral symbols $\omega$ belonging to the $L^{1}$ class. Finally, numerical validations are provided to the necessary condition for a uniform spectral approximation.

In Subsection 4.2 and Subsection 4.3 we generalize what observed in the previous subsection to the case of central FD and IgA methods of higher order. By means of a suitable non-uniform grid suggested by the spectral symbol, we observe by numerical experiments that the discretization methods produce a uniform relative approximation of the spectrum of the continuous differential operator, as the mesh fineness parameter $n$ and the order of the methods increase.

• •

In Section 5 we collect the proof of the results discussed previously. In particular:

• –

it is proved that in general a uniform sampling of the spectral symbol $\omega$ does not produce an accurate approximation nor of the eigenvalues of $\mathcal{L}$, nor of the eigenvalues of $\hat{\mathcal{L}}^{(n)}$;

• –

it is provided a necessary condition to a numerical matrix method for the relative uniform approximation of the spectrum of $\mathcal{L}$;

All the results can be extended to different matrix methods and differential operators, once there are known the spectral symbol $\omega$ of the matrix method and the asymptotic expansion of the eigenvalues of the differential operator.

1.1 About the notations

We will use the notation $A^{(n)},B^{(n)},C^{(n)},T^{(n)}$ to indicate general $n\times n$ square matrices, and we will use the notation $\mathcal{L}^{(n,\eta)}$ to indicate a $n\times n$ square matrix which is the discretization of a differential operator $\mathcal{L}$ by means of a numerical method of order of approximation $\eta$. In the case that the approximation order $\eta$ is implied by the context, then we will omit it. If the discretized operator $\mathcal{L}^{(n)}$ is weighted by a constant depending by the finesses mesh parameter $n$, then we will denote it with $\hat{\mathcal{L}}^{(n)}$. We will use the subscripts dir and BCs to indicate a (discretized) differential operator characterized by Dirichlet or generic BCs, respectively. When it will be necessary to highlight the dependency of the differential operator with respect to the coefficients $p(x),q(x),w(x)$ and the endpoints $a,b$, we will write them as right and left subscripts/supscripts. So, for example, the weighted discretization of a Sturm-Liouville operator with Dirichlet BCs by means of the IgA method of order $\eta$ will be denoted by

[TABLE]

In the special case of the (negative) Laplace operator we will use the symbol $-\Delta$, and all the other previous notations will apply.

For all the constants we will use the letter $c$, making explicit the dependency to other parameters if needed.

Finally, if not stated differently, we will always consider a non-decreasing ordering of the eigenvalues $\lambda_{k}$ of a given operator.

2 The model equation

Our model equation will be the following Sturm-Liouville equation with separated boundary conditions (BCs),

[TABLE]

and its associated weighted Sturm-Liouville eigenvalue problem (SLP)

[TABLE]

There is an extensive literature concerning the Sturm-Liouville eigenvalue boundary problem (2.2), see as a (not exhaustive) collection of references [22, 43, 54]. We will write $\partial_{x}f$ and $f^{\prime}$ equivalently to indicate the first derivative of a function $f$ with respect to its argument.

If not otherwise explicitly stated, we will suppose that $p,p^{\prime},w,w^{\prime},q,(pw)^{\prime},(pw)^{\prime\prime}\in C([a,b])$, $p,w>0$, $q\geq 0$ and $\sigma_{1}^{2}+\sigma_{2}^{2}>0,\zeta_{1}^{2}+\zeta_{2}^{2}>0$. Under these conditions, problem (2.2) is regular and the differential operator

[TABLE]

defined over an appropriate domain $\textnormal{dom}\left(\mathcal{L}\right)\subset L^{1}([a,b])$, depending on the boundary conditions, is self-adjoint and has an increasing sequence of positive real eigenvalues $0\leq\lambda_{1}<\lambda_{2}<\cdots<\lambda_{n}$ of multiplicity one and such that

[TABLE]

By the Liouville transformation

[TABLE]

the regular SLP (2.2) converts to the Liouville normal form

[TABLE]

where $\Delta$ is the Laplace operator and

[TABLE]

Indicating with $v(y,\lambda)$ a solution of the differential equation in (2.5) and normalizing it by the initial conditions at $y=0$, namely

[TABLE]

then the corresponding solution $u(x,\lambda)=(r(x)p(x))^{-1/4}v(y(x),\lambda)$ of (1.2) is normalized by

[TABLE]

Therefore, the eigenvalues of the SLP (2.2) are the same of the eigenvalues of the SLP (2.5), being the zeros of the entire function

[TABLE]

In particular, the eigenvalues and eigenfunction norms remain invariant under the change from (2.2) to (2.5).

3 Generalized Locally Toeplitz sequences and spectral symbol

In this section we are going to provide a brief summary about the theory of GLT sequences, starting from the spectral symbol, see Definition 3.1.1. For a detailed treatment of the theory of the GLT sequences we invite the reader to look at [24, 25] and all the references therein. In Subsections 3.2 and 3.3 we present the spectral symbols concerning the discretization of (2.2) by means of $(2\eta+1)$-points central FD approximation and Isogeometric Galerkin approximation (IgA) based on B-splines of degree $\eta$ and smoothness $C^{\eta-1}$. We will not discuss directly here the case of spectral symbols concerning Finite Elements method or IgA with more general smoothness assumptions, in order to avoid that this paper becomes unnecessarily long.

In Subsection 3.4, we present a natural way to approximate the monotone rearrangement of a spectral symbol $\omega$.

3.1 Preliminaries on GLT sequences and spectral symbol

A matrix $T^{(n)}$ is said to be Toeplitz if it is a matrix with constant coefficients along its diagonals, i.e., if it is of the form $T^{(n)}_{i,j}=t_{i-j}$ for all $i,j=1,\ldots,n$, and with $\mathbf{t}=\left[t_{-n+1},\ldots,t_{0},\ldots,t_{n-1}\right]\in\mathbb{C}^{2n-1}$. If $t_{k}$ is the $k$-th Fourier coefficient of a complex integrable function $f$ defined over the interval $[-\pi,\pi]$, then $T^{(n)}=T^{(n)}(f)$ is said to be the Toeplitz matrix generated by $f$. Matrices with a Toeplitz-related structure naturally arise when discretizing, over a uniform grid, problems which have a translation invariance property, such as differential operators with constant coefficients.

From [24] we have the following definitions.

Definition 3.1.1** (Spectral symbol)**

Let $\{T^{(n)}\}_{n\in\mathbb{N}}$ be a sequence of matrices tending to infinity as $n\to\infty$. Let $\omega:D\subset\mathbb{R}^{m}\to\mathbb{C}$ be a measurable function defined on a set $D$ with $0<\mu_{m}\left(D\right)<\infty$, where $\mu_{m}$ is the Lebesgue measure on $\mathbb{R}^{m}$, and

[TABLE]

We say that $\{T^{(n)}\}_{n\in\mathbb{N}}$ has a spectral (or eigenvalue) distribution described by $\omega$, and we write

[TABLE]

if for all functions $F\in C_{c}(\mathbb{C})$ (i.e., continuous with compact support) we have

[TABLE]

In this case, $\omega$ is called spectral (or eigenvalue) symbol of $\{T^{(n)}\}_{n\in\mathbb{N}}$. See Subsection 3.4 for its monotone rearrangement.

Relation (3.2) is satisfied for example by Hermitian Toeplitz matrices generated by real-valued functions $\omega\in L^{1}([-\pi,\pi])$, i.e., $T^{(n)}(\omega)\sim_{\lambda}\omega$, see [28, 51, 52]. For a general overview on Toeplitz operators and spectral symbol, see [9].

Remark 1

In particular, if $D$ is compact, $\omega$ real and continuous, and $\lambda_{k}\in\left[\min\omega,\max\omega\right]$ for every $k\in\mathbb{N}$, then taking $F(\lambda)=\lambda\chi_{D}(\lambda)$, with $\chi_{D}$ a $C^{\infty}$ cut-off of $D$ such that $\chi_{D}(t)\equiv 1$ on $D$, then

[TABLE]

Because the Riemannian sum over equispaced points converges to the integral of the right hand side of the above formula, then (3.3) could suggest us that the eigenvalues $\lambda_{k}(T^{(n)})$ can be approximated by a pointwise evaluation of the symbol $\omega({\bf x})$ over an equispaced grid of $D$, for $n\to\infty$, expect for at most an $o(n)$ of outliers, see Definition 3.1.5. This is mostly the content of [24, Remark 3.2] and [26, Section 2.2].

We remark that in the special cases of Toeplitz matrices generated by enough regular functions $\omega:[0,\pi]\to\mathbb{R}$, then

[TABLE]

see [53] and [8]. In some sense, this justifies the informal meaning given above to the spectral distribution $\omega$.

An important consequence of (3.1) is contained in the next Theorem 3.1.1, but we need a couple of definitions more.

Definition 3.1.2** (Clustering)**

Given a subset $\Omega\subset\mathbb{C}$ and a positive number $\epsilon>0$, define the $\epsilon$-expansion of $\Omega$ as

[TABLE]

We say that a matrix sequence $\left\{T^{(n)}\right\}_{n\in\mathbb{N}}$ is weakly clustered at $\Omega$ if for every $\epsilon>0$ it holds that

[TABLE]

Definition 3.1.3** (Spectral attraction)**

Let $z\in\mathbb{C}$. Given a matrix sequence $\left\{T^{(n)}\right\}_{n\in\mathbb{N}}$, let us order the eigenvalues of $T^{(n)}$ according to their distance from $z$, i.e.,

[TABLE]

We say that $z$ strongly attracts the spectrum of $T^{(n)}$ with infinite order if

[TABLE]

Definition 3.1.4** (Essential range)**

Let $\omega$ be a $\mu_{m}$-measurable function and define the set $R_{\omega}\subset\mathbb{C}$ as

[TABLE]

We call $R_{\omega}$ the essential range of $\omega$. $R_{\omega}$ is closed.

Theorem 3.1.1

Let $\{T^{(n)}\}$ be a matrix sequence such that

[TABLE]

Then $\{T^{(n)}\}_{n}$ is weakly clustered at $R_{\omega}$ and every point of $R_{\omega}$ strongly attracts the spectrum of $T^{(n)}$ with infinite order.

Proof. See [24, Theorem 3.1].

Definition 3.1.5** (Outliers)**

Given a matrix sequence $\{T^{(n)}\}$ such that $\left\{T^{(n)}\right\}_{n\in\mathbb{N}}\sim_{\lambda}\omega$, if $\lambda_{k}\left(T^{(n)}\right)\notin R_{\omega}$, we call it an outlier.

The matrices which represent the discrete version of a general differential operator do not always own a Toeplitz structure and therefore we cannot know beforehand whether they posses a spectral symbol or not. Hereafter we begin to introduce the concepts which will bring us to the definition of GLT sequences.

Definition 3.1.6** (Approximating class of sequences)**

Let $\{A^{(n)}\}_{n\in\mathbb{N}}$ be a matrix-sequence and let $\{\{B^{(n,m)}\}_{n\in\mathbb{N}}\}_{m\in\mathbb{N}}$ be a sequence of matrix-sequences. We say that $\{\{B^{(n,m)}\}_{n}\}_{m}$ is an approximating class of sequences (a.c.s.) for $\{A^{(n)}\}_{n}$, and we write

[TABLE]

if the following condition is met: for every $m$ there exists $n_{m}$ such that, for $n>n_{m}$,

[TABLE]

where $n_{m},c_{1}(m),c_{2}(m)$ depend only on $m$, and

[TABLE]

Roughly speaking, $\{\{B^{(n,m)}\}_{n}\}_{m}$ is an a.c.s. for $\{A^{(n)}\}_{n}$ if, for all sufficiently large $m$, the sequence $\{B^{(n,m)}\}_{n}$ approximates $\{A^{(n)}\}_{n}$ in the sense that $A^{(n)}$ is eventually equal to $B^{(n,m)}$ plus a small-rank matrix (with respect to the matrix size $n$) plus a small-norm matrix.

Definition 3.1.7** (Locally Toeplitz operators and Locally Toeplitz sequences)**

Let $m,n\in\mathbb{N}$, let $a:[0,1]\to\mathbb{C}$ and let $f\in L^{1}[-\pi,\pi]$. The locally Toeplitz (LT) operator is defined as the following $n\times n$ matrix:

[TABLE]

where $T^{\left(\lfloor n/m\rfloor\right)}(f)$ is a Toeplitz matrix of size $\lfloor n/m\rfloor$ generated by $f$ and $O_{n\bmod m}$ is the zero matrix of size $n\bmod m$.

Let $\{A^{(n)}\}_{n\in\mathbb{N}}$ be a matrix-sequence. We say that $\{A^{(n)}\}_{n}$ is locally Toeplitz (LT) sequence with symbol $a\otimes f$, and we write $\{A^{(n)}\}_{n}\sim_{\textnormal{LT}}a\otimes f$, if

[TABLE]

The functions $a$ and $f$ are, respectively, the weight function and the generating function of $\{A^{(n)}\}_{n}$.

We can finally give the definition of GLT sequence.

Definition 3.1.8** (GLT sequence)**

Let $\{A^{(n)}\}_{n\in\mathbb{N}}$ be a matrix-sequence and let $\omega:[0,1]\times[-\pi,\pi]\to\mathbb{C}$ be a measurable function. We say that $\{A^{(n)}\}_{n}$ is a GLT sequence with symbol $\omega$, and we write $\{A^{(n)}\}_{n}\sim_{\textnormal{GLT}}\omega$, if the following condition is met: for every $\epsilon>0$ there exists a finite number of LT sequences $\{A^{(n)}_{(i,\epsilon)}\}_{n}\sim_{\textnormal{LT}}a_{i,\epsilon}\otimes f_{i,\epsilon}$, $i=1,\cdots,N_{\epsilon}$, such that

1. (1)

$\sum_{i=1}^{N_{\epsilon}}a_{i,\epsilon}\otimes f_{i,\epsilon}\to\omega$* in measure as $\epsilon\to 0$;* 2. (2)

$\left\{\sum_{i=1}^{N_{\epsilon}}A_{(i,\epsilon)}^{(n)}\right\}\to\{A_{n}\}_{n}$* a.c.s. as $\epsilon\to 0$.*

We have the following main property (see [24, Property GLT 1 p. 170]) which connects the GLT symbol with the spectral symbol of Definition 3.1.1.

Proposition 3.1.1

If $\{A^{(n)}\}_{n}\sim_{\textnormal{GLT}}\omega$ and $A^{(n)}$ are Hermitian, then $\{A^{(n)}\}_{n}\sim_{\lambda}\omega$.

3.2 $(2\eta+1)$-points central FD discretization of the SLP (2.2)

The (central) Finite Difference method basically consists in the approximation of the $j\textsuperscript{th}$-derivative $u^{(j)}(x_{0})$ by means of $2\eta$ Taylor expansions, centered at $x_{0}\in(a,b)$, at $2\eta$ points $\{x_{-\eta},\ldots,x_{-1},x_{1},\ldots,x_{\eta}\}$ with $2\eta+1\geq j$ and such that $x_{-\eta}<\ldots<x_{0}<\ldots<x_{\eta}$. Given a standard equispaced grid ${\bf x}=\left\{x_{j}\right\}_{j=1-\eta}^{n+\eta}\subseteq\left[\bar{a},\bar{b}\right]$, with

[TABLE]

let us consider a $C^{1}$-diffeomorphism $\tau:[a,b]\to[a,b]$ such that $\tau^{\prime}(x)\neq 0$, $\tau(a)=a,\tau(b)=b$ and let us consider its piecewise $C^{1}$-extension $\bar{\tau}:\left[\bar{a},\bar{b}\right]\to\left[\bar{a},\bar{b}\right]$ such that

[TABLE]

By means of the piecewise $C^{1}$-diffeomorphism $\bar{\tau}$ we have a new grid $\bar{{\bf x}}=\left\{\bar{x}_{j}\right\}_{j=1-\eta}^{n+\eta}\subset[\bar{a},\bar{b}]$, non necessarily uniformly equispaced, with $\bar{x}_{j}=\bar{\tau}(x_{j})$. Combining together the high-order central FD schemes in [3, 4, 31], it is not difficult to obtain the following general matrix eigenvalue problem to approximate the SLP (2.2) in the case of Dirichlet BCs:

[TABLE]

with

[TABLE]

[TABLE]

and finally

[TABLE]

where, if we define the $C^{0}$-extension of $p(x)$ to $[\bar{a},\bar{b}]$ as

[TABLE]

and the element $l_{i,j}$ as

[TABLE]

then the generic matrix element of $\prescript{a}{b}{\mathcal{L}^{(n,\eta)}_{\textnormal{dir},p(\bar{x})}}$ is given by

[TABLE]

With abuse of notation, we will call $\eta$ the order of approximation of the central FD method. We have the following results.

Theorem 3.2.1

In the above assumptions, for $\eta\geq 1$, defining

[TABLE]

it holds that

1. (i)

the eigenvalues $\lambda_{k}\left(\prescript{a}{b}{\mathcal{L}^{(n,\eta)}_{\textnormal{dir},p(\bar{x}),q(\bar{x}),w(\bar{x})}}\right)$ are real for every $k$ and

[TABLE] 2. (ii)

[TABLE]

where

[TABLE]

and

[TABLE]

Proof. See Appendix A.

Corollary 3.2.1

The function $f_{\eta}$ from Theorem 3.2.1 is differentiable, nonnegative, monotone increasing on $[0,\pi]$ and it holds that

[TABLE]

Proof. See Appendix A.

3.3 Isogeometric Galerkin discretization of the SLP (2.2) by B-slpine of degree $\eta$ and smoothness $C^{\eta-1}$

The Galerkin discretization approach deals with the SLP in its weak formulation. For simplicity, let $q(x)\equiv 0$ and $\sigma_{2}=\zeta_{2}=0$ in (2.2). In the standard Galerkin method, fix a set of basis functions $\{\vartheta_{1},\ldots,\vartheta_{N}\}\subset H^{1}_{0}((a,b))$ defined on the reference domain $[a,b]$ and vanishing on the boundary, and look for approximations of the exact eigenpairs $\left\{\lambda_{k};u_{k}\right\}$ of Problem (2.2) by solving the following matrix eigenvalue problem:

[TABLE]

where

[TABLE]

are the stiffness and mass matrices, respectively. In the isogeometric Galerkin method (IgA), the physical domain $[a,b]$ is described by a global geometry map $\tau:[a,b]\to[a,b]$ such that $\tau\in C^{1}([a,b])$, $\tau(a)=a,\tau(b)=b$ and $\tau^{\prime}(x)\neq 0$. Within this reparametrization of the domain, the stiffness and mass matrices assume the following form

[TABLE]

The numerical eigenvalues are simply the eigenvalues of $L_{N}=M_{N}^{-1}K_{N}$. Due to the assumption of strictly positivity of $p$ and $w$, then both $K_{N}$ and $M_{N}$ are symmetric positive definite matrices for any basis functions and diffeomorphism $\tau$. In the next theorem we collect the results of [24, Theorem 10.15] in a shortened way and according to our notations. For a detailed treatment of the IgA, see [14].

Theorem 3.3.1

Let $(a,b)$ be discretized by a uniform mesh $\left\{x_{j}\right\}_{j=1}^{n}$ of step-size $(n+1)^{-1}$ (i.e., it divides the interval $(a,b)$ into $n+1$ equispaced subintervals) and let $\tau:[a,b]\to[a,b]$ a $C^{1}$-diffeomorphism such that $\tau(a)=a,\tau(b)=b$ and $\tau^{\prime}(x)\neq 0$ for every $x\in[a,b]$. Let $\eta\geq 1$ and let the basis functions $\{\vartheta_{1},\ldots,\vartheta_{N}\}$ be taken as $\{B_{1,[\eta,\eta-1]},\ldots,B_{n+\eta-2,[\eta,\eta-1]}\}$, where $B_{j,[\eta,\eta-1]}$ for $j=1,\cdots,(n+1)+\eta-2$ are the B-splines of degree $\eta$ and smoothness $C^{\eta-1}([0,1])$ defined on the knot sequence

[TABLE]

Indicating with $\psi_{[s]}$ the cardinal B-spline of degree $s$, then

1. (i)

for every fixed $k$

[TABLE] 2. (ii)

it holds the following spectral distribution

[TABLE]

where

[TABLE]

and $g_{\eta}(\theta)$ and $h_{\eta}(\theta)$ are given by

[TABLE]

For $\eta=1,2,3,4$, equation (3.7) gives

[TABLE]

Proof. For item (i) see [6, 41]. For item (ii), see [24, Theorem 10.15].

We have an analogue of Corollary 3.2.1.

Corollary 3.3.1

The function $f_{\eta}$ is differentiable, nonnegative, monotone increasing on $[0,\pi]$ and it holds that

[TABLE]

Proof. See [20, Theorem 1, Theorem 2 and Lemma 1].

With abuse of notation, we will call $\eta$ the order of approximation of the IgA method.

3.4 The monotone rearrangement of $\omega$

From here on, if not differently stated, we will consider $\omega:D\subset\mathbb{R}^{m}\to\mathbb{R}$ to be a real valued measurable function in $L^{1}\left(D,\mu_{m}\right)$.

It is always nice to deal with an univariate and monotone spectral symbol $\omega({\bf x})=\omega(\theta^{1})$. Unfortunately, in general $\omega$ is multivariate or not monotone. Nevertheless, in such cases, it is possible to consider a rearrangement $\tilde{\omega}:[0,1]\to[\inf R_{\omega},\sup R_{\omega}]$, where $R_{\omega}$ is the essential range of $\omega$ as in Definition 3.1.4, such that $\tilde{\omega}$ is univariate, monotone nondecreasing and

[TABLE]

i.e.,

[TABLE]

This can be achieved by defining

[TABLE]

where

[TABLE]

Clearly, $\tilde{\omega}$ is a.e. unique, univariate, monotone strictly increasing and right-continuous, which make it a good choice for an equispaced sampling. On the other hand, $\tilde{\omega}$ could not have an analytical expression or it could be not feasible to calculate, therefore it is often needed an approximation $\tilde{\omega}_{r}$. Hereafter we summarize the steps presented in [24, Example 10.2] and [26, Section 3] in order to approximate the eigenvalues $\lambda_{k}\left(T^{(n)}\right)$ by means of an equispaced sampling of the rearranged symbol $\tilde{\omega}$ (or its approximated version $\tilde{\omega}_{r}$). For the sake of clarity and without loss of generality, we suppose that $D=[0,1]\times[0,\pi]$ and $\omega$ continuous.

Algorithm I

Fix $r\in\mathbb{N}$ such that $r=r(n)\geq n$, and fix the equispaced grid $\{(x_{i},\theta_{j})\}$ over $[0,1]\times[0,\pi]$, where $x_{i}=\frac{i}{r+1}$, $\theta_{j}=\frac{j\pi}{r+1}$ for $i,j=1,\cdots,r$; 2. 2)

Get the set of samplings $\left\{\omega(x_{i},\theta_{j})\right\}_{i,j=1}^{r}$ and form a nondecreasing sequence $\{\omega_{1}\leq\omega_{2}\leq\cdots\leq\omega_{r^{2}}\}$; 3. 3)

Define $\tilde{\omega}_{r}:[0,1]\to[\min\omega,\max\omega]$ as the piecewise linear nondecreasing function which interpolates the samples $\{\min\omega=\omega_{0}\leq\omega_{1}\leq\omega_{2}\leq\cdots\leq\omega_{r^{2}}\leq\omega_{r^{2}+1}=\max\omega\}$ over the nodes $\{0,\frac{1}{r^{2}+1},\frac{2}{r^{2}+1},\cdots,\frac{r^{2}}{r^{2}+1},1\}$; 4. 4)

Sample $\tilde{\omega}_{r}$ over the set $\left\{\frac{k}{n+1}\right\}_{j=1}^{n}$ and define

[TABLE] 5. 5)

Finally, approximate the eigenvalues of $T_{n}$ by $\tilde{\omega}_{r,k}^{(n)}$, i.e.,

[TABLE]

Obviously, if $\tilde{\omega}$ is available then use it instead of $\tilde{\omega}_{r}$ and define

[TABLE]

As standard result in approximations of monotone rearrangements, it holds that $\|\tilde{\omega}_{r(n)}-\tilde{\omega}\|_{\infty}\to 0$ as $n\to\infty$, see [12, 49]. See moreover [16, Definition 3.1 and Theorem 3.3] and [44], where the monotone rearrangement were first introduced in the context of spectral symbol. We have the following limit relation.

Theorem 3.4.1** (Discrete Weyl’s law)**

Let $\tilde{\omega}:[0,1]\to[\min R_{\omega},\max R_{\omega}]$ be the monotone rearrangement 3.9 of a spectral symbol $\omega$ of the matrix sequence $\left\{T^{(n)}\right\}_{n}$. Let $\tilde{\omega}$ be piecewise Lipschitz continuous. Then

[TABLE]

In particular, let $k=k(n)$ be such that $k(n)/n\sim x$ as $n\to\infty$ for a fixed $x\in[0,1]$ and $\lambda_{k(n)}\left(T^{(n)}\right)\in\tilde{\omega}\left([0,1]\right)=R_{\omega}$ definitely. Then

[TABLE]

Proof. By our assumptions, $\tilde{\omega}\left([0,1]\right)=R_{\omega}=\bigcup_{j=1}^{m}[a_{j},b_{j}]$ such that $a_{1}=\min R_{\omega}$, $b_{m}=\max R_{\omega}$, and since $\tilde{\omega}$ is rectifiable then it induces a positive Borel measure on $\tilde{\omega}\left([0,1]\right)$, which we call $d\tilde{\omega}(t)$. If the subsets $[a_{j},b_{j}]$ are disjoint, then $\tilde{\omega}^{-1}$ can be extended on all over $[\min R_{\omega},\max R_{\omega}]$ by defining $\tilde{\omega}^{-1}(t)=\tilde{\omega}^{-1}(a_{j})$ for every $t\in[b_{j-1},a_{j})$. The same applies to $d\tilde{\omega}(t)$, which can be extended to a positive measure defined on $[\min R_{\omega},\max R_{\omega}]$, which we keep calling $d\tilde{\omega}(t)$, such that $d\tilde{\omega}\left((b_{j-1},a_{j})\right)=0$. It holds that

[TABLE]

Then equation (3.11) follows immediately from [38, Theorem 7.1] and the definition of asymptotic distribution function mod $1$, [38, Definition 7.1], and recalling that $\left\{T^{(n)}\right\}_{n}$ is weakly clustered at $R_{\omega}$ by Theorem 3.1.1. Equation (3.12) follows instead from (3.11) and again by Theorem 3.1.1, taking $t=\lambda_{k(n)}\left(T^{(n)}\right)$ and applying $\tilde{\omega}$ to both sides.

Corollary 3.4.1

In the same hypothesis of Theorem 3.4.1 and in presence of no outliers as in Definition 3.1.5, then the absolute error between a uniform sampling of $\tilde{\omega}$ and the eigenvalues of $T^{(n)}$ converges to zero, namely

[TABLE]

Proof. Let us suppose that there exists a sequence $\left\{k(n)\right\}_{n}$ such that $\left|\lambda_{k(n)}\left(T^{(n)}\right)-\tilde{\omega}\left(\frac{k(n)}{n+1}\right)\right|\geq c>0$ for every $n$. Since $\left\{\frac{k(n)}{n}\right\}\subset[0,1]$ is bounded then there exists a convergent subsequence, which we keep calling $k(n)$, such that $k(n)/n\to x\in[0,1]$. Therefore $\left|\lambda_{k(n)}\left(T^{(n)}\right)-\tilde{\omega}\left(x\right)\right|\geq c>0$, which contradicts (3.12).

In some sense, $\tilde{\omega}$ can be intended as the inverse cumulative distribution function of the eigenvalue distribution of $\left\{T^{(n)}\right\}_{n}$, and the limit relation (3.8) as the strong law for large numbers for specially chosen sequences of dependent complex-valued random variables. See [32] as a recent survey about equidistributions from a probabilistic point of view.

4 Application to the Euler-Cauchy differential operator

We are now in position to begin our analysis with respect to a toy-model example. The main tasks of this section are:

• •

to disprove that, in general, a uniform sampling of the spectral symbol $\omega({\bf x})$ can provide an accurate approximation of the eigenvalues of the weighted and un-weighted discrete operators $\hat{\mathcal{L}}^{(n)}$ and $\mathcal{L}^{(n)}$, respectively;

• •

to show numerical evidences of Theorem 5.1.1, i.e., that the spectral symbol $\omega({\bf x})$ measures how far a matrix discretization method is from a uniform approximation of the spectrum of the continuous operator $\mathcal{L}$, in the sense of the relative error.

• •

to show that a matrix discretization method, if paired with a suitable (non-uniform) grid can produce a uniform relative approximation of $\mathcal{L}$, as the order of approximation of the method increases.

Let us fix $\alpha>0$ and let us consider the following SLP with Dirichlet BCs,

[TABLE]

that is (2.2) with $p(x)=\alpha x^{2}$, $q(x)\equiv 0$, $(a,b)=(1,\textrm{e}^{\sqrt{\alpha}})$ and $\sigma_{1}=\zeta_{1}=1$,$\sigma_{2}=\zeta_{2}=0$. It is an Euler-Cauchy differential equation and by means of the transformation (2.4), i.e., $y(x)=\int_{1}^{x}\left(\sqrt{\alpha}t\right)^{-1}dt$, Problem (4.1) is (spectrally) equivalent to

[TABLE]

which is the normal form of a Sturm-Liouville equation, with constant reaction term (or potential) $V(y)\equiv-\alpha/4$. If we write

[TABLE]

for the self-adjoint operator in Equation (4.1) and

[TABLE]

for the self-adjoint operator in the form of Equation (4.2), then it is clear that

[TABLE]

For later reference, notice that

[TABLE]

namely, the diffusion coefficient $p(x)=\alpha x^{2}$ produces a constant shift of $\alpha/4$ to the eigenvalues of the unperturbed Laplacian operator with Dirichlet BCs, i.e., $-\Delta_{\textnormal{dir}}$.

We introduce the following definition.

Definition 4.0.1** (Numerical and analytic relative error)**

Let $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n)}}$ be the discrete differential operator obtained from (4.3) by means of a generic numerical discretization method, and suppose that

[TABLE]

Fix $n,n^{\prime},r\in\mathbb{N}$, with $n^{\prime}>>n$ and $r=r(n)\geq n$, and compute the following quantities

[TABLE]

where $\prescript{}{\alpha}{\tilde{\omega}_{r,k}^{(n)}}=\prescript{}{\alpha}{\tilde{\omega}_{r}}\left(\frac{k}{n+1}\right)$ and $\prescript{}{\alpha}{\tilde{\omega}_{r}}$ is the (approximated) monotone rearrangement of the spectral symbol $\prescript{}{\alpha}{\omega}$ obtained by the procedure described in Algorithm I. We call $\prescript{}{\alpha}{\textnormal{{err}}_{k}^{(n)}}$ the numerical relative error and $\prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{r,k}^{(n)}}$ the analytic relative error. We say that $\prescript{}{\alpha}{\omega}$ spectrally approximates the discrete differential operator $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n)}}$ if

[TABLE]

4.1 Approximation by $3$-points central FD method on uniform grid

In our example, if we apply the standard central $3$-point FD scheme as in (3.4) with $\eta=1$ and $\tau(x)=x$, then the weighted discretization matrix $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n)}}$ of Problem (4.1) has spectral symbol

[TABLE]

Working with this toy-model problem in the $3$-points central FD scheme provides us a further advantage, since we can calculate the exact monotone rearrangement $\prescript{}{\alpha}{\tilde{\omega}}$, or at least a finer approximation than $\prescript{}{\alpha}{\tilde{\omega}_{r}}$ which does not depend on the extra parameter $r$ and which is less computationally expensive. Indeed, from equation (3.10) we have that

[TABLE]

where

[TABLE]

and where

[TABLE]

Clearly, $\prescript{}{\alpha}{\tilde{\omega}}=\prescript{}{\alpha}{\phi}^{-1}$, and having an analytic expression for $\prescript{}{\alpha}{\phi}(t)$, it is then possible to compute a numerical approximation of its inverse $\prescript{}{\alpha}{\phi}^{-1}:[0,1]\to\left[0,\frac{4\alpha\textrm{e}^{2\sqrt{\alpha}}}{\left(\textrm{e}^{\sqrt{\alpha}}-1\right)^{2}}\right]$ over the uniform grid $\left\{\frac{k}{n+1}\right\}_{k=1}^{n}$, for example by means of a Newton method. This approximation of $\prescript{}{\alpha}{\phi}^{-1}$ does not depend on the extra parameter $r$ and with abuse of notation we will call it $\prescript{}{\alpha}{\tilde{\omega}}$. Henceforth, we will work with both $\prescript{}{\alpha}{\tilde{\omega}_{r}}$ and $\prescript{}{\alpha}{\tilde{\omega}}$, and when we will compute the analytical relative error with respect to $\prescript{}{\alpha}{\tilde{\omega}}$, we will write $\prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{k}^{(n)}}$ without the subscript $r$.

Finally, let us begin our analysis. In the above Figure 1 it is possible to check how an equispaced sampling of $(n+1)^{2}\prescript{}{\alpha}{\tilde{\omega}_{r}}$ seems to distribute exactly as the eigenvalues of the unweighted discrete operator $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n)}}$.

Moreover, according with equation (4.4), we observe that

[TABLE]

and so

[TABLE]

which means that the monotone rearrangement $\prescript{}{\alpha}{\tilde{\omega}}$ converges to the spectral symbol $\omega$ as $\alpha\to 0$, with

[TABLE]

that is the spectral symbol which characterizes the differential operator $-\prescript{1}{0}{\Delta_{\textnormal{dir}}}$ discretized by means of a 3-points FD scheme. In Figure 2 and Table 1 it is visually and numerically summarized this last observation. The eigenvalues of $(n+1)^{-2}\left(-\prescript{1}{0}{\Delta_{\textnormal{dir}}^{(n)}}\right)$ are the exact sampling of $\omega(\theta)=2-2\cos(\theta\pi)$ over the uniform grid $\left\{\frac{k}{n+1}\right\}_{k=1}^{n}$, see [48, p. 154].

All these remarks would suggest that $\prescript{}{\alpha}{\tilde{\omega}}$, or equivalently $\prescript{}{\alpha}{\tilde{\omega}_{r}}$, spectrally approximates the weighted discrete operator $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\hat{\mathcal{L}}_{\textnormal{dir},\alpha x^{2}}^{(n)}}$.

Unfortunately, this conjecture looks to be partially proven wrong by Figure 3. There it is shown the comparison between the graphs of the numerical relative error $\prescript{}{\alpha}{\textnormal{{err}}_{k}^{(n)}}$ and the analytic relative error $\prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{r,k}^{(n)}}$, for several different increasing values of the parameter $r$. We observe a discrepancy in the analytical prediction of the eigenvalue error $\prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{r,k}^{(n)}}$, for small $k<<n$, with respect to the numerical relative error $\prescript{}{\alpha}{\textnormal{{err}}_{k}^{(n)}}$.

In particular, the maximum discrepancy is achieved at the first eigenvalue approximation $k=1$, for every $r$. The discrepancy apparently decreases as the number of grid points $r$ increases, as well observed in [26, Figure 48] for some test-problems in the setting of Galerkin discretization by linear $C^{0}$ B-spline. In that same paper, some plausible hypothesis and suggestions were advanced:

• •

the discrepancy could depend on the fact that it has been used $\prescript{}{\alpha}{\tilde{\omega}_{r}}$ instead of $\prescript{}{\alpha}{\tilde{\omega}}$, and then that discrepancy should tend to zero in the limit $r\to\infty$, since $\prescript{}{\alpha}{\tilde{\omega}_{r}}\to\prescript{}{\alpha}{\tilde{\omega}}$.

• •

numerical instability of the analytic relative error $\prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{r,k}^{(n)}}$ for small eigenvalues,[26, Remark 3.1].

• •

Change the sampling grid into an “almost uniform”grid: see [26, Rermark 3.2] for details.

The problem is that these hypothesis, which stem from numerical observations, cannot be validated: the descent to zero of the observed discrepancy as $r$ increases is only apparent. Indeed, as we are going to show below, for every fixed $k$ it holds that

[TABLE]

with $c_{\alpha,k}$ independent of $n$ and $n^{\prime}$. Let us explain here the reason for this mismatch between the numerical relative error $\prescript{}{\alpha}{\textnormal{{err}}_{k}^{(n)}}$ and the analytic relative error $\prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{r,k}^{(n)}}$.

First of all, as we already observed in Subsection 3.4, for $r\to\infty$ it holds that $\|\tilde{\omega}_{r}-\tilde{\omega}\|_{\infty}\to 0$. Then we can work directly with $\tilde{\omega}$, avoiding us to pass through its approximation $\tilde{\omega}_{r}$. Moreover, for fixed $n<<n^{\prime}$ and $n^{\prime}$ large enough, it holds that

[TABLE]

see for example [27, Lemma 7.1 and Lemma 7.2]. So, from here after we will consider $n^{\prime}$ large enough such that the above approximation holds. We can then rewrite $\prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{r,k}^{(n)}}$ as

[TABLE]

From (4.5), for $0\leq t\leq\frac{4\alpha}{\left(\textrm{e}^{\sqrt{\alpha}}-1\right)^{2}}$, it holds that

[TABLE]

with $\Theta(t,\alpha)$ a monotone increasing function with respect to $t$ and such that $\Theta(t,\alpha)=o(t)$, for every fixed $\alpha$. In particular, for $0<t<<1$ it holds that

[TABLE]

Therefore, for every $0\leq x<1$ sufficiently small, from (3.9) we have that

[TABLE]

and then, provided that $\frac{N}{n}<<1$, for every $k=1,\ldots,N$

[TABLE]

If now we let $r=r(n)\to\infty$ as $n\to\infty$, then for every fixed $k\in\mathbb{N}$ it holds that

[TABLE]

We then have a lower bound for the analytic relative error which can not be avoided by refining the grid points. Of course, as $n\to\infty$, then $c_{\alpha,k}\to 0$ as $k$ increases. This is summarized below by Figure 4 and Table 2.

The problem lies on the wrong informal interpretation given to the limit relation in Definition 3.1.1, and suggested by Remark 1. Indeed, that limit relation tells us that

[TABLE]

or equivalently

[TABLE]

for every $k(n)$ such that $k(n)/n\to x\in[0,1]$, see Theorem 3.4.1. Therefore, since $\prescript{}{\alpha}{\tilde{\omega}}\in C^{1}([0,1])$, it follows that

[TABLE]

as observed for example in [24, Example 10.2 p. 198]. On the contrary, a uniform sampling of the symbol $\omega$ does not necessarily provide an accurate approximation of the eigenvalues of the weighted operator $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\hat{\mathcal{L}}_{\textnormal{dir},\alpha x^{2}}^{(n)}}$, in the sense of the relative error. The uniform sampling of the symbol works perfectly only on certain subclasses of symbol functions $\omega$, but it fails in general. We loose convergence even for the absolute error as soon as we relax the hypothesis on $\omega$, requiring it to be just Lebesgue integrable over $D$, see Subsection 4.4.

As a last remark, in general there does not exist an “almost”uniform grid as well, nor in an asymptotic sense as described in [26, Rermark 3.2]. Knowing the exact sampling grid $\left\{\tau\left(\frac{j}{n+1}\right)\right\}_{j=1}^{n}$ which guarantees $\prescript{}{\alpha}{\tilde{\omega}}$ to spectrally approximate the discrete differential operator $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}^{(n)}}$ is equivalent to know the eigenvalue distribution of the continuous differential operator $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}$.

What we can say instead is that every $t\in[\min\prescript{}{\alpha}{\omega},\max\prescript{}{\alpha}{\omega}]$ strongly attracts the spectrum of $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\hat{\mathcal{L}}_{\textnormal{dir},\alpha x^{2}}^{(n)}}$ with infinite order, see Definition 3.1.3 and Theorem 3.1.1. In particular, in presence of no outliers, it holds that

[TABLE]

and so

[TABLE]

In Table 3 we can see that (4.10) is confirmed.

Therefore, any numerical scheme which is characterized by a spectral symbol of the form $\prescript{}{\alpha}{\omega(x,\theta)}=\frac{p(x)}{(\textrm{e}^{\sqrt{\alpha}}-1)^{2}}f(\theta)$ such that

[TABLE]

will not provide a relative uniform approximation of the eigenvalues $\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right)$ of the continuous differential operator of Problem 4.1, namely

[TABLE]

see Corollary 5.1.3. This is the case of the method we implemented in this subsection, i.e., the 3-points uniform FD scheme, as it can be easily checked by Table 3 and Figure 5.

Actually, Theorem 5.1.1 says more: it says that if $\prescript{}{\alpha}{\tilde{\omega}}(x)\neq x^{2}\pi^{2}$ and in presence of no outliers, it holds that

[TABLE]

In Figure 5 and Table 4 it is validated numerically (4.13), i.e., the thesis of Theorem 5.1.1 set in this specific case.

4.2 Discretization by $(2\eta+1)$-points central FD method on non-uniform grid

Clearly, everything said in the preceding Subsection 4.1 remains valid even if we increase the order of accuracy of the FD method, namely, the spectral symbol $\omega_{\eta}$ of equation (3.4) does not spectrally approximate the discrete differential operator $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha\tau(x)^{2}}^{(n,\eta)}}$, in the sense of the relative error, for any $\eta\geq 1$. See Figure 6 in relation with Figure 3.

What is interesting instead is to change the sampling grid and to increase the order of accuracy $\eta$ of the FD discretization method. Indeed, as it was observed in the relations (4.11), (4.12) and (4.13), it is not possible to achieve a relative uniform approximation of the eigenvalues $\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right)$ if $\prescript{}{\alpha}{\tilde{\omega}_{\eta}}(x)\neq x^{2}\pi^{2}$. From (3.4), for every $\eta\geq 1$ it is easy to check that $\max_{(x,\theta)\in[a,b]\times[0,\pi]}\prescript{}{\alpha}{\omega_{\eta}}(x,\theta)=\prescript{}{\alpha}{\tilde{\omega}_{\eta}}(1)\neq\pi^{2}$, and so we do not have any improvement by just increasing the order of accuracy $\eta$. But let us observe that:

• •

from Corollary 3.2.1 and equation (3.4),

[TABLE]

• •

$\tau(y)=\textrm{e}^{\sqrt{\alpha}y}$ is such that $\frac{\alpha\tau(y)^{2}}{\left(\tau^{\prime}(y)\right)^{2}}\equiv 1$.

In some sense, the spectral symbol $\omega_{\eta}$ suggests us to change the uniform grid

[TABLE]

by means of the diffeomorphism induced by the Liouville transformation. Indeed, from (2.4) we have that

[TABLE]

and therefore we can construct a $C^{1}$-diffeomorphism $\tau:[1,\textrm{e}^{\sqrt{\alpha}}]\to[1,\textrm{e}^{\sqrt{\alpha}}]$ such that

[TABLE]

[TABLE]

The new non-uniform grid is then given by

[TABLE]

and it holds that

[TABLE]

and therefore

[TABLE]

Of course, this last equality is not enough to guarantee that

[TABLE]

since

[TABLE]

is only a necessary condition, see Theorem 5.1.1. Nevertheless, from figures 7, 8 and Table 5 we can see that (4.16) seems to be validated numerically. This would suggest that condition (4.15) becomes sufficient if paired by a suitable rearranged (non-uniform) grid and an increasing refinement of the method order. This phenomenon is not confined to this specific case but it occurs in several many different other cases we tested. More investigations will be carried on in future works.

Finally, in Table 6 we check numerically the validity of Theorem 5.1.1.

4.3 IgA discretization by B-slpine of degree $\eta$ and smoothness $C^{\eta-1}$

In this subsection we continue our analysis in the IgA framework. We just collect all the numerical results of the tests, which confirm again what observed in subsections 4.1 and 4.2. The only difference relies on the fact that we took out the largest eigenvalues of the discrete operator $\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha\tau(x)^{2}}^{(n+\eta-1,\eta)}}$. This is due to the fact that the IgA discretization suffers of a fixed number of outliers which depends on the degree $\eta$ and it is independent of $n$, see [14, Chapter 5.1.2 p. 153]. So, we consider only the eigenvalues $\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha\tau(x)^{2}}^{(n+\eta-1,\eta)}}\right)$ such that

[TABLE]

We stress out the fact that the number of outliers is fixed for every $n$, in accordance with Theorem 3.1.1.

In Figure 9 we observe again the discrepancy between the analytic relative errors and the numerical relative error.

In Figure 10 and Figure 11 we compare the graphs of the eigenvalue distributions and the relative errors, respectively, between the discrete eigenvalues $\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha\tau(x)^{2}}^{(n+\eta-1,\eta)}}\right)$ and the exact eigenvalues $\lambda_{k}\left(\prescript{\textrm{e}^{\sqrt{\alpha}}}{1}{\mathcal{L}_{\textnormal{dir},\alpha x^{2}}}\right)$, for different values of $\eta$ on uniform and non-uniform grids. They line up with the numerics of Table 7: if the sampling grid is given by (4.14), the maximum relative error decreases as the order of approximation increases.

4.4 The $L^{1}$ case

We leave for a moment the case of regular SLPs. In [23] it was addressed the issue of extending the spectral symbol analysis to the case of $\textnormal{L}^{1}$ coefficients. We show in the next example that the sampling of the spectral symbol does not provide accurate approximation of the eigenvalues of the weighted matrix discretization operator even in the sense of the absolute error.

Let us fix $p(x)=x^{-1/2}\in\textnormal{L}^{1}(0,1)$, $q(x)\equiv 0$ in (2.1) with Dirichlet BCs, namely

[TABLE]

Then the spectral symbol given by a $3$-points central FD scheme is

[TABLE]

see [24, Theorem 10.5] (the proof works fine even if $p(x)\in C([0,1])$ is replaced by $p(x)\in L^{1}([0,1])$). It is not difficult to prove that the monotone rearrangement $\tilde{\omega}:[0,1)\to[0,+\infty)$ is such that

[TABLE]

and that

[TABLE]

where the upper bound on the largest eigenvalue is proven by the Gershgorin theorems. The limit relation in Remark 1 still holds,

[TABLE]

but on the other hand,

[TABLE]

In Figure 12 and Table 9 it is possible to see a summary of these last considerations.

Remark 2

Let us observe that Corollary 3.4.1 is still valid on every compact subset $K\Subset[0,1)$. For example, let $x_{0}\in(0,1)$: then

[TABLE]

Indeed, $\tilde{\omega}$ is absolute continuous on every compact subset which does not contain $x=1$.

5 Theoretical results

5.1 Proofs of results presented in Section 4

We summarize and generalize the results of Section 4.

Proposition 5.1.1

Let us consider the Sturm-Liouville eigenvalue problem

[TABLE]

such that

1. (i)

$p,p^{\prime},w,w^{\prime},q,(pw)^{\prime},(pw)^{\prime\prime}\in C([a,b])$; 2. (ii)

$p,w>0$; 3. (iii)

$\sigma_{1}^{2}+\sigma_{2}^{2}>0$, $\zeta_{1}^{2}+\zeta_{2}^{2}>0$.

Discretize the above Problem (5.1) by means of a numerical matrix method, and let $\prescript{b}{a}{\mathcal{L}_{\textnormal{BCs},p,q,w}^{(n,\eta)}}$ be the correspondent discrete operator, where $n$ is the mesh finesse parameter and $\eta$ is the order of approximation of the numerical method. Define $B=\int_{a}^{b}\sqrt{\frac{w(x)}{p(x)}}dx$. If:

1. (a)

for every fixed $k\in\mathbb{N}$,

[TABLE] 2. (b)

there exists $\omega:[a,b]\times[0,\pi]\to\mathbb{R}$, $\omega\in L^{1}\left([a,b]\times[0,\pi]\right)$ such that

[TABLE] 3. (c)

for every $\eta$, the monotone rearrangement $\tilde{\omega}_{\eta}$, as defined in (3.9), is such that

[TABLE]

Then, for every fixed $k\in\mathbb{N}$

[TABLE]

with $c_{k}$ independent of $\eta$. The inequality is strict, i.e., $c_{k}>0$ for every $k$ such that

[TABLE]

where $\prescript{B}{0}{\mathcal{L}_{\textnormal{BCs},V}}$ is the differential operator associated to the normal form (2.5) of Problem (5.1) by the Liouville transform (2.4), and $-\prescript{B}{0}{\Delta_{\textnormal{dir}}}$ is the (negative) Laplacian operator with Dirichlet BCs over the domain $[0,B]$.

In particular,

[TABLE]

where $\prescript{}{\alpha}{\tilde{\textnormal{{err}}}_{r,k}^{(n)}}$ is the analytical error defined in Definition 4.0.1.

Proof. By hypothesis, for every fixed $k$

[TABLE]

where $\lambda_{k}\left(\prescript{b}{a}{\mathcal{L}_{\textnormal{BCs},p(x),q(x),w(x)}}\right)$ are the eigenvalues of the continuous differential operator, and

[TABLE]

where $\prescript{B}{0}{\mathcal{L}_{\textnormal{BCs},V(y)}}$ is the differential operator associated to the normal form (2.5) of Problem (5.1) through the Liouville transform. Let us observe that for $V(y)\equiv 0$ and Dirichlet BCs, then

[TABLE]

Therefore, from item (c)

[TABLE]

Then it is immediate to prove that if

[TABLE]

then

[TABLE]

Moreover,

[TABLE]

and then $c_{k}\to 0$ as $k\to\infty$.

Corollary 5.1.2

If

[TABLE]

with $f_{\eta}(\theta)$ nonnegative, nondecreasing and such that $f_{\eta}(\theta)\sim\theta^{2}$ as $\theta\to 0$, then item (c) of Proposition 5.1.1 is satisfied.

Proof. From (3.10) and (3.9), for all $t\in\left[0,t_{0}\right]$, with

[TABLE]

we have that

[TABLE]

By the monotonicity of $f_{\eta}$, it holds that $\theta\to 0$ as $t\to 0$. For every $\epsilon>0$ there exists $\delta_{\epsilon}>0$ such that for every $t\in[0,\min\{t_{0};t_{0}/\delta_{\epsilon}\}]$ then $(1-\epsilon)\theta^{2}<f_{\eta}(\theta)<(1+\epsilon)\theta^{2}$, and so

[TABLE]

with

[TABLE]

So we have that

[TABLE]

By definition (3.9), $t\to 0$ as $x\to 0$ and then it holds that

[TABLE]

and the thesis follows.

Remark 3

The matrix methods of subsections 3.2, 3.3 satisfy the hypothesis of Proposition 5.1.1, see theorems 3.2.1, 3.3.1 and corollaries 3.2.1, 3.3.1. Therefore, in general, a uniform sampling of their spectral symbols does not provide an accurate approximation of the eigenvalues $\lambda_{k}\left(\prescript{b}{a}{\mathcal{L}_{\textnormal{BCs},p(x),q(x),w(x)}^{(n,\eta)}}\right)$ and $\lambda_{k}\left(\prescript{b}{a}{\mathcal{L}_{\textnormal{BCs},p(x),q(x),w(x)}}\right)$, in the sense of the relative error. See subsections 4.1,4.2 and 4.3 for numerical examples. On the other hand,

[TABLE]

if we exclude the outliers, see Corollary 3.4.1 and figures (1) and (12).

Theorem 5.1.1

Let us consider Problem 5.1 such that items (i)-(iv) of Proposition 5.1.1 are satisfied. Let $\prescript{b}{a}{\mathcal{L}_{\textnormal{BCs},p(x),q(x),w(x)}^{(n)}}$ be the discrete operator of Problem 5.1 obtained by means of any matrix method such that

1. (a)

[TABLE] 2. (b)

[TABLE] 3. (c)

[TABLE]

with $\omega\in\textnormal{L}^{\infty}\left([a,b]\times[0,\pi]\right)$ and real; 4. (d)

the monotone rearrangement $\tilde{\omega}:[0,1]\to[\min R_{\omega},\max R_{\omega}]$, defined as in equation (3.9), is piecewise Lipschitz.

Then

[TABLE]

Moreover, if for $n$ large enough there are not outliers as in Definition 3.1.5, i.e., if for every $n\geq N$

[TABLE]

then

[TABLE]

Proof. Since $\omega$ is bounded then $R_{\omega}$ is compact and $\min R_{\omega},\max R_{\omega}$ are well-defined. Let $x\in[0,1]$.

Without loss of generality, let us suppose moreover that $\lambda_{k}\left(\prescript{b}{a}{\mathcal{L}_{\textnormal{BCs},p(x),q(x),w(x)}^{(n)}}\right)$ are distinct for every $k$,

[TABLE]

so we can uniquely define the permutation index function $\sigma_{n}:\left\{1,\ldots,n\right\}\to\left\{1,\ldots,n\right\}$ such that the new index $\sigma_{n}(k)$ is reordered according to the distance of $(n+1)^{-2}\lambda_{k}\left(\prescript{b}{a}{\mathcal{L}_{\textnormal{BCs},p(x),q(x),w(x)}^{(n)}}\right)$ from $\tilde{\omega}(x)$, in ascending order. Namely, $\sigma_{n}(k)=j_{n}$ such that

[TABLE]

where

[TABLE]

and

[TABLE]

For every $n$, choose $k=k(n)$ such that $\sigma_{n}(k)=1$. By Theorem 3.1.1 we have that

[TABLE]

and by Theorem 3.4.1 it holds that

[TABLE]

Therefore,

[TABLE]

and then for every $\epsilon>0$

[TABLE]

The thesis follows at once.

Corollary 5.1.3

In the same hypothesis of Theorem 5.1.1, if

[TABLE]

then the numerical method used for the discretization of Problem 5.1 can not provide a relative uniform approximation of the eigenvalues of the continuous operator $\prescript{b}{a}{\mathcal{L}_{\textnormal{BCs},p,q,w}}$. In particular, it is enough that

[TABLE]

Proof. Immediate from Theorem 5.1.1.

Corollary 5.1.4

In the same hypothesis of Theorem 5.1.1, for every diffeomorphism $\tau:[a,b]\to[a,b]$, the numerical methods in subsections 3.2, 3.3 verify condition (5.3), for any fixed order of approximation $\eta$.

Proof. For every fixed $\eta$, it is easy to check that both the methods satisfy the hypothesis of Theorem 5.1.1, in particular $\tilde{\omega}_{\eta}\in C^{1}([0,1])$. By the regularity of the functions $f_{\eta}$, it holds that $\tilde{\omega}_{\eta}$ is twice differentiable almost everywhere and $\tilde{\omega}_{\eta}^{-1}=\phi_{\eta}$, with $\phi_{\eta}$ as in (3.10). At the point $x_{0}=\phi_{\eta}\left(t_{0}\right)$, with

[TABLE]

there exists $\epsilon>0$ such that

[TABLE]

It can be proved by direct computation, using the same approach as in Corollary 5.1.2, we skip the details.

Therefore, even if $\tilde{\omega}_{\eta}(x)=\frac{x^{2}\pi^{2}}{B^{2}}$ for every $x\in[0,x_{0}]$, there exists a non-negligible compact set $I\subseteq[x_{0},1]$ such that $\tilde{\omega}_{\eta}(x)\neq\frac{x^{2}\pi^{2}}{B^{2}}$.

On the contrary, if the order of approximation is let free to increase, then the maximum of the relative error seems to decrease to zero, provided that the discretization is made on an appropriate non-uniform grid. See subsections 4.2 and 4.3.

5.2 Further generalizations

As we stated at the beginning of the Introduction, the results of this section can be generalized to different kind of differential operators in any dimension. In particular, Theorem 5.1.1 basically relies just on Theorem 3.4.1 and the Weyl’s law for Sturm-Liouville operators. Therefore, if we know the asymptotic distribution of the discrete spectrum of the operator $\mathcal{L}$, we can measure the lower bound of the maximum relative approximation error by means of the spectral symbol which characterizes the matrix method used for the discretization.

As a plain example, consider the following Dirichlet boundary value problem on the unit square,

[TABLE]

The eigenvalues are given by $\lambda\left(-\Delta_{\textnormal{dir}}\right)=\pi^{2}\left(i^{2}+j^{2}\right)$ for $i,j\in\mathbb{N}$. If we indicate with $\tilde{\Omega}(x)$ the inverse cumulative distribution function for the eigenvalue distribution of $\left\{n^{-2}\lambda_{k}\left(-\Delta_{\textnormal{dir}}\right)\,:\,k=1,\ldots,n^{2}\right\}_{n}$, even if there is not a closed formula as in the one dimensional case for all $x\in[0,1]$, clearly we have that $\tilde{\Omega}(1)=2\pi^{2}$.

On the other hand, discretizing problem (5.4) by means of $3$-points FD as described in [25, Chapter 7.3], we get

[TABLE]

Even in the discrete case there is not a closed formula for $\tilde{\omega}$, but nevertheless we have that

[TABLE]

6 Conclusions

Although in the present paper for simplicity all the examples and the theory were developed mostly among the setting of regular Sturm-Liouville problems, a generalization to a wider class of differential operators in dimension $d\geq 1$ is feasible by means of the techniques presented in this paper.

Given a differential operator $\mathcal{L}$ discretized by means of a numerical scheme, the sampling of the spectral symbol calculated by the theory of GLT sequences does not provide in general an accurate approximation of the eigenvalues $\lambda_{k}\left(\mathcal{L}\right)$. Nevertheless, by Theorem 5.1.1 the knowledge of the spectral symbol provides to a numerical discretization scheme a necessary condition for the uniform spectral approximation of $\mathcal{L}$, in the sense of the relative error. It can measure how far the discretization method is from a uniform approximation of all the modes of the differential operator and this can be useful for engineering applications, see for example [29]. Moreover, the condition seems to become sufficient if the discretization method is paired with a suitable (non-uniform) grid and an increasing refinement of the order of approximation of the method. In light of this, it becomes a priority to devise new specific discretization schemes with mesh-dependent order of approximation which guarantee a good balance between convergence to zero of the relative spectral error and computational costs.

Finally, in reference with Theorem 3.4.1, corollaries 3.2.1, 3.3.1, and [35, Remark 15], since the spectral symbol is deeply related to the Weyl’s asymptotic distribution of the eigenvalues of the differential operator, this connection can be exploited to give better estimates of the Weyl function of generic elliptic operators on manifolds with bounded geometry, through a smart discretization of the elliptic operator itself and the analysis of the associated spectral symbol generated by the discretization scheme.

Appendix A Proofs of Subsection 3.2

Theorem 3.2.1 Proof. The proof of item (i) is long and technical, and we avoid to present it here. Let us just mention that it can be proved by a straightforward generalization of standard techniques, see [11, Theorem 1] and [27, 15]. About item (ii), let us preliminarily observe that in case of $p(x)\equiv 1$ and $\tau(x)=x$, then $\prescript{a}{b}{\mathcal{L}^{(n,\eta)}_{\textnormal{dir},p(x)\equiv 1}}$ is a symmetric Toeplitz matrix defined by

[TABLE]

where $d_{\eta,k}$ are the coefficients defined in (3.5), see [31, Corollary 2.2] and [34, Equation (27)]. Therefore,

[TABLE]

Let us now define an approximation of $\prescript{a}{b}{\mathcal{L}^{(n,\eta)}_{\textnormal{dir},p(\bar{x})}}$, namely

[TABLE]

$\prescript{a}{b}{\mathfrak{L}^{(n,\eta)}_{\textnormal{dir},p(\bar{x})}}$ is symmetric and $\left\|(n+1)^{-2}\prescript{a}{b}{\mathfrak{L}^{(n,\eta)}_{\textnormal{dir},p(\bar{x})}}\right\|\leq\max_{[a,b]}\left\{p(\tau(x))/\tau^{\prime}(x)^{2}\right\}\max_{[0,\pi]}\left\{|f_{\eta}(\theta)|\right\}$. Moreover, by [24, GLT 3-4 p. 160] and Proposition 3.1.1, it holds that

[TABLE]

By the regularity of $p(x)$ and $\tau(x)$, it is not difficult to prove that

[TABLE]

such that

[TABLE]

Combining now (A.1) and (A.2), by [24, Property S 4 p. 156] we get that

[TABLE]

Finally, it is immediate to check that $\left\|(n+1)^{-2}Q^{(n)}\right\|\leq cn^{-2}$ and again, by [24, properties GLT 3-4 p. 160] and [24, Property S 4 p. 156] we conclude that

[TABLE]

Corollary 3.2.1. Proof. $f_{\eta}(\theta)$ is obviously $C^{\infty}([0,\pi])$. Let us begin to prove that $f_{\eta}(\theta)\sim\theta^{2}$ as $\theta\to 0$ and that it is monotone nonnegative on $[0,\pi]$. By the Taylor expansion at $\theta=0$ we get

[TABLE]

Moreover, let us observe that

[TABLE]

Define then

[TABLE]

It is immediate to check that $a_{1}\geq a_{2}\geq\ldots\geq a_{n}>0$ and that

[TABLE]

By [5, Theorem 1] we can conclude that $g(\psi)>0$ on $(0,\pi)$ and then $f_{\eta}^{\prime}(\theta)>0$ on $(0,\pi)$. Since $f_{\eta}(0)=0$, we deduce that $f_{\eta}(\theta)\geq 0$ on $[0,\pi]$.

The second part of the thesis is an immediate consequence of identities (3.5). Indeed, for every fixed $k\geq 1$ it holds that

[TABLE]

and

[TABLE]

Since $\sum_{k=1}^{\infty}\frac{(-1)^{k}}{k^{2}}=-\frac{\pi^{2}}{12}$, we conclude that

[TABLE]

Therefore, for every fixed $\theta\in[0,\pi]$,

[TABLE]

being $\left\{\pi^{2}/3\right\}\cup\left\{(-1)^{k}4/k^{2}\right\}_{k\geq 1}$ the Fourier coefficients of $\theta^{2}$ on $[0,\pi]$, and the convergence is uniform.