Loading [MathJax]/extensions/TeX/boldsymbol.js
Research article Special Issues

Statistical inference of an α-quantile past lifetime function with applications

  • Received: 16 February 2024 Revised: 04 April 2024 Accepted: 23 April 2024 Published: 28 April 2024
  • MSC : 62N01, 62N05

  • In reliability engineering and survival analysis, quantile functions are fundamental and often the most natural way to represent probability distributions and data samples. In this paper, the α-quantile function of past lifetime was estimated for right-censored data by applying the Kaplan-Meier survival estimator. The weak convergence of the proposed estimator to a Gaussian process was investigated. A confidence interval for the α-quantile of the past life function that does not depend on the density function was proposed. The strong convergence of the estimator to a Gaussian process was also discussed. The properties of the estimator and the confidence interval were investigated in a simulation study. Finally, two real datasets were analyzed.

    Citation: Mohamed Kayid. Statistical inference of an α-quantile past lifetime function with applications[J]. AIMS Mathematics, 2024, 9(6): 15346-15360. doi: 10.3934/math.2024745

    Related Papers:

    [1] Attaullah, Mansour F. Yassen, Sultan Alyobi, Fuad S. Al-Duais, Wajaree Weera . On the comparative performance of fourth order Runge-Kutta and the Galerkin-Petrov time discretization methods for solving nonlinear ordinary differential equations with application to some mathematical models in epidemiology. AIMS Mathematics, 2023, 8(2): 3699-3729. doi: 10.3934/math.2023185
    [2] Yayun Fu, Mengyue Shi . A conservative exponential integrators method for fractional conservative differential equations. AIMS Mathematics, 2023, 8(8): 19067-19082. doi: 10.3934/math.2023973
    [3] Yu He, Jianing Yang, Theodore E. Simos, Charalampos Tsitouras . A novel class of Runge-Kutta-Nyström pairs sharing orders 8(6). AIMS Mathematics, 2024, 9(2): 4882-4895. doi: 10.3934/math.2024237
    [4] Sara S. Alzaid, Pawan Kumar Shaw, Sunil Kumar . A numerical study of fractional population growth and nuclear decay model. AIMS Mathematics, 2022, 7(6): 11417-11442. doi: 10.3934/math.2022637
    [5] Muhammad Bilal Riaz, Syeda Sarwat Kazmi, Adil Jhangeer, Jan Martinovic . Unveiling solitons and dynamic patterns for a (3+1)-dimensional model describing nonlinear wave motion. AIMS Mathematics, 2024, 9(8): 20390-20412. doi: 10.3934/math.2024992
    [6] Xintian Pan . A high-accuracy conservative numerical scheme for the generalized nonlinear Schrödinger equation with wave operator. AIMS Mathematics, 2024, 9(10): 27388-27402. doi: 10.3934/math.20241330
    [7] M. Adams, J. Finden, P. Phoncharon, P. H. Muir . Superconvergent interpolants for Gaussian collocation solutions of mixed order BVODE systems. AIMS Mathematics, 2022, 7(4): 5634-5661. doi: 10.3934/math.2022312
    [8] Shumoua F. Alrzqi, Fatimah A. Alrawajeh, Hany N. Hassan . An efficient numerical technique for investigating the generalized Rosenau–KDV–RLW equation by using the Fourier spectral method. AIMS Mathematics, 2024, 9(4): 8661-8688. doi: 10.3934/math.2024420
    [9] Nurain Zulaikha Husin, Muhammad Zaini Ahmad . Hybridization of the shooting and Runge-Kutta Cash-Karp methods for solving Fuzzy Boundary Value Problems. AIMS Mathematics, 2024, 9(11): 31806-31847. doi: 10.3934/math.20241529
    [10] Dan-Dan Dai, Wei Zhang, Yu-Lan Wang . Numerical simulation of the space fractional (3+1)-dimensional Gray-Scott models with the Riesz fractional derivative. AIMS Mathematics, 2022, 7(6): 10234-10244. doi: 10.3934/math.2022569
  • In reliability engineering and survival analysis, quantile functions are fundamental and often the most natural way to represent probability distributions and data samples. In this paper, the α-quantile function of past lifetime was estimated for right-censored data by applying the Kaplan-Meier survival estimator. The weak convergence of the proposed estimator to a Gaussian process was investigated. A confidence interval for the α-quantile of the past life function that does not depend on the density function was proposed. The strong convergence of the estimator to a Gaussian process was also discussed. The properties of the estimator and the confidence interval were investigated in a simulation study. Finally, two real datasets were analyzed.



    The nonlinear partial differential equations have significant roles in a variety of fields in engineering and science (see, e.g., [4,39]), including quantum field theory, nonlinear optics, propagation of dislocations in crystals, nucleation, and solid state physics. In this paper, we consider the following nonlinear wave equation in two dimensions:

    2ut2κ2(2ux2+2uy2)=f(u(x,y,t)),(x,y)(0,L1)×(0,L2),t[t0,T], (1.1)

    subject to the initial conditions

    u(x,y,t0)=φ0(x,y),ut(x,y,t0)=φ1(x,y),x[0,L1]×[0,L2], (1.2)

    and the periodic boundary conditions

    u(x,y,t)=u(x+L1,y,t),u(x,y,t)=u(x,y+L2,t),(x,y)ˉΩ,t[t0,T], (1.3)

    where κ2 is a dimensionless positive parameter, φ0(x,y) and φ1(x,y) are the given (L1,L2)- periodic functions, and L1 and L2 are the basic positive periods. In the literature, many works have been made to explore the analytical solution for the nonlinear wave equations (see, e.g., [1,39]). However, it is difficult to obtain the general exact solutions for all the nonlinear wave equations. Therefore, the development of efficient and high-precision numerical methods for solving the two-dimensional nonlinear wave equations has became much more important. A great number of excellent numerical strategies have been proposed to study the nonlinear wave equations, including the finite difference methods (see, e.g., [12,17,18]), the finite element methods (see, e.g., [2,3,33]), the spectral methods [25], the domain decomposition methods [19], and the radial basis functions methods [11].

    If the nonlinear function f(u) is the negative derivative of a nonnegative function F(u), i.e., f(u)=dF(u)du, and the solution of (1.1) satisfies (u,ut)H1(Ω)×L2(Ω), then the nonlinear wave Eqs (1.1)–(1.3) could conserve the energy

    E(t):=12Ω(u2t(x,y,t)+κ2|u(x,y,t)|2+2F(u(x,y,t)))dxdy12Ω(φ21(x,y)+κ2|φ0(x,y)|2+2F(φ0(x,y)))dxdy=E(t0),tt0. (1.4)

    The energy conservation (1.4) is a significant property of the nonlinear wave equations, and plays prominent roles in investigating soliton theory. Under this case, the nonlinear wave equations like (1.1) are called nonlinear Hamiltonian wave equations. We know that the energy conservation along the exact flow is one most characteristic properties of the nonlinear Hamiltonian wave Eq (1.1). The energy-conserving numerical schemes usually yield correct physical phenomenons and numerical stability (see, e.g., [24,29,32]). Therefore, it will be meaningful to design suitable numerical schemes which could exactly preserve the discrete energy and symmetry of the nonlinear Hamiltonian wave Eq (1.1).

    The development of energy-preserving numerical schemes for nonlinear Hamiltonian wave equations has garnered significant attention across various fields of mechanics. For example, Li et al. [21] proposed several finite difference schemes that preserve specific algebraic invariants of the nonlinear Klein-Gordon equations. Moreover, based on the concept of the discrete line integral method (see, e.g., [5,6]), L. Brugnano et al. developed the energy-preserving Hamiltonian boundary value methods (HBVMs) to solve the nonlinear Hamiltonian PDEs (see, e.g., [7,8]). The energy-preserving average vector field (AVF) method was initially developed for solving Hamiltonian ordinary differential equations (ODEs). Recently, the AVF method, when combined with appropriate spatial semi-discretization techniques, has been utilized to numerically investigate nonlinear Hamiltonian wave equations, thereby attracting significant attention from researchers. For instance, AVF finite element methods have been introduced to solve one-dimensional Hamiltonian wave equations (see [10,33]). In [31,32], combining the AVF method with the spatial fourth-order finite difference semidiscretisation, the authors developed energy-preserving methods for one- and two- dimensional Hamiltonian wave equations with Neumann boundary conditions. However, the previous schemes have only second-order accuracy in time. To enhance temporal accuracy, Hou et al. [17,18] integrated the fourth-order AVF temporal integrator with spatial compact finite difference (CFD) discretization to construct and analyze high-order energy-preserving schemes for solving one- and two-dimensional nonlinear wave equations with variable coefficients. This represents a significant advancement in the field of energy-preserving methods for nonlinear Hamiltonian wave equations. Building on these contributions, we aim to develop and analyze a high-order energy-preserving and symmetric scheme for two-dimensional nonlinear Hamiltonian wave equations by combining the continuous-stage Runge-Kutta-Nyström time integrator with Fourier pseudo-spectral spatial discretization.

    The rest of the paper is organized as follows: In Section 2, the two-dimensional nonlinear wave Eq (1.1) will first be reformulated as an abstract infinite-dimensional separable Hamiltonian ODE system in an appropriate function space. Then, the application of a continuous-stage Runge-Kutta-Nyström time integrator to the yielded ODE's system to derive the time-stepping scheme is presented. The energy preservation and symmetry of the proposed time-stepping scheme will be investigated. Furthermore, by approximating the spatial differential operator using the two-dimensional Fourier pseudo-spectral method, we derive a fully discrete scheme. A rigorous analysis of the energy conservation properties of this scheme is then conducted. The error analysis demonstrates that the proposed scheme achieves sixth-order accuracy in the relatively low regularity function space C2([t0,T],B). Numerical experiments are presented in Section 4. Lastly, a concise conclusion is provided in Section 5.

    In this section, we will first represent the two-dimensional nonlinear wave Eqs (1.1)–(1.3) as an abstract nonlinear ODE on an appropriate infinite-dimensional Hilbert space. Then, we will develop and analyze a novel energy-conserving time-stepping scheme for the abstract ODE.

    According to the analysis in references [30,37], by defining the mapping

    u(t):=[(x,y)u(x,y,t)],

    we can express the nonlinear wave Eqs (1.1)–(1.3) as the following abstract ODE (e.g., [29,30,38]):

    {u(t)=Au(t)+f(u(t)):≜g(u(t)),t[t0,T],u(t0)=φ0(x,y),u(t0)=φ1(x,y),(x,y)ˉΩ, (2.1)

    where A is the linear differential operator

    Au(t)=κ2Δu(t),u(t)B,

    and B is the infinite-dimensional Hilbert space

    B={uH2(Ω):u(x,y)=u(x+L1,y),u(x,y)=u(x,y+L2)}. (2.2)

    For any ϕ(x,y),ψ(x,y)B, we introduce the inner product and the norms

    (ϕ(x,y),ψ(x,y))=Ωϕ(x,y)ψ(x,y)dxdy,ϕ=(ϕ(x,y),ϕ(x,y)),|ϕ|1=(Δϕ(x,y),ϕ(x,y)).

    Then, by taking the inner product of the abstract ODE's system (2.1) with u(t), we are able to find that system (2.1) can preserve the separable energy

    H[u(t),u(t)]:=H1[u(t)]+H2[u(t)]H1[u(t0)]+H2[u(t0)]=H[u(t0),u(t0)], (2.3)

    where the kinetic energy part H1[u(t)] and the potential energy part H2[u(t)] are

    H1[u(t)]=12u(t)2andH2[u(t)]=κ22|u(t)|21+(F(u(t)),1), (2.4)

    respectively. Obviously, the energy H[u(t),u(t)] of the abstract ODE (2.1) is the same as the energy E(t) of the two-dimensional nonlinear wave Eqs (1.1)–(1.3). Moreover, the abstract ODE's system (2.1) is actually a Hamiltonian system

    ddt[u(t)v(t)]=S[δH2[u(t)]δuδH1[v(t)]δv], (2.5)

    where v(t)=u(t) and S is a skew-adjoint operator

    S=[0110].

    In light of the definition of the variational derivatives, we are able to check that

    δH1[v(t)]δv=v(t)andδH2[u(t)]δu=Au(t)f(u(t))=g(u(t)). (2.6)

    The main purpose of this work is to design a suitable time-stepping scheme for the two-dimensional nonlinear wave Eqs (1.1)–(1.3), which could exactly preserve the separable energy H[u(t),u(t)] or E(t). To achieve this purpose, the temporal discretization strategy will be first considered for the abstract ODE (2.1) in the infinite-dimensional function space.

    For any positive integer N, we define the temporal mesh grid as

    ΩN:={tn|tn=t0+nΔt,n=0,1,,N} (2.7)

    with time step size Δt=(Tt0)/N, and introduce the following approximations:

    unu(tn),vnu(tn),Unτu(tn+τΔt),τ[0,1].

    Then, applying the energy-preserving integrators, which are proposed for the second-order Hamiltonian ordinary differential systems (see [22,26]), to the abstract ODEs (2.1), we can establish the time-stepping scheme for the two-dimensional nonlinear wave Eqs (1.1)–(1.3).

    Definition 2.1. For any one temporal single step tn to tn+1, a continuous-stage Runge-Kutta-Nyström (RKN) time-stepping scheme for the abstract ODE (2.1) is defined as

    {Unτ=un+τΔtvn+Δt210P3,2(τ,σ)g(Unσ)dσ,un+1=un+Δtvn+Δt210(1τ)g(Unτ)dτ,vn+1=vn+Δt10g(Unτ)dτ, (2.8)

    where the weight function P3,2(τ,σ) is a cubic binary polynomial of the form

    P3,2(τ,σ)=τ2(16σ+6σ2+3τ6σ2τ2τ2+4στ2),(τ,σ)[0,1]×[0,1].

    Remark 2.1. In [22,26], the authors introduced a framework for energy-preserving continuous-stage RKN methods for solving second-order Hamiltonian ODEs. Drawing upon the methodologies proposed in [22,26], we develop a novel energy-preserving time-stepping scheme utilizing the weight function P3,2(τ,σ), and extend this scheme to the two-dimensional nonlinear wave Eqs (1.1)(1.3). Furthermore, it is important to emphasize that the selection of the weight function P3,2(τ,σ) is not unique. Different choices of weight functions can result in numerical methods exhibiting varying accuracy.

    Now, we focus on verifying the energy conservation of the continuous-stage RKN time-stepping scheme defined in Definition 2.1 for the two-dimensional nonlinear wave Eqs (1.1)–(1.3).

    Theorem 2.1. The continuous-stage Runge-Kutta-Nyström time-stepping scheme (2.8) can exactly preserve the energy H[u(t),v(t)] of the two-dimensional nonlinear wave Eqs (1.1)(1.3) or the infinite-dimensional abstract ODE's system (2.1), that is,

    H[un+1,vn+1]H[un,vn],n=0,1,2,,N1. (2.9)

    Proof. Noticing the form of the separable energy (2.3) and (2.4), we have

    H[un+1,vn+1]=H1[vn+1]+H2[un+1]=12(vn+1,vn+1)+H2[un+1]. (2.10)

    It follows from inserting the expression of vn+1 into (2.10) and after a careful calculation that

    H[un+1,vn+1]=12(vn+Δt10g(Unτ)dτ,vn+Δt10g(Unτ)dτ)+H2[un+1]=H[un,vn]+Δt(vn,10g(Unτ)dτ)+Δt22(10g(Unτ)dτ,10g(Unτ)dτ)+H2[un+1]H2[un]. (2.11)

    Moreover, it is evident that

    H2[un+1]H2[un]=10dH2[Unτ]=10(δH2[Unτ]δu,dUnτdτ)dτ=10(g(Unτ),dUnτdτ)dτ. (2.12)

    Substituting the expressions of Unτ into (2.12) leads to

    H2[un+1]H2[un]=10(g(Unτ),Δtvn+Δt210P3,2(τ,σ)τg(Unσ)dσ)dτ=Δt(vn,10g(Unτ)dτ)Δt210(g(Unτ),10P3,2(τ,σ)τg(Unσ)dσ)dτ. (2.13)

    Noticing the form of the weight function P3,2(τ,σ), we have

    P3,2(τ,σ)τ=123σ+3τ+3σ23τ26σ2τ+6στ2,(τ,σ)[0,1]×[0,1].

    Therefore, Eq (2.13) can be simplified as

    H2[un+1]H2[un]=Δt(vn,10g(Unτ)dτ)Δt22(10g(Unτ)dτ,10g(Unσ)dσ). (2.14)

    Comparing (2.11) with (2.14), we obtain

    H[un+1,un+1]H[un,un],n=0,1,2,,N1.

    The conclusion of the theorem is confirmed.

    The symmetric time integration method usually exhibits superior long time computational behavior along the numerical flows (see Chapter V in [16]). We know that the two-dimensional nonlinear wave Eqs (1.1)–(1.3) are temporal reversible (see, e.g., [24,28,29]). Therefore, it will be significant to investigate the symmetry of the energy-preserving continuous stage RKN time-stepping scheme.

    Theorem 2.2. The energy-preserving continuous-stage Runge-Kutta-Nyström time-stepping scheme (2.1) for solving the two-dimensional nonlinear wave Eqs (1.1)(1.3) is temporal reversible.

    Proof. According to the concept of the time reversible integration method (see Chapter V in [16]), and applying the following transformations

    ΔtΔt,un+1un,vn+1vn,τ=1τ

    to the time-stepping scheme (2.1), we obtain the adjoint scheme

    {Un1τ=un+1(1τ)Δtvn+1+Δt210P3,2(1τ,σ)g(Unσ)dσ,un=un+1Δtvn+1+Δt210σg(Unσ)dσ,vn=vn+1Δt10g(Unσ)dσ. (2.15)

    After a careful calculation, the last two equations of the adjoint scheme (2.15) can be rewritten as

    {un+1=un+Δtvn+Δt210σg(Unσ)dσ,vn+1=vn+Δt10g(Unσ)dσ. (2.16)

    Inserting Eq (2.16) into the first equation of (2.15), we obtain

    Un1τ=un+τΔtvn+Δt210(σ(1τ)+P3,2(1τ,σ))g(Unσ)dσ. (2.17)

    The integral transformation τ=1ξ yields that

    1ξ(1τ)+P3,2(1τ,1ξ)=P3,2(τ,ξ). (2.18)

    Therefore, we see that the adjoint scheme (2.15) is the same as scheme (2.8). That means the energy-preserving continuous stage RKN time-stepping scheme is symmetric or temporal reversible.

    Utilizing the variation-of-constants formula on the infinite-dimensional abstract ODE system (2.1), the exact solution of the abstract system (2.1) can be expressed as

    u(tn+τΔt)=u(tn)+τΔtv(tn)+Δt2τ0(τσ)g(u(tn+σΔt))dσ,τ[0,1]. (2.19)

    Furthermore, it is easy to obtain from Eq (2.19) that

    {u(tn+1)=u(tn)+Δtv(tn)+Δt210(1τ)g(u(tn+τΔt))dτ,v(tn+1)=v(tn)+Δt10g(u(tn+τΔt))dτ. (2.20)

    Inserting the exact solution u(t) of the infinite-dimensional abstract ODE system (2.1) into the time-stepping scheme (2.8), we have

    {u(tn+τΔt)=u(tn)+τΔtv(tn)+Δt210P3,2(τ,σ)g(u(tn+σΔt))dσ+Rn(τ),u(tn+1)=u(tn)+Δtv(tn)+Δt210(1τ)g(u(tn+τΔt))dτ,v(tn+1)=v(tn)+Δt10g(u(tn+τΔt))dτ, (2.21)

    where the residual Rn(τ) is a function of τ[0,1]. Applying the Taylor expansion with integral remainder

    u(tn+σΔt)=u(tn)+σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz (2.22)

    to the nonlinear integrands appearing in (2.19) and the first equation of (2.21) results in

    u(tn+τΔt)=u(tn)+τΔtv(tn)+Δt2τ0(τσ)g(u(tn)+σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz)dσ, (2.23)

    and

    u(tn+τΔt)=u(tn)+τΔtv(tn)+Δt210P3,2(τ,σ)g(u(tn)+σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz)dσ+Rn(τ). (2.24)

    Comparing the formulae (2.23) with (2.24), and noticing g(u)=Au+f(u), we can approximate the local residuals Rn(τ),0τ1, in the following theorem.

    Theorem 2.3. Suppose that the exact solution u of the abstract ODE's system (2.1) satisfies uC2([t0,T],B) and the nonlinear function fL([t0,T],B). Then, the remainder Rn(τ) satisfies the estimations

    Rn(τ)CΔt4,0τ1, (2.25)

    where C is a constant and independent of Δt.

    Proof. Subtracting (2.24) from (2.23) and noticing g(u)=Au+f(u), we obtain

    Rn(τ)=Θn1(τ)+Θn2(τ), (2.26)

    where

    Θn1(τ)=Δt2τ0(τσ)A(u(tn)+σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz)dσ+Δt210P3,2(τ,σ)A(u(tn)+σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz)dσ, (2.27)

    and

    Θn2(τ)=Δt2τ0(τσ)f(u(tn)+σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz)dσΔt210P3,2(τ,σ)f(u(tn)+σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz)dσ. (2.28)

    It follows from the definition of the bilinear polynomial weight function P3,2(τ,σ) that

    τ0(τσ)dσ=10P3,2(τ,σ)dσ=τ22andτ0σ(τσ)dσ=10σP3,2(τ,σ)dσ=τ36. (2.29)

    Therefore, it is easy to check that

    Θn1(τ)=Δt4τ0σ0(τσ)(σz)Au(tn+zΔt)dzdσ+Δt410σ0P3,2(τ,σ)(σz)Au(tn+zΔt)dzdσ. (2.30)

    Utilizing the Taylor expansion of f(), i.e.,

    f(u(tn)+σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz)=f(u(tn))+f(u(tn))(σΔtv(tn)+Δt2σ0(σz)u(tn+zΔt)dz)+

    in (2.28) and recalling Eq (2.29), we have

    Θn2(τ)=Δt4f(u(tn))τ0σ0(τσ)(σz)u(tn+zΔt)dzdσΔt4f(u(tn))10σ0P3,2(τ,σ)(σz)u(tn+zΔt)dzdσ+O(Δt5). (2.31)

    Inserting the results (2.30) and (2.31) into (2.26), and taking the L2 norms on both sides of (2.29), it is easy to verify the estimated result of the theorem.

    In this section, by combing the Fourier pseudo-spectral spatial approximation with the continuous-stage RKN time-stepping scheme (2.8), we will construct the energy-preserving fully discrete scheme for the two-dimensional nonlinear wave Eqs (1.1) and (1.2).

    Choose M1 and M2 to be any even integers, and define Δx:=L1M1 and Δy:=L2M2 as the spatial steps. Then, the temporal-spatial grid points are denoted as ΩNM=ΩM×ΩN, where the temporal grid ΩN is given by (2.7), and the spatial grid ΩM is defined as

    ΩM:={(xj,xk)|xj=jΔx,j=0,1,,M11,yk=kΔy,k=0,1,,M21}. (3.1)

    The grid function space VM defined on ΩM is given by

    VM={u|u=(uj,k)withuj,k=u(xj,yk),(xj,yk)ΩM}.

    For any u=(uj,k)VM, we can reformulate it as the vector form

    u=(u0,0,,uM11,0,u0,1,,uM11,1,,u0,M21,,uM11,M21).

    Therefore, the vector space of the grid functions VM, which is identical to VM, can be presented as

    VM={u|u=(u0,0,,uM11,0,u0,1,,uM11,1,,u0,M21,,uM11,M21)withu=(uj,k)VM}.

    In addition, the corresponding discrete inner product and norm are defined as

    u,v=ΔxΔyM11j=0M21k=0unj,kvnj,k,u=u,u,u,vVM.

    Define the interpolation space as

    SpM:=span{gj(x)gk(y),0jM11,0kM21},

    where gj(x) and gk(y) are trigonometric polynomials

    gj(x)=1M1M1/2k1=M1/21ck1eik1μ1(xxj),gk(y)=1M2M2/2k2=M2/21ck2eik2μ2(yyk)

    with μ1=2πL1,μ2=2πL2, and

    ck1={1,|k1|<M1/2,2,|k1|=M1/2,ck2={1,|k2|<M2/2,2,|k2|=M2/2.

    Therefore, for any periodic function u(x,y)L2p(Ω), the interpolation operator IM:L2p(Ω)SpM is

    IMu(x,y)=M11j=0M21k=0u(xj,yk)gj(x)gk(x)=M1/2k1=M1/2M2/2k2=M2/2˜uk1,k2eik1μ1xeik2μ2y, (3.2)

    where the Fourier coefficients ˜uk1,k2 are

    ˜uk1,k2=1M1ck11M2ck2M11l=0M21k=0u(xl,yk)eik1μ1xleik2μ2yk. (3.3)

    Moreover, it is simple to check that

    ˜uM1/2,=˜uM1/2,and˜u,M2/2=˜u,M2/2.

    It follows from applying the differential operator A to the interpolation IMu(x,y) that

    AIMu(x,y)|x=xi,y=yj=M1/2k1=M1/2M2/2k2=M2/2κ2[(μ1k1)2+(μ2k2)2]˜uk1,k2eik1μ1xieik2μ2yj=M1/2k2=M1/2(M2/2k1=M2/2κ2(μ1k1)2˜uk1,k2eik1μ1xi)eik2μ2yj+M2/2k1=M2/2(M1/2k2=M1/2κ2(μ2k2)2˜uk1,k2eik2μ2yj)eik1μ1xi=((IM2Dx2+Dy2IM1)u)i,j, (3.4)

    where IMi,i=1,2 are the unity matrices, and Dx2=FHM1Λ1FM1 and Dy2=FHM2Λ2FM2 are the spectral differential matrices. Here, we should point out that FM is the discrete Fourier transform matrix with elements (FM)j,k=1Me2πi(j)(k)M,j,k=0,1,,M1, FHM is the conjugate transformation matrix of FM, and Λ1,Λ2 are the frequency matrices with entries

    Λ1=diag(λDx2,0,λDx2,1,,λDx2,M11),λDx2,j={κ2(μ1j)2,0jM1/2,κ2(μ1(jM1))2,M1/2<jM11,Λ2=diag(λDy2,0,λDy2,1,,λDy2,M21),λDy2,j={κ2(μ2j)2,0jM2/2,κ2(μ2(jM2))2,M2/2<jM21. (3.5)

    Thus, the spectral differential matrix A for approximating the 2D differential operator A can be expressed as (see, e.g., [13,14,40])

    Au=(IM2Dx2+Dy2IM1)u=(FHM2IM2FM2FHM1Λ1FM1+FHM2Λ2FM2FHM1IM1FM1)u=((FM2FM1)H(IM2Λ1)(FM2FM1)+(FM2FM1)H(Λ2IM1)(FM2FM1))u=((FM2FM1)H(IM2Λ1+Λ2IM1)(FM2FM1))u. (3.6)

    Actually, the spectral differential matrix A is a symmetric and semi-positive matrix, and Au can be fast computed by the two-dimensional FFT function ifft2(Λ.fft2(u)) built in MATLAB.

    Using the spectral differential matrix A to approximate the differential operator A, the two-dimensional nonlinear wave Eqs (1.1) and (1.2) or the abstract ODE system (2.1) can be converted into the semi-discrete system

    {u(t)=g(u(t)),t[t0,T],u(t0)=φ0,u(t0)=φ1, (3.7)

    where u(t)VM and g(u(t))=Au(t)+f(u(t)) with

    f(u(t))=(f(u0,0(t)),,f(uM11,0(t)),f(u0,1(t)),,f(uM11,1(t)),,f(u0,M21(t)),,f(uM11,M21(t))).

    Taking the discrete inner product on both sides of the semi-discrete system (3.7) with u(t), we have

    u(t),u(t)+Au(t),u(t)f(u(t)),u(t)=0,

    which means

    ddt(12u(t),u(t)+12Au(t),u(t)+F(u(t)),1)=0.

    Therefore, we can conclude that the semi-discrete system (3.7) is energy-conserving.

    Theorem 3.1. Suppose that u(t)VM is the solution of the semi-discrete system (3.7). Then, the semi-discrete system (3.7) can conserve the discrete energy

    H[u(t),u(t)]:=H1[u(t)]+H2[u(t)], (3.8)

    where the discrete kinetic energy H1[u(t)] and discrete potential energy H2[u(t)] are given by

    H1[u(t)]=12u(t),u(t)andH2[u(t)]=12Au(t),u(t)+F(u(t)),1. (3.9)

    Proof. The proof process is straightforward along the above analysis. Therefore, we omit the details.

    Remark 3.1. Actually, the energies (3.8) and (3.9) the discrete versions of the energies (2.3) and (2.4) of the the two-dimensional nonlinear wave Eqs (1.1) and (1.2) or the abstract system (2.1). Therefore, we will explore the fully discrete scheme which can exactly preserve the discrete energies (3.8) and (3.9) in this work.

    The main strategy of the construction of the energy-preserving fully discrete scheme is to approximate the differential operator A in the continuous-stage RKN time-stepping scheme (2.8) by the spectral differential matrix A. Therefore, the following theorem will show that the fully discrete scheme could preserve the discrete energy (3.8) and (3.9) exactly.

    Theorem 3.2. By following the continuous stage RKN Fourier pseudo-spectral scheme

    {Unτ=un+τΔtvn+Δt210P3,2(τ,σ)g(Unσ)dσ,un+1=un+Δtvn+Δt210(1τ)g(Unτ)dτ,vn+1=vn+Δt10g(Unτ)dτ, (3.10)

    with g(u)=Au+f(u) and Δt the time step size, the discrete energies (3.8) and (3.9) are conserved, i.e.,

    H[un+1,vn+1]=H[un,vn]. (3.11)

    Proof. We calculate the separable energy

    H[un+1,vn+1]=H1[vn+1]+H2[un+1]. (3.12)

    Inserting vn+1 into H1[vn+1] and keeping that A is a symmetric matrix in mind gives

    H1[vn+1]=12vn+1,vn+1=12vn+Δt10g(Unτ)dτ,vn+Δt10g(Unτ)dτ=12vn,vn+Δtvn,10g(Unτ)dτ+Δt2210g(Unτ)dτ,10g(Unτ)dτ=H1[vn]+Δtvn,10g(Unτ)dτ+Δt2210g(Unτ)dτ,10g(Unτ)dτ. (3.13)

    On the other hand, we have

    H2[un+1]H2[un]=10dH2[Unτ]=10g(Unτ),dUnτdτdτ=10g(Unτ),Δtvn+Δt210P3,2(τ,σ)τg(Unσ)dσdτ. (3.14)

    It follows from inserting

    P3,2(τ,σ)τ=123σ+3τ+3σ23τ26σ2τ+6στ2

    into (3.14) that

    H2[un+1]H2[un]=Δt10g(Unτ)dτ,vnΔt2210g(Unτ)dτ,10g(Unσ)dσ. (3.15)

    Combining the results of (3.12), (3.13), and (3.15), we have the desired result.

    Remark 3.2. Similar as the proof process of Theorem 2.2, it can be verified that the continuous stage RKN Fourier pseudo-spectral scheme (3.10) is time reversible. Moreover, we have noticed that the authors in [40] considered the integrating factor technique and the 4th-order (2-stage) Gauss-Legendre Runge-Kutta scheme to propose a symplectic time integration method for three-dimensional nonlinear water waves. This method can sufficiently use the oscillation generated by the spatial discretisation. Moreover, we have observed that the authors in [40] employed the integrating factor technique and the fourth-order (two-stage) Gauss-Legendre Runge-Kutta scheme to develop a symplectic time integration method for three-dimensional nonlinear water waves. This approach efficiently utilizes the oscillations resulting from spatial discretization, and typically yields accurate results at a reasonable computational cost. Perhaps, the combination of our proposed energy-preserving time integrator with the integrating factor technique could lead to a more efficient energy-preserving scheme. This will be considered in our future research.

    Remark 3.3. In practice, the integrals in the fully discrete scheme usually cannot be easily calculated. Therefore, the s-point Gauss-Legendre formula (bi,ci)si=1 will be used to evaluate the integrals

    {Unci=un+ciΔtvn+Δt2sj=1bjP3,2(ci,cj)g(Uncj),i=1,2,,s,un+1=un+Δtvn+Δt2si=1bi(1ci)g(Unci),vn+1=vn+Δtsi=1big(Unci). (3.16)

    Since the s-point GL quadrature formula is symmetric, the formula (3.16) is also symmetric.

    To date, we have developed an energy-preserving fully discrete scheme for solving the two-dimensional nonlinear wave Eqs (1.1) and (1.2). This was achieved by initially semidiscretizing the temporal derivatives using a continuous-stage RKN method, followed by applying the Fourier spectral differential matrix A to approximate the spatial differential operator A. It has been observed that the Fourier pseudo-spectral method for approximating spatial derivatives can achieve spectral precision of order O(Mr) provided that the spatial regularity conditions are adequately satisfied.

    Assume that u(t) and v(t) represent the exact solution and its derivative of the abstract ODE's system (2.1), while u(t) and v(t) denote the exact solution and its derivative of the semi-discrete system (3.7). Additionally, un and vn signify the numerical solutions obtained from the continuous stage RKN Fourier pseudo-spectral scheme (3.10). It follows from inserting the exact solution u(t) of the semi-discrete system (3.7) into the continuous stage RKN Fourier pseudo-spectral scheme (3.10) that

    {u(tn+τΔt)=u(tn)+τΔtv(tn)+Δt210P3,2(τ,σ)g(u(tn+σΔt))dσ+Rn(τ),u(tn+1)=u(tn)+Δtv(tn)+Δt210(1τ)g(u(tn+τΔt))dτ,v(tn+1)=v(tn)+Δt10g(u(tn+τΔt))dτ, (3.17)

    where Rn(τ)VM is the temporal local truncation error. Similar to the analysis of Theorem 2.3, we obtain the estimation for the residual Rn(τ) in the following theorem.

    Theorem 3.3. Suppose that the semi-discrete system (3.7) is well-posed and satisfies u(t)C2([t0,T]) and the nonlinear function fL([t0,T]). Then, the local truncation error Rn(τ) could be estimated as

    Rn(τ)˜CΔt4,0τ1, (3.18)

    where ˜C is a constant and independent of Δt.

    Proof. The details of the proof are similar to the process of Theorem 2.3, so we omit the details.

    Letting

    \begin{equation} e^n = u(t_n)-\mathbf{u}^n, \quad \eta^n = v(t_n)-\mathbf{v}^n, \quad \mathbf{e}^n = \mathbf{u}(t_n)-\mathbf{u}^n, \quad \mathit{\boldsymbol{\eta}}^n = \mathbf{v}(t_n)-\mathbf{v}^n, \quad \mathbf{E}_{\tau}^n = \mathbf{u}(t_n+\tau\Delta t)-\mathbf{U}_{\tau}^n, \end{equation} (3.19)

    and subtracting (3.10) from (3.17), we obtain

    \begin{equation} \left\{\begin{aligned} &\mathbf{E}_\tau^n = \mathbf{e}^n+\tau \Delta t \mathit{\boldsymbol{\eta}}^n+\Delta t^2\int_0^1P_{3, 2}(\tau, \sigma)\Big(g\big(\mathbf{u}(t_n+\sigma\Delta t\big)-g\big(\mathbf{U}_\sigma^n)\big)\Big){\rm d}\sigma +\mathbf{R}^n(\tau), \\ &\mathbf{e}^{n+1} = \mathbf{e}^n+\Delta t \mathit{\boldsymbol{\eta}}^n+\Delta t^2\int_0^1(1-\tau)\Big(g\big(\mathbf{u}(t_n+\tau\Delta t)\big)-g\big(\mathbf{U}_{\tau}^n\big)\Big){\rm d}\tau, \\ &\mathit{\boldsymbol{\eta}}^{n+1} = \mathit{\boldsymbol{\eta}}^n+\Delta t\int_0^1\Big(g\big(\mathbf{u}(t_n+\tau\Delta t)\big)-g\big(\mathbf{U}_{\tau}^n\big)\Big){\rm d}\tau. \end{aligned}\right. \end{equation} (3.20)

    We suppose the two-dimensional nonlinear wave Eqs (1.1) and (1.2) are well-posed. Subsequently, we present the error estimation for the fully discrete scheme (3.10) as detailed in the following theorem.

    Theorem 3.4. If the exact solution u(x, y, t) of the two-dimensional nonlinear wave Eqs (1.1) and (1.2) satisfies u(x, y, t)\in C^2\big([t_0, T], \mathcal{B}\big) , and the nonlinear function f(\cdot) satisfies f'\in L^{\infty}\big([t_0, T], \mathcal{B}\big) , then under the limitation of 0 < \Delta t\leq h_0 with a sufficiently small h_0 such that h_0\|A\| < 1 , we obtain the error bounds

    \begin{equation} \|e^n\|+\Delta t\|\eta^n\|\lesssim M^{-r}+\Delta t^4. \end{equation} (3.21)

    Here, we should point out that A\lesssim B means there is a constant C such that A\leq CB , and M = M_1 = M_2 is the spatial grid scale. Moreover, the constant C depends on T , but is independent of M, \|A\| , and \Delta t .

    Proof. The concept of the temporal-spatial error splitting method suggests that

    \begin{equation} \begin{aligned} \|e^n\|+\Delta t\|\eta^n\|\leq\; & \|u(t_n)-\mathbf{u}(t_n)\|+\Delta t\|v(t_n)-\mathbf{v}(t_n)\| +\|\mathbf{u}(t_n)-\mathbf{u}^n\|+\Delta t\|\mathbf{v}(t_n)-\mathbf{v}^n\| \\ \leq\; &\mathcal{O}(M^{-r})+\|\mathbf{e}^n\|+\Delta t\|\mathit{\boldsymbol{\eta}}^n\|. \end{aligned} \end{equation} (3.22)

    Therefore, to obtain the accuracy of the fully discrete scheme, it is essential to concentrate on the analysis of temporal accuracy. Taking norms on both sides of (3.20) leads to

    \begin{equation} \left\{\begin{aligned} &\|\mathbf{E}_\tau^n\|\lesssim\|\mathbf{e}^n\|+\tau \Delta t \|\mathit{\boldsymbol{\eta}}^n\|+\Delta t^2\|A\|\int_0^1\|\mathbf{E}_\sigma^n\|{\rm d}\sigma +\mathcal{O}(\Delta t^4), \\ &\|\mathbf{e}^{n+1}\|\lesssim\|\mathbf{e}^n\|+\Delta t \|\mathit{\boldsymbol{\eta}}^n\|+\Delta t^2\|A\|\int_0^1(1-\tau)\|\mathbf{E}_\tau^n\|{\rm d}\tau, \\ &\|\mathit{\boldsymbol{\eta}}^{n+1}\|\lesssim\|\mathit{\boldsymbol{\eta}}^n\|+\Delta t\|A\|\int_0^1\|\mathbf{E}_\tau^n\|{\rm d}\tau. \end{aligned}\right. \end{equation} (3.23)

    Then, under the restriction of the time step size \Delta t\leq h_0 with sufficiently small h_0 satisfying h_0\|A\| < 1 , the first inequality of (3.23) results in

    \begin{equation} \int_0^1\|\mathbf{E}_\tau^n\|{\rm d}\tau\lesssim\|\mathbf{e}^n\|+\Delta t \|\mathit{\boldsymbol{\eta}}^n\|+\mathcal{O}(\Delta t^4). \end{equation} (3.24)

    Summing up the last two inequalities of (3.23), we have

    \begin{equation} \|\mathbf{e}^{n+1}\|+\Delta t\|\mathit{\boldsymbol{\eta}}^{n+1}\|\lesssim\; \|\mathbf{e}^{n}\|+2\Delta t\|\mathit{\boldsymbol{\eta}}^{n}\|+\Delta t^2\|A\|\int_0^1\|\mathbf{E}_\tau^n\|{\rm d}\tau. \end{equation} (3.25)

    Moreover, the third inequality of (3.23) results in

    \begin{equation} \|\mathit{\boldsymbol{\eta}}^n\|\lesssim\; \Delta t\|A\|\sum\limits_{k = 0}^{n-1}\int_0^1\|\mathbf{E}_\tau^k\|{\rm d}\tau. \end{equation} (3.26)

    Combining (3.25) and (3.26), we obtain

    \begin{equation} \|\mathbf{e}^{n+1}\|+\Delta t\|\mathit{\boldsymbol{\eta}}^{n+1}\|\lesssim\; \|\mathbf{e}^{n}\|+\Delta t\|\mathit{\boldsymbol{\eta}}^{n}\|+\Delta t^2\|A\|\sum\limits_{k = 0}^{n}\int_0^1\|\mathbf{E}_\tau^k\|{\rm d}\tau. \end{equation} (3.27)

    The mathematical induction will be an efficient approach to prove the result of the theorem.

    Step Ⅰ. Letting n = 0 in (3.24) and noticing \mathbf{e}^0 = 0, \mathit{\boldsymbol{\eta}}^0 = 0 , we have

    \begin{equation*} \int_0^1\|\mathbf{E}_\tau^0\|{\rm d}\tau = \mathcal{O}(\Delta t^4). \end{equation*}

    Furthermore, noticing the limitation of the time step size again, the inequality (3.25) leads to

    \begin{equation*} \|\mathbf{e}^{1}\|+\Delta t\|\mathit{\boldsymbol{\eta}}^{1}\|\lesssim\Delta t^2\|A\|\int_0^1\|\mathbf{E}_{\tau}^0\|{\rm d}\tau = \mathcal{O}(\Delta t^4). \end{equation*}

    Step Ⅱ. Now, we assume that the estimation (3.21) is valid for 1\leq n\leq m-1 . That is,

    \begin{equation*} \|\mathbf{e}^{n}\|+\Delta t\|\mathit{\boldsymbol{\eta}}^{n}\| = \mathcal{O}(\Delta t^4), \qquad n = 1, 2, \dots, m-1. \end{equation*}

    Then, by mathematical induction, we only need to verify that the estimation (3.21) is still valid for n = m . Setting n = m-1 in (3.27) and using the above assumptions leads to

    \begin{equation*} \begin{aligned} \|\mathbf{e}^{m}\|+\Delta t\|\mathit{\boldsymbol{\eta}}^{m}\|\lesssim\; & \|\mathbf{e}^{m-1}\|+\Delta t\|\mathit{\boldsymbol{\eta}}^{m-1}\|+\Delta t^2\|A\|\sum\limits_{k = 0}^{m-1}\int_0^1\|\mathbf{E}_\tau^k\|{\rm d}\tau \\ = \; & \mathcal{O}(\Delta t^4)+\Delta t^2\|A\|\sum\limits_{k = 0}^{m-1}\mathcal{O}(\Delta t^4) \lesssim\Delta t^4. \end{aligned} \end{equation*}

    Therefore, the proof of (3.21) is completed.

    Remark 3.4. The conclusion of Theorem 3.3 indicates that the continuous-stage RKN Fourier pseudo-spectral scheme has at least fourth-order accuracy in the temporal domain. Owing to the temporal reversibility of the scheme (3.10), the forthcoming numerical experiments demonstrate that the proposed energy-preserving continuous-stage RKN Fourier pseudo-spectral scheme (3.10) can achieve sixth-order convergence in time.

    Remark 3.5. In general, the preservation of energy typically ensures the stability of the fully discrete scheme. The analysis process of Theorem 3.2 demonstrates that the proposed energy-preserving scheme (3.10) is unconditionally stable. However, according to the result presented in Theorem 3.4, it can be concluded that the continuous stage RKN Fourier pseudo-spectral scheme (3.10) exhibits convergence under the condition 0 < \Delta t \leq h_0 with h_0\|A\| < 1 . In fact, the constraint h_0\|A\| < 1 corresponds to the CFL condition, as the differential matrix is intrinsically linked to the spatial discretization scale.

    We observe that the weight function P_{3, 2}(\tau, \sigma) of the the continuous stage RKN Fourier pseudo-spectral scheme defined in Definition 2.1 is a cubic binary polynomial. Hereafter, the continuous stage RKN Fourier pseudo-spectral scheme will be denoted as EP3-FP. In this section, we will calculate the two-dimensional Klein-Gordon equation and the two-dimensional sine-Gordon equation to verify the precision, the efficiency, and the energy preservation of the derived EP3-FP scheme. Additionally, the following energy-preserving time integrators are chosen for comparison:

    AVF: the energy-preserving second-order averaged vector field method (see, e.g., [15,16]);

    HEP3: the symmetric sixth-order energy-preserving integrator constructed by Hairer in [15];

    SRKN3: the continuous-stage symplectic RKN method of order six (see, e.g., [34,35,36]).

    The fully discrete scheme is obtained after discretizing the spatial derivatives with the Fourier pseudo-spectral method. We compute the temporal convergence rate by the following formula:

    \begin{equation} {\rm Rate} = \log_2\frac{{\rm GE}(h)}{{\rm GE}(h/2)}\qquad {\rm with}\qquad {\rm GE}(h) = \|U(T;h)-u(T;h)\|, \end{equation} (4.1)

    where the global error {\rm GE}(h) is the difference of the exact solution U(T; h) with the numerical solution u(T; h) at time T with step h . Moreover, it is known that the exact solution of the two-dimensional sine-Gordon equation could not be explicitly represented. Therefore, we will use the posterior error (see, e.g., [9,20]) to calculate the convergence rate, i.e.,

    \begin{equation} {\rm Rate} = \log_2\frac{{\rm PE}(h)}{{\rm PE}(h/2)}\qquad {\rm with}\qquad {\rm PE}(h) = \|u(T;h)-u(T;h/2)\|. \end{equation} (4.2)

    Furthermore, it is important to emphasize that the energy-preserving continuous-stage RKN time-stepping method introduced in this work, as well as the numerical methods selected for comparison, are closely associated with nonlinear integrals. To approximate these nonlinear integrals, the four-point Gauss-Legendre quadrature formula will be employed in numerical simulation.

    Problem 1. Consider the two-dimensional nonlinear periodic Klein-Gordon equation (see, e.g., [18])

    \begin{equation} \left\{ \begin{aligned} &\frac{\partial^2u}{\partial t^2}-c^2\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right)+au+bu^3 = f(x, y, t), \quad (x, y)\in (-2, 2)\times(-2, 2), \quad t\in(0,100], \\ &u(x, y, 0) = \cos\big(\pi(x+y)\big), \quad \frac{\partial u(x, y, 0)}{\partial t} = \pi\sin\big(\pi(x+y)\big), \quad (x, y)\in [-2, 2]\times[-2, 2], \end{aligned} \right. \end{equation} (4.3)

    with the right-hand function f(x, y, t) = \cos\big(\pi(x+y-t)\big)+\cos^3\big(\pi(x+y-t)\big) . The exact solution is given by

    u(x, y, t) = \cos\big(\pi(x+y-t)\big).

    We choose the parameters as a = 1, b = 1, and c^2 = \frac{1}{2} . In Table 1, we list the errors and the corresponding convergence rates of the proposed EP3-FP scheme by varying the spatial and temporal step sizes. In Figure 1, we set the spatial grid scales as M = M_x = M_y = 64 , and the logarithms of the global errors \log_{10}(GE) against different time steps and the CPU times are plotted in Figure 1(a) and Figure 1(b), respectively. The logarithms of the energy errors of the EP3-FP scheme are plotted in Figure 1(b), which show that the proposed scheme is energy-preserving. The numerical results in Table 1 and Figure 1(a) illustrate that the proposed EP3-FP scheme achieves sixth-order temporal accuracy. Figure 1(c) shows that the EP3 time-stepping scheme has better computational efficiency than the chosen numerical methods.

    Table 1.  The global errors and temporal convergence rates of the "EP3-FP" scheme for solving Problem 1.
    Error \Delta t=0.08 \Delta t=0.04 \Delta t=0.02 \Delta t=0.01
    M=8 4.7423{\rm E}-08 7.4247{\rm E}-10 1.1636{\rm E}-11 2.0639{\rm E}-13
    \textbf{Rate} \textbf{5.9971} \textbf{5.9957} \textbf{5.8170}
    M=16 5.5603{\rm E}-08 8.7012{\rm E}-10 1.3624{\rm E}-11 2.4547{\rm E}-13
    \textbf{Rate} \textbf{5.9978} \textbf{5.9970} \textbf{5.7944}
    M=32 5.9563{\rm E}-08 9.3049{\rm E}-10 1.4581{\rm E}-11 2.4236{\rm E}-13
    \textbf{Rate} \textbf{6.0003} \textbf{5.9958} \textbf{5.9108}
    M=64 5.9670{\rm E}-08 9.3484{\rm E}-10 1.4611{\rm E}-11 2.4836{\rm E}-13
    \textbf{Rate} \textbf{5.9962} \textbf{5.9996} \textbf{5.8785}

     | Show Table
    DownLoad: CSV
    Figure 1.  Results for Problem 1 with M_x = M_y = 64 . (a) the log-log plot of the global error against different steps h . (b) the log-log plot of the global error against CPU time. (c) the log plot of the relative energy error against integrate time with H(0) = 123.9352528130722 .

    Problem 2. Consider the two-dimensional sine-Gordon equations (see, e.g., [26,27])

    \begin{equation} \frac{\partial^2 u}{\partial t^2}-\kappa^2\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)+\sin\big(u(x, y, t)\big) = 0, \quad (x, y)\in[-1, 1]\times[-1, 1], \; \; \; \; t\in(0,100], \end{equation} (4.4)

    with the dimensionless parameter \kappa = 1/20 , and the initial conditions

    u(x, y, 0) = 4\arctan\Big(\exp\big(3-\sqrt{x^2+y^2}\big/\kappa^2\big)\Big), \qquad \frac{\partial u(x, y, 0)}{\partial t} = 0.

    Suppose that the two-dimensional sine-Gordon Eq (4.4) is equipped with periodic boundary conditions. Some snapshots of the numerical solution by the EP3-FP scheme are shown in Figure 2. These results demonstrate that the proposed EP3-EP scheme can efficiently simulate the two-dimensional sine-Gordon Eq (4.4) in a relatively long time domain. Moreover, it can be clearly observed from Figure 2 that the ring soliton shrinks during the initial stage, and oscillations and radiations have emerged by t = 34.2 . Furthermore, the graphs also illustrate that the pulse simulated by the 2D sine-Gordon equation exhibits periodic oscillation. These phenomena are indeed valid, as other numerical methods have been employed to simulate this problem and exhibit similar phenomena. Here, we do not display the graphs obtained by other numerical methods. The numerical data listed in Table 2 demonstrates the convergence rate of the proposed EP3-FP scheme by varying the spatial and temporal step sizes. In Figure 3, after discretizing the spatial derivatives by the Fourier pseudo-spectral method with the fixed spatial mesh grid scales M = M_x = M_y = 64 , the problem is calculated by a different time-stepping scheme. These phenomena further validate the accuracy, efficiency, and long-term energy conservation of the EP3-FP scheme.

    Figure 2.  Snapshots of the numerical solution of the proposed EP3-FP scheme for solving Problem 2 at different times with the spatial mesh grid scales M = M_x = M_y = 256 and time step size \Delta t = 0.01 .
    Table 2.  The posterior errors and temporal convergence rates of the "EP3-FP" scheme for solving Problem 2.
    Error \Delta t=0.4 \Delta t=0.2 \Delta t=0.1 \Delta t=0.05
    M=8 5.0772{\rm E}-06 8.1347{\rm E}-08 1.2794{\rm E}-09 1.9459{\rm E}-11
    \textbf{Rate} \textbf{5.9638} \textbf{5.9906} \textbf{6.0388}
    M=16 3.4833{\rm E}-05 5.6543{\rm E}-07 8.9029{\rm E}-09 1.4501{\rm E}-10
    \textbf{Rate} \textbf{5.9449} \textbf{5.9889} \textbf{5.9401}
    M=32 3.0811{\rm E}-04 5.3090{\rm E}-06 8.5278{\rm E}-08 1.3415{\rm E}-09
    \textbf{Rate} \textbf{5.8588} \textbf{5.9601} \textbf{5.9902}
    M=64 1.8629{\rm E}-03 2.9562{\rm E}-05 4.7003{\rm E}-07 7.3794{\rm E}-09
    \textbf{Rate} \textbf{5.9777} \textbf{5.9748} \textbf{5.9931}

     | Show Table
    DownLoad: CSV
    Figure 3.  Results for Problem 2 with M_x = M_y = 64 . (a) the log-log plot of the posterior error against different steps h . (b) the log-log plot of the posterior error against CPU time. (c) the log plot of the relative energy error against integrate time with H(0) = 0.377193865036316 .

    In this paper, based on the blend of the energy-preserving continuous-stage RKN time integrator with the Fourier pseudo-spectral spatial discretization, we presented a novel energy-preserving and symmetric fully discrete scheme for solving the two-dimensional nonlinear wave equations. The discrete energy of the two-dimensional nonlinear wave Eqs (1.1) and (1.2) is well conserved by the proposed scheme. Meanwhile, another significant discovery is that the derived EP3-FP scheme can achieve sixth-order temporal accuracy under the low regularity assumption u\in C^2\big([t_0, T], \mathcal{B}\big) . Numerical experiments are presented to illustrate the accuracy, efficiency, and long-term energy conservation of the derived EP3-FP scheme.

    In light of a similar process, the derived EP3-FP scheme could be generalized to investigate other Hamiltonian PDEs, including the fractional nonlinear Hamiltonian wave equation, the Klein-Gordon equation with weak nonlinearity, the Klein-Gordon equation in the nonrelativistic limit regime, and the Klein-Gordon-Zakharov system.

    Dongjie Gao wrote the main manuscript. Peiguo Zhang prepared numerical experiments. Longqin Wang, Zhenlong Dai and Yonglei Fang made the corrections. All the authors have agreed and given their consent for the publication of this research paper.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    The work is supported by the Natural Science Foundation of China under Grant 11801280 and 12071419, the Natural Science Foundation of the Jiangsu Higher Education Institutions under Grant 21KJD110002, the Natural Science Foundation of Shandong Province under Grant ZR2024MA056, the foundation of innovative science and technology for youth in universities of Shandong Province under Grant 2023KJ278.

    The authors declare no conflict of interest.



    [1] N. Unnikrishnan, B. Vineshkumar, Reversed percentile residual life and related concepts, J. Korean Stat. Soc., 40 (2011), 85–92. https://doi.org/10.1016/j.jkss.2010.06.001 doi: 10.1016/j.jkss.2010.06.001
    [2] M. Shafaei Noughabi, S. Izadkhah, On the quantile past lifetime of the components of the parallel systems, Commun. Stat. Theory Meth., 45 (2016), 2130–2136. https://doi.org/10.1080/03610926.2013.875573 doi: 10.1080/03610926.2013.875573
    [3] M. Shafaei Noughabi, Solving a functional equation and characterizing distributions by quantile past lifetime functions, Econ. Qual. Control, 31 (2016), 55–58. https://doi.org/10.1515/eqc-2015-0017 doi: 10.1515/eqc-2015-0017
    [4] M. Mahdy, M. Further results involving percentile inactivity time order and its inference, Metron, 72 (2014), 269–282. https://doi.org/10.1007/s40300-013-0017-9 doi: 10.1007/s40300-013-0017-9
    [5] L. Balmert, J. H. Jeong, Nonparametric inference on quantile lost lifespan, Biometrics, 73 (2017), 252–259. https://doi.org/10.1111/biom.12555 doi: 10.1111/biom.12555
    [6] L. C. Balmert, R. Li, L. Peng, J. H. Jeong, Quantile regression on inactivity time, Stat. Meth. Med. Res., 30 (2021), 1332–1346. https://doi.org/10.1177/0962280221995977 doi: 10.1177/0962280221995977
    [7] A. M. Teamah Abd-Elmonem, A. Elbanna, A. M. Gemeay, Frechet-Weibull mixture distribution: properties and applications, Appl. Math. Sci., 14 (2020), 75–86. https://doi.org/10.12988/ams.2020.912165 doi: 10.12988/ams.2020.912165
    [8] I. Elbatal, S. M. Alghamdi, F. Jamal, S. Khan, E. M. Almetwally, M. Elgarhy, Kavya-manoharan weibull-g family of distributions: statistical inference under progressive type-ii censoring scheme, Adv. Appl. Stat., 87 (2023), 191–223. https://doi.org/10.17654/0972361723034 doi: 10.17654/0972361723034
    [9] C. F. Chung, Confidence bands for percentile residual lifetime under random censorship model, J. Multivariate Anal., 29 (1989), 94–126. https://doi.org/10.1016/0047-259X(89)90079-1 doi: 10.1016/0047-259X(89)90079-1
    [10] M. D. Burke, S. Csorgo, L. Horvath, Strong approximations of some biometric estimates under random censorship, Z. Wahrsch. Verw. Gebiete, 56 (1981), 87–112.
    [11] S. Csorgo, L. Horvath, The rate of strong uniform consistency for the product-limit estimator, Z. Wahrsch. Verw. Gebiete, 62 (1983), 411–426.
    [12] E. E. Aly, M. Csorgo, L. Horvath, Strong apprroximations of the quantile process of the product-limit estimator, J. Multivariate Anal., 16 (1985), 185–210. https://doi.org/10.1016/0047-259X(85)90033-8 doi: 10.1016/0047-259X(85)90033-8
    [13] M. Csorgo, S. Scorgo, Estimation of percentile residual life, Oper. Res., 35 (1987), 598–606. https://doi.org/10.1287/opre.35.4.598 doi: 10.1287/opre.35.4.598
    [14] J. R. Blum, V. Susarla, Maximal deviation theory of density and failure rate function estimates based on censored data, Multivariate Anal., 5 (1980), 213–222.
    [15] M. D. Burke, L. Horvath, Density and failure rate estimation in competing risks model, Sankhya A., 46 (1982), 135–154. https://doi.org/10.1016/j.ress.2010.04.006 doi: 10.1016/j.ress.2010.04.006
    [16] J. Mielniczuk, Some asymptotic properties of kernel estimators of a density function in case of censored data, Ann. Statist., 14 (1986), 766–773. https://doi.org/10.1214/aos/1176349954 doi: 10.1214/aos/1176349954
    [17] J. S. Marron, W. J. Padgett, Asymptotically optimal bandwidth selection for kernel density estimators from randomly right-censored samples, Ann. Statist., 15 (1987), 1520–1535. https://doi.org/10.1214/aos/1176350607 doi: 10.1214/aos/1176350607
    [18] S. H. Lo, Y. P. Mack, J. L. Wang, Density and hazard rate estimation for censored data via strong representation of the Kaplan-Meier estimator, Probab. Theory Related Fields, 80 (1989), 461–473.
    [19] T. R. Fleming, D. P. Harrington, Counting processes and survival analysis, New York: Wiley, 1991. https://doi.org/10.1002/9781118150672
  • Reader Comments
  • © 2024 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(974) PDF downloads(46) Cited by(4)

Figures and Tables

Figures(4)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog