Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
Research article

Rooftop solar PV for urban residential buildings of Nigeria: A preliminary attempt towards potential estimation

  • To decarbonise the Nigerian electricity sector and ensure stable power supply, rooftop solar PV will play a major role. However, studies aimed at estimating the technical potential of rooftop solar PV in Nigeria are limited. Here, a preliminary attempt has been made using a computationally logical methodology to estimate the technical potential of rooftop solar PV in urban residential buildings of Nigeria. We use the PVSyst® software to estimate the annual energy yield of rooftop solar PV in selected cities across the country. The paper also heads on to estimate the levelized unit cost of electricity and the break-even capital cost of rooftop solar PV. The available roof area for solar PV in urban residential buildings of Nigeria is estimated at 796 km2 and the technical potential at around 124 GWp. Annual energy yield and levelized unit cost of electricity analysis shows that Kano and Port Harcourt cities have the highest and lowest potentials of rooftop solar PV respectively. Break-even capital cost analysis suggest that rooftop solar PV is not financially attractive for consumers in R1 and R2s categories in the country presently. The enabling conditions and policy implications for deployment of rooftop solar PV in the country are also highlighted.

    Citation: Michael O. Dioha, Atul Kumar. Rooftop solar PV for urban residential buildings of Nigeria: A preliminary attempt towards potential estimation[J]. AIMS Energy, 2018, 6(5): 710-734. doi: 10.3934/energy.2018.5.710

    Related Papers:

    [1] Junjie Wang, Yaping Zhang, Liangliang Zhai . Structure-preserving scheme for one dimension and two dimension fractional KGS equations. Networks and Heterogeneous Media, 2023, 18(1): 463-493. doi: 10.3934/nhm.2023019
    [2] Tingting Ma, Yayun Fu, Yuehua He, Wenjie Yang . A linearly implicit energy-preserving exponential time differencing scheme for the fractional nonlinear Schrödinger equation. Networks and Heterogeneous Media, 2023, 18(3): 1105-1117. doi: 10.3934/nhm.2023048
    [3] Fengli Yin, Dongliang Xu, Wenjie Yang . High-order schemes for the fractional coupled nonlinear Schrödinger equation. Networks and Heterogeneous Media, 2023, 18(4): 1434-1453. doi: 10.3934/nhm.2023063
    [4] Yves Achdou, Victor Perez . Iterative strategies for solving linearized discrete mean field games systems. Networks and Heterogeneous Media, 2012, 7(2): 197-217. doi: 10.3934/nhm.2012.7.197
    [5] Wen Dong, Dongling Wang . Mittag-Leffler stability of numerical solutions to linear homogeneous time fractional parabolic equations. Networks and Heterogeneous Media, 2023, 18(3): 946-956. doi: 10.3934/nhm.2023041
    [6] Farman Ali Shah, Kamran, Dania Santina, Nabil Mlaiki, Salma Aljawi . Application of a hybrid pseudospectral method to a new two-dimensional multi-term mixed sub-diffusion and wave-diffusion equation of fractional order. Networks and Heterogeneous Media, 2024, 19(1): 44-85. doi: 10.3934/nhm.2024003
    [7] Jin Cui, Yayun Fu . A high-order linearly implicit energy-preserving Partitioned Runge-Kutta scheme for a class of nonlinear dispersive equations. Networks and Heterogeneous Media, 2023, 18(1): 399-411. doi: 10.3934/nhm.2023016
    [8] Yaxin Hou, Cao Wen, Yang Liu, Hong Li . A two-grid ADI finite element approximation for a nonlinear distributed-order fractional sub-diffusion equation. Networks and Heterogeneous Media, 2023, 18(2): 855-876. doi: 10.3934/nhm.2023037
    [9] Jinhu Zhao . Natural convection flow and heat transfer of generalized Maxwell fluid with distributed order time fractional derivatives embedded in the porous medium. Networks and Heterogeneous Media, 2024, 19(2): 753-770. doi: 10.3934/nhm.2024034
    [10] Caihong Gu, Yanbin Tang . Global solution to the Cauchy problem of fractional drift diffusion system with power-law nonlinearity. Networks and Heterogeneous Media, 2023, 18(1): 109-139. doi: 10.3934/nhm.2023005
  • To decarbonise the Nigerian electricity sector and ensure stable power supply, rooftop solar PV will play a major role. However, studies aimed at estimating the technical potential of rooftop solar PV in Nigeria are limited. Here, a preliminary attempt has been made using a computationally logical methodology to estimate the technical potential of rooftop solar PV in urban residential buildings of Nigeria. We use the PVSyst® software to estimate the annual energy yield of rooftop solar PV in selected cities across the country. The paper also heads on to estimate the levelized unit cost of electricity and the break-even capital cost of rooftop solar PV. The available roof area for solar PV in urban residential buildings of Nigeria is estimated at 796 km2 and the technical potential at around 124 GWp. Annual energy yield and levelized unit cost of electricity analysis shows that Kano and Port Harcourt cities have the highest and lowest potentials of rooftop solar PV respectively. Break-even capital cost analysis suggest that rooftop solar PV is not financially attractive for consumers in R1 and R2s categories in the country presently. The enabling conditions and policy implications for deployment of rooftop solar PV in the country are also highlighted.


    This paper mainly focuses on constructing and analyzing an efficient energy-preserving finite difference method (EP-FDM) for solving the nonlinear coupled space-fractional Klein-Gordon (KG) equations:

    uttκ2dk=1αkxku+a1u+b1u3+c1uv2=0, (1.1)
    vttκ2dk=1αkxkv+a2v+b2v3+c2u2v=0, (1.2)

    with (x,t)Ω×[0,T] and the following widely used boundary and initial conditions

    (u(x,t),v(x,t))=(0,0),(x,t)Ω×[0,T], (1.3)
    (u(x,0),v(x,0))=(ϕ1(x),ϕ2(x)),xˉΩ, (1.4)
    (ut(x,0),vt(x,0))=(φ1(x),φ2(x)),xˉΩ. (1.5)

    Here, x=(x1,...,xd)T(d=1,2,3)ΩRd, Ω is the boundary of Ω, ˉΩ=ΩΩ, κ is a constant and ai,bi,ci are all positive constants. ϕ1,ϕ2,φ1,φ2 are all known sufficiently smooth functions. u(x,t), v(x,t) are interacting relativistic fields of masses, αkxku and αkxkv stand for the Riesz fractional operator with 1<αk2,(k=1,...,d) in xk directions, which are well defined as follows

    αkxku(x,t)=12cos(αkπ/2)[Dαkxku(x,t)+xkDαk+u(x,t)], (1.6)

    where Dαkxku(x,t) and xkDαk+u(x,t) are the left and right Riemann-Liouville fractional derivative.

    Plenty of physical phenomena, such as the long-wave dynamics of two waves, are represented by the system (1.2). For example, these equations are used to study a number of issues in solid state physics, relativistic mechanics, quantum mechanics, and classical mechanics [1,2,3,4].

    Especially, when αk tends to 2, the fractional derivative αkxk would converge to the second-order Laplace operator, and thus Eqs (1.1) and (1.2) reduce to the classical system of multi-dimensional coupled KG equations [5,6,7]. The system has the following conserved energy, which is mentioned in detail in [11],

    E(t)=12Ω[1c1(ut)2+κ2c1|u|2+1c2(vt)2+κ2c2|v|2+2G(u,v)]dΩ=E(0),

    where

    G(u,v)=b14c1u4+b24c2v4+a12c1u2+a22c2v2+12u2v2.

    The coupled KG equations is initially introduced in [8] and is applied to model the usual motion of charged mesons within a magnetic field. There have been many works for solving the classical KG equations. Tsutsumi [9] considered nonrelativistic approximation of nonlinear KG equations and proved the convergence of solutions rigorously. Joseph [10] obtained some exact solutions for these systems. Deng [11] developed two kinds of energy-preserving finite difference methods for the systems of coupled sine-Gordon (SG) equations or coupled KG equations in two dimensions. He [12] analyzed two kinds of energy-preserving finite element approximation schemes for a class of nonlinear wave equation. Zhu [13] developed the finite element method and the mesh-free deep neural network approach in a comparative fashion for solving two types of coupled nonlinear hyperbolic/wave partial differential equations. Deng [14] proposed a two-level linearized compact ADI method for solving the nonlinear coupled wave equations. More relevant and significant references can be found in [15,16,17].

    However, it has been found that fractional derivatives can be used to describe some physical problems with the spatial non-locality of anomalous diffusion. Therefore, more attentions have been paid to fractional KG equations. There are also some related numerical methods for the related models. These methods may be applied to solve the fractional KG systems. For example, Cheng [18] constructed a linearized compact ADI scheme for the two-dimensional Riesz space fractional nonlinear reaction-diffusion equations. Wang [19] proposed Fourier spectral method to solve space fractional KG equations with periodic boundary condition. Liu [20] presented an implicit finite difference scheme for the nonlinear time-space-fractional Schrödinger equation. Cheng [21] constructed an energy-conserving and linearly implicit scheme by combining the scalar auxiliary variable approach for the nonlinear space-fractional Schrödinger equations. Similar scalar auxiliary variable approach can also be found in [22,23]. Wang et al. [24,25] developed some energy-conserving schemes for space-fractional Schrödinger equations. Meanwhile, the equations are also investigated by some analytical techniques, such as the Fourier transform method [26], the Mellin transform method [27] and so on. Besides, the spatial disccretization of the KG equations usually gives a system of conservative ordinary differential equations. There are also some energy-conserving time discretizations, such as the implicit midpoint method [28], some Runge–Kutta methods [28,29], relaxation methods [30,31,32] and so on [33,34]. To the best of our knowledge, there exist few reports on numerical methods for coupled space-fractional KG equations. Most references focus on the KG equations rather than the coupled systems.

    The main purpose of this paper is to develop an EP-FDM for the system of nonlinear coupled space-fractional KG equations. Firstly, we transform the coupled systems of KG equations into an equivalent general form and provide energy conservation for the new system. Secondly, we propose a second-order consistent implicit three-level scheme by using the finite difference method to solve problems (1.1) and (1.2). Thirdly, we give the proof of the discrete energy conservation, boundedness of numerical solutions and convergence analysis in discrete L2 norm. More specifically, the results show that the proposed schemes are energy-conserving. And the schemes have second-order accuracy in both the temporal and spatial directions. Finally, numerical experiments are presented to show the performance of our proposed scheme in one and two dimensions. They confirm our obtained theoretical results very well.

    The rest of the paper is organized as follows. Some denotations and preliminaries are given in Section 2. An energy-preserving scheme is constructed in Section 3. The discrete conservation law and boundedness of numerical solutions are given in Section 4. The convergence results are given in Section 5. Several numerical tests are offered to validate our theoretical results in Section 6. Finally, some conclusions are given in Section 7.

    Throughout the paper, we set C as a general positive constant that is independent of mesh sizes, which may be changed under different circumstances.

    We first rewrite Eqs (1.1) and (1.2) into an equivalent form

    αuttβdk=1αkxku+Gu(u,v)=0, (2.1)
    γvttσdk=1αkxkv+Gv(u,v)=0, (2.2)

    with the widely used boundary and initial conditions

    (u(x,t),v(x,t))=(0,0),(x,t)Ω×[0,T], (2.3)
    (u(x,0),v(x,0))=(ϕ1(x),ϕ2(x)),xˉΩ, (2.4)
    (ut(x,0),vt(x,0))=(φ1(x),φ2(x)),xˉΩ, (2.5)

    where G(u,v)=b14c1u4+b24c2v4+a12c1u2+a22c2v2+12u2v2, and α=1/c1, β=κ2/c1, γ=1/c2, σ=κ2/c2. A similar treatment is mentioned in [11]. The definition of operator αkxk is already presented in Eq (1.6), where the left and right Riemann-Liouville fractional derivatives in space of order α are defined as

    Dαxu(x,t)=1Γ(2α)2x2xu(ξ,t)(xξ)α1dξ,(x,t)Ω,xDα+u(x,t)=1Γ(2α)2x2xu(ξ,t)(ξx)α1dξ,(x,t)Ω.

    Theorem 1. Let u(x,t),v(x,t) be the solutions of this systems (2.1)–(2.5), the energy conservation law is defined by

    E(t)=12[αut2L2+βdk=1αk/2xku2L2+γvt2+σdk=1αk/2xkv2L2+2G(u,v),1]. (2.6)

    Namely, E(t)=E(0), where u(,t)2L2=Ω|u(x,t)|2dx and G(u,v),1=ΩG(u,v)dx.

    Proof. Taking inner product of Eqs (2.1) and (2.2) with ut and vt, then summing the obtained equations, and finally applying a integration over the time interval [0,t], it yields the required result.

    The finite difference method is used to achieve spatial and temporal discretization in this paper. We now denote temporal step size by τ, let τ=T/N, tn=nτ. For a list of functions {wn}, we define

    wˉn=wn+1+wn12,δtwn=wn+1wnτ,μtwn=wn+1+wn2,Dtwn=wn+1wn12τ=δtwn+δtwn12,δ2twn=wn+12wn+wn1τ2=δtwnδtwn1τ.

    Let Ω=(a1,b1)×(ad,bd), with the given positive integers M1,,Md, for the convenience of subsequent proofs, we have set it uniformly to M, so we get hk=(bkak)/Mk=h (k=1,,d) be the spatial stepsizes in xk-direction, then the spatial mesh is defined as ˉΩh={(xk1,xk2,,xkd)0ksMs,s=1,,d}, where xks=as+kshs.

    Moreover, we define the space V0h as follows by using the grid function on Ωh,

    V0h:={v=vnk1kdvnk1kd=0for(k1,,kd)Ωh},

    where 1ksMs1, s=1,,d, 0nN. Then we write δx1uk1kd=uk1+1kduk1kdh. Notations δxsuk1kd (s=2,,d) are defined similarly.

    We then introduce the discrete norm, respectively. For u,vV0h, denote

    (u,v)=hdM11k1=1Md1kd=1uk1kdvk1kd,u=(u,u),|U|2H1=ds=1δxsU2,us=[hdM11k1=1Md1kd=1(uk1kd)s]1s.

    Based on the definitions, we give the following lemmas which are important for this paper.

    Lemma 1. ([35]) Suppose p(x)L1(R) and

    p(x)C2+α(R):={p(x)+(1+|k|)2+α|ˆp(k)|dk<},

    where ˆp(k) is the Fourier transformation of p(x), then for a given h, it holds that

    Dαxp(x)=1hα+k=0w(α)kp(x(k1)h)+O(h2),xDα+p(x)=1hα+k=0w(α)kp(x+(k1)h)+O(h2),

    where w(α)k are defined by

    {w(α)0=λ1g(α)0,w(α)1=λ1g(α)1+λ0g(α)0,w(α)k=λ1g(α)k+λ0g(α)k1+λ1g(α)k2,k2, (2.7)

    where λ1=(α2+3α+2)/12,λ0=(4α2)/6,λ1=(α23α+2)/12 and g(α)k=(1)k(αk).

    In addition, we arrange in this section some of the lemmas that are necessary for the demonstration of later theorems in this paper.

    Lemma 2. ([36]) For any two grid functions u,vV0h, there exists a linear operator Λα such that (δαxu,v)=(Λα2u,Λα2v), where the difference operator Λα2 is defined by Λα2u=Lu, and matrix L satisfying C=LTL is the cholesky factor of matrix C=1/(2hαcos(απ/2))(P+PT) with

    P=[w(α)1w(α)0w(α)2w(α)1w(α)0w(α)2w(α)1w(α)M2w(α)0w(α)M1w(α)M2w(α)2w(α)1](M1)×(M1).

    While for multi-dimensional case, we give a further lemma.

    Lemma 3. ([18]) For any two grid functions u,vV0h, there exists a linear operator Λαk2k such that (δαkxku,v)=(Λαk2ku,Λαk2kv), k=1,,d, where Λαk2k is defined by Λαk2ku=[2cos(αkπ/2)hαk]1/2Lku, and matrix Lk is given by IDαkI= [2cos(αkπ/2)hαk]1LTkLk. I is a unit matrix and matrix Dαk is given by Dαk=1/(2cos(αkπ/2)hαk) (Pk+PTk), Pk is the matrix P in the case α=αk as defined in Lemma 2.

    Lemma 4. ([11]) Let g(x)C4(I), then x0I,x0+ΔxI, we have

    g(x0+Δx)2g(x0)+g(x0Δx)Δx2=g(x0)+Δx2610[g(4)(x0+λΔx)+g(4)(x0λΔx)](1λ)3dλ,g(x0+Δx)+g(x0Δx)2=g(x0)+Δx210[g(x0+λΔx)+g(x0λΔx)](1λ)dλ.

    Lemma 5. ([11]) Let u(x,t),v(x,t)C4,4(Ω×[0,T]), and G(u,v)C4,4(R1×R1). Then we have

    G(un+1,vn)G(un1,vn)un+1un1=Gu(u(x,tn),v(x,tn))+O(τ2),G(un,vn+1)G(un,vn1)vn+1vn1=Gv(u(x,tn),v(x,tn))+O(τ2).

    Lemma 6. For any grid function uV0h, it holds that

    upCuCp1(Cp2|u|H1+1lu)Cp3,2p<,

    where Cp1,Cp2,Cp3 are constants related to p, l=min, and d is the dimension of space \mathcal{V}_h^0 .

    Specially, for two-dimensional case, the parameters C_{p_1} = \frac{2}{p} , C_{p_2} = \max \left\{2 \sqrt{2}, \frac{p}{\sqrt{2}}\right\} and C_{p_3} = 1-\frac{2}{p} are shown in [37,38].

    While in the case of three dimensions, C_{p_1} = \frac{p+6}{4p} , C_{p_2} = \max \left\{2 \sqrt{3}, \frac{p}{\sqrt{3}}\right\} and C_{p_3} = \frac{3p-6}{4p} , the proof is given in Appendix.

    Lemma 7. ([39]) For M \geq 5, 1 \leq \alpha \leq 2 and any v \in \mathcal{V}_h^0 , there exists a positive constant C_1 , such that

    \|v\|^2 < \frac{\cos (\alpha \pi / 2)}{C_1 \ln 2}\left(\delta_x^\alpha v, v\right) = -\frac{\cos (\alpha \pi / 2)}{C_1 \ln 2}\left\|\Lambda^{\frac{\alpha}{2}} v\right\|^2 .

    Specially, for multi-dimensional case, it can be written as \|v\|^2 < C \sum_{k = 1}^d\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}} v\right\|^2 , where C is a positive constant.

    Lemma 8. ([40]) Assume that \left\{g^n \mid n \geq 0\right\} is a nonnegative sequence, \psi^0 \geq 0 , and the nonnegative sequence \left\{G^n \mid n \geq 0\right\} satisfies

    G^n \leq \psi^0+\tau \sum\limits_{l = 0}^{n-1} G^l+\tau \sum\limits_{l = 0}^n g^l, \quad n \geq 0 .

    Then it holds that

    G^n \leq e^{n \tau}\left(\psi^0+\tau \sum\limits_{l = 0}^n g^l\right), \quad n \geq 0.

    Lemma 9. For any grid function u \in V_h^0 , V_h^0 is defined in Section 2 for the the three-dimensional case, let p \leq r \leq q, \alpha \in[0, 1] satisfying \frac{1}{r} = \frac{\alpha}{p}+\frac{1-\alpha}{q} , then

    \|u\|_r \leq\|u\|_p^\alpha \cdot\|u\|_q^{1-\alpha} .

    Proof. By using Hölder inequality, we have

    \begin{aligned} h_1 h_2 h_3 \sum\limits_{i = 1}^{M_1-1} \sum\limits_{j = 1}^{M_2-1} \sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^r & = h_1 h_2 h_3 \sum\limits_{i = 1}^{M_1-1} \sum\limits_{j = 1}^{M_2-1} \sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^{\alpha r+(1-\alpha) r} \\ & \leq\left(h_1 h_2 h_3 \sum\limits_{i = 1}^{M_1-1} \sum\limits_{j = 1}^{M_2-1} \sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^{\alpha r \frac{p}{\alpha r}}\right)^{\frac{\alpha r}{p}}\left(h_1 h_2 h_3 \sum\limits_{i = 1}^{M_1-1} \sum\limits_{j = 1}^{M_2-1} \sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^{(1-\alpha) r \frac{q}{(1-\alpha) r}}\right)^{\frac{(1-\alpha) r}{q}} \\ & = \|u\|_p^{r \alpha} \cdot\|u\|_q^{r(1-\alpha)} . \end{aligned}

    This completes the proof.

    Now we are ready to construct the fully-discrete numerical scheme for systems (2.1) and (2.2).

    With the help of Lemma 1 and for clarity of description, we will denote the space fractional operator under one-dimensional case firstly.

    \begin{aligned} & \delta_{x, +}^\alpha v_j^n = \frac{1}{h^\alpha} \sum\limits_{k = 0}^j w_k^{(\alpha)} v_{j-k+1}^n, \quad \delta_{x, -}^\alpha v_j^n = \frac{1}{h^\alpha} \sum\limits_{k = 0}^{M-j} w_k^{(\alpha)} v_{j+k-1}^n, \\ & \delta_x^\alpha v_j^n = -1 /(2 \cos (\alpha \pi / 2))\left(\delta_{x, +}^\alpha v_j^n+\delta_{x, -}^\alpha v_j^n\right), \end{aligned}

    where w_{k}^{(\alpha)} is given in Eq (2.7). In the multi-dimensional case, the definitions of \delta_{x_{k}}^{\alpha_{k}} are similar to it.

    For numerically solving systems (2.1)–(2.5), we propose a three-level scheme. We firstly define the following approximations.

    Let u^{n}_{k_{1}\cdots k_{d}} = u(\boldsymbol{x}, t_{n}) and v^{n}_{k_{1}\cdots k_{d}} = v(\boldsymbol{x}, t_{n}) , for ease of presentation, we shall henceforth write u^{n}_{k_{1}\cdots k_{d}} for u^{n} . Denote numerical solutions of u^{n} and v^{n} by U^{n} and V^{n} , respectively.

    With the definition of G(u, v) in systems (2.1) and (2.2) and by using Lemma 5, then we have

    \begin{align} & \frac{G\left(u^{n+1}, v^n\right)-G\left(u^{n-1}, v^n\right)}{u^{n+1}-u^{n-1}} = \frac{\partial G}{\partial u}\left(u\left(\boldsymbol{x}, t_n\right), v\left(\boldsymbol{x}, t_n\right)\right)+\mathcal{O}\left(\tau^2\right), \end{align} (3.1)
    \begin{align} & \frac{G\left(u^n, v^{n+1}\right)-G\left(u^n, v^{n-1}\right)}{v^{n+1}-v^{n-1}} = \frac{\partial G}{\partial v}\left(u\left(\boldsymbol{x}, t_n\right), v\left(\boldsymbol{x}, t_n\right)\right)+\mathcal{O}\left(\tau^2\right), \end{align} (3.2)

    which is given in [11]. Further, using the space fractional operator which is already introduced above and second-order centered finite difference operator to approximate at node (\boldsymbol{x}, t_{n}) , it holds that

    \begin{align} & \alpha \delta_t^2 u^n-\beta\sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} u^{\bar{n}}+\frac{G\left(u^{n+1}, v^n\right)-G\left(u^{n-1}, v^n\right)}{u^{n+1}-u^{n-1}} = R_{1}^n , \quad 2\leq n\leq N-1 , \end{align} (3.3)
    \begin{align} & \gamma \delta_t^2 v^n-\sigma\sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} v^{\bar{n}}+\frac{G\left(u^n, v^{n+1}\right)-G\left(u^n, v^{n-1}\right)}{v^{n+1}-v^{n-1}} = R_{2}^n , \quad 2\leq n\leq N-1, \end{align} (3.4)

    and

    \begin{align} & u^1 = \phi_1\left(\boldsymbol{x}\right)+\tau \varphi_1\left(\boldsymbol{x}\right)+\frac{\tau^2}{2 \alpha}\left[\beta \sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} \phi_1\left(\boldsymbol{x}\right)-\frac{\partial G}{\partial u}\left(\phi_1\left(\boldsymbol{x}\right), \phi_2\left(\boldsymbol{x}\right)\right)\right]+R_{1}^1, \end{align} (3.5)
    \begin{align} & v^1 = \phi_2\left(\boldsymbol{x}\right)+\tau \varphi_2\left(\boldsymbol{x}\right)+\frac{\tau^2}{2 \gamma}\left[\sigma \sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} \phi_2\left(\boldsymbol{x}\right)-\frac{\partial G}{\partial v}\left(\phi_1\left(\boldsymbol{x}\right), \phi_2\left(\boldsymbol{x}\right)\right)\right]+R_2^1, \end{align} (3.6)

    where R_1^n and R_2^n are the truncation errors.

    Let u(\boldsymbol{x}, t), v(\boldsymbol{x}, t) \in \mathcal{C}^{4, 4}(\Omega \times[0, T]) . Combining Lemma 4 with Eqs (3.1) and (3.2), the truncation errors can be estimated as follows.

    \begin{align} \max _{1 \leq n \leq N-1}\left\{\left\|R_1^n\right\|^2, \left\|R_2^n\right\|^2\right\} \leq C\left(\tau^2+h_1^2+\cdots +h_d^2\right)^2 , \end{align} (3.7)

    where C is a positive constant and d means the dimension of space.

    Omitting the truncation errors in Eqs (3.3)–(3.6), we can get the three-level EP-FDM:

    \begin{align} & \alpha \delta_t^2 U^n-\beta\sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} U^{\bar{n}}+\frac{G\left(U^{n+1}, V^n\right)-G\left(U^{n-1}, V^n\right)}{U^{n+1}-U^{n-1}} = 0, \end{align} (3.8)
    \begin{align} & \gamma \delta_t^2 V^n-\sigma\sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} V^{\bar{n}}+\frac{G\left(U^n, V^{n+1}\right)-G\left(U^n, V^{n-1}\right)}{V^{n+1}-V^{n-1}} = 0 , \end{align} (3.9)

    and

    \begin{align} & U^n = V^n = 0, \ \ \boldsymbol{x}\in \partial\Omega_{h}, \ \ 0\leq n\leq N, \end{align} (3.10)
    \begin{align} & U^1 = \phi_1\left(\boldsymbol{x}\right)+\tau \varphi_1\left(\boldsymbol{x}\right)+\frac{\tau^2}{2 \alpha}\left[\beta \sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} \phi_1\left(\boldsymbol{x}\right)-\frac{\partial G}{\partial u}\left(\phi_1\left(\boldsymbol{x}\right), \phi_2\left(\boldsymbol{x}\right)\right)\right], \end{align} (3.11)
    \begin{align} & V^1 = \phi_2\left(\boldsymbol{x}\right)+\tau \varphi_2\left(\boldsymbol{x}\right)+\frac{\tau^2}{2 \gamma}\left[\sigma \sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} \phi_2\left(\boldsymbol{x}\right)-\frac{\partial G}{\partial v}\left(\phi_1\left(\boldsymbol{x}\right), \phi_2\left(\boldsymbol{x}\right)\right)\right], \end{align} (3.12)

    where U^{1} and V^{1} are obtained by applying Taylor expansion to expand u(\boldsymbol{x}, \tau) and v(\boldsymbol{x}, \tau) at (\boldsymbol{x}, 0) , and by Eq (2.4) we know that U^0 = \phi_1\left(\boldsymbol{x}\right) , V^0 = \phi_2\left(\boldsymbol{x}\right) .

    For contrast, by doing explicit treatment of nonlinear terms \frac{\partial G}{\partial u} and \frac{\partial G}{\partial v} , we introduce an explicit scheme as follows

    \begin{align} & \alpha \delta_t^2 U^n-\beta\sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} U^{\bar{n}}+\frac{\partial G}{\partial u}\left(U^n, V^n\right) = 0, \end{align} (3.13)
    \begin{align} & \gamma \delta_t^2 V^n-\sigma\sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} V^{\bar{n}}+\frac{\partial G}{\partial v}\left(U^n, V^n\right) = 0 , \end{align} (3.14)
    \begin{align} & U^n = V^n = 0, \ \ \boldsymbol{x}\in \partial\Omega_{h}, \ \ 0\leq n\leq N, \end{align} (3.15)
    \begin{align} & U^1 = \phi_1\left(\boldsymbol{x}\right)+\tau \varphi_1\left(\boldsymbol{x}\right)+\frac{\tau^2}{2 \alpha}\left[\beta \sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} \phi_1\left(\boldsymbol{x}\right)-\frac{\partial G}{\partial u}\left(\phi_1\left(\boldsymbol{x}\right), \phi_2\left(\boldsymbol{x}\right)\right)\right], \end{align} (3.16)
    \begin{align} & V^1 = \phi_2\left(\boldsymbol{x}\right)+\tau \varphi_2\left(\boldsymbol{x}\right)+\frac{\tau^2}{2 \gamma}\left[\sigma \sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} \phi_2\left(\boldsymbol{x}\right)-\frac{\partial G}{\partial v}\left(\phi_1\left(\boldsymbol{x}\right), \phi_2\left(\boldsymbol{x}\right)\right)\right], \end{align} (3.17)

    which will be used in Section 6 later.

    In this section, we give the energy conservation of the fully-discrete schemes (3.8)–(3.12) and boundedness of numerical solutions. Here, the lemmas given in Section 2 are applied.

    Now, we present the energy conservation of the EP-FDMs (3.8)–(3.12).

    Theorem 2. Let U^n, V^n\in \mathcal{V}_{h}^{0} be numerical solutions of the three-level FDMs (3.8)–(3.12). Then, the energy, which is defined by

    \begin{align} E^n = &\frac{\alpha}{2}\|\delta_{t}U^n\|^{2}+\frac{\beta}{2}\sum\limits_{k = 1}^{d}\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^n\|^2+ \frac{\gamma}{2}\|\delta_{t}V^n\|^{2}+\frac{\sigma}{2}\sum\limits_{k = 1}^{d}\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}V^n\|^2\\ &+ \frac{1}{2}h^d \sum\limits_{k_1 = 1}^{M_1-1} \cdots \sum\limits_{k_d = 1}^{M_d-1}\left[G(U_{k_{1}\cdots k_{d}}^{n+1}, V_{k_{1}\cdots k_{d}}^{n})+G(U^{n}_{k_{1}\cdots k_{d}}, V_{k_{1}\cdots k_{d}}^{n+1})\right] \end{align} (4.1)

    is conservative. Namely, E^n = E^0 , for n = 1, \cdots, N-1, where \Lambda_{k}^{\frac{\alpha_{k}}{2}} is already introduced by Lemma 3.

    Proof. Multiplying h^dD_{t}U_{k_{1}\cdots k_{d}}^n to both sides of Eq (3.8), summing them over \Omega_{h} , by using Lemma 3, we obtain

    \begin{align} &\frac{\alpha}{2\tau}\left( \|\delta_{t}U^n\|^2-\|\delta_{t}U^{n-1}\|^2\right)+\frac{\beta}{4\tau}\sum\limits_{k = 1}^d\left( \|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^{n+1}\|^2-\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^{n-1}\|^2\right)\\ &+\frac{1}{2\tau}h^d \sum\limits_{k_1 = 1}^{M_1-1} \cdots \sum\limits_{k_d = 1}^{M_d-1}\left[G(U_{k_{1}\cdots k_{d}}^{n+1}, V_{k_{1}\cdots k_{d}}^{n})-G(U^{n-1}_{k_{1}\cdots k_{d}}, V_{k_{1}\cdots k_{d}}^{n})\right] = 0, \end{align} (4.2)

    where the second term can be reduced to

    \|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^{n+1}\|^2-\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^{n-1}\|^2 = 2\left(\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^{n}\|^2-\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^{n-1}\|^2\right),

    then Eq (4.2) turned into

    \begin{align} &\frac{\alpha}{2\tau}\left( \|\delta_{t}U^n\|^2-\|\delta_{t}U^{n-1}\|^2\right)+\frac{\beta}{2\tau}\sum\limits_{k = 1}^d\left( \mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^{n}\|^2-\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^{n-1}\|^2\right)\\ &+\frac{1}{2\tau}h^d \sum\limits_{k_1 = 1}^{M_1-1} \cdots \sum\limits_{k_d = 1}^{M_d-1}\left[G(U_{k_{1}\cdots k_{d}}^{n+1}, V_{k_{1}\cdots k_{d}}^{n})-G(U^{n-1}_{k_{1}\cdots k_{d}}, V_{k_{1}\cdots k_{d}}^{n})\right] = 0. \end{align} (4.3)

    Similarly, multiplying h^dD_{t}V_{k_{1}\cdots k_{d}}^n to both sides of Eq (3.9), summing them over \Omega_{h} , by using Lemma 3, we obtain

    \begin{align} &\frac{\gamma}{2\tau}\left( \|\delta_{t}V^n\|^2-\|\delta_{t}V^{n-1}\|^2\right)+\frac{\sigma}{2\tau}\sum\limits_{k = 1}^d\left( \mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}V^{n}\|^2-\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}V^{n-1}\|^2\right)\\ &+\frac{1}{2\tau}h^d \sum\limits_{k_1 = 1}^{M_1-1} \cdots \sum\limits_{k_d = 1}^{M_d-1}\left[G(U_{k_{1}\cdots k_{d}}^{n}, V_{k_{1}\cdots k_{d}}^{n+1})-G(U^{n}_{k_{1}\cdots k_{d}}, V_{k_{1}\cdots k_{d}}^{n-1})\right] = 0. \end{align} (4.4)

    Adding up Eqs (4.3) and (4.4) yields that (E^n-E^{n-1})/\tau = 0 , which infers that E^n = E^{n-1} .

    By Theorem 2, we present the following estimation.

    Theorem 3. Let U^n, V^n\in \mathcal{V}_{h}^{0} be numerical solutions of the EP-FDMs (3.8)–(3.12). Then, the following estimates hold:

    \begin{align} \max _{1 \leq n \leq N}\left\{\left\|\delta_t U^{n}\right\|, \left\|\delta_t V^{n}\right\|, \left\|U^n\right\|, \left\|V^n\right\|, \left\| \Lambda_{k}^{\frac{\alpha_{k}}{2}}U^n\right\|, \left\| \Lambda_{k}^{\frac{\alpha_{k}}{2}}V^n\right\|\right\} \leq C, \end{align} (4.5)

    where C is a positive constant independent of \tau and h and 1\leq \alpha_{k}\leq 2. Specially, when \alpha_{k} = 2, it holds that |U^n|_{H_{1}}\leq C , |V^n|_{H_{1}}\leq C .

    Proof. It follows from Theorem 2, there exists a constant C such that

    \begin{align} E^n = &\frac{\alpha}{2}\|\delta_{t}U^n\|^{2}+\frac{\beta}{2}\sum\limits_{k = 1}^{d}\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^n\|^2+ \frac{\gamma}{2}\|\delta_{t}V^n\|^{2}+\frac{\sigma}{2}\sum\limits_{k = 1}^{d}\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}V^n\|^2\\ &+ \frac{1}{2}h^d \sum\limits_{k_1 = 1}^{M_1-1} \cdots \sum\limits_{k_d = 1}^{M_d-1}\left[G(U_{k_{1}\cdots k_{d}}^{n+1}, V_{k_{1}\cdots k_{d}}^{n})+G(U^{n}_{k_{1}\cdots k_{d}}, V_{k_{1}\cdots k_{d}}^{n+1})\right] = E^0 = C, \end{align}

    then, we obtain

    \|\delta_{t}U^n\|\leq C, \quad \|\delta_{t}V^n\|\leq C, \quad \|\Lambda_{k}^{\frac{\alpha_{k}}{2}}U^n\|\leq C, \quad\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}V^n\|\leq C.

    By \|\delta_{t}U^n\|\leq C , we have \|U^{n+1}-U^n\|\leq C\tau , then it is easy to check that

    \|U^n\| = \|U^0+\tau\sum\limits_{i = 0}^{n-1}\delta_{t}U^i\|\leq\|U^0\|+\tau\sum\limits_{i = 0}^{n-1}\|\delta_{t}U^i\|\leq C.

    This completes the proof.

    In this section, the convergence analysis of the proposed scheme is given, which is based on some important lemmas presented in Section 2.

    We first give the error equations of the EP-FDMs (3.8) and (3.9). Let e^n = u^n-U^n , \theta^n = v^n-V^n and for more readability we denote

    \begin{align} & \varepsilon_1\left(u^{n+1}, U^{n+1}\right) = \frac{G\left(u^{n+1}, v^n\right)-G\left(u^{n-1}, v^n\right)}{u^{n+1}-u^{n-1}}-\frac{G\left(U^{n+1}, V^n\right)-G\left(U^{n-1}, V^n\right)}{U^{n+1}-U^{n-1}}, \end{align} (5.1)
    \begin{align} & \varepsilon_2\left(v^{n+1}, V^{n+1}\right) = \frac{G\left(u^n, v^{n+1}\right)-G\left(u^n, v^{n-1}\right)}{v^{n+1}-v^{n-1}}-\frac{G\left(U^n, V^{n+1}\right)-G\left(U^n, V^{n-1}\right)}{V^{n+1}-V^{n-1}}. \end{align} (5.2)

    By deducting Eqs (3.8) and (3.9) from Eqs (3.3) and (3.4), we have

    \begin{align} & \alpha \delta_t^2 e^n-\beta\sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} e^{\bar{n}}+\varepsilon_1\left(u^{n+1}, U^{n+1}\right) = R_{1}^n , \quad 1\leq n\leq N-1 , \end{align} (5.3)
    \begin{align} & \gamma \delta_t^2 \theta^n-\sigma\sum\limits_{i = 1}^d \delta_{x_i}^{\alpha_i} \theta^{\bar{n}}+\varepsilon_2\left(v^{n+1}, V^{n+1}\right) = R_{2}^n , \quad 1\leq n\leq N-1, \end{align} (5.4)
    \begin{align} & e^n = \theta^n = 0, \boldsymbol{x}\in \partial\Omega_h, \ 1\leq n\leq N-1, \end{align} (5.5)
    \begin{align} & e^0 = \theta^0 = 0, \boldsymbol{x}\in \bar{\Omega}_h, \end{align} (5.6)
    \begin{align} &\|e^1\|\leq c_1\tau^3, \ \boldsymbol{x}\in \Omega_h, \end{align} (5.7)
    \begin{align} &\|\theta^1\|\leq c_2\tau^3, \ \boldsymbol{x}\in \Omega_h. \end{align} (5.8)

    Before giving a proof of convergence, we provide the following estimates for Eqs (5.1)–(5.2).

    Lemma 10. On \bar{\Omega}_h , we have

    \begin{align} & \left(\varepsilon_1\!\left(u^{n\!+\!1}\!, \! U^{n\!+\!1}\right)\!, \! D_t e^n\right) \!\leq \! C\!\left(\sum\limits_{k = 1}^d\!\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\!\theta^n\right\|^2\!+\!\sum\limits_{k = 1}^d\!\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\!e^{n\!+\!1}\right\|^2 \!+\!\sum\limits_{k = 1}^d\!\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\!e^{n-1}\right\|^2\!+\!\left\|\delta_t e^{n}\right\|^2\!+\!\left\|\delta_t e^{n-1}\right\|^2\right)\!, \! \end{align} (5.9)
    \begin{align} & \left(\varepsilon_2\!\left(v^{n\!+\!1}\!, \! V^{n\!+\!1}\right)\!, \! D_t \theta^n\right) \!\leq\! C\!\left(\sum\limits_{k = 1}^d\!\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\!e^n\right\|^2\!+\!\sum\limits_{k = 1}^d\!\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\!\theta^{n\!+\!1}\right\|^2 \!+\!\sum\limits_{k = 1}^d\!\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\!\theta^{n-1}\right\|^2\!+\!\left\|\delta_t \theta^{n}\right\|^2\!+\!\left\|\delta_t \theta^{n-1}\right\|^2\right)\!, \! \end{align} (5.10)

    where C > 0 is a constant, independent of grid parameters \tau, h_1, \cdots, h_d .

    Proof. Recalling the definition of G(u, v) , we can obtain

    \begin{align*} \varepsilon_1\left(u^{n+1}, U^{n+1}\right) & = \frac{b_1}{2 c_1}\left\{\left[\left(u^{n+1}\right)^2+\left(u^{n-1}\right)^2\right] u^{\bar{n}}-\left[\left(U^{n+1}\right)^2+\left(U^{n-1}\right)^2\right] U^{\bar{n}}\right\} \\ & +\left[\left(v^n\right)^2\left(u^{\bar{n}}\right)-\left(V^n\right)^2 U^{\bar{n}}\right]+\frac{a_1}{c_1} e^{\bar{n}} = \sum\limits_{k = 1}^3 Q_k . \end{align*}

    Noting that U^k = u^k-e^k and V^k = v^k-\theta^k\ (k = n-1, n, n+1) , then we get

    \begin{align} Q_1 = & \frac{b_1}{2 c_1}\left[2 u^{n+1} e^{n+1}-\left(e^{n+1}\right)^2+2 u^{n-1} e^{n-1}-\left(e^{n-1}\right)^2\right] u^{\bar{n}} \\ & +\frac{b_1}{2 c_1}\left[\left(u^{n+1}\right)^2-2 u^{n+1} e^{n+1}+\left(e^{n+1}\right)^2+\left(u^{n-1}\right)^2-2 u^{n-1} e^{n-1}+\left(e^{n-1}\right)^2\right] e^{\bar{n}}, \end{align} (5.11)
    \begin{align} Q_2 = & 2 u^{\bar{n}} v^n \theta^n-u^{\bar{n}}\left(\theta^n\right)^2+\left(V^n\right)^2 e^{\bar{n}} . \end{align} (5.12)

    When d = 2 , combining Theorem 3, Lemma 6 with Lemma 7, we can get the estimation of \|e^m\|_{4}^{4} , \|e^m\|_{6}^{6} , \|e^m\|_{8}^{8} , that is

    \begin{align} & \left\|e^m\right\|_4^4 \leq\left\|e^m\right\|^2\left(2\left|e^m\right|_{H^1}+\frac{1}{l}\left\|e^m\right\|\right)^2 \\ & \leq\left\|e^m\right\|^2\left[8\left(\left|u^m\right|_{H^1}^2+\left|U^m\right|_{H^1}^2\right)+\frac{2}{l^2}\left(\left\|u^m\right\|^2+\left\|U^m\right\|^2\right)\right] \\ & \leq C\left\|e^m\right\|^2 \leq C\sum\limits_{k = 1}^d\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}} e^m\right\|^2. \end{align} (5.13)

    The same reasoning can be used to prove that

    \begin{align} \left\|e^m\right\|_6^6 \leq C\sum\limits_{k = 1}^d\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}} e^m\right\|^2, \ \left\|e^m\right\|_8^8 \leq C\sum\limits_{k = 1}^d\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}} e^m\right\|^2, \end{align} (5.14)

    Smilarly, when d = 3 the results can be found in the same way.

    By using Cauchy-Schwarz inequality and the widely used inequality [(a+b)/2]^s\leq(a^s+b^s)/2 \ (a\geq0, b\geq0, s\geq1) , multiplying both sides of Eq (5.11) by h^dD_{t}e^n , then summing it on whole \Omega_h , it follows that

    \begin{align} \left(Q_1, D_t e^n\right) \leq &\frac{b_1}{4 c_1}\left[\frac{5 M^2}{2}\left(\left\|e^{n+1}\right\|^2+\left\|e^{n-1}\right\|^2\right)+\left(3 M+\frac{1}{4}\right)\left(\left\|e^{n+1}\right\|_4^4\right.\right. \\ & \left.\left.+\left\|e^{n-1}\right\|_4^4\right)+\frac{1}{2}\left(\left\|e^{n+1}\right\|_6^6+\left\|e^{n-1}\right\|_6^6\right)+\frac{1}{8}\left(\left\|e^{n+1}\right\|_8^8+\left\|e^{n-1}\right\|_8^8\right)\right] \\ & +\frac{b_1}{4 c_1}\left(5 M^2+6 M+1\right)\left\|D_t e^n\right\|^2 \\ \leq & C\sum\limits_{k = 1}^d\left(\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^{n+1}\right\|^2+\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^{n-1}\right\|^2\right) +\frac{b_1}{8 c_1}\left(5 M^2+6 M+1\right)\left(\left\|\delta_t e^{n}\right\|^2+\left\|\delta_t e^{n-1}\right\|^2\right), \end{align} (5.15)

    the last inequality is derived by inequalities (5.13) and (5.14), similarly, we can also obtain

    \begin{align} \left(Q_2, D_t e^n\right) & \leq M^2\left\|\theta^n\right\|^2+\frac{M}{2}\left\|\theta^n\right\|_4^4+\frac{M^2}{4}\left(\left\|e^{n+1}\right\|^2+\left\|e^{n-1}\right\|^2\right)+\left(\frac{3 M^2}{4}+\frac{M}{4}\right)\left\|D_t e^n\right\|^2 \\ & \leq C\sum\limits_{k = 1}^d\left(\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\theta^{n}\right\|^2+\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^{n+1}\right\|^2 +\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^{n-1}\right\|^2\right)\\ &\quad +\left(\frac{3 M^2}{4}+\frac{M}{4}\right)\left(\left\|\delta_t e^{n}\right\|^2+\left\|\delta_t e^{n-1}\right\|^2\right), \end{align} (5.16)
    \begin{align} \left(Q_3, D_t e^n\right) & \leq \frac{a_1^2 C}{4 c_1^2}\sum\limits_{k = 1}^d\left(\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^{n+1}\right\|^2+\left\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^{n-1}\right\|^2\right)+\frac{1}{4}\left(\left\|\delta_t e^{n}\right\|^2+\left\|\delta_t e^{n-1}\right\|^2\right), \end{align} (5.17)

    combine inequalities (5.15)–(5.17), then we get inequality (5.9) is proved. We can demonstrate that inequality (5.10) is likewise true using techniques similar to inequality (5.9). This completes the proof.

    Now we further investigate the accuracy of the proposed scheme with the help of the above lemmas, see Theorem 4.

    Theorem 4. Assume that u(\boldsymbol{x}, t), v(\boldsymbol{x}, t) \in \mathcal{C}^{4, 4}(\Omega \times[0, T]) are exact solutions of systems (2.1)–(2.5), let u_{k_{1}\cdots k_{d}}^n = u(\boldsymbol{x}, t) and v_{k_{1}\cdots k_{d}}^n = v(\boldsymbol{x}, t) , denote numerical solutions by U_{k_{1}\cdots k_{d}}^n and V_{k_{1}\cdots k_{d}}^n , define e^n = u^n-U^n , \theta^n = v^n-V^n (1\leq n\leq N) . Then suppose that \tau is sufficiently small. The error estimates of the EP-FDM are

    \begin{align} &\sum\limits_{k = 1}^d\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^n\|^2\leq C(\tau^2+h_1^2+\cdots+h_d^2)^2, \quad \|e^n\|\leq C(\tau^2+h_1^2+\cdots+h_d^2), \\ &\sum\limits_{k = 1}^d\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\theta^n\|^2\leq C(\tau^2+h_1^2+\cdots+h_d^2)^2, \quad \|\theta^n\|\leq C(\tau^2+h_1^2+\cdots+h_d^2), \end{align}

    where C is a positive constant, independent of grid parameters \tau, h_1, \cdots, h_d .

    Proof. Noting that at every time level, the systems defined in Eqs (3.8) and (3.9) is a linear PDE. Obviously, the existence and uniqueness of the solution can be obtained.

    For ease of expression, we write

    I^n = \alpha\|\delta_{t}e^n\|^2+\beta\sum\limits_{k = 1}^d\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^n\|^2 +\gamma\|\delta_{t}\theta^n\|^2+\sigma\sum\limits_{k = 1}^d\mu_{t}\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\theta^n\|^2.

    Apparently, we have that I^1\leq C(\tau^2+h_1^2+\cdots+h_d^2)^2 .

    Multiplying h^dD_{t}e^n and h^dD_{t}\theta^n to both sides of Eqs (5.3) and (5.4), then summing it over the whole \Omega_h respectively. Then adding up the obtained results, it follows that

    \begin{align} \frac{I^n-I^{n-1}}{2\tau}+(\varepsilon_1\left(u^{n+1}, U^{n+1}\right), D_{t}e^n)+(\varepsilon_2\left(v^{n+1}, V^{n+1}\right), D_{t}\theta^n) = (R_1^n, D_te^n)+(R_2^n, D_t\theta^n), \end{align} (5.18)

    by using Cauchy-Schwarz inequality, we have

    \begin{align} \frac{I^n-I^{n-1}}{2\tau}\leq&|(\varepsilon_1\left(u^{n+1}, U^{n+1}\right), D_{t}e^n)|+|(\varepsilon_2\left(v^{n+1}, V^{n+1}\right), D_{t}\theta^n)|\\ &+\frac{1}{2}\|R_1^n\|^2+\frac{1}{4}(\|\delta_te^n\|^2+\|\delta_te^{n-1}\|^2)\\ &+\frac{1}{2}\|R_2^n\|^2+\frac{1}{4}(\|\delta_t\theta^n\|^2+\|\delta_t\theta^{n-1}\|^2), \end{align} (5.19)

    multiplying 2\tau to both sides of inequality (5.19), and using Lemma 10, then we get

    \begin{align} I^{n}-I^{n-1}\leq2C\tau(I^n+I^{n-1})+\tau\|R_1^n\|^2+\tau\|R_2^n\|^2. \end{align} (5.20)

    Thus, \forall K (2\leq n\leq K\leq N-1) , summing n from 2 to K , we get

    \begin{align} (1-2C\tau)I^K\leq I^1+4C\tau\sum\limits_{n = 1}^{K-1}I^n+\sum\limits_{n = 2}^{K}\tau(\|R_1^n\|^2+\|R_2^n\|^2), \end{align} (5.21)

    when C\tau\leq\frac{1}{3} , inequality (5.21) is turned into

    \begin{align} I^K\leq 3I^1+12C\tau\sum\limits_{n = 1}^{K-1}I^n+3\tau\sum\limits_{n = 2}^{K}(\|R_1^n\|^2+\|R_2^n\|^2), \end{align} (5.22)

    then by using Lemma 8 and inequality (3.7), we obtain

    \begin{align} I^K&\leq e^{n\tau}(3I^1+3\tau\sum\limits_{n = 2}^{K}(\|R_1^n\|^2+\|R_2^n\|^2))\\ &\leq C(\tau^2+h_1^2+\cdots+h_d^2)^2. \end{align} (5.23)

    By the definition of I , it is easy to conclude that

    \begin{align} &\sum\limits_{k = 1}^d\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}e^n\|^2\leq C(\tau^2+h_1^2+\cdots+h_d^2)^2, \quad \|\delta_{t}e^n\|\leq C(\tau^2+h_1^2+\cdots+h_d^2), \\ &\sum\limits_{k = 1}^d\|\Lambda_{k}^{\frac{\alpha_{k}}{2}}\theta^n\|^2\leq C(\tau^2+h_1^2+\cdots+h_d^2)^2, \quad \|\delta_{t}\theta^n\|\leq C(\tau^2+h_1^2+\cdots+h_d^2), \end{align}

    furthermore, we have

    \|e^n\| = \|e^0+\tau\sum\limits_{i = 0}^{n-1}\delta_{t}e^i\|\leq\tau\sum\limits_{i = 0}^{n-1}\|\delta_{t}e^i\|\leq C(\tau^2+h_1^2+\cdots+h_d^2).

    Similarly, \|\theta^n\|\leq C(\tau^2+h_1^2+\cdots+h_d^2). This completes the proof.

    We carry out several numerical examples to support the theoretical results in this section. All computations are performed with Matlab. Throughout the experiments, the spatial domain is divided into M parts in every direction uniformly, that is, in the 1D case, we set M_{1} = M , while in the 2D case, we set M_{1} = M_{2} = M , and the time interval [0, T] is also divided uniformly into N parts. Then we use the discrete L^\infty -norm to measure the global error of the scheme, namely,

    \begin{align*} E_u(M, N) = \| U^{N} - u(T)\|_{\infty}, \quad E_v(M, N) = \| V^{N} - v(T)\|_{\infty}, \end{align*}

    Example 1. Consider the following one-dimensional coupled KG model

    \begin{align*} u_{t t}-\kappa^2 \partial_x^{\alpha} u+a_1 u+b_1 u^3+c_1 u v^2 = g, \quad& (x, t) \in \Omega \times[0, T], \\ v_{t t}-\kappa^2 \partial_x^{\alpha} v+a_2 v+b_2 v^3+c_2 u^2 v = g, \quad& (x, t) \in \Omega \times[0, T], \end{align*}

    with \Omega = [0, 1] . The initial and boundary conditions are determined by the exact solutions

    \begin{align*} u(x, t) = x^4(1-x)^4e^{-t}, \quad v(x, t) = x^5(1-x)^5\cos(1+t), \end{align*}

    as well as the source term g . Here, we take a_1 = a_2 = 1 , b_1 = -1 , b_2 = -2 , c_1 = 1 , c_2 = 0.5 and \kappa = 1 .

    The precision of the scheme in spatial direction is first tested by fixing N = 1000 . We compute the global errors at T = 1 with different mesh sizes, and the numerical results with \alpha = 1.2, 1.5, 1.8 are listed in Table 1 and Table 2. As can be seen in the table, the proposed scheme can have second order convergence in space, which confirms the results of theoretical analysis in Theorem 4. To track the evolution of the discrete energy, we preserve the initial value condition in this case and set the source term to g = 0 . Additionally, for the terminal time T = 50 , we fix h = 0.05 and \tau = 0.05 . The evolutionary trend image for scheme1 (3.8)–(3.12) and explicit scheme2 (3.13)–(3.17) with various \alpha are displayed in Figure 1. Then we further verify that the proposed scheme1 (3.8)–(3.12) preserves the discrete energy very well but scheme2 (3.13)–(3.17) does not.

    Table 1.  L^\infty error and spatial convergence rates of scheme1 (3.8)–(3.12) for Example 1.
    \alpha =1.2 \alpha =1.5 \alpha =1.8
    M E_u(M, N) order(u) E_u(M, N) order(u) E_u(M, N) order(u)
    32 1.51e-05 * 5.86e-06 * 4.30e-06 *
    64 3.69e-06 2.03 1.46e-06 2.01 1.09e-06 1.98
    128 9.18e-07 2.01 3.66e-07 2.00 2.74e-07 1.99
    256 2.28e-07 2.01 9.08e-08 2.01 6.75e-08 2.02

     | Show Table
    DownLoad: CSV
    Table 2.  L^\infty error and spatial convergence rates of scheme1 (3.8)–(3.12) for Example 1.
    \alpha =1.2 \alpha =1.5 \alpha =1.8
    M E_v(M, N) order(v) E_v(M, N) order(v) E_v(M, N) order(v)
    32 1.61e-06 * 3.33e-06 * 2.21e-06 *
    64 4.11e-07 1.97 8.14e-07 2.03 5.30e-07 2.06
    128 1.04e-07 1.99 2.02e-07 2.01 1.31e-07 2.01
    256 2.60e-08 2.00 5.06e-08 2.00 3.28e-08 2.00

     | Show Table
    DownLoad: CSV
    Figure 1.  The long time discrete energy of Example 1 with h = 0.05 , \tau = 0.05 for scheme1 (3.8)–(3.12) and explicit scheme2 (3.13)–(3.17).

    Example 2. Consider the following two-dimensional coupled KG model

    \begin{align*} u_{t t}-\kappa^2 \partial_{x}^{\alpha_1} u-\kappa^2 \partial_{y}^{\alpha_2} u+a_1 u+b_1 u^3+c_1 u v^2 = g, \quad& (x, y, t) \in \Omega \times[0, T], \\ v_{t t}-\kappa^2 \partial_{x}^{\alpha_1} v-\kappa^2 \partial_{y}^{\alpha_2} v+a_2 v+b_2 v^3+c_2 u^2 v = g, \quad& (x, y, t) \in \Omega \times[0, T], \end{align*}

    with \Omega = [0, 2] \times[0, 2] . The initial and boundary conditions are determined by the exact solutions

    \begin{align*} u(x, y, t) = x^2(2-x)^2y^2(2-y)^2e^{-t}, \quad v(x, y, t) = x^4(2-x)^4y^4(2-y)^4\sin(1+t), \end{align*}

    as well as the source term g . Here, we take a_1 = a_2 = 1 , b_1 = -1 , b_2 = -2 , c_1 = 1 , c_2 = 0.5 and \kappa = 1 .

    Similar to Example 1, we verify the convergence orders of the scheme in spatial direction at T = 1 . For spatial convergence order, we still set N = 1000 and thus the temporal error of the scheme can be negligible. The numerical results are presented in Table 3 and Table 4 with different values of \alpha_1 and \alpha_2 which are in the x and y directions, respectively. The second-order accuracy of the scheme is achieved. Moreover, for the terminal time T = 100 , Figure 2 shows the evolution of discrete energy for scheme1 (3.8)–(3.12) and explicit scheme2 (3.13)–(3.17) when g(x, y, t) = 0 . The figure indicate that the discrete conservation law holds very well if the proposed scheme1 (3.8)–(3.12) are used. In contrast, scheme2 (3.13)–(3.17) cannot preserve the discrete energy. Both tables and figure further confirm the theoretical results.

    Table 3.  L^\infty error and spatial convergence rates of scheme1 (3.8)–(3.12) for Example 2.
    \alpha_{1} =1.3, \alpha_{2} =1.6 \alpha_1 =1.5, \alpha_2 =1.5 \alpha_1 =1.7, \alpha_2 =1.2
    M E_u(M, N) order(u) E_u(M, N) order(u) E_u(M, N) order(u)
    8 3.76e-02 * 3.76e-02 * 3.82e-02 *
    16 9.44e-03 2.00 9.28e-03 2.02 9.63e-03 1.99
    32 2.32e-03 2.03 2.30e-03 2.01 2.37e-03 2.02
    64 5.70e-04 2.02 5.65e-04 2.03 5.85e-04 2.02

     | Show Table
    DownLoad: CSV
    Table 4.  L^\infty error and spatial convergence rates of scheme1 (3.8)–(3.12) for Example 2.
    \alpha_{1} =1.3, \alpha_{2} =1.6 \alpha_1 =1.5, \alpha_2 =1.5 \alpha_1 =1.7, \alpha_2 =1.2
    M E_v(M, N) order(v) E_v(M, N) order(v) E_v(M, N) order(v)
    8 2.98e-01 * 3.02e-01 * 2.77e-01 *
    16 6.16e-02 2.28 6.13e-02 2.30 5.72e-02 2.28
    32 1.46e-02 2.08 1.45e-02 2.08 1.36e-02 2.08
    64 3.60e-03 2.02 3.56e-03 2.02 3.34e-03 2.02

     | Show Table
    DownLoad: CSV
    Figure 2.  The long time discrete energy of Example 2 with h = 0.1 , \tau = 0.05 for scheme1 (3.8)–(3.12) and explicit scheme2 (3.13)–(3.17).

    Example 3. Consider the following two-dimensional coupled KG model

    \begin{align*} u_{t t}-\kappa^2 \partial_{x}^{\alpha_1} u-\kappa^2 \partial_{y}^{\alpha_2} u+a_1 u+b_1 u^3+c_1 u v^2 = 0, \quad& (x, y, t) \in \Omega \times[0, T], \\ v_{t t}-\kappa^2 \partial_{x}^{\alpha_1} v-\kappa^2 \partial_{y}^{\alpha_2} v+a_2 v+b_2 v^3+c_2 u^2 v = 0, \quad& (x, y, t) \in \Omega \times[0, T], \end{align*}

    and

    \begin{align*} & (u(x, y, t), v(x, y, t)) = (0, 0), \quad(x, y, t) \in \partial \Omega \times[0, T], \\ & (u(x, y, 0), v(x, y, 0)) = \left(u_0(x, y), v_0(x, y)\right), \quad (x, y) \in \bar{\Omega}, \\ & \left(u_t(x, y, 0), v_t(x, y, 0)\right) = (0, 0), \quad (x, y) \in \bar{\Omega}, \end{align*}

    with \Omega = [0, 1] \times[0, 1] .

    Here, we take

    u_0(x, y) = 2[1-\cos (2 \pi x)][1-\cos (2 \pi y)] \operatorname{sech}(x+y),
    v_0(x, y) = 4 \sin (\pi x) \sin (\pi y) \tanh (x+y)

    and

    a_1 = 10, a_2 = 4, b_1 = 6, b_2 = 5, c_1 = 2, c_2 = 3, \kappa = 1.

    The scheme1 (3.8)–(3.12) with

    \tau = h = 0.05, \alpha_1 = \alpha_2 = 1.5

    are used to Example 3. Figure 3 and Figure 4 show the surfaces of U_{ij}^n and V_{ij}^n at different times, respectively. The significant dynamical evolutionary features of the numerical solutions U_{ij}^n and V_{ij}^n , such as radiation and oscillation, can be found in Figure 3 and Figure 4.

    Figure 3.  Surfaces of U_{ij}^n at different times of Example 3 with \alpha_1 = \alpha_2 = 1.5 for scheme1 (3.8)–(3.12).
    Figure 4.  Surfaces of V_{ij}^n at different times of Example 3 with \alpha_1 = \alpha_2 = 1.5 for scheme1 (3.8)–(3.12).

    In this paper, the three-level energy-preserving scheme is proposed for the space-fractional coupled KG systems. The scheme is derived by using the finite difference method. The discrete conservation law, boundedness of numerical solutions and the global error of the scheme are further discussed. It is shown that the scheme can have second order convergence in both temporal direction and spatial direction. Several numerical examples are performed to support the theoretical results in the paper. Moreover, due to the nonlocal derivative operator and considering that the implicit methods involve Toeplitz matrices, fast methods are fairly meaningful to reduce the computational cost of the proposed scheme; refer to the recent work [41,42] for this issue.

    This work is supported by NSFC (Grant Nos. 11971010, 12001067).

    The authors declared that they have no conflicts of interest to this work.

    In the following, we present the proof of Lemma 6.

    Proof of Lemma 6: Obviously, the result holds for p = 2 . We prove the conclusion for p > 2 .

    For any m, s = 1, 2, \ldots, M_1-1 , and m > s , using mean value theorem, we have

    \left|u_{m j k}\right|^{\frac{p}{3}}-\left|u_{s j k}\right|^{\frac{p}{3}} = \sum\limits_{i = s}^{m-1}\left(\left|u_{i+1, j k}\right|^{\frac{p}{3}}-\left|u_{i j k}\right|^{\frac{p}{3}}\right) = \frac{p}{3} \sum\limits_{i = s}^{m-1}\left(\left|u_{i+1, j k}\right|-\left|u_{i j k}\right|\right) \xi_{i j k}^{\frac{p}{3}-1},

    where

    \xi_{i j k} \in\left(\min \left\{\left|u_{i j k}\right|, \left|u_{i+1, j k}\right|\right\}, \max \left\{\left|u_{i j k}\right|, \left|u_{i+1, j k}\right|\right\}\right) .

    Then,

    \begin{aligned} \left|u_{m j k}\right|^{\frac{p}{3}}-\left|u_{s j k}\right|^{\frac{p}{3}} & \leq \frac{p}{3} \sum\limits_{i = s}^{m-1}\left|u_{i+1, j k}-u_{i j k}\right|\left(\left|u_{i j k}\right|^{\frac{p}{3}-1}+\left|u_{i+1, j k}\right|^{\frac{p}{3}-1}\right) \\ & = p h_1 \sum\limits_{i = s}^{m-1}\left|\delta_{x_1} u_{i j k}\right| \frac{\left|u_{i j k}\right|^{\frac{p}{3}-1}+\left|u_{i+1, j k}\right|^{\frac{p}{3}-1}}{2} \\ & \leq p\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|u_{i j k}\right|^{\frac{2p}{3}-2}\right)^{\frac{1}{2}}\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|\delta_{x_1} u_{i j k}\right|^2\right)^{\frac{1}{2}} . \end{aligned}

    It is easy to verify the above inequality also holds for m \leq s . Thus, we have

    \left|u_{m j k}\right|^{\frac{p}{3}} \leq p\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|u_{i j k}\right|^{\frac{2p}{3}-2}\right)^{\frac{1}{2}}\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|\delta_{x_1} u_{i j k}\right|^2\right)^{\frac{1}{2}}+\left|u_{s j k}\right|^{\frac{p}{3}}, \quad \forall 1 \leq m, s \leq M_1-1.

    Multiplying the above inequality by h_1 and summing up for s from 1 to M_1-1 , we have

    l_1\left|u_{m j k}\right|^{\frac{p}{3}} \leq l_1 p\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|u_{i j k}\right|^{\frac{2p}{3}-2}\right)^{\frac{1}{2}}\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|\delta_{x_1} u_{i j k}\right|^2\right)^{\frac{1}{2}}+h_1 \sum\limits_{i = 1}^{M_1-1}\left|u_{i j k}\right|^{\frac{p}{3}} .

    Dividing the result by l_1 , and noticing that the above inequality holds for m = 1, 2, \ldots, M_1-1 , we have

    \max _{1 \leq m \leq M_1-1}\left|u_{m j k}\right|^{\frac{p}{3}} \leq p\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|u_{i j k}\right|^{\frac{2p}{3}-2}\right)^{\frac{1}{2}}\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|\delta_{x_1} u_{i j k}\right|^2\right)^{\frac{1}{2}}+\frac{1}{l_1} h_1 \sum\limits_{i = 1}^{M_1-1}\left|u_{i j k}\right|^{\frac{p}{3}} .

    Multiplying the above inequality by h_2h_3 and summing over j, k , then applying the Cauchy-Schwarz inequality, we obtain

    \begin{align} &\quad h_2h_3 \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1} \max _{1 \leq i \leq M_1-1}\left|u_{i j k}\right|^{\frac{p}{3}} \\ & \leq p h_2h_3 \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1}\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|u_{i j k}\right|^{\frac{2p}{3}-2}\right)^{\frac{1}{2}}\left(h_1 \sum\limits_{i = 1}^{M_1-1}\left|\delta_{x_1} u_{i j k}\right|^2\right)^{\frac{1}{2}}+\frac{1}{l_1}\left(\|u\|_{\frac{p}{3}}\right)^{\frac{p}{3}} \\ & \leq p\left(h_2h_3 \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1} h_1 \sum\limits_{i = 1}^{M_1-1}\left|u_{i j k}\right|^{\frac{2p}{3}-2}\right)^{\frac{1}{2}}\left(h_2h_3 \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1} h_1 \sum\limits_{i = 1}^{M_1-1}\left|\delta_{x_1} u_{i j k}\right|^2\right)^{\frac{1}{2}}+\frac{1}{l_1}\left(\|u\|_{\frac{p}{3}}\right)^{\frac{p}{3}} \\ & = p\left(\|u\|_{\frac{2p}{3}-2}\right)^{\frac{p}{3}-1} \cdot\left\|\delta_{x_1} u\right\|+\frac{1}{l_1}\left(\|u\|_{\frac{p}{3}}\right)^{\frac{p}{3}}. \end{align} (A.1)

    Multiply both sides of inequality (A.1) by (h_2h_3)^{\frac{1}{2}} , it follows easily that there exists a constant C such that (h_2h_3)^{\frac{1}{2}}\leq C , we obtain

    \begin{align} (h_2h_3)^{\frac{1}{2}} \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1} \max _{1 \leq i \leq M_1-1}\left|u_{i j k}\right|^{\frac{p}{3}}\leq Cp\left(\|u\|_{\frac{2p}{3}-2}\right)^{\frac{p}{3}-1} \cdot\left\|\delta_{x_1} u\right\|+\frac{C}{l_1}\left(\|u\|_{\frac{p}{3}}\right)^{\frac{p}{3}}. \end{align} (A.2)

    Similarly to the previous analysis, we have

    \begin{align} &(h_1h_3)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{k = 1}^{M_3-1} \max _{1 \leq j \leq M_2-1}\left|u_{i j k}\right|^{\frac{p}{3}}\leq Cp\left(\|u\|_{\frac{2p}{3}-2}\right)^{\frac{p}{3}-1} \cdot\left\|\delta_{x_2} u\right\|+\frac{C}{l_2}\left(\|u\|_{\frac{p}{3}}\right)^{\frac{p}{3}}. \end{align} (A.3)
    \begin{align} &(h_1h_2)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1} \max _{1 \leq k \leq M_3-1}\left|u_{i j k}\right|^{\frac{p}{3}}\leq Cp\left(\|u\|_{\frac{2p}{3}-2}\right)^{\frac{p}{3}-1} \cdot\left\|\delta_{x_3} u\right\|+\frac{C}{l_3}\left(\|u\|_{\frac{p}{3}}\right)^{\frac{p}{3}}. \end{align} (A.4)

    Using the Cauchy-Schwarz inequality, we have

    \left(\|u\|_{\frac{p}{3}}\right)^{\frac{p}{3}} = h_1 h_2 h_3 \sum\limits_{i = 1}^{M_1-1} \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^{\frac{p}{3}} = h_1 h_2 h_3\sum\limits_{i = 1}^{M_1-1} \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right| \cdot\left|u_{i j k}\right|^{\frac{p}{3}-1} \leq\|u\| \cdot\left(\|u\|_{\frac{2p}{3}-2}\right)^{\frac{p}{3}-1} .

    Substituting the above inequality into inequalities (A.2)–(A.4), we have

    \begin{align} &(h_2h_3)^{\frac{1}{2}} \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1} \max _{1 \leq i \leq M_1-1}\left|u_{i j k}\right|^{\frac{p}{3}}\leq C\left(\|u\|_{\frac{2p}{3}-2}\right)^{\frac{p}{3}-1} \cdot\left(p\left\|\delta_{x_1} u\right\|+\frac{1}{l_1}\|u\|\right). \end{align} (A.5)
    \begin{align} &(h_1h_3)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{k = 1}^{M_3-1} \max _{1 \leq j \leq M_2-1}\left|u_{i j k}\right|^{\frac{p}{3}}\leq C\left(\|u\|_{\frac{2p}{3}-2}\right)^{\frac{p}{3}-1} \cdot\left(p\left\|\delta_{x_2} u\right\|+\frac{1}{l_2}\|u\|\right). \end{align} (A.6)
    \begin{align} &(h_1h_2)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1} \max _{1 \leq k \leq M_3-1}\left|u_{i j k}\right|^{\frac{p}{3}}\leq C\left(\|u\|_{\frac{2p}{3}-2}\right)^{\frac{p}{3}-1} \cdot\left(p\left\|\delta_{x_3} u\right\|+\frac{1}{l_3}\|u\|\right). \end{align} (A.7)

    We now estimate \|u\|_p^p ,

    \begin{align} &\quad \|u\|_p^p = h_1 h_2 h_3 \sum\limits_{i = 1}^{M_1-1} \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^p \\ & = (h_1h_2)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1}\left((h_1h_3)^{\frac{1}{2}} (h_2h_3)^{\frac{1}{2}}\sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^{\frac{2p}{3}}\left|u_{i j k}\right|^{\frac{p}{3}}\right) \\ & \leq (h_1h_2)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1}\left(\max _{1 \leq k \leq M_3-1}\left|u_{i j k}\right|^{\frac{p}{3}}\cdot(h_1h_3)^{\frac{1}{2}} (h_2h_3)^{\frac{1}{2}}\sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^{\frac{2p}{3}}\right) \\ & \leq\left((h_1h_2)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1}\left(\max _{1 \leq k \leq M_3-1}\left|u_{i j k}\right|^{\frac{p}{3}}\right)\right)\cdot\left(\sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1}(h_1h_3)^{\frac{1}{2}} (h_2h_3)^{\frac{1}{2}}\sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^{\frac{2p}{3}}\right) \\ & \leq\left((h_1h_2)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1}\left(\max _{1 \leq k \leq M_3-1}\left|u_{i j k}\right|^{\frac{p}{3}}\right)\right)\cdot\left((h_2h_3)^{\frac{1}{2}} \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1}\left(\max _{1 \leq i \leq M_1-1}\left|u_{i j k}\right|^{\frac{p}{3}}\right)\right) \\ & \ \ \ \cdot\left((h_1h_3)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1}\left|u_{i j k}\right|^{\frac{p}{3}}\right) \\ & \leq\left((h_1h_2)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{j = 1}^{M_2-1}\left(\max _{1 \leq k \leq M_3-1}\left|u_{i j k}\right|^{\frac{p}{3}}\right)\right)\cdot\left((h_2h_3)^{\frac{1}{2}} \sum\limits_{j = 1}^{M_2-1}\sum\limits_{k = 1}^{M_3-1}\left(\max _{1 \leq i \leq M_1-1}\left|u_{i j k}\right|^{\frac{p}{3}}\right)\right) \\ & \ \ \ \cdot\left((h_1h_3)^{\frac{1}{2}} \sum\limits_{i = 1}^{M_1-1}\sum\limits_{k = 1}^{M_3-1}\left(\max _{1 \leq j \leq M_2-1}\left|u_{i j k}\right|^{\frac{p}{3}}\right)\right) \\ & \leq C^3\left(\|u\|_{\frac{2p}{3}-2}\right)^{p-3} \cdot\left(p\left\|\delta_{x_1} u\right\|+\frac{1}{l_1}\|u\|\right) \cdot\left(p\left\|\delta_{x_2} u\right\|+\frac{1}{l_2}\|u\|\right)\cdot\left(p\left\|\delta_{x_3} u\right\|+\frac{1}{l_3}\|u\|\right), \end{align} (A.8)

    the last inequality is obtained by inequalities (A.5)–(A.7).

    In addition, we set l = \min \left\{l_1, l_2, l_3\right\} , by using mean value inequality then we have

    \begin{align} &\quad \left(p\left\|\delta_{x_1} u\right\|+\frac{1}{l_1}\|u\|\right) \cdot\left(p\left\|\delta_{x_2} u\right\|+\frac{1}{l_2}\|u\|\right)\cdot\left(p\left\|\delta_{x_3} u\right\|+\frac{1}{l_3}\|u\|\right)\\ &\leq \left(p\left\|\delta_{x_1} u\right\|+\frac{1}{l}\|u\|\right) \cdot\left(p\left\|\delta_{x_2} u\right\|+\frac{1}{l}\|u\|\right)\cdot\left(p\left\|\delta_{x_3} u\right\|+\frac{1}{l}\|u\|\right)\\ &\leq p^3\left\|\delta_{x_1} u\right\| \cdot\left\|\delta_{x_2} u\right\| \cdot\left\|\delta_{x_3} u\right\|+\frac{p}{l^2}\|u\|^2 \cdot\left(\left\|\delta_{x_1} u\right\|+\left\|\delta_{x_2} u\right\|+\left\|\delta_{x_3} u\right\|\right) \\ &\quad +\frac{p^2}{l}\|u\| \cdot\left(\left\|\delta_{x_1} u\right\| \cdot\left\|\delta_{x_3} u\right\|+\left\|\delta_{x_2} u\right\| \cdot\left\|\delta_{x_3} u\right\|+\left\|\delta_{x_1} u\right\| \cdot\left\|\delta_{x_2} u\right\|\right)+\frac{1}{l^3}\|u\|^3 \\ &\leq \left(\frac{p}{\sqrt{3}}\right)^3 \cdot\left(\left\|\delta_{x_1} u\right\|^2+\left\|\delta_{x_2} u\right\|^2+\left\|\delta_{x_3} u\right\|^2\right)^{\frac{3}{2}}+\frac{\sqrt{3} p}{l^2}\|u\|^2 \cdot\left(\left\|\delta_{x_1} u\right\|^2+\left\|\delta_{x_2} u\right\|^2+\left\|\delta_{x_3} u\right\|^2\right)^{\frac{1}{2}} \\ &\quad +\frac{p^2}{l}\|u\| \cdot\left(\left\|\delta_{x_1} u\right\|^2+\left\|\delta_{x_2} u\right\|^2+\left\|\delta_{x_3} u\right\|^2\right)+\frac{1}{l^3}\|u\|^3 \\ & = \left(\frac{p}{\sqrt{3}}\right)^3 |u|_{H^{1}}^3+\frac{\sqrt{3} p}{l^2}|u|_{H^{1}} \cdot \| u\|^2+\frac{p^2}{l}|u|_{H^{1}}^2 \cdot\| u\|+\frac{1}{l_3}\| u \|^3 \\ & = \left(\frac{p}{\sqrt{3}}|u|_{H^{1}}+\frac{1}{l} \| u\|\right)^3. \end{align} (A.9)

    Combining inequalities (A.8) and (A.9) yields

    \begin{align} \left(\|u\|_p\right)^p \leq C^3\left(\|u\|_{\frac{2p}{3}-2}\right)^{p-3} \cdot\left(\frac{p}{\sqrt{3}}|u|_{H^1}+\frac{1}{l}\|u\|\right)^3 . \end{align} (A.10)

    We consider the case p \geq 6 , applying Lemma 9 for p \geq 6 , it holds

    \left(\|u\|_{\frac{2p}{3}-2}\right)^{p-3} \leq\|u\|^{\frac{p+6}{p-2}}\left(\|u\|_p\right)^{\frac{p(p-6)}{p-2}} .

    Substituting the above inequality into inequality (A.10), we get

    \left(\|u\|_p\right)^{\frac{4 p}{p-2}} \leq C^3\|u\|^{\frac{p+6}{p-2}} \cdot\left(\frac{p}{\sqrt{3}}|u|_{H^1}+\frac{1}{l}\|u\|\right)^3 ,

    that is

    \begin{align} \|u\|_p \leq C^3\|u\|^{\frac{p+6}{4p}} \cdot\left(\frac{p}{\sqrt{3}}|u|_{H^1}+\frac{1}{l}\|u\|\right)^{\frac{3 p-6}{4p}} . \end{align} (A.11)

    Thus, we have proved the result for p \geq 6 . Taking p = 6 in inequality (A.11) yields

    \begin{align} \|u\|_6 \leq C^3\|u\|^{\frac{1}{2}} \cdot\left(2\sqrt{3}|u|_{H^1}+\frac{1}{l}\|u\|\right)^{\frac{1}{2}} . \end{align} (A.12)

    When 2 < p < 6 , using Lemma 9 and inequality (A.12), we have

    \begin{aligned} \|u\|_p \leq\|u\|^{\frac{6-p}{2p}}\|u\|_6^{\frac{3(p-2)}{4p}} & \leq C^3\|u\|^{\frac{6-p}{2p}}\left[\|u\|^{\frac{1}{2}} \cdot\left(2\sqrt{3}|u|_{H^1}+\frac{1}{l}\|u\|\right)^{\frac{1}{2}}\right]^{\frac{3p-6}{4p}} \\ & = C^3\|u\|^{\frac{p+6}{4p}}\left(2 \sqrt{3}|u|_{H^1}+\frac{1}{l}\|u\|\right)^{\frac{3p-6}{4p}} . \end{aligned}

    This completes the proof.

    [1] Akorede MF, Hizam H, Pouresmaeil E (2010) Distributed energy resources and benefits to the environment. Renew Sust Energ Rev 14: 724–734. doi: 10.1016/j.rser.2009.10.025
    [2] Arnold W, Reitze J (2010) Electric power in a carbon constrained world. William Mary Environ Law Policy Rev 34: 821.
    [3] Kumar A, Kandpal TC (2007) Potential and cost of CO2 emissions mitigation by using solar photovoltaic pumps in India. Int J Sol Energy 26: 159–166.
    [4] Kumar A, Kandpal TC (2007) CO2 emissions mitigation potential of some renewable energy technologies in India. Energ Source 29: 1203–1214. doi: 10.1080/009083190965343
    [5] Purohit P, Kumar A, Kandpal TC (2002) Renewable energy technologies for domestic cooking in India: Estimation of CO2 emissions mitigation potential. Int J Amb Energ 23: 127–135. doi: 10.1080/01430750.2002.9674881
    [6] Kumar A, Kandpal TC (2005) Solar drying and CO2 emissions mitigation: Potential for selected cash crops in India. Sol Energy 78: 321–329. doi: 10.1016/j.solener.2004.10.001
    [7] Harmsen R, Moth L, Kumar A (2014) Applicability of energy saving obligations to Indian electricity efficiency efforts. Energy Strateg Rev 2: 298–306. doi: 10.1016/j.esr.2013.07.001
    [8] REN 21. Renewables 2017: Global Status Report. (Paris: REN21 Secretariat). 2017.
    [9] REN 21. Renewables 2014: Global Status Report. (Paris: REN21 Secretariat). 2014.
    [10] CIA, The World Factbook-Central Intelligence Agency, 2017. Available from: https://www.cia.gov/library/publications/the-world-factbook/geos/ni.html (accessed July 17, 2017).
    [11] ECN, Nigeria Energy Calculator 2050 (NECAL 2050) report 2015:105. Available from: http://www.energy.gov.ng/index.php?option=com_docman&task=doc_view&gid=123&Itemid=49 (accessed October 14, 2017).
    [12] IEA, World Energy Outlook 2017. Special Report: Energy Access Outlook. Paris: International Energy Agency (IEA) and the Organisation of Economic Co-Operation and Development (OECD); 2017.
    [13] Boz MB, Calvert K, Brownson JRS (2015) An automated model for rooftop PV systems assessment in ArcGIS using LIDAR. AIMS Energy 3: 401–420. doi: 10.3934/energy.2015.3.401
    [14] Brito MC, Gomes N, Santos T, et al. (2012) Photovoltaic potential in a Lisbon suburb using LiDAR data. Sol Energy 86: 283–288. doi: 10.1016/j.solener.2011.09.031
    [15] Wiginton LK, Nguyen HT, Pearce JM (2010) Quantifying rooftop solar photovoltaic potential for regional renewable energy policy. Comput Environ Urban Syst 34: 345–357. doi: 10.1016/j.compenvurbsys.2010.01.001
    [16] Izquierdo S, Rodrigues M, Fueyo N (2008) A method for estimating the geographical distribution of the available roof surface area for large-scale photovoltaic energy-potential evaluations. Sol Energy 82: 929–939. doi: 10.1016/j.solener.2008.03.007
    [17] Kumar D, Shekhar S (2014) Computing building rooftop solar potential for photovoltaics. ISPRS TC VIII Int. Symp. Oper. Remote Sens. Appl. Oppor. Prog. Challenges Hyderabad, India, December 9–12, 2014.
    [18] Johnson DO, Ogunseye AA (2017) Grid-connected photovoltaic system design for local government offices in Nigeria. Niger J Technol 36: 571–581. doi: 10.4314/njt.v36i2.33
    [19] Adaramola MS (2014) Viability of grid-connected solar PV energy system in Jos, Nigeria. Int J Electr Power Energ Syst 61: 64–69. doi: 10.1016/j.ijepes.2014.03.015
    [20] Akinyele DO, Rayudu RK (2014) Distributed Photovoltaic Power Generation for Energy-Poor Households: The Nigerian Perspective. Power Energy Eng Conf 2014: 1–6
    [21] Udoh SP, Umoren AM, Okpura NI (2017) Techno-economic analysis of building rooftop photovoltaic power system for lecture hall at Imo State. Renew Energ Res 1: 8–16.
    [22] Adaramola MS, Paul SS (2017) Economic analysis and potential feed-in tariff of grid-connected PV systems in Nigeria. Environ Prog Sustain 36: 305–314. doi: 10.1002/ep.12502
    [23] NPBR, Nigeria Power Baseline Report, 2015. Available from: http://www.nesistats.org/uploads/3/6/3/6/3636925/20150916_nigeria_energy_power_report_final.pdf (accessed October 12, 2017).
    [24] NESP. The Nigerian energy sector: An overview with a special emphasis on renewable energy, energy efficiency and rural electrification. 2015.
    [25] GENI, Global Energy Network Institute: Solar Energy in Nigeria n.d. Available from: http://www.geni.org/globalenergy/library/renewable-energy-resources/world/africa/solar-africa/solar-nigeria.shtml (accessed July 26, 2017).
    [26] Akorede MF, Ibrahim O, Amuda SA, et al. (2017) Current status and outlook of renewable energy development in Nigeria. Niger J Technol 36: 196–212.
    [27] Chineke TC, Okoro UK (2010) Application of Sayigh "Universal Formula" for global solar radiation estimation in the Niger Delta region of Nigeria. Renew Energ 35: 734–749. doi: 10.1016/j.renene.2009.08.010
    [28] Shaaban M, Petinrin J (2014) Renewable energy potentials in Nigeria: Meeting rural energy needs. Renew Sust Energ Rev 29: 72–84. doi: 10.1016/j.rser.2013.08.078
    [29] Ohunakina OS, Adaramola MS, Oyewola OM (2014) Solar energy applications and development in Nigeria: Drivers and barriers. Renew Sust Energ Rev 32: 294–301. doi: 10.1016/j.rser.2014.01.014
    [30] ECN, Research Centres, 2012. Available from: http://www.energy.gov.ng/index.php?option=com_content&view=article&id=54 (accessed July 12, 2017).
    [31] NREEEP, Approved National Renewable Energy and Energy Efficiency Policy (NREEEP) for the power sector, 2015. Available from: http://www.power.gov.ng/download/NREEE POLICY 2015- FEC APPROVED COPY.pdf (accessed October 15, 2017).
    [32] NIPC, Nigeria Investment Promotion Commission, Investment Incentives, 2017. Available from: http://www.invest-nigeria.com/investment-incentives/ (accessed July 19, 2017).
    [33] NTA, Acting President Launches Solar Home Systems In Wuna Village 2017. Available from: http://www.nta.ng/news/20170131-acting-president-solar-home-systems-wuna/ (accessed July 19, 2017).
    [34] ESI AFRICA, Nigerian minister unveils rooftop solar energy project 2017. Available from: https://www.esi-africa.com/news/nigeria-unveils-rooftop-solar-energy-project/ (accessed July 19, 2017).
    [35] Voivontas D, Assimacopoulos D, Mourelatos A, et al. (1998) Evaluation of renewable energy potential using a GIS decision support system. Renew Energ 13: 333–344. doi: 10.1016/S0960-1481(98)00006-8
    [36] Nguyen QK (2005) Long term optimization of energy supply and demand in Vietnam with special reference to the potential of renewable energy. Von der Carl von Ossietzky Universität Oldenburg, 2005.
    [37] TERI, Reaching the sun with rooftop solar 2014. Available from: http://mnre.gov.in/file-manager/UserFiles/Rooftop-SPV-White-Paper-low.pdf (accessed July 19, 2017).
    [38] Sun YW, Hof A, Wang R, et al. (2013) GIS-based approach for potential analysis of solar PV generation at the regional scale: A case study of Fujian Province. Energ Policy 58: 248–259. doi: 10.1016/j.enpol.2013.03.002
    [39] Lopez A, Roberts B, Heimiller D, et al. (2012) U.S. renewable energy technical potentials: A GIS-based analysis. Contract 82: 40
    [40] Gagnon P, Margolis R, Melius J, et al. (2016) Rooftop solar photovolatic technical potential in the united states: A detailed assessment.
    [41] Zhang X, Walker R, Salisbury M, et al. (2009) Creating a solar city: Determining the potential of solar rooftop systems in the city of Newark. Newark, DE: University of Delaware, Center for Energy and Environmental Policy.
    [42] Bright J, Burman K (2010) Portland Ecodistrict Solar Site Assessment. Golden, CO: National Renewable Energy Laboratory (for Solar Cities America).
    [43] Ntsoane M (2017) Rooftop solar PV potential assessment in the city of Johannesburg. Stellenbosch University.
    [44] NPC, Population and Housing Census of the Federal Republic of Nigeria 2006. Available from: http://www.population.gov.ng/ (accessed July 12, 2017).
    [45] UN, World Population Prospects-Population Division-United Nations, 2015. Available from: https://esa.un.org/unpd/wpp/ (accessed April 25, 2017).
    [46] NBS, Harmonized Nigeria living standard survey 2009/10: Core welfare indicator questinnaire survey 2009 (Part A), 2010. Available from: http://www.nigerianstat.gov.ng/pdfuploads/HARMONIZED NIGERIA LIVING STANDARD SURVEY 2009 Part A.pdf (accessed May 2, 2017).
    [47] NBC, National Building Code, 2006. Available from: http://sdngnet.com/Files/Lectures/FUTA-ARC-807-Professional_Practice_and_Procedure/CD 2013-2014/National Building Code of Nigeria 2006.pdf (accessed July 6, 2017).
    [48] TEDA, Guidelines for grid-connected small scale (Rooftop) solar PV systems for Tamil Nadu. TEDA: Tamil Nadu Energy Development Agency, 2014.
    [49] Quaschning V, Hanitsch R (1996) Numerical simulation of current-voltage characteristics of photovoltaic systems with shaded solar cells. Sol Energy 56: 513–520. doi: 10.1016/0038-092X(96)00006-0
    [50] Anigstein PA, Pena RSS (1998) Analysis of solar panel orientation in low altitude satellites. IEEE T Aero Elec Sys 34: 569–578. doi: 10.1109/7.670337
    [51] Quaschning V, Hanitsch R (1998) Irradiance calculation on shaded surfaces. Sol Energy 62: 369–375. doi: 10.1016/S0038-092X(98)00018-8
    [52] Abubakar IR (2014) Abuja city profile. Cities 41: 81–91. doi: 10.1016/j.cities.2014.05.008
    [53] Nguyen HT, Pearce JM (2012) Incorporating shading losses in solar photovoltaic potential assessment at the municipal scale. Sol Energy 86: 1245–1260. doi: 10.1016/j.solener.2012.01.017
    [54] Pillai IR, Banerjee R (2007) Methodology for estimation of potential for solar water heating in a target area. Sol Energy 81: 162–172. doi: 10.1016/j.solener.2006.04.009
    [55] Real Estate and Property in Nigeria for Sale and Rent-Nigerian Real Estate &amp; Property-ToLet.com.ng n.d. Available from: https://www.tolet.com.ng/ (accessed July 28, 2017).
    [56] NREL, PVWatts Calculator 2017. Available from: http://pvwatts.nrel.gov/pvwatts.php (accessed August 17, 2017).
    [57] SolarWorld, Sunmodule Plus SW 260 poly n.d. Available from: http://www.tehnosat.ro/pdf/PVmodules/SW_poly_260.pdf (accessed July 19, 2017).
    [58] Mermoud A, Wittmer B. PVSYST user's manual. Switzerland, January 2014.
    [59] Nouni MR, Mullick SC, Kandpal TC (2006) Photovoltaic projects for decentralized power supply in India: A financial evaluation. Energ Policy 34: 3727–3738. doi: 10.1016/j.enpol.2005.08.015
    [60] Purohit I, Purohit P, Shekhar S (2013) Evaluating the potential of concentrating solar power generation in Northwestern India. Energ Policy 62: 157–175. doi: 10.1016/j.enpol.2013.06.069
    [61] NERC, Multi­year tariff order for the determination of the cost of electricity generation for the period 1 June 2012 to 31 May 2017 2012. Available from: http://www.ecowrex.org/system/files/documents/2012_multiyear-tariff-order-generation_nerc.pdf (accessed August 17, 2017).
    [62] MNRE, Guidelines for smooth implementation of grid connected solar rooftop projects by SNAs/Discoms 2016. Available from: http://mnre.gov.in/file-manager/UserFiles/gcrt-guide-080916.pdf (accessed August 19, 2017).
    [63] CBN, Central Bank of Nigeria | Home 2017. Available from: https://www.cbn.gov.ng/ (accessed August 19, 2017).
    [64] TERI, Technical manual for banks & FIs on grid-connected rooftop solar power 2015. Available from: http://mnre.gov.in/file-manager/UserFiles/TERI-Technical-Manual-Banks-FIs.pdf (accessed August 11, 2017).
    [65] FGN, Electric power sector reform act, 2005. Available from: http://www.power.gov.ng/download/Electric Power Sector Reform Act 2005.pdf (accessed August 19, 2017).
    [66] NERC, Nigerian Electricity Regulatory Commission 2017. Available from: http://www.nercng.org/ (accessed September 8, 2017).
  • This article has been cited by:

    1. Dongdong Hu, Fully decoupled and high-order linearly implicit energy-preserving RK-GSAV methods for the coupled nonlinear wave equation, 2024, 445, 03770427, 115836, 10.1016/j.cam.2024.115836
  • Reader Comments
  • © 2018 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(6828) PDF downloads(1666) Cited by(17)

Figures and Tables

Figures(7)  /  Tables(11)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog