1 Introduction↩︎

The biharmonic equation plays a fundamental role in various fields of mathematical physics and engineering. Over the past few decades, extensive research has been devoted to developing efficient and accurate numerical schemes for this typical fourth-order partial differential equation. Conventionally, finite element methods and finite difference schemes [1] have been widely employed to solve such problems. In recent years, various stabilized and efficient mixed methods have been proposed to improve computational stability and convergence. To alleviate the stringent \(C^1\)-continuity requirement of standard conforming elements, the mixed finite element method was introduced and subsequently analyzed in [2], [3]. For more details, we refer the readers to [4]–[7] and the references cited therein. On the other hand, spectral methods are widely adopted in scientific computing for their capability to achieve exponential convergence rates for biharmonic equations. While extensive studies exist for rectangular or polygonal domains, employing Legendre and Chebyshev polynomials[8], [9], spectral collocation methods[10], as well as Legendre spectral-element approaches [11], research on canonical domains such as the three-dimensional unit ball requires special consideration. The theoretical foundation for using orthogonal polynomials of several variables was established in [12], with subsequent spectral approximation theories on the unit ball in [13]. However, efficient direct spectral solvers and rigorous a-priori error analysis for the biharmonic equations remain largely unexplored compared to their second-order counterparts.

In this paper, we propose and analyze a novel mixed spectral-Galerkin method for solving the biharmonic equation with simply supported boundary conditions: \[\label{eq:biharmonic95model} \left\{ \begin{align} \Delta^2 u &= f, \quad \text{in}\;\mathbb{B}^3, \\ u & = 0, \quad \text{on}\;\mathbb{S}^2,\\ \Delta u & = 0, \quad \text{on}\;\mathbb{S}^2, \end{align} \right.\tag{1}\] where \(\mathbb{B}^3\) denotes the three-dimensional unit ball and \(\mathbb{S}^2\) denotes the unit sphere of \(\mathbb{B}^3\). We construct an efficient direct spectral-Galerkin scheme by utilizing generalized ball polynomials as basis functions. A significant advantage of this approach is that it is directly formulated on the ball, naturally avoiding the coordinate singularities associated with standard spherical transformations. By introducing an auxiliary variable, the discrete scheme yields a strictly diagonal stiffness matrix, which greatly reduces the computational complexity. Beyond the solver construction, a primary focus of this work is the establishment of a rigorous theoretical framework for the a-priori error analysis in both \(L^2\)- and \(H^1\)-norms. Furthermore, we conduct extensive numerical experiments to validate our theoretical results and demonstrate the efficiency of the proposed scheme. We demonstrate that the numerical method achieves spectral accuracy for analytic and transcendental smooth solutions. The numerical results highly align with our theoretical error estimates, confirming the robustness and high precision of the mixed spectral formulation for biharmonic problems on the unit ball.

The remainder of this paper is organized as follows. In Section 2, we introduce the necessary notations and fundamental properties of Jacobi polynomials and generalized ball polynomials. Section 3 details the mixed spectral-Galerkin formulation and rigorously derives the corresponding a-priori error estimates in both the \(L^2\)- and \(H^1\)-norms. Section 4 provides numerical experiments to validate the theoretical results and demonstrate the efficiency of our proposed method. Finally, we draw a short conclusion in Section 5.

2 Preliminaries and Notations↩︎

We first introduce some standard notations that will be used throughout this paper. Let \(\mathbb{N}\) and \(\mathbb{N}_0\) denote the set of positive and nonnegative integers, respectively. For any vector \(x \in \mathbb{R}^3\), we write \(x = (x_1, x_2, x_3)^{\top}\) as a column vector, where \((\cdot)^{\top}\) denotes the transpose operation. Let \(\mathbb{B}^3 = \{x \in \mathbb{R}^3: \|x\| \leq 1\}\) be unit ball centered at the origin. We adopt the standard Sobolev space \(H^m(\mathbb{B}^3)\) with \(m \in \mathbb{N}_0\) on \(\mathbb{B}^3\), which is equipped with the norm \(\|\cdot\|_{m,\mathbb{B}^3}\) and semi-norm \(|\cdot|_{m,\mathbb{B}^3}\). The space \(H_{0}^m(\mathbb{B}^3)\) is defined as the closure of \(C_0^\infty(\mathbb{B}^3)\) within \(H^m(\mathbb{B}^3)\). Specially, \(H^0(\mathbb{B}^3) = L^2(\mathbb{B}^3)\), and the corresponding inner-product reads: \[\label{eq:inner95form} (u, v)_{L^2(\mathbb{B}^3)} = \int_{\mathbb{B}^3} u(x) v(x) \, \mathrm{d}x,\tag{2}\] and the corresponding norm is given by \(\|u\| = \sqrt{(u, u)_{L^2(\mathbb{B}^3)}}\).

Throughout this paper, we use \(C\) and \(c\) to denote generic positive constants that are independent of the polynomial degree and the domain, which may take different values at different occurrences.

2.1 Jacobi polynomials and their Properties↩︎

Let \(I = (-1, 1)\) and \(w^{\alpha, \beta}(t) = (1-t)^\alpha (1+t)^\beta\) be a weight function with \(\alpha, \beta > -1\). The Jacobi polynomials \(\{P_n^{\alpha, \beta}(t)\}_{n \ge 0}\) form an orthogonal basis set for \(L^2_{w^{\alpha, \beta}}(I)\) [12], which satisfy the following orthogonality identity: \[\label{eq:Jacobi-orth} \int_I P_n^{\alpha, \beta}(t) P_m^{\alpha, \beta}(t) w^{\alpha, \beta}(t) \, \mathrm{d}t = \hat{h}_n^{\alpha, \beta} \delta_{nm},\tag{3}\] where \[\label{eq:Jacobi-norm} \hat{h}_n^{\alpha, \beta} = \frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1} \frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{n!\Gamma(n+\alpha+\beta+1)}.\tag{4}\]

Note that \(P_n^{\alpha, \beta}\) satisfies the Sturm-Liouville equation: \[\label{eq:Jacobi-SL} -\partial_t \left( w^{\alpha+1, \beta+1}(t) \partial_t P_n^{\alpha, \beta}(t) \right) = \lambda_n^{\alpha, \beta} w^{\alpha, \beta}(t) P_n^{\alpha, \beta}(t)\tag{5}\] with the corresponding eigenvalue \(\lambda_n^{\alpha, \beta} = n(n+\alpha+\beta+1)\). Moreover, the derivative of \(P_n^{\alpha, \beta}(t)\) relates to a parameter-shifted Jacobi polynomial as follows: \[\label{eq:Jacobi-diff} \partial_t P_n^{\alpha, \beta}(t) = \frac{1}{2}(n+\alpha+\beta+1) P_{n-1}^{\alpha+1, \beta+1}(t).\tag{6}\] This index-shifting property is crucial for constructing Sobolev-orthogonal basis functions that are compatible with the mixed weak formulation involving second-order derivatives.

2.2 Spherical harmonics and ball polynomials↩︎

For any non-zero vector \({x} \in \mathbb{B}^3\), we adopt the polar decomposition \({x} = r{\xi}\), where \(r = \|{x}\|\) represents the radial distance and \({\xi} = \frac{x}{r} \in \mathbb{S}^2\) denotes the angular direction. The standard \(L^2\)-inner product on the unit sphere is denoted by \[\label{eq:SH-inner} (u, v)_{\mathbb{S}^2} = \int_{\mathbb{S}^2} u({\xi}) v({\xi}) \, \mathrm{d}\omega({\xi}),\tag{7}\] where \(\mathrm{d}\omega\) denotes the spherical surface measure.

We denote by \(\mathcal{H}_n(\mathbb{S}^2)\) the finite dimension space with \(n\)-order spherical harmonics, which is defined as the restriction of homogeneous harmonic polynomials of degree \(n\) to the unit sphere [12]. The dimension of this space is given by \[\label{eq:dim-Hn} d_n := \dim \mathcal{H}_n(\mathbb{S}^2) = 2n+1.\tag{8}\] For each \(n \ge 0\), we fix an orthogonal basis set \(\{Y_l^n\}_{l=1}^{2n+1}\) for \(\mathcal{H}_n(\mathbb{S}^2)\), which satisfies \[\label{eq:SH-orth} (Y_l^n, Y_{l'}^{n'})_{\mathbb{S}^2} = \delta_{nn'} \delta_{ll'}, \quad (n, l), (n', l') \in \Upsilon^3,\tag{9}\] where the index sets are defined as \[\Upsilon^3 = \left\{(n, l) : 1 \le l \le 2n+1, \, n \in \mathbb{N}_0 \right\}, \quad \Upsilon_N^3 = \left\{(n, l) \in \Upsilon^3 : n \le N \right\}.\]

Recalling the spherical coordinates \(x = (\sin\theta \cos\phi, \sin\theta \sin\phi, \cos\theta)^{\top}\) on the unit sphere \(\mathbb{S}^2\), i.e., \(\|x\|=1\), the orthonormal bases \(\{Y_l^n\}_{l=1}^{2n+1}\) can be explicitly constructed by Jacobi polynomials as follows [12]: \[\begin{align} Y_1^n(\theta, \phi) &= \sqrt{\frac{2n+1}{4\pi}} P_n^{0,0}(\cos\theta), \\ Y_{2j}^n(\theta, \phi) &= C_{n,j} (\sin\theta)^j \cos(j\phi) P_{n-j}^{j,j}(\cos\theta), \quad Y_{2j+1}^n(\theta, \phi) = C_{n,j} (\sin\theta)^j \sin(j\phi) P_{n-j}^{j,j}(\cos\theta) , \end{align}\] where \(1 \le j \le n\) and the normalization constant \[C_{n,j} = \frac{1}{2^{j}} \sqrt{\frac{(2n+1) \Gamma(n-j+1) \Gamma(n+j+1)}{2\pi \Gamma(n+1)^2}}.\]

For \(\alpha > -1\), we define the ball polynomials as \[\label{eq:ball95poly95def} B_{k,l}^{\alpha, n}(x) = \frac{(n+k+\frac{3}{2})_k}{(n+k+\frac{3}{2}+\alpha)_k} P_k^{\alpha, n+\frac{1}{2}}(2\|x\|^2 - 1) Y_l^n(x),\tag{10}\] where \(x \in \mathbb{B}^3\), \(n\) is the total degree of the spherical harmonic part, and \(k \in \mathbb{N}_0\) is the radial degree. Note that the total degree of \(B_{k,l}^{\alpha, n}(x)\) is \(n+2k\).

The ball polynomials are mutually orthogonal with respect to the weight function \(\omega_\alpha({x}) = (1 - \|x\|^2)^\alpha\) for \(\alpha > -1\) [13]: \[\label{eq:ball95orth} \int_{\mathbb{B}^3} B_{k,l}^{\alpha, n}(x) B_{\iota,\ell}^{\alpha, m}(x) \omega_{\alpha}(x) \, \mathrm{d}x = h_k^{\alpha, n} \delta_{k,\iota} \delta_{l,\ell} \delta_{n,m},\tag{11}\] where the normalization constant is given by \[h_k^{\alpha, n} = \frac{\Gamma(k+\alpha+1)\Gamma(n+2k+\frac{3}{2}) (n+k+\frac{3}{2})_k}{2 \Gamma(k+1) \Gamma(n+2k+\alpha+\frac{5}{2}) (n+k+\alpha+\frac{3}{2})_k}.\]

Specifically, we highlight two special cases \[h_k^{0,n} = \frac{1}{2n+4k+3}, \quad h_k^{1,n} = \frac{2(k+1)(2n+2k+3)}{(2n+4k+3)^2(2n+4k+5)}.\] The orthogonal ball polynomials \(\{B_{k,l}^{\alpha, n}(x)\}\) are the eigenfunctions of the following second-order differential operator: \[-(1-\|x\|^2)^{-\alpha}\nabla \cdot \left[(\mathrm{I}-xx^{\mathrm{t}})(1-\|x\|^2)^{\alpha}\nabla B_{k,l}^{\alpha, n}(x)\right] = (n+2k)(n+2k+2\alpha+3) B_{k,l}^{\alpha, n}(x),\] where \(\rm{I}\) denotes the identity operator.

2.3 Model Problem and Weak Formulation↩︎

In this section, we define the standard bilinear form \[\label{eq:bilinear95form} a(u,v) = \int_{\mathbb{B}^3} \nabla u \cdot \nabla v \, \mathrm{d}x.\tag{12}\]

By introducing an auxiliary variable \(\sigma = -\Delta u\), we propose the following mixed weak formulation of 1 : find \((\sigma, u) \in H^1_0(\mathbb{B}^3) \times H^1_0(\mathbb{B}^3)\) such that \[\label{eq:weak95mixed} \begin{cases} a(\sigma, v) = (f, v)_{L^2(\mathbb{B}^3)}, & \forall v \in H^1_0(\mathbb{B}^3), \\ a(u, \tau) = (\sigma, \tau)_{L^2(\mathbb{B}^3)}, & \forall \tau \in H^1_0(\mathbb{B}^3). \end{cases}\tag{13}\] This splitting framework is standard for the simply supported biharmonic equation [2], [9]. The well-posedness of 13 is guaranteed by the Lax-Milgram theorem applied sequentially to the two symmetric positive-definite subproblems. This relies on the continuity and coercivity of \(a(\cdot, \cdot)\) on \(H^1_0(\mathbb{B}^3)\), i.e., there hold: \[\begin{align} |a(u,v)| \le C \|u\|_{H^1(\mathbb{B}^3)} \|v\|_{H^1(\mathbb{B}^3)}, \quad a(v,v) \ge c \|v\|_{H^1(\mathbb{B}^3)}^2, \quad \forall u,v \in H^1_0(\mathbb{B}^3). \label{eq:coercivity} \end{align}\tag{14}\] For comprehensive details on these fundamental properties and the associated variational theory, we refer the reader to classical monographs such as [14].

3 Novel Spectral-Galerkin Approximation↩︎

To enforce the homogeneous boundary conditions, we employ generalized ball polynomials to construct an approximation space [12], [13]. The generalized ball polynomials with \(\alpha=-1\), which ensure the gradient orthogonality of the basis functions, as follows: \[B_{0,l}^{-1,n}(x) := B_{0,l}^{0,n}(x), \qquad B_{k,l}^{-1,n}(x) := B_{k,l}^{0,n}(x) - B_{k-1,l}^{0,n}(x), \quad k \in \mathbb{N},\quad (n,l)\in \Upsilon^3.\] Here, \(k \in \mathbb{N}\) represents the radial degree, which guarantees that the basis functions vanish on the boundary \(\mathbb{S}^2\).

3.1 Basis function and mixed formulation↩︎

Hence, we define finite-dimensional approximation spaces for our mixed spectral-Galerkin method. Let \(V_N \subset H^1_{0}(\mathbb{B}^3)\) be the space spanned by basis polynomials of total degree at most \(N\): \[\label{eq:VN} V_N := \text{span} \left\{ B_{k,l}^{-1, n}({x}) : 2 \le 2k+n \le N, \, k \in \mathbb{N}, \, (n, l) \in \Upsilon_N^3 \right\}.\tag{15}\]

A fundamental advantage of this approximation space is that it completely diagonalizes the stiffness matrix [13], whose entries are defined by 12 . In fact, the following gradient orthogonality relation holds: \[\label{grad-orth} a(B_{k,l}^{-1,n}, B_{k',l'}^{-1,n'}) = \left( \nabla B_{k,l}^{-1,n}, \nabla B_{k',l'}^{-1,n'} \right)_{L^2(\mathbb{B}^3)} = \lambda_{k,n} \, \delta_{k,k'} \, \delta_{n,n'} \, \delta_{l,l'},\tag{16}\] where \[\lambda_{k,n} = 2n + 4k + 1.\]

Using the approximation space defined in 15 , we now state the discrete mixed formulation of 13 : find \((\sigma_N, u_N) \in V_N \times V_N\) such that: \[\label{eq:discrete95problem} \begin{cases} a(\sigma_N,v_N) = (f,v_N)_{L^2(\mathbb{B}^3)}, & \forall v_N \in V_N,\\ a(u_N,\tau_N) = (\sigma_N,\tau_N)_{L^2(\mathbb{B}^3)}, & \forall \tau_N \in V_N. \end{cases}\tag{17}\] It is straightforward to conclude that the discrete system 17 is well-posed, which is guaranteed by the judicious choice of the approximation spaces. For more details, we refer the readers to [9], [14].

3.2 A-priori error estimates↩︎

In this subsection, we investigate the a-priori error estimates for \(u\) and \(\sigma\) in our proposed mixed spectral-Galerkin formulation 17 , respectively. To facilitate our analysis, we recall a standard orthogonal projection (i.e., the Ritz projection), \(\Pi_N: H^1_0(\mathbb{B}^3) \mapsto V_N\), satisfying \[\label{grad-dis} a(v - \Pi_N v, w_N) = 0, \quad \forall w_N \in V_N.\tag{18}\]

Subsequently, we derive a-priori error estimates in \(L^2\)- and \(H^1\)-norms by utilizing the aforementioned approximation properties, respectively.

4 Numerical Experiments↩︎

In this section, we present two numerical examples to validate the accuracy and convergence rates of the proposed mixed spectral-Galerkin method. Throughout this section, the maximum polynomial degree is denoted by \(N\). To quantify the exponential convergence, we calculate the experimental convergence rate: \[\text{Rate}_v=\frac{\ln \left( E_{i-1}(v) / E_i(v) \right)}{N_{i} - N_{i-1}} ,\] where \(E_i(v)\) denotes the corresponding \(L^2\)-error of \(v\) at \(N_i\). We employ \(r = \sqrt{x_1^2 + x_2^2 + x_3^2}\) to denote the radial distance from the origin.

4.1 Case 1: Analytic function↩︎

To verify the accuracy, we construct an exact solution: \[\label{eq:exact95u} u(r) = \begin{cases} \frac{\sin(\pi r)}{r}, & r > 0, \\ \pi, & r = 0. \end{cases}\tag{29}\] Note that this function is smooth. We evaluate the proposed mixed spectral-Galerkin method by measuring the numerical errors in both the \(L^2\)- and the \(H^1\)-norms. Table 1 presents the numerical results for a sequence of polynomial degrees \(N = 4, 8, 12, 16\).

The quantitative results in Table 1 clearly demonstrate spectral accuracy of our proposed scheme. As the polynomial degree \(N\) increases, the errors decay exponentially, which is consistent with the high smoothness of the exact solution. To further visualize this convergence behavior, Figure 3 presents the semi-logarithmic profiles of the \(L^2\)-errors for both \(u\) and \(\sigma\) comparing with a red dashed line segment equipped slope \(-1\). The strictly linear trajectory observed in these plots is a hallmark of spectral convergence. This confirms that our mixed formulation effectively captures the smooth features of the solution with high precision and stability.

4.2 Case 2: Transcendental function↩︎

To further assess the accuracy of the proposed method for smooth but non-polynomial functions. We construct the exact solution as: \[u(r) = \mathrm{e}^{r^2} - \frac{5\mathrm{e}}{3}r^2 + \frac{2\mathrm{e}}{3}.\] This solution is analytic but not a polynomial and exhibits monotonic behavior in the radial direction. Table 2 presents the convergence results for the primal and auxiliary variables. Consistent with the previous case, the method exhibits excellent spectral accuracy, with the \(L^2\)-error of \(u\) decaying rapidly from \(\mathcal{O}(10^{-2})\) at \(N=4\) to \(\mathcal{O}(10^{-11})\) at \(N=16\).

The exponential decay of the numerical errors is further illustrated in Figure 6, which is compared with a red dashed line segment equipped slope \(-1\). The distinct linear trajectories observed in the semi-logarithmic profiles for \(u\) and \(\sigma\) strongly confirm the spectral convergence. These results reinforce the high precision of the proposed mixed spectral-Galerkin method, even when applied to non-oscillatory transcendental solutions.

5 Conclusions↩︎

In this paper, we have proposed and analyzed a mixed spectral-Galerkin method for solving the biharmonic equation with simply supported boundary conditions on the unit ball. By introducing an auxiliary variable, the original fourth-order problem is successfully decoupled into a system of second-order equations. A key feature of our method is the basis construction using generalized ball polynomials, which yields a strictly diagonal stiffness matrix, thereby significantly reducing the computational cost. Furthermore, we rigorously proved the a-priori error estimates in both the \(L^2\)- and \(H^1\)-norms for the primal variable and auxiliary variable, respectively. The extensive numerical experiments highly align with our theoretical analysis, further confirming the spectral accuracy and efficiency of the proposed mixed method.

References↩︎