Conservative Finite Element Methods for Kirchhoff-Type Equations

1.1 Problem Formulation

• The time-dependent spatially nonlocal coefficient depends on total Dirichlet energy, unlike the predominantly local or damped structures in related quasilinear-wave analyses.In two dimensions, the model describes membrane vibrations with effective tension depending on total deformation energy.Source: §1.1 · p. 2, §1.2 · p. 3

1.2 Related work and contribution

• The fully discrete method combines symmetric three-layer time stepping, conforming finite elements of arbitrary degree p, and Picard iteration for time-dependent-coefficient Kirchhoff equations.The nonlinear scheme has a unique solution at each time level, and the Picard iteration converges.Source: §1.2 · p. 3
• For the homogeneous problem with time-independent coefficients, the symmetric three-layer scheme conserves a discrete energy exactly for each choice of its free parameter.The companion paper [19] established temporal energy conservation and local second-order convergence; this paper supplies the spatial discretization.Source: §1.2 · p. 3
• The full error is O(h^p+1 + τ^2) in L2 and O(h^p + τ^2) in H1, matching the rates exhibited by the numerical experiments.The half-order gap in the estimates comes from the norms used in the proof, not from the scheme.Source: §1.2 · p. 3
• Compared with Dautov and Ivanova, the paper treats time-dependent coefficients using a symmetric three-layer scheme for the second-order formulation and arbitrary-degree finite elements.Their cited scheme is an implicit two-layer method for a first-order formulation and reports O(h2 + τ2) L2 accuracy in numerical tests.Source: §1.2 · p. 3

2 Setting and the Semi-Discrete Scheme

• The fractional powers of the Dirichlet Laplacian are defined spectrally, and the analysis uses operator domains for temporal regularity and Sobolev spaces for spatial approximation.The interpolation embedding is used to apply literature results stated in the Sobolev scale under operator-domain hypotheses.Source: §2 · pp. 3–4
• The semi-discrete time scheme is the starting point for the paper’s spatial discretization, and its solution converges locally with second order in time under stated regularity hypotheses.The central-difference approximation of the first time derivative has the same order.Source: §2 · p. 4
• The companion-paper results and existence of a time-discrete trajectory are not assumed; the fully discrete solution’s existence is proved in this paper.The continuous-problem assumptions are stated with each result, and higher regularity is required for convergence.Source: §2 · p. 4
• For time-dependent coefficients, existence of the continuous solution is taken as a hypothesis because the nonautonomous coefficient prevents applying the cited constant-coefficient theorem.The authors also do not claim that the cited theorem’s domain regularity hypotheses hold for every admitted domain.Source: §2 · p. 4

3.1 Preliminaries

• The analysis uses solution-centered balls to establish that fully discrete solutions remain bounded, supporting unique existence and convergence.The section states that the discrete solution lies inside these balls at all times.Source: §3.1 · p. 4
• On a convex polygonal domain with a quasi-uniform triangular mesh, the conforming finite element space uses piecewise polynomials of arbitrary degree p.The discrete operator Lh plays the role of the continuous operator L and is linked to the Ritz projection.Source: §3.1 · p. 5

3.2 Fully discrete scheme and energy conservation

• For homogeneous problems, the exact nonlinear scheme conserves its defined fully discrete energy at every time level.The conservation statement applies when the nonlinear scheme is solved exactly.Source: §3.2 · p. 6
• The exact-solution conservation identity does not hold exactly for Picard iterates; their energy conservation is limited by the nonlinear residual.In implementation, conservation is also limited by the nonlinear-solve tolerance and floating-point accuracy of the linear solves.Source: §3.2 · p. 7
• For θ in [1/4, 1/2), the discrete energy is not conserved; the stated estimate is formal and is not used in the analysis.The formal estimate relies on bounds stronger than those used for the discrete solution.Source: §3.2 · p. 7
• The deviation is bounded by O(τ^2) uniformly in n and does not accumulate over [0, T].The passage reports numerical confirmation in Section 5.Source: §3.2 · pp. 7–8

4.1 Local existence of a solution

• For sufficiently small τ, the nonlinear equation has a unique solution at each time level, and its Picard iteration converges, provided the discrete solutions satisfy the stated bounds.The contraction threshold is independent of the time-step index and mesh size; the local argument alone does not prevent recursively growing solution norms.Source: §4.1 · pp. 8, 10

4.2 Convergence of the local solutions

• Convergence to the Ritz projection supplies the induction step that keeps discrete solutions bounded and allows the local existence result to be applied recursively.The convergence theorem assumes θ ∈ [1/4, 1/2], with τ ≤ κh at θ = 1/4, as well as stronger spatial and temporal regularity and boundedness at previous time steps.Source: §4.2 · pp. 10, 12
• Summation by parts avoids a τ^-1 factor on the defect term, allowing the order-τ^2 defects to suffice in the local-solution analysis.The resulting second asserted estimate follows using uniform positive definiteness of the discrete Laplacian.Source: §4.2 · pp. 13–14

4.3 Initial data and closure of the induction

• The Taylor-projected starting values give purely temporal O(τ^2) error relative to the Ritz projection, with no spatial contribution.The starting values mirror those of the time-discrete scheme.Source: §4.3 · p. 14
• Under Theorem 4.4’s assumptions, including τ ≤ κh at θ = 1/4, the induction closes uniformly across all time levels for sufficiently small τ and h.The argument uses constants independent of n, τ, and h, so a single pair of smallness bounds applies at every step.Source: §4.3 · pp. 14–15
• The full error is O(h^p+1 + τ^2) in L2 and O(h^p + τ^2) in H1, both optimal in h for degree-p continuous piecewise polynomials.The bounds hold uniformly over all time levels.Source: §4.3 · p. 15

Numerical setup

• The experiments use manufactured solutions on Ω=(0,1)^2 with continuous Lagrange elements of degrees p=1, 2, and 3.The source is constructed from the prescribed exact solution and coefficients; λ controls spatial oscillation.Source: §5 · pp. 15–16

Example 1: Optimal spatial and temporal convergence studies

• For λ=1 and 3, observed spatial and temporal orders are independent of λ, although larger λ increases error; Picard stopping succeeds at every time step.The iteration never reaches Lmax in the reported runs.Source: §5 · p. 16

Example 2: Conservation of the discrete energy

• With α=β=1, zero source, and θ=1/2, energy is conserved over 100 periods for γ=0, 1/4, and 1/2, up to nonlinear-solve error.For γ=0 the linearly implicit scheme has round-off-level drift; the other cases' drift grows linearly at the Picard-tolerance rate.Source: §5 · pp. 16–17

Example 1: Optimal spatial and temporal convergence studies

• Spatial rates p+1 in L2 and p in H1, and second-order temporal convergence, are observed until spatial error dominates.The final p=3 spatial point at h=2^-6 is already at the temporal plateau, so its reduced rate is expected.Source: §5 · pp. 16–17, Figure 5.1 · p. 17

Example 2: Conservation of the discrete energy

• For γ=1/2, energy spread is independent of τ when θ=1/2; when θ≠1/2, energy is approximated to second order in τ.The plotted spread is the largest relative difference between any two energy values over the run.Source: §5 · p. 17, Figure 5.2 · p. 18

Conclusion

• The scheme combines symmetric three-layer time stepping, conforming finite elements of arbitrary degree p, and Picard iteration for the nonlinearity.The nonlinear equation at each time level has a unique solution, and the Picard iteration converges to it.Source: §6 · p. 18
• For the homogeneous problem with constant coefficients and θ=1/2, the exact scheme solution conserves the discrete energy identically.The authors identify this as the fully discrete counterpart of temporal invariants and attribute its survival under spatial discretization to conformity.Source: §6 · p. 18
• The full error converges at order O(h^(p+1)+τ^2) in L2 and O(h^p+τ^2) in H1; experiments exhibit these rates.The temporal second-order rate is observed until spatial error becomes dominant.Source: §6 · p. 18
• The estimates hold only on the local interval [0,T] with T<T**, require strong regularity, and do not establish unconditional stability for time-dependent coefficients.Extending the convergence results to blow-up time would require a different argument; optimal rates under weaker regularity remain open.Source: §6 · p. 18
• The authors leave open whether non-coercivity at θ=1/4 is necessary for the scheme or only for the proof, and whether energy-identity verification with solver tolerances applies to convergence-study spaces.A systematic study of sensitivity to the Picard and conjugate-gradient tolerances is also outside scope.Source: §6 · pp. 18–19

Declarations

• The work reports support from the DFG, including SFB 1173 and Excellence Strategy funding, and from the Shota Rustaveli National Science Foundation of Georgia.The DFG project identifiers and grant number are specified in the declarations.Source: Declarations · p. 19
• The authors declare no competing interests.Source: Declarations · p. 19
• Raw error data, figure sources, and the Python/FEniCSx implementation are available in the accompanying repository.Source: Declarations · p. 19
• The authors disclose using Claude to improve manuscript formulation and assist with code scripts and software presentation, and state that they reviewed and edited the content.Source: Declarations · p. 19

A.1 Proof of the full-order convergence

• The defect estimate obtained from the appendix argument does not depend on the parameter a.The argument uses boundedness of the fully discrete solution at prior levels and the uniform bounds supplied by Lemma 4.1 at the top index.Source: §A.1 · pp. 21–22
• The first source-defect term is rewritten as P_h(f^m−J_hf^m), since J_hf^m∈V_h, then estimated using (3.5) at a∈{0,1/4} followed by (3.4).Source: §A.1 · p. 22
• The central second-difference consistency term is written in integral form so that the proof requires only the stated fourth-time-derivative regularity.The passage identifies ∂^4_t u∈L∞(0,T;D(L^(1/4))) as the regularity used for this term.Source: §A.1 · p. 22
• Ritz orthogonality makes the projection contribution ρ^m equal to a negative squared H1 projection error, bounded by O(h^(2p)).Nonnegative weights summing to one preserve this bound for their weighted combination.Source: §A.1 · p. 23

• P_h and J_h are L^2- and H^1_0-stable, supporting fractional-power bounds by interpolation between those endpoint norms.P_h has L^2 stability with constant 1; its H^1_0 stability holds on quasi-uniform meshes, or on shape-regular families under a local mesh condition.Source: §A.2 · p. 24
• The discrete Laplacian's eigenvalues are bounded below independently of h, while its largest eigenvalue is bounded by C_inv^2 h^-2.These spectral bounds support the stated estimates for fractional powers of L_h.Source: §A.2 · p. 24
• R_h is H^1_0-stable with constant 1, and its weaker-norm estimates use the approximation bounds and commutation identity.Source: §A.2 · pp. 24–25
• The Ritz projection satisfies L_h R_h v = P_h L v, yielding a constant-one bound in the discrete Laplacian norm.This follows from the Ritz projection's H^1_0 orthogonality and the L^2 stability of P_h.Source: §A.2 · pp. 24–25