A two-dimensional discrete delay-differential system model of viremia

1.   Introduction

In this paper, we aim to develop a novel weak Galerkin (WG) finite element method for the biharmonic equation that is applicable to non-convex polytopal meshes and eliminates the need for traditional stabilizers. To this aim, we consider the biharmonic equation with Dirichlet and Neumann boundary conditions. The goal is to find an unknown function $u$ satisfying

where $\Omega\subset \mathbb R^d$ is an open bounded domain with a Lipschitz continuous boundary $\partial\Omega$. The domain $\Omega$ considered in this paper can be of any dimension $d\geq 2$. For the sake of clarity in presentation, we will focus on the case where $d = 2$ throughout this paper. However, the analysis presented here can be readily extended to higher dimensions ($d\geq 3$) without significant modifications.

The variational formulation of the model problem (1.1) is as follows: Find an unknown function $u\in H^2(\Omega)$ satisfying $u|_{\partial\Omega} = \xi$ and $\frac{\partial u}{\partial {{\mathbf{n}}}}|_{\partial\Omega} = \nu$, and the following equation

where $\partial^2_{ij}$ denotes the second order partial derivative with respect to $x_i$ and $x_j$, and $H_0^2(\Omega) = \{v\in H^2(\Omega): v|_{\partial\Omega} = 0, \nabla v |_{\partial\Omega} = 0\}$.

The WG finite element method offers an innovative framework for the numerical solution of partial differential equations (PDEs). This approach approximates differential operators within a structure inspired by the theory of distributions, particularly for piecewise polynomial functions. Unlike traditional methods, WG reduces the regularity requirements on approximating functions through the use of carefully designed stabilizers. Extensive studies have highlighted the versatility and effectiveness of WG methods across a wide range of model PDEs, as demonstrated by numerous references ^{[1,2,3,4,5,6]} for an incomplete list, establishing WG as a powerful tool in scientific computing. The defining feature of WG methods lies in their innovative use of weak derivatives and weak continuities to construct numerical schemes based on the weak forms of the underlying PDEs. This unique structure provides WG methods with exceptional flexibility, enabling them to address a wide variety of PDEs while ensuring both stability and accuracy in their numerical solutions.

This paper presents a simplified formulation of the WG finite element method, capable of handling both convex and non-convex elements in finite element partitions. A key innovation of our method is the elimination of stabilizers through the use of higher-degree polynomials for computing weak second-order partial derivatives. This design preserves the size and global sparsity of the stiffness matrix while substantially reducing the programming complexity associated with traditional stabilizer-dependent methods. The method leverages bubble functions as a critical analytical tool, representing a significant improvement over existing stabilizer-free WG methods ^[7], which are limited to convex polytopal meshes. Our approach is versatile, accommodating arbitrary dimensions and polynomial degrees in the discretization process. In contrast, prior stabilizer-free WG methods ^[7] often require specific polynomial degree combinations and are restricted to 2D or 3D settings. Theoretical analysis establishes optimal error estimates for the WG approximations in both the discrete $H^2$ norm and an $L^2$ norm.

This paper is organized as follows. Section 2 provides a brief review of the definition of the weak-second order partial derivative and its discrete counterpart. In Section 3, we introduce an efficient WG scheme that eliminates the need for stabilization terms. Section 4 establishes the existence and uniqueness of the solution. The error equation for the proposed WG scheme is derived in Section 5. Section 6 focuses on obtaining the error estimate for the numerical approximation in the discrete $H^2$ norm, while Section 7 extends the analysis to derive the error estimate in the $L^2$ norm.

Throughout this paper, we adopt standard notations. Let $D$ be any open, bounded domain with a Lipschitz continuous boundary in $\mathbb{R}^d$. The inner product, semi-norm, and norm in the Sobolev space $H^s(D)$ for any integer $s\geq0$ are denoted by $(\cdot, \cdot)_{s, D}$, $|\cdot|_{s, D}$ and $\|\cdot\|_{s, D}$ respectively. For simplicity, when the domain $D$ is $\Omega$, the subscript $D$ is omitted from these notations. In the case $s = 0$, the notations $(\cdot, \cdot)_{0, D}$, $|\cdot|_{0, D}$ and $\|\cdot\|_{0, D}$ are further simplified as $(\cdot, \cdot)_D$, $|\cdot|_D$ and $\|\cdot\|_D$, respectively.

2.   Discrete weak second order partial derivatives

This section provides a brief review of the definition of weak weak-second partial derivatives and their discrete counterparts, as introduced in ^[5].

Let $T$ be a polygonal element with boundary $\partial T$. A weak function on $T$ is represented as $v = \{v_0, v_b, {{\mathbf{v}}}_g\}$, where $v_0\in L^2(T)$, $v_b\in L^{2}(\partial T)$ and ${{\mathbf{v}}}_g\in [L^2(\partial T)]^2$. The first component, $v_0$, denotes the value of $v$ within the interior of $T$, while the second component, $v_b$, represents the value of $v$ on the boundary of $T$. The third component ${{\mathbf{v}}}_g\in \mathbb R^2$ with components $v_{gi}$ ($i = 1, 2$) approximates the gradient $\nabla v$ on the boundary $\partial T$. In general, $v_b$ and ${{\mathbf{v}}}_g$ are treated as independent of the traces of $v_0$ and $\nabla v_0$, respectively.

The space of all weak functions on $T$, denote by $W(T)$, is defined as

The weak second order partial derivative, $\partial^2_{ij, w}$, is a linear operator mapping $W(T)$ to the dual space of $H^2(T)$. For any $v\in W(T)$, $\partial^2_{ij, w} v$ is defined as a bounded linear functional on $H^2(T)$, given by:

where ${{\mathbf{n}}}$, with components $n_i\; (i = 1, 2)$, represents the unit outward normal vector to $\partial T$.

For any non-negative integer $r\ge 0$, let $P_r(T)$ denote the space of polynomials on $T$ with total degree at most $r$. A discrete weak second order partial derivative, $\partial^2_{ij, w, r, T}$, is a linear operator mapping $W(T)$ to $P_r(T)$. For any $v\in W(T)$, $\partial^2_{ij, w, r, T}v$ is the unique polynomial in $P_r(T)$ satisfying

For a smooth $v_0\in H^2(T)$, applying standard integration by parts to the first term on the right-hand side of (2.1) yields:

for any $\varphi \in P_r(T)$.

3.   Weak Galerkin algorithms without stabilization terms

Let ${\mathcal{T}}_h$ be a finite element partition of the domain $\Omega\subset \mathbb R^2$ into polygons. Assume that ${\mathcal{T}}_h$ satisfies the shape regularity condition ^[8]. Let ${\mathcal E}_h$ represent the set of all edges in ${\mathcal{T}}_h$, and denote the set of interior edges by ${\mathcal E}_h^0 = {\mathcal E}_h \setminus \partial\Omega$. For any element $T\in {\mathcal{T}}_h$, let $h_T$ be its diameter, and define the mesh size as $h = \max_{T\in {\mathcal{T}}_h}h_T$.

Let $k$, $p$ and $q$ be integers such that $k\geq p\geq q\geq 1$. For any element $T\in{{\mathcal{T}}}_h$, the local weak finite element space is defined as:

By combining the local spaces $V(k, p, q, T)$ across all elements $T\in {\mathcal{T}}_h$ and ensuring continuity of $v_b$ and ${{\mathbf{v}}}_g$ along the interior edges ${{\mathcal{E}}}_h^0$, we obtain the global weak finite element space:

The subspace of $V_h$ consisting of functions with vanishing boundary values on $\partial\Omega$ is defined as:

For simplicity, the discrete weak second order partial derivative $\partial^2_{ij, w}v$ is used to denote the operator $\partial^2_{ij, w, r, T} v$ defined in (2.1) on each element $T\in {\mathcal{T}}_h$, as:

On each element $T\in{{\mathcal{T}}}_h$, let $Q_0$ denote the $L^2$ projection onto $P_k(T)$. On each edge $e\subset\partial T$, let $Q_b$ and $Q_n$ denote the $L^2$ projection operators onto $P_{p}(e)$ and $P_{q}(e)$, respectively. For any $w\in H^2(\Omega)$, the $L^2$ projection into the weak finite element space $V_h$ is denoted by $Q_h w$, defined as:

The simplified WG numerical scheme, free from stabilization terms, for solving the biharmonic equation (1.1) is formulated as follows:

Weak Galerkin Algorithm 3.1. Find $u_h = \{u_0, u_b, {{\mathbf{u}}}_g\} \in V_h$ such that $u_b = Q_b\xi$, ${{\mathbf{u}}}_g\cdot{{\mathbf{n}}} = Q_n\nu$ and ${{\mathbf{u}}}_g\cdot{\boldsymbol{\tau}} = Q_n(\nabla\xi\cdot{\boldsymbol{\tau}})$ on $\partial\Omega$, and satisfy:

where ${\boldsymbol{\tau}}\in \mathbb R^2$ is the tangential direction along $\partial\Omega$, and the terms are defined as:

4.   Solution existence and uniqueness

Recall that ${\mathcal{T}}_h$ is a shape-regular finite element partition of the domain $\Omega$. Consequently, for any $T\in {\mathcal{T}}_h$ and $\phi\in H^1(T)$, the following trace inequality holds ^[8]:

If $\phi$ is a polynomial on the element $T\in {\mathcal{T}}_h$, a simpler form of the trace inequality applies ^[8]:

For any $v = \{v_0, v_b, {{\mathbf{v}}}_g\}\in V_h$, define the norm:

and introduce the discrete $H^2$- semi-norm:

Lemma 4.1. For $v = \{v_0, v_b, {{\mathbf{v}}}_g\}\in V_h$, there exists a constant $C$ such that for $i, j = 1, 2$,

Proof. Let $T\in {\mathcal{T}}_h$ be a polytopal element with $N$ edges denoted as $e_1, \cdots$, $e_N$. Importantly, $T$ can be non-convex. On each edge $e_i$, construct a linear function $l_i(x)$ satisfying $l_i(x) = 0$ on $e_i$ as:

where $A$ is a fixed point on $e_i$, $X$ is any point on $e_i$, ${{\mathbf{n}}}_i$ is the normal vector to $e_i$, and $h_T$ is the diameter of $T$.

Define the bubble function for $T$ as:

It is straightforward to verify that $\Phi_B = 0$ on $\partial T$. The function $\Phi_B$ can be scaled such that $\Phi_B(M) = 1$ where $M$ is the barycenter of $T$. Additionally, there exists a subdomain $\hat{T}\subset T$ such that $\Phi_B\geq \rho_0$ for some constant $\rho_0 > 0$.

For $v = \{v_0, v_b, {{\mathbf{v}}}_g\}\in V_h$, let $r = 2N+k-2$ and choose $\varphi = \Phi_B \partial^2_{ij} v_0\in P_r(T)$ in (2.2). This yields:

where we applied $\Phi_B = 0$ on $\partial T$.

Using the domain inverse inequality ^[8], there exists a constant $C$ such that

By applying the Cauchy-Schwarz inequality to (4.5) and (4.6), we obtain

which implies:

This completes the proof.

Remark 4.1. If the polytopal element $T$ is convex, the bubble function in Lemma 4.1 can be simplified to:

This simplified bubble function satisfies 1) $\Phi_B = 0$ on $\partial T$, 2) there exists a subdomain $\hat{T}\subset T$ such that $\Phi_B\geq \rho_0$ for some constant $\rho_0 > 0$. The proof of Lemma 4.1 follows the same approach, using this simplified bubble function. In this case, we set $r = N+k-2$.

By constructing an edge-based bubble function,

it can be easily verified that 1) $\varphi_{e_k} = 0$ on each edge $e_i$ for $i \neq k$, and 2) there exists a subdomain $\widehat{e_k}\subset e_k$ such that $\varphi_{e_k} \geq \rho_1$ for some constant $\rho_1 > 0$. Let $\varphi = (v_b-v_0)l_k \varphi_{e_k}$. It is straightforward to verify the following properties: 1) $\varphi = 0$ on each edge $e_i$ for $i = 1, \cdots, N$, 2) $\nabla \varphi = 0$ on each edge $e_i$ for $i \neq k$, and 3) $\nabla \varphi = (v_0-v_b)(\nabla l_k) \varphi_{e_k} = \mathcal{O}(\frac{ (v_0-v_b)\varphi_{e_k}}{h_T}\textbf{C})$ on $e_k$, where $\textbf{C}$ is a constant vector.

Lemma 4.2. ^[9] For $\{v_0, v_b, {{\mathbf{v}}}_g\}\in V_h$, let $\varphi = (v_b-v_0)l_k \varphi_{e_k}$. The following inequality holds:

Lemma 4.3. For $\{v_0, v_b, {{\mathbf{v}}}_g\}\in V_h$, let $\varphi = (v_{gi} -\partial_i v_0) \varphi_{e_k}$. The following inequality holds:

Proof. Define the extension of ${{\mathbf{v}}}_g$, originally defined on the edge $e_k$, to the entire polytopal element $T$ as:

where $X = (x_1, x_2)$ is any point in $T$, $Proj_{e_k} (X)$ denotes the orthogonal projection of $X$ onto the plane $H\subset\mathbb R^2$ containing $e_k$. If $Proj_{e_k} (X)$ is not on $e_k$, ${{\mathbf{v}}}_g(Proj_{e_k} (X))$ is defined as the extension of ${{\mathbf{v}}}_g$ from $e_k$ to $H$. The extension preserves the polynomial nature of ${{\mathbf{v}}}_g$ as demonstrated in ^[9].

Let $v_{trace}$ denote the trace of $v_0$ on $e_k$. Define its extension to $T$ as:

This extension is also polynomial, as demonstrated in ^[9].

Let $\varphi = (v_{gi} -\partial_i v_0) \varphi_{e_k}$. Then,

where we used the facts that 1) $\varphi_{e_k} = 0$ on each edge $e_i$ for $i \neq k$, 2) there exists a subdomain $\widehat{e_k}\subset e_k$ such that $\varphi_{e_k} \geq \rho_1$ for some constant $\rho_1 > 0$, and applied the properties of the projection.

This completes the proof of the lemma.

Lemma 4.4. There exist two positive constants, $C_1$ and $C_2$, such that for any $v = \{v_0, v_b, {{\mathbf{v}}}_g\} \in V_h$, the following equivalence holds:

Proof. Consider the edge-based bubble function defined as

First, extend $v_b$ from the edge $e_k$ to the element $T$. Similarly, let $v_{trace}$ denote the trace of $v_0$ on the edge $e_k$ and extend $v_{trace}$ to the element $T$. For simplicity, we continue to use $v_b$ and $v_0$ to represent their respective extensions. Details of these extensions can be found in Lemma 4.3. Substituting $\varphi = (v_b-v_0)l_k\varphi_{e_k}$ into (2.2), we obtain

where we used 1) $\varphi = 0$ on each edge $e_i$ for $i = 1$, $\cdots$, $N$, 2) $\nabla \varphi = 0$ on each edge $e_i$ for $i \neq k$, and 3) $\nabla \varphi = (v_0-v_b)(\nabla l_k) \varphi_{e_k} = \mathcal{O}(\frac{ (v_0-v_b)\varphi_{e_k}}{h_T}\textbf{C})$ on $e_k$, where $\textbf{C}$ is a constant vector.

Recall that 1) $\varphi_{e_k} = 0$ on each edge $e_i$ for $i \neq k$, and 2) there exists a subdomain $\widehat{e_k}\subset e_k$ such that $\varphi_{e_k} \geq \rho_1$ for some constant $\rho_1 > 0$. Using Cauchy-Schwarz inequality, the domain inverse inequality ^[8], (4.10) and Lemma 4.2, we deduce:

which, from Lemma 4.1, leads to:

Next, extend ${{\mathbf{v}}}_g$ from the edge $e_k$ to the element $T$, denoting the extension by the same symbol for simplicity. Details of this extension are in Lemma 4.3. Substituting $\varphi = (v_{gi}-\partial_i v_0)\varphi_{e_k}$ into (2.2), we obtain:

where we used $\varphi_{e_k} = 0$ on edge $e_i$ for $i \neq k$, and the fact that there exists a sub-domain $\widehat{e_k}\subset e_k$ such that $\varphi_{e_k} \geq \rho_1$ for some constant $\rho_1 > 0$. This, together with Cauchy-Schwarz inequality, the domain inverse inequality ^[8], the inverse inequality, the trace inequality (4.2), (4.11) and Lemma 4.3, gives

Applying Lemma 4.1, gives

Using Lemma 4.1, Eqs (4.11), (4.13), (4.3) and (4.4), we deduce:

Finally, using the Cauchy-Schwarz inequality, inverse inequalities, and the trace inequality (4.2) in (2.2), we derive:

which gives:

and further gives

This completes the proof.

Theorem 4.5. The WG scheme 3.1 admits a unique solution.

Proof. Assume that $u_h^{(1)}\in V_h$ and $u_h^{(2)}\in V_h$ are two distinct solutions of the WG scheme 3.1. Define $\eta_h = u_h^{(1)}-u_h^{(2)}\in V_h^0$. Then, $\eta_h$ satisfies

Choosing $v = \eta_h$ yields $||| \eta_h||| = 0$. From the equivalence of norms in (4.9), it follows that $\|\eta_h\|_{2, h} = 0$, which yields $\partial^2_{ij} \eta_0 = 0$ for $i, j = 1, 2$ on each $T$, $\eta_0 = \eta_b$ and $\nabla \eta_0 = {\boldsymbol{\eta}}_g$ on each $\partial T$. Consequently, $\eta_0$ is a linear function on each element $T$ and $\nabla \eta_0 = C$ on each $T$.

Since $\nabla \eta_0 = {\boldsymbol{\eta}}_g$ on each $\partial T$, it follows that $\nabla \eta_0$ is continuous across the entire domain $\Omega$. Thus, $\nabla \eta_0 = C$ throughout $\Omega$. Furthermore, the condition ${\boldsymbol{\eta}}_g = 0$ on $\partial\Omega$ implies $\nabla \eta_0 = 0$ in $\Omega$ and ${\boldsymbol{\eta}}_g = 0$ on each $\partial T$. Therefore, $\eta_0$ is a constant on each element $T$.

Since $\eta_0 = \eta_b$ on $\partial T$, the continuity of $\eta_0$ over $\Omega$ implies $\eta_0$ is globally constant. From $\eta_b = 0$ on $\partial\Omega$, we conclude $\eta_0 = 0$ throughout $\Omega$. Consequently, $\eta_b = \eta_0 = 0$ on each $\partial T$, which implies $\eta_h\equiv 0$ in $\Omega$. Thus, $u_h^{(1)}\equiv u_h^{(2)}$, proving the uniqueness of the solution.

5.   Error equations

Let $Q_r$ denote the $L^2$ projection operator onto the finite element space of piecewise polynomials of degree at most $r$.

Lemma 5.1. The following property holds:

Proof. For any $u\in H^2(T)$, using (2.2), we have

for all $\varphi\in P_r(T)$. This completes the proof.

Let $u$ be the exact solution of the biharmonic equation (1.1), and $u_h \in V_{h}$ its numerical approximation obtained from the WG scheme 3.1. The error function, denoted by $e_h$, is defined as

Lemma 5.2. The error function $e_h$ defined in (5.2) satisfies the following error equation:

where

Proof. Using (5.1), standard integration by parts, and substituting $\varphi = Q_r \partial_{ij}^2 u$ into (2.2), we obtain

where we used (1.1), $\partial_{ij}^2 v_0\in P_{k-2}(T)$ and $r = 2N+k-2\geq k-2$, $\sum_{T\in {\mathcal{T}}_h} \sum_{i, j = 1}^2 \langle \partial_{ij}^2 u, v_{gi}\cdot n_j\rangle_{\partial T} = \sum_{T\in {\mathcal{T}}_h} \sum_{i, j = 1}^2 \langle \partial_{ij}^2 u, v_{gi}\cdot n_j\rangle_{\partial \Omega} = 0$ since $v_{gi} = 0$ on $\partial \Omega$, and $\sum_{T\in {\mathcal{T}}_h} \sum_{i, j = 1}^2 \langle \partial_j(\partial_{ij}^2u)\cdot n_i, v_b\rangle_{\partial T} = \sum_{T\in {\mathcal{T}}_h} \sum_{i, j = 1}^2 \langle \partial_j(\partial_{ij}^2u)\cdot n_i, v_b\rangle_{\partial \Omega} = 0$ since $v_{b} = 0$ on $\partial \Omega$.

Subtracting (3.1) from (5.4) yields

This concludes the proof.

6.   Error estimates in the $H^2$ norm

Lemma 6.1. ^[5] Let ${\mathcal{T}}_h$ be a finite element partition of the domain $\Omega$ satisfying the shape regularity assumption specified in ^[8]. For any $0\leq s \leq 2$ and $1\leq m \leq k$, the following estimates hold:

Lemma 6.2. If $u\in H^{k+1}(\Omega)$, then there exists a constant $C$ such that

Proof. Utilizing (2.2), the trace inequalities (4.1) and (4.2), the inverse inequality, and the estimate (6.2) for $m = k$ and $s = 0, 1, 2$, we analyze the following summation for any $\varphi\in P_r(T)$:

Letting $\varphi = \partial^2_{ij, w}(u-Q_hu)$ gives

This completes the proof.

Theorem 6.3. Suppose the exact solution $u$ of the biharmonic equation (1.1) satisfies $u\in H^{k+1}(\Omega)$. Then, the error estimate satisfies:

Proof. Note that $r\geq 1$. For the first term on the right-hand side of the error equation (5.3), using Cauchy-Schwarz inequality, the trace inequality (4.1), the estimate (6.1) with $m = k$ and $s = 1, 2$, and (4.9), we have

For the second term on the right-hand side of the error equation (5.3), using the Cauchy-Schwarz inequality, the trace inequality (4.1), the estimate (6.1) with $m = k$ and $s = 0, 1$, and (4.9), we have

Substituting (6.5) and (6.6) into (5.3) gives

Using Cauchy-Schwarz inequality, letting $v = Q_hu-u_h$ in (6.7), the error estimate (6.3) gives

which further gives

This completes the proof.

7.   $L^2$ error estimate

To derive the error estimate in the $L^2$ norm, we use the standard duality argument. The error is expressed as $e_h = u-u_h = \{e_0, e_b, {{\mathbf{e}}}_g\}$, and we define $\zeta_h = Q_hu - u_h = \{\zeta_0, \zeta_b, {{\boldsymbol{\zeta}}}_g\}\in V_h^0$. Consider the dual problem associated with the biharmonic equation (1.1), which seeks a function $w \in H_0^2(\Omega)$ satisfying:

We assume the following regularity condition for the dual problem:

Theorem 7.1. Let $u\in H^{k+1}(\Omega)$ be the exact solution of the biharmonic equation (1.1), and let $u_h\in V_h$ denote the numerical solution obtained using the weak Galerkin scheme 3.1. Assume that the $H^4$-regularity condition (7.2) holds. Then, there exists a constant $C$ such that

Proof. Testing the dual problem (7.1) with $\zeta_0$ and applying integration by parts, we derive:

where we used $\sum_{T\in {\mathcal{T}}_h} \sum_{i, j = 1}^2 \langle \partial^2_{ij} w, \zeta_{gi} \cdot n_j \rangle_{\partial T} = \sum_{i, j = 1}^2\langle \partial^2_{ij} w, \zeta_{gi} \cdot n_j \rangle_{\partial \Omega} = 0$ due to ${{\boldsymbol{\zeta}}}_g = 0$ on $\partial\Omega$, and $\sum_{T\in {\mathcal{T}}_h} \sum_{i, j = 1}^2 \langle \partial_j(\partial^2_{ij} w)\cdot n_i, \zeta_b \rangle_{\partial T} = \sum_{i, j = 1}^2\langle \partial_j(\partial^2_{ij} w)\cdot n_i, \zeta_b \rangle_{\partial \Omega} = 0$ due to $\zeta_b = 0$ on $\partial\Omega$.

Letting $u = w$ and $v = \zeta_h$ in (5.4) gives

which is equivalent to

Substituting the above equation into (7.3) and using (5.3) gives

We will estimate the four terms $J_i \; (i = 1$, $\cdots$, $4)$ on the last line of (7.4) individually.

For $J_1$, using the Cauchy-Schwarz inequality, the trace inequality (4.1), the inverse inequality, the estimate (6.1) with $m = k$ and $s = 0, 1, 2$, the estimate (6.2) with $m = 3$ and $s = 0, 1, 2$, gives

For $J_2$, using Cauchy-Schwarz inequality, (6.3) with $k = 3$ and (6.4) gives

For $J_3$, denote by $Q^1$ a $L^2$ projection onto $P_1(T)$. Using (2.1) gives

where we used $\partial^2_{ji} (Q^1\partial^2_{ij, w} w) = 0$, $\partial_j (Q^1\partial^2_{ij, w} w) = C$ and the property of the projection operators $Q_b$ and $Q_n$ and $p\geq q\geq 1$.

Using (7.7), Cauchy-Schwarz inequality, (5.1) and (6.3), gives

For $J_4$, using Cauchy-Schwarz inequality, the trace inequality (4.1), Lemma 4.4, the estimate (6.1) with $m = 3$ and $s = 0, 1$, (6.3), (6.4) gives

Substituting (7.5), (7.6), (7.8) and (7.9) into (7.4), and using (7.2), gives

This gives

which, using the triangle inequality and (6.2) with $m = k$ and $s = 0$, gives

This completes the proof of the theorem.

Use of AI tools declaration

The author declares she has not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research of Chunmei Wang was partially supported by National Science Foundation Grant DMS-2136380.

Conflict of interest

The author declares there is no conflict of interest.