A combinational efficient algorithm for elliptic eigenvalue problems

Abstract

In this paper, we present a combinational efficient algorithm for solving second-order elliptic ordinary differential equation eigenvalue problems on four kinds of meshes. Firstly, a combinational scheme is proposed by leveraging the linear finite element scheme and the second-order finite difference scheme based on mass lumping method. An efficient algorithm is constructed by using the combinational scheme with a constant combined parameter. Then, we improve the efficient algorithm by computing the quasi-optimal combined parameters for different eigenvalues. Finally, numerical experiments tested on both uniform and nonuniform meshes are shown to illustrate the efficiency of our algorithms.

Introduction

Studies on the eigenvalue problem have been a hot topic for many years. Problems of this kind have many applications in both applied mathematics and engineering, such as quantum mechanics, eigenfaces in image processing, vibration analysis, etc. There are many methods, such as the linear finite element method or the second-order finite difference method which always involve inevitably massive meshes and large-scale computation.^1–3

To overcome the challenge of large-scale computation, one often exploits distributed or parallel algorithms. Although parallel algorithm is a good choice, it may not avoid expensive computational cost. Alternatively, one can design efficient algorithm with high accuracy for specific problem which can avoid large-scale computation. This is just our theme in this paper. At present, there are many researches on constructing high-accuracy algorithms for solving eigenvalue problems.^4–11 In this paper, we will give an efficient algorithm for the ODE eigenvalue problems (1) with variable coefficient on four kinds of meshes by leveraging the linear finite element method and the second-order finite difference method. Especially, we design a second-order finite difference scheme with mass lumping method in the text and then compute the quasi-optimal combined parameters for different eigenvalues to improve our efficient algorithm. The second-order elliptic ODE eigenvalue problem can be described as follows where $p\in C1(I)$ , $p(x)\ge \underset{x\in I}{min}p(x)>0$ , and $q\in C1(I)\overline{I}=[a,b]$ .

Throughout this paper, ( f, g) denotes ${\int }_0^{1}f(t)g(t)\mathrm{dt}$ .

As is well known, the Rayleigh quotient of ODE eigen-equation $\mathrm{Lu}=\lambda u$ is defined as

$R[v]:=\frac{(\mathrm{Lv},v)}{(v,v)}=\frac{a(v,v)}{(v,v)}$

(2)

where a(v, v) is a weak form of the original energy form. For $(\mathrm{Lv},v)>0$ , eigenvalues satisfy the following min–max principle

$\lambda 1=\underset{v\in H}{min}R[v],\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\lambda k=\underset{Ek}{min}\underset{v\in Ek}{max}R[v]$

(3)

$0<\lambda 1\le \lambda 2\le \dots \le \lambda k\to \infty$

(4)

where H is called admissible space and E_k is a k-dimensional subspace of H.

For a given finite-dimensional approximation subspace V^h in H, the Petrov-Galerkin discrete Rayleigh quotient can be expressed as

$R[vh]:=\frac{a(vh,vh)}{(vh,vh)}.$

(5)

The min–max principle also holds for the discrete form

${\lambda }_k^h=\underset{Vk}{min}\underset{vh\in Vk}{max}R[vh],$

(6)

$0<{\lambda }_1^h\le {\lambda }_2^h\le \dots \le {\lambda }_k^h\le \dots \le {\lambda }_{\mathrm{DOF}}^h,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{DOF}=\dim (Vh),$

(7)

Here, V_k ranges over all k-dimensional subspace in V^h.

For a given basis $\{\phi 1,\dots ,\phi N\}$ in V^h and $vh=\sum qj\phi j$ , the discrete Rayleigh quotient becomes

$R[vh]=\frac{qTKhq}{qTMhq},$

(8)

where K^h and M^h are stiffness matrix and mass matrix, respectively. K^h is symmetric positive semidefinite and M^h is symmetric positive definite. Then, this discrete Rayleigh quotient form (8) is equivalent to the following generalized matrix eigenvalue problem

$KhQh=\lambda hMhQh,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}Kh=(a(\phi j,\phi k)),\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}Mh=((\phi j,\phi k)).$

(9)

Sun ⁴ presented two fourth-order accuracy schemes for solving the elliptic eigenvalue problems on the uniform mesh. Furthermore, Sun ⁵ gave some sixth-order schemes by linear combination of the fourth-order accuracy schemes with proper combination coefficients. But all those high-order accuracy schemes are just applicable to uniform meshes. In Zhang et al., ⁷ they generalized Sun’s idea to nonuniform meshes, gave two fourth-order accuracy schemes and two sixth-order accuracy schemes. But all these high-order accuracy schemes only aim at eigenvalue problem (1) with constant coefficient in the eigen-equation, namely p and q are both constants.

In this paper, we aim at the case with variable coefficient in the eigen-equation (1). By using the techniques of combinational method and mass lumping method,^2,12 we construct an efficient algorithm. It is applicable to general eigenvalue problem (1) with variable coefficient q or p or both in the eigen-equation. That is,

Case I

Case II

$\mathrm{Lu}=(-\mathrm{pu}\text{'})\text{'}=\lambda u\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}x\in (a,b),(u2)\text{'}|\partial \Omega =0$

(11)

Case III

Here, $(u2)\text{'}|\partial \Omega =0$ stands for the boundary conditions of the following cases:

• ♦ Direchlet condition: $u(a)=u(b)=0$ ;

• ♦ Niemann condition: $u\text{'}(a)=u\text{'}(b)=0$ ;

• ♦ Mixed condition: $u(a)=u\text{'}(b)=0$ or $u\text{'}(a)=u(b)=0$ .

Although our efficient algorithm and its improved algorithm may not guarantee high-order accuracy, they both converge much faster for currently numerical experiments than the single linear finite element method or the corresponding finite difference method on the same mesh. Some numerical experiments tested on both uniform meshes and nonuniform meshes are given to illustrate the computational cost measured by the CPU time and grid size. For efficiency, all the matrices are assembled with sparse storage in the algorithm implementation. Typically, all the numerical experiments about the eigenvalue problems (10) to (12) are conducted when the boundary condition is $u(0)=u(1)=0$ . All the numerical experiments are tested on both the uniform and nonuniform meshes.

The rest of the paper is organised as follows. In An efficient algorithm for ODE eigenvalue problems, we propose an efficient algorithm by combining the linear finite element method and the corresponding second-order finite difference method with a constant combined parameter. In An improved efficient algorithm, we give an improved algorithm for the efficient algorithm given in the above section by computing the quasi-optimal combined parameters for different eigenvalues. In the following section, we give some numerical experiments to test our algorithms and in finally we give some conclusions.

An efficient algorithm for ODE eigenvalue problems

In this section, an efficient algorithm is constructed by taking a weighted mean combination of the linear finite element method and corresponding second-order finite difference method. Especially, we design the second-order finite difference scheme with mass lumping method.

Firstly, we suppose u is the eigenfunction of problem (1), s is the linear interpolation of u on the mesh

$a=x0<x1<\dots xn=b.$

(13)

Henceforth, h_i denotes $xi-xi-1,\mathrm{mesh}[h]$ denotes such a mesh where h is the average step length of the mesh, namely, $h=\frac{\sum _{i}hi}{n}$ . We define vector U as $U(i)=u(xi)$ . Then, the eigenvalue expression with finite element method is as below

$\lambda F:=\frac{(-\mathrm{ps}\text{'},s\text{'})+(\mathrm{qs},s)}{(s,s)}.$

(14)

It is equivalent to

$\lambda F=\frac{UTKhU}{UTMhU}.$

(15)

where $(-\mathrm{ps}\text{'},s\text{'})+(\mathrm{qs},s)=(UTAhU)+(UTChU)=(UTKhU)$ , $Kh=Ah+Ch$ , and $(s,s)=UTMhU$ . We call equation (15) as finite element method (FEM) scheme.

Just like in literature,^4,7 we can define the second-order finite difference scheme corresponding to the above linear finite element scheme as below.

$\stackrel{\wedge }{\lambda }D:=\frac{(-\mathrm{ps}\text{'},s\text{'})+(\mathrm{qu},u)}{(u,u)}$

(16)

However, we can not directly use the Rayleigh-quotient forms (14) and (17) to compute the numerical eigenvalues just like in literature^4,7 because we do not know the genuine eigenfunctions. It is practically difficult to find some kind of meshes satisfying $UTU=(u,u)$ . In this case, we can only compute the numerical eigenvalues $(\lambda D)$ by solving a generalized matrix eigenvalue problem in the form of equation (9). So, we propose a modified second-order finite difference scheme corresponding to the above linear finite element scheme as below

$\lambda D:=\frac{(-\mathrm{ps}\text{'},s\text{'})+〈\mathrm{qu},u〉}{〈u,u〉}$

(17)

where $〈·,·〉$ denote the approximate inner product of $(·,·)$ .

It is equivalent to

$\lambda D=\frac{UT\tilde{K}hU}{UTDhU}$

(18)

where $(-\mathrm{ps}\text{'},s\text{'})+〈\mathrm{qu},u〉=(UTAhU)+(UTEhU)=(UT\tilde{K}hU)$ , $\tilde{K}h=Ah+Eh$ ; here, E^h is the mass lumping matrix obtained from C^h in the FEM scheme (15). $〈u,u〉=UTDhU$ , where D^h is the mass lumping matrix obtained from M^h in the FEM scheme (15). We call equation (18) as FDM scheme.

In the two schemes of equations (14) and (16), the discrete form for p(x) is $pi=p((xi+xi-1)/2)$ .

From Strang and Fix, ¹ we know the numerical eigenvalues computed with the finite element scheme are not smaller than the genuine eigenvalues of (1). Conversely, according to George and Wolgang, ¹³ we know the FDM scheme on the uniform mesh leads to a lower bound of genuine eigenvalues under some constraints. Based on these two facts, we know any combinational scheme in the form of

$\lambda C=\alpha \lambda D+(1-\alpha )\lambda F,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{for}\alpha >0$

(19)

can obtain a better approximate solution for equation (1) than individual linear finite element scheme λ_F or the corresponding second-order finite difference scheme λ_D on the uniform mesh under some constraints.

How is it for the combinational scheme on the nonuniform mesh? Next, we choose proper combination coefficient α to construct efficient algorithm for solving the eigenproblems (10) to (12).

Theoretical analysis

We first take the arithmetic mean of equations (14) and (16) to construct our efficient algorithm.

We take a mean combination of the FEM scheme λ_F and our modified FDM scheme λ_D for the eigenvalue problems with variable coefficient q or p or both on the uniform and nonuniform meshes.

Next, we present our algorithm.

Algorithm structure

In our algorithm, we assemble K^h with two parts A^h and C^h, assemble $\tilde{K}h$ with A^h and E^h, where E^h is the lumping matrix of C^h with proper boundary treatment. For the eigenproblem of equation (11) (Case I), the mass part C^h in K^h and E^h in $\tilde{K}h$ is equal to zero. Specifically, $Kh=Ah+Ch,Ah=(ai,j)=(\phi \text{'}j,\phi \text{'}i),\hspace{0.5em}Ch=(ci,j)=(q\phi j,\phi i)$ , where $\phi i,i=1,\dots$ is the basis of s. $\tilde{K}h=Ah+Eh$ . $Mh=(bi,j)=(\phi j,\phi i)$ , D^h is the lumping matrix of M^h with proper boundary treatment. For efficiency and memory limit, we assemble all the matrices with sparse storage. Then we use ARPACK software to solve the corresponding generalized eigenvalue problem.

n denotes grid size of the interval $[a,b]$ , m denotes the number of eigenvalues to compute, and sp denotes the boundary condition. Our algorithm structure is described as Algorithm 1.

Algorithm 1 MCM (n, m)

Require: Two integer $n\ge 0,m\ge 0$ , Real a, b, Bool sp.

Ensure: The numerical eigenvalues of $\lambda 1\dots \lambda m$ .

1: compute $\mathrm{ld}=\mathrm{FDM}(n,m,\mathrm{sp})$ and $\mathrm{lf}=\mathrm{FEM}(n,m,\mathrm{sp})$ concurrently;

2: $\mathrm{la}=(\mathrm{lf}+\mathrm{ld})/2$ ;

Return la.

MCM: method of computing eigenvalues by using combinational method.

For a given initial value n = n₀, one can iterate the MCM algorithm to compute eigenvalues by updating n = 2n until converge. The complete iterative structure is described as Algorithm 2.

Algorithm 2 ITMCM(n, m, eps,  $\lambda 0$ )

Require: Two integers $n0\ge 0,m\ge 0$ , Real number eps and Real vector $\lambda 0$ .

Ensure: The numerical eigenvalues of $\lambda 1\dots \lambda m$ , Meshes number n.

1: $\mathrm{la}=\mathrm{MCM}(n,m);$

2: $\mathrm{er}=\mathrm{la}-\lambda 0;$ $\mathrm{error}=\mathrm{norm}(\mathrm{er},1);$

3: while (error > eps)

  $\lambda 0=\mathrm{la};$

  $n=2n;$

  $\mathrm{la}=\mathrm{MCM}(n,m);$

  $\mathrm{er}=\mathrm{la}-\lambda 0;$

  $\mathrm{error}=\mathrm{norm}(\mathrm{er},1);$

Return la, n.

ITMCM: Iterate the MCM algorithm until converge.

An improved efficient algorithm

In this section, we give an improved algorithm for the algorithm given in the above section by computing the quasi-optimal combined parameters for every eigenvalue.

As we have known, eigenvalues obtained by the FEM scheme are not smaller than the exact ones while the FDM scheme are not bigger than those. The errors of the two schemes on $\mathrm{mesh}[h]$ and $\mathrm{mesh}[2h]$ can be expressed as

${\lambda }_F^h-\lambda =C1h2+o(h2),{\lambda }_F^{2h}-\lambda =C1(2h)2+o(h2),{\lambda }_D^h-\lambda =C2h2+o(h2),{\lambda }_D^{2h}-\lambda =C2(2h)2+o(h2)$

Then we can compute the quasi-optimal combined parameters as below

$\mu =-({\lambda }_D^h-{\lambda }_D^{2h})/({\lambda }_F^h-{\lambda }_F^{2h}),\alpha =1/(1+\mu )$

Here, μ and α depend on both the eigenvalues and the mesh. In this way, we can improve the error

$\lambda C=\alpha \lambda D+(1-\alpha )\lambda F-\lambda =o(h2)$

(higher than second order) theoretically. This improved scheme λ_C means higher accuracy approximation of the true eigenvalues.

For a given initial value n = n₀, one can iterate this algorithm to compute eigenvalues by updating n = 2n until converge. Firstly, we compute $\mathrm{ld}1=\mathrm{FDM}(n,m),\mathrm{lf}1=\mathrm{FEM}(n,m),\mathrm{ld}2=\mathrm{FDM}(2n,m),\mathrm{lf}2=\mathrm{FEM}(2n,m)$ concurrently. Then, define

$\mu :=-(\mathrm{ld}2-\mathrm{ld}1)./(\mathrm{lf}2-\mathrm{lf}1),\alpha :=1./(1+\mu )$

Here, α and μ are both vectors, α is just the vector of quasi-optimal combined parameters for every eigenvalue of equation (19). Lastly, we iterate the above steps and update α until converge. The improved algorithm structure is described as Algorithm 3.

Algorithm 3 QOCM (n, m, eps,  $\lambda 0$ )

Require: Two integers $n0\ge 0,m\ge 0,\mathrm{sp}$ , Real eps, Vector λ₀.

Ensure: The numerical eigenvalues of $\lambda 1\dots \lambda m$ , Meshes number n.

1: compute concurrently: $\mathrm{ld}1=\mathrm{FDM}(n,m,\mathrm{sp})$ , $\mathrm{lf}1=\mathrm{FEM}(n,m,\mathrm{sp})$ , $\mathrm{ld}2=\mathrm{FDM}(2n,m,\mathrm{sp})$ ,

$\mathrm{lf}2=\mathrm{FEM}(2n,m,\mathrm{sp})$ ;

2: $\mathrm{alf}:=-(\mathrm{ld}2-\mathrm{ld}1)./(\mathrm{lf}2-\mathrm{lf}1);$

  $a=1./(1+\mathrm{alf});$

3: For $i=1,\dots ,m$

  $\mathrm{lc}(i)=(1-a(i))*\mathrm{lf}2(i)+a(i)*\mathrm{ld}2(i);$

4: $\mathrm{er}=\mathrm{lc}-\lambda 0;$

  $\mathrm{error}=\mathrm{norm}(\mathrm{er},1);$

5: while ( $\mathrm{error}>\mathrm{eps}$ )

  $\lambda 0=\mathrm{lc};$

  $\mathrm{lf}1=\mathrm{lf}2;$

  $\mathrm{ld}1=\mathrm{ld}2;$

  $n=2n;$

 compute concurrently: $\mathrm{ld}2=\mathrm{FDM}(n,m);$ $\mathrm{lf}2=\mathrm{FEM}(n,m)$

  $\mathrm{alf}:=-(\mathrm{ld}2-\mathrm{ld}1)./(\mathrm{lf}2-\mathrm{lf}1);$

  $a=1./(1+\mathrm{alf});$

 For $i=1,\dots ,m$

  $\mathrm{lc}(i)=(1-a(i))*\mathrm{lf}2(i)+a(i)*\mathrm{ld}2(i);$

  $\mathrm{er}=\mathrm{lc}-\lambda 0;$

  $\mathrm{error}=\mathrm{norm}(\mathrm{er},1);$

Return $\mathrm{lc},$ n.

QOCM:Algorithm by using quai-optimal combinational method.

Algorithm 4 FEM (n,m,sp)

Require: Two integers $n\ge 0,m\ge 0$ , Real $x1,x2$ , Bool $\mathrm{sp}$ .

Ensure: The numerical eigenvalues of $\lambda 1\dots \lambda m$ .

1: If $\mathrm{sp}=0$

  $n=n-1;$

 If $\mathrm{sp}=1$

  $n=n+1;$

 2: $\mathrm{col}=\mathrm{zeros}(3*n-2,1)$ ; $\mathrm{lin}=\mathrm{col};$

   $A=\mathrm{double}(\mathrm{zeros}(3*n-2,1));$

   $B=A;C=A;$

 3: Obtain the Interval subdivision of $[a,b].$

   $[x,H]=\mathrm{mesh}(n,x1,x2);$

 4: Assign $\mathrm{col},\mathrm{lin};$

  Left boundary treatment;

 5: For $i=1,\dots ,n-2$ .

  Assign C by computing $ci,i-1,ci,i,ci+1,i$ ,

  Assign A by computing $ai,i-1,ai,i,ai+1,i$ ,

  Assign B by computing $bi,i-1,bi,i,bi+1,i$ ,

 6: Right boundary treatment;

 7: $\mathrm{AF}=A+C;$

  K = sparse $(\mathrm{lin},\mathrm{col},\mathrm{AF},n-1,n-1);$

  M = sparse $(\mathrm{lin},\mathrm{col},B,n-1,n-1);$

 8: $\mathrm{lf}$  = eigs $(K,M,m,\text{'}\mathrm{sm}\text{'});$

   $\mathrm{lf}$  = sort $(\mathrm{lf})$ ;

Return $\mathrm{lf}$ .

Algorithm to compute eigenvalues by using the FEM scheme.

Algorithm 5 FDM(n,m,sp)

Require: Two integers $n\ge 0,m\ge 0$ , Real $x1,x2$ , Bool $\mathrm{sp}$ .

Ensure: The numerical eigenvalues of $\lambda 1\dots \lambda m$ .

1: Update $n\mathbf{according}\mathrm{sp};$

2: Initialize vectors $A,E,D;\mathrm{col},\mathrm{lin};$

3: Obtain the Interval subdivision of $[x1,x2]$ .

  $[x,H]=\mathrm{mesh}(n,x1,x2);$

4: Assign $\mathrm{col},\mathrm{lin};$

 Left boundary treatment;

5: For $i=1,\dots ,n-2$

 Assign E by computing $\mathrm{cd}i,i,$

 Assign A by computing $ai,i-1,ai,i,ai+1,i,$

 Assign D by computing $di,i,$

6: Right boundary treatment;

7: $\tilde{K}=A+E;$

  $\tilde{K}$  = sparse $(\mathrm{lin},\mathrm{col},\tilde{K},n-1,n-1);$

  $\tilde{M}$  = sparse $(\mathrm{lin},\mathrm{col},D,n-1,n-1);$

8: $\mathrm{ld}$  = eigs $(\tilde{K},\tilde{M},m,$ ‘sm’ );

  $\mathrm{ld}$  = sort $(\mathrm{ld});$

Return $\mathrm{ld}$ .

Algorithm to compute eigenvalues by using the FDM scheme

Numerical experiments

In this section, we give some numerical experiments to test the efficient algorithm $\mathrm{MCM}(n,m)$ and its improved algorithm $\mathrm{QOCM}(n,m)$ , respectively, for eigenvalue problems (10) to (12). Then, compare with the algorithms $\mathrm{FEM}(n,m),\mathrm{FDM}(n,m)$ computed with finite element scheme and the corresponding finite difference scheme. Typically, we conduct all the numerical experiments about the eigenvalue problems (10) to (12) when the boundary condition is $u(0)=u(1)=0$ . Calling the algorithms $\mathrm{FEM}(n,m),\mathrm{FDM}(n,m),\hspace{1em}\mathrm{MCM}(n,m)\hspace{1em}\text{and}\hspace{1em}\mathrm{QOCM}(n,m)$ , we respectively compute the smallest four eigenvalues $\mathrm{ld},\mathrm{lf},\mathrm{la},\mathrm{lc}$ for case I: $q(x)=2\cos (2x)$ , case II: $p(x)=1+ex$ , and case III: $q(x)=2\cos (2x)\hspace{1em}\text{and}\hspace{1em}p(x)=1+ex$ , on both uniform mesh and nonuniform mesh. We take $x\in [0,\pi /2]$ . Give the symbol description as follows:

$\Lambda 0$ initial guess of the eigenvalues vector for the smallest four eigenvalues: $\Lambda 0=\overrightarrow{0}$ ;

$\mathrm{er}d$ the iterative error vector: $\mathrm{er}d[2n]=\mathrm{ld}[2n]-\mathrm{ld}[n]$ , specially, $\mathrm{er}d[n]=\mathrm{ld}[n]-\mathrm{ld}[0],\mathrm{ld}[0]=\Lambda 0,\mathrm{ld}[n]$ is eigenvalues vector of $\mathrm{ld}$ when the mesh size is n. In the following figures referred to iterative error, we express $\mathrm{er}d[2n]$ as the iterative error from $n\_=n\hspace{1em}\text{to}\hspace{1em}n\_=2n,\mathrm{er}d$ as the iterative error when converge, so it is with $\mathrm{er}f,\mathrm{er}a,\mathrm{er}c$ ;

$\Vert \mathrm{ld}\Vert$ 1-norm of the iterate error $\mathrm{er}d$ : $\Vert \mathrm{ld}\Vert \hspace{1em}\text{denote}{\max }_{i=1}^4|\mathrm{er}d(i)|$ of the smallest four eigenvalues, the iteration terminate when $\Vert \mathrm{ld}\Vert$ is less than $\mathrm{eps}=1.0e-4$ . So, it is with $\Vert \mathrm{lf}\Vert ,\Vert \mathrm{la}\Vert ,\Vert \mathrm{lc}\Vert$ .

In all the experiments, we set grid size starter $n0=20$ and refine the mesh with grid size $n=2i*n0$ in the $i-\mathrm{th}$ iteration.

With the same accuracy requirement as above, we test the gird sizes n and time cost of the four algorithms when converge. All the example are tested on four kinds of meshes as follows

No.1:  $xj=\pi (1-\cos (j\pi h))/4,h=1/n,j=1,\dots ,n$

No.2:  $xj=\pi \tan (j\pi h/4)/2,h=1/n,j=1,\dots ,n$

No.3:  $xj=\pi (\mathrm{jh})2/2,h=1/n,j=1,\dots ,n$

No.4:  $xj=\pi (\mathrm{jh})/2,h=1/n,j=1,\dots ,n$

Case I: $q(x)=2\cos (2x)$

In this case, the eigen-equation is just the Mathieu equation. The smallest four numerical eigenvalues are $\Lambda =(3.9170,16.0330,36.0143,64.0079)$ .

Firstly, we compare the efficient algorithm $\mathrm{MCM}$ to algorithms FEM and FDM by observing $\mathrm{er}d[80],\mathrm{er}f[80],\mathrm{er}a[80]$ , namely, the iterative errors from $n\_=40$ to $n\_=80$ for the four algorithms tested on No.1–No.4 meshes as shown in Figure 1 and Table 1.

Figure 1.

The iterative error vectors $\mathrm{er}d[80],\mathrm{er}f[80],\mathrm{er}a[80]$ for $q(x)=2\cos (2x)$ from $n\_=40$ to $n\_=80$ , tested on No.1–No.4 meshes, x-axis denote the eigenvalues indices.

Table 1.

The iterative error vectors $\mathrm{er}f[80],\mathrm{er}d[80],\mathrm{er}a[80]$ for $q(x)=2\cos (2x)$ from $n\_=40$ to $n\_=80$ , tested on No.1–No.4 meshes, $\mathrm{er}x(i)$ denote the $i-\mathrm{th}$ component of vector $\mathrm{er}x[80]$ .

Mesh no.	No.1	No.2	No.3	No.4
$\mathrm{er}f(1)$	−0.0032	−0.0017	−0.0031	−0.0017
$\mathrm{er}f(2)$	−0.430	−0.0282	−0.0494	−0.0254
$\mathrm{er}f(3)$	−0.2107	−0.1435	−0.2497	−0.1267
$\mathrm{er}f(4)$	−0.6593	−0.4543	−0.7885	−0.3991
$\mathrm{er}d(1)$	0.0007	0.0032	0.0044	0.0022
$\mathrm{er}d(2)$	0.0298	0.0318	0.0504	0.0252
$\mathrm{er}d(3)$	0.1816	0.1501	0.2507	0.1251
$\mathrm{er}d(4)$	0.6066	0.4642	0.7905	0.3937
$\mathrm{er}a(1)$	−0.0013	0.0008	0.0006	0.0003
$\mathrm{er}a(2)$	−0.0066	0.0020	0.0005	0.0002
$\mathrm{er}a(3)$	−0.0145	0.0037	0.0005	−0.0000
$\mathrm{er}a(4)$	−0.0263	0.0056	0.0010	−0.0014

Figure 1 and Table 1 show that when the mesh is refined to $n=80$ , the algorithm $\mathrm{MCM}(n,m)$ has much higher accuracy than algorithms $\mathrm{FEM}(n,m)$ and $\mathrm{FDM}(n,m)$ on all the four kinds of meshes.

When all the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ , the error vectors are shown in Figure 2, time cost is shown in Figure 3, and the final refined grid sizes n are shown in Table 2. All experiments are tested on No.1–No.4 meshes.

Figure 2.

The iterative errors $\mathrm{er}d,\mathrm{er}f,\mathrm{er}a,\mathrm{er}c$ when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $q(x)=2\cos (2x)$ , tested on No.1–No.4 meshes, x-axis denote the eigenvalues indices.

Figure 3.

The time cost when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $q(x)=2\cos (2x)$ , x-axis denote the algorithms: 1 denotes $\mathrm{QOCM}$ , 2 denotes $\mathrm{MCM}$ , 3 denotes $\mathrm{FEM}$ , and 4 denotes $\mathrm{FDM}$ .

Table 2.

The final refined grid sizes of when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $q(x)=2\cos (2x)$ tested on No.1–No.4 meshes.

Mesh no.	QOCM	MCM	FEM	FDM
1	320	2560	10,240	10,240
2	320	1280	10,240	10,240
3	160	640	10,240	10,240
4	320	320	10,240	10,240

Case II: $p(x)=1+ex$

In this case, the smallest four numerical eigenvalues are $\Lambda =(13.18,51.45,115.2,204.4)$ . We still firstly compare the efficient algorithm $\mathrm{MCM}$ to algorithms $\mathrm{FEM}$ and $\mathrm{FDM}$ by observing $\mathrm{er}d[80],\mathrm{er}f[80],\mathrm{er}a[80]$ , namely, the iterative errors from $n\_=40$ to $n\_=80$ for the four algorithms tested on No.1–No.4 meshes as shown in Figure 4 and Table 3.

Figure 4.

The iterative error vectors $\mathrm{er}d[80],\mathrm{er}f[80],\mathrm{er}a[80]$ for $p(x)=1+ex$ from $n\_=40$ to $n\_=80$ , tested on NO.1–NO.4 meshes, x-axis denote the eigenvalues indices.

Table 3.

The iterative error vectors $\mathrm{er}f[80],\mathrm{er}d[80],\mathrm{er}a[80]$ for $p(x)=1+ex$ from $n\_=40$ to $n\_=80$ , tested on No.1–No.4 meshes, $\mathrm{er}x(i)$ denote the $i-\mathrm{th}$ component of vector $\mathrm{er}x[80]$ .

Mesh no.	No.1	No.2	No.3	No.4
$\mathrm{er}f(1)$	−0.0110	−0.0043	−0.0071	−0.0055
$\mathrm{er}f(2)$	−0.1458	−0.0763	−0.1229	−0.0852
$\mathrm{er}f(3)$	−0.7078	−0.3957	−0.6315	−0.4290
$\mathrm{er}f(4)$	−2.2039	−1.2631	−2.0079	−1.3545
$\mathrm{er}d(1)$	0.0014	0.0069	0.0087	0.0055
$\mathrm{er}d(2)$	0.1067	0.0880	0.1297	0.0849
$\mathrm{er}d(3)$	0.6187	0.4212	0.6459	0.4271
$\mathrm{er}d(4)$	2.0409	1.3019	2.0275	1.3454
$\mathrm{er}a(1)$	−0.0048	0.0004	0.0018	0.0001
$\mathrm{er}a(2)$	−0.0196	−0.0201	0.0034	−0.0524
$\mathrm{er}a(3)$	−0.0448	−0.1139	0.0072	−0.2506
$\mathrm{er}a(4)$	−0.0823	−0.3643	0.0098	−0.7411

Figure 4 and Table 3 show that when mesh is refined to $n=80$ , the algorithm $\mathrm{MCM}(n,m)$ has much higher accuracy than algorithms $\mathrm{FEM}(n,m)$ and $\mathrm{FDM}(n,m)$ on all the four kinds of meshes. Moreover, the uniform mesh does not always have advantage over nonuniform mesh. Case II with variable coefficient of $p(x)$ in the eigen-equation is just a counter-example.

When all the four algorithms iteration terminate within tolerance $\mathrm{eps}=1.0e-4$ , the error vectors are as shown in Figure 5, time cost is shown in Figure 6, and the final refined grid sizes n are shown in Table 4. All the experiments are tested on No.1–No.4 meshes.

Figure 5.

The iterative errors $\mathrm{er}d,\mathrm{er}f,\mathrm{er}a,\mathrm{er}c$ when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $p(x)=1+ex$ , tested on No.1–No.4 meshes, x-axis denote the eigenvalues indices.

Figure 6.

The time cost when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $p(x)=1+ex$ , x-axis denote the algorithms: 1 denotes $\mathrm{QOCM}$ , 2 denotes $\mathrm{MCM}$ , 3 denotes $\mathrm{FEM}$ , and 4 denotes $\mathrm{FDM}$ .

Table 4.

The final refined grid sizes when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $p(x)=1+ex$ .

Mesh no.	QOCM	MCM	FEM	FDM
1	320	5120	20,480	20,480
2	2560	1280	20,480	20,480
3	320	2560	20,480	20,480
4	2560	2560	20,480	20,480

Case III: $q(x)=2\cos (2x)$ and $p(x)=1+ex$

In this case, The smallest four numerical eigenvalues are $\Lambda =(13.4780,51.6827,115.4283,204.6768)$ .

We firstly compare the efficient algorithm $\mathrm{MCM}$ to algorithms $\mathrm{FEM}$ and $\mathrm{FDM}$ by observing $\mathrm{er}d[80],\mathrm{er}f[80],\mathrm{er}a[80]$ , namely, the iterative errors from $n\_=40$ to $n\_=80$ for the four algorithms tested on No.1–No.4 meshes as shown in Figure 7 and Table 5.

Figure 7.

The iterative error vectors $\mathrm{er}d[80],\mathrm{er}f[80],\mathrm{er}a[80]$ for $q(x)=2\cos (2x)$ and $p(x)=1+ex$ from $n\_=40$ to $n\_=80$ , tested on No.1–No.4 meshes, x-axis denote the eigenvalues indices.

Table 5.

The iterative error vectors $\mathrm{er}f[80],\mathrm{er}d[80],\mathrm{er}a[80]$ , for $q(x)=2\cos (2x)$ and $p(x)=1+ex$ from $n\_=40$ to $n\_=80$ , tested on No.1–No.4 meshes, $\mathrm{er}x(i)$ denote the $i-\mathrm{th}$ component of vector $\mathrm{er}x[80]$ .

Mesh no.	No.1	No.2	No. 3	No.4
$\mathrm{er}f(1)$	−0.0106	−0.0042	−0.0072	−0.0052
$\mathrm{er}f(2)$	−0.1449	−0.0764	−0.1245	−0.0841
$\mathrm{er}f(3)$	−0.7060	−0.3963	−0.6355	−0.4267
$\mathrm{er}f(4)$	−2.2008	−1.2643	−2.0152	−1.3506
$\mathrm{er}d(1)$	0.0018	0.0077	0.0108	0.0054
$\mathrm{er}d(2)$	0.1062	0.0888	0.1326	0.0841
$\mathrm{er}d(3)$	0.6172	0.4224	0.6513	0.4251
$\mathrm{er}d(4)$	2.0381	1.3038	2.0364	1.3415
$\mathrm{er}a(1)$	−0.0044	0.0018	0.0018	0.0001
$\mathrm{er}a(2)$	−0.0194	0.0062	0.0041	−0.0000
$\mathrm{er}a(3)$	−0.0444	0.0131	0.0079	−0.0008
$\mathrm{er}a(4)$	−0.0814	0.0198	0.0106	−0.0046

Figure 7 and Table 5 show that when mesh is refined to $n=80$ , the algorithm $\mathrm{MCM}(n,m)$ has much higher accuracy than algorithms $\mathrm{FEM}(n,m)$ and $\mathrm{FDM}(n,m)$ on all the four kinds of meshes.

When all the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ , the error vectors are as shown in Figure 8, time cost is shown in Figure 9, and the final refined grid sizes n are shown in Table 6. All of experiments are tested on No.1–No.4 meshes.

Figure 8.

The iterative errors $\mathrm{er}d,\mathrm{er}f,\mathrm{er}a,\mathrm{er}c$ when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $q(x)=2\cos (2x)$ and $p(x)=1+ex$ , tested on No.1–No.4 meshes, x-axis denote the eigenvalues indices.

Figure 9.

The time cost when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $q(x)=2\cos (2x)$ and $p(x)=1+ex$ , x-axis denote the algorithms: 1 denotes $\mathrm{QOCM}$ , 2 denotes $\mathrm{MCM}$ , 3 denotes $\mathrm{FEM}$ , and 4 denotes $\mathrm{FDM}$ .

Table 6.

The final refined grid sizes when the four algorithms converge within tolerance $\mathrm{eps}=1.0e-4$ for $q(x)=2\cos (2x)$ and $p(x)=1+ex$ .

Mesh no.	QOCM	MCM	FEM	FDM
1	320	5120	20,480	20,480
2	320	2560	20,480	20,480
3	320	2560	20,480	20,480
4	320	320	20,480	20,480

From Figures 1, 4, and 7 and Table 1, 3, and 5, we can know that the iterative errors of the $\mathrm{MCM}$ algorithms have similar variation trend for different variable coefficients $q(x),p(x)$ on all the four kinds of meshes. So, it is with $\mathrm{FDM}$ and $\mathrm{FEM}$ . As $\mathrm{FEM}$ and $\mathrm{FDM}$ still change in opposite trends (just as the case in $q(x),p(x)$ are constant), the combinational scheme of (19) can obtain a better approximate solution for the variable coefficients eigenproblems of cases (10) to (12) than individual linear finite element scheme $\lambda F$ or the corresponding second finite difference scheme $\lambda D$ . As the iterative error of the $\mathrm{MCM}$ varies almost as a horizontal line comparing to $\mathrm{FEM},\mathrm{FDM}$ , we know $\alpha =0.5$ is a good combined parameter for (19), particular for case I. As Figures 2, 3, 5, 6, 9, and 8 and Tables 2, 4, and 6 show, with the same accuracy requirement, for different variable coefficient $q(x),p(x)$ and different kinds of meshes, the efficient algorithm $\mathrm{QOCM}$ is always much better than the algorithms $\mathrm{FEM}$ and $\mathrm{FDM}$ , especially for case I. Algorithm $\mathrm{QOCM}$ costs much less (CPU) time and much less mesh refinement than $\mathrm{FEM}$ and $\mathrm{FDM}$ algorithms under the same accuracy requirement. Moreover, the grid size that algorithms $\mathrm{FEM}$ and $\mathrm{FDM}$ need are nearly 30 times of $\mathrm{QOCM}$ for case I and even 60 times for case III. Although the algorithm $\mathrm{MCM}$ is no better than the $\mathrm{QOCM}$ , it is also more efficient than the algorithms $\mathrm{FEM}$ and $\mathrm{FDM}$ . Particularly, on the uniform mesh, $\mathrm{MCM}$ has the same efficiency with $\mathrm{QOCM}$ . While for case I and case III the eigenvalues obtained with uniform mesh have better accuracy, the uniform mesh does not always have advantage over nonuniform mesh. Case II with variable coefficient of $p(x)$ in the eigen-equation is just a counter-example.

Conclusions

By leveraging the linear finite element method and second-order finite difference method, the combinational scheme we propose for solving eigenvalue problems with variable coefficient is effective and quite efficient. Although our efficient algorithm and its improved algorithm may not guarantee high-order accuracy, they both converge much faster than the single linear finite element method or the corresponding finite difference method on the same mesh. The numerical experiments in section IV show that, with the same accuracy requirement, for different variable coefficient $q(x),p(x)$ , the algorithm $\mathrm{QOCM}$ and $\mathrm{MCM}$ take less time and need much smaller grid size than the algorithms $\mathrm{FEM}$ and $\mathrm{FDM}$ on all the four kinds of meshes. The efficient algorithm may be more applicable for eigenproblem (11), while its improved algorithm for eigenproblems (10) and (12). The grid size that Algorithms $\mathrm{FEM}$ and $\mathrm{FDM}$ need are nearly 30 times of $\mathrm{QOCM}$ for case I and even 60 times for case III. Namely, our two algorithms are really more efficient than algorithms $\mathrm{FEM}$ and $\mathrm{FDM}$ and are applicable to compute several eigenvalues at the same time. Moreover, both our efficient algorithm and its improved algorithm are easy to parallelize.