B-Spline Solutions of General Euler-Lagrange Equations

Abstract

The Euler-Lagrange equations are useful for solving optimization problems in mechanics. In this paper, we study the B-spline solutions of the Euler-Lagrange equations associated with the general functionals. The existing conditions of B-spline solutions to general Euler-Lagrange equations are given. As part of this work, we present a general method for generating B-spline solutions of the second- and fourth-order Euler-Lagrange equations. Furthermore, we show that some existing techniques for surface design, such as Coons patches, are exactly the special cases of the generalized Partial differential equations (PDE) surfaces with appropriate choices of the constants.

1. Introduction

Computer-aided geometric design (CAGD) is an area where the improvement of surface generation techniques is an ongoing demand. Partial differential equations (PDEs) were introduced as a valuable tool for geometric modeling, since they offer several features from which these areas can benefit [1,2,3]. The use of PDEs for shape design is conceptually different to conventional methods such as splines. The basic philosophy behind this method is that shape design is effectively treated as a mathematical boundary-value problem, i.e., shapes are produced by finding the solutions to a suitably chosen PDE that satisfies certain boundary conditions [4,5].

The importance of such a mathematical tool is that almost every physical phenomenon is modeled by a PDE. Now, PDEs have been introduced to many areas, such as computer graphics and animation and manufacturing designs, where they have been capable of solving a variety of problems efficiently [6]. It is well known that both the harmonic and biharmonic operators associated with Laplacian and bi-Laplacian equations are widely used in many application areas [7]. For example, the harmonic operator is widely used in physical problems, such as gravity, electromagnetism, and fluid flows. Similarly, the biharmonic operator is also associated with a variety of physical problems, such as tension in elastic membranes and the study of stress and strain in physical structures.

A Lagrangian functional defined on the space of smooth patches $\overrightarrow{\mathbf{x}}:\Omega \to {\mathbb{R}}^3,\Omega \subset {\mathbb{R}}^2$ is

$$\begin{array}{c}\hfill L(\overrightarrow{\mathbf{x}})=L(\overrightarrow{\mathbf{x}},{\overrightarrow{\mathbf{x}}}_u,{\overrightarrow{\mathbf{x}}}_v,{\overrightarrow{\mathbf{x}}}_{uu},{\overrightarrow{\mathbf{x}}}_{uv},{\overrightarrow{\mathbf{x}}}_{vv}).\end{array}$$

We take a functional to be

$$\begin{array}{c}\hfill I(\overrightarrow{\mathbf{x}})={\int }_{\Omega }L(\overrightarrow{\mathbf{x}})dudv.\end{array}$$

Minimizing the functional I is equivalent to requiring that the first variation of I is zero which then gives rise to the corresponding Euler-Lagrange equation [8]. For instance, the general quadratic functional [9] is

$$\begin{array}{c}\hfill L_{a,b,c}(\overrightarrow{\mathbf{x}})={\displaystyle \frac{1}{2}}{\int }_{\Omega }(a{\parallel {\overrightarrow{\mathbf{x}}}_u\parallel }^2+b\langle {\overrightarrow{\mathbf{x}}}_u,{\overrightarrow{\mathbf{x}}}_v\rangle +c{\parallel {\overrightarrow{\mathbf{x}}}_v\parallel }^2)dudv,\end{array}$$

(1)

where $a,b,c$ are constants and the Euler-Lagrange equation associated with the above general functional is

$$\begin{array}{c}\hfill a{\overrightarrow{\mathbf{x}}}_{uu}+b{\overrightarrow{\mathbf{x}}}_{uv}+c{\overrightarrow{\mathbf{x}}}_{vv}=0,\end{array}$$

which is a second-order partial differential equation, thereby a general form of the second-order PDEs is

$$\begin{array}{c}\hfill a{\overrightarrow{\mathbf{x}}}_{uu}+b{\overrightarrow{\mathbf{x}}}_{uv}+c{\overrightarrow{\mathbf{x}}}_{vv}+d{\overrightarrow{\mathbf{x}}}_u+e{\overrightarrow{\mathbf{x}}}_v+f\overrightarrow{\mathbf{x}}=g(u,v),\end{array}$$

where $a,b,c,d,e,f$ are constants and subscripts denote derivatives. The general Euler-Lagrange equations is an effective tool in geometrical modeling. Some of the existing techniques for surface design, such as the Coons patches [10] and harmonic surfaces [11], are the particular cases of the generalized PDE surfaces with appropriate choices of the constants in above equations.

In 1989, a method was presented for approximating surfaces which are the solutions of PDEs [2]. Monterde [11] focused on harmonic and biharmonic Bézier surfaces and presented a method of surface generation from prescribed boundaries based on the elliptic partial differential operators. In 2006, Monterde [12] also presented a method for generating Bézier surfaces from the boundary information based on the general 4th-order Euler-Lagrange equations which are a generalization of harmonic and biharmonic Bézier surfaces corresponding to Laplace and the standard biharmonic equations. A sixth-order PDE was presented by Zhang and You [13], which provided enough degrees of freedom not only to accommodate tangent, but also curvature boundary conditions and offered more shape control parameters to serve as user controls for the manipulation of surface shapes. Arnal [14] presented methods for designing triangular Bézier PDE surfaces by solving the second-order and fourth-order PDEs. Wang [15] presented a general 8th-order PDE method to generate Bézier surfaces from the boundary with positions and tangent vector information and extended Monterde’s results [11] to tri-harmonic Bézier surfaces. In 2015, Beltran [16] studied polynomial solutions in the Bézier form of the wave equation in dimensions one and two and determined which control points of the Bézier solution at two different times fix the solution. In various literature [17,18,19,20], they have indicated the ways in which, by a suitable choice of PDEs and boundary conditions imposed upon its solutions, surfaces satisfying a wide range of functional requirements can be created.

The work presented in this paper is geared to show the B-spline solutions of the general Euler-Lagrange equations. For this purpose, we study the existing conditions of B-spline solutions and then give a general method to compute the B-spline solutions of the general PDEs. Furthermore, we show that some of the existing techniques for surface design, such as the Coons patches and harmonic surfaces, are the special cases of the generalized PDE surfaces we present here.

The paper is organized as follows. In Section 2, we recall some notations and backgrounds about the Euler-Lagrange equations. In Section 3, we study the existing conditions of B-spline solutions to the Euler-Lagrange equations and the general method to compute the B-spline solutions are developed. As part of this work, the B-spline solutions to the 4th PDEs are discussed in Section 4. In Section 5, we indicate some of the existing methods for surface design are the special cases of the general PDE surfaces. The last Section concludes this paper.

2. Preliminary

For notational purposes, we assume that a parametric surface patch $\overrightarrow{\mathbf{x}}:\Omega \to {\mathbb{R}}^3,\Omega \subset {\mathbb{R}}^2$. We also represent the usual derivatives within the parametric domain such that ${\overrightarrow{\mathbf{x}}}_u=\frac{\partial \overrightarrow{\mathbf{x}}}{\partial u}$ and ${\overrightarrow{\mathbf{x}}}_v=\frac{\partial \overrightarrow{\mathbf{x}}}{\partial v}$. We impose the standard harmonic and biharmonic PDE operators respectively acting on these surfaces to be, $∆\overrightarrow{\mathbf{x}}=0$ and ${∆}^2\overrightarrow{\mathbf{x}}=0$ where

$$∆=(\frac{{\partial }^2}{\partial u^2}+\frac{{\partial }^2}{\partial v^2}),and{∆}^2={(\frac{{\partial }^2}{\partial u^2}+\frac{{\partial }^2}{\partial v^2})}^2.$$

We study functionals defined on the space of smooth patches $\overrightarrow{\mathbf{x}}:\Omega \to {\mathbb{R}}^3,\Omega \subset {\mathbb{R}}^2$. Given a Lagrangian functional as

$$\begin{array}{c}\hfill L(\overrightarrow{\mathbf{x}})=L(\overrightarrow{\mathbf{x}},{\overrightarrow{\mathbf{x}}}_u,{\overrightarrow{\mathbf{x}}}_v,{\overrightarrow{\mathbf{x}}}_{uu},{\overrightarrow{\mathbf{x}}}_{uv},{\overrightarrow{\mathbf{x}}}_{vv}),\end{array}$$

Equation (1) displayed in the Introduction is the general first-order quadratic form of the above functional. The general second-order quadratic form is

$$\begin{array}{c}\hfill L_{a,b,c,d,e,f}(\overrightarrow{\mathbf{x}})={\displaystyle \frac{1}{2}}{\int }_{\Omega }(a{\parallel {\overrightarrow{\mathbf{x}}}_{uu}\parallel }^2+b\langle {\overrightarrow{\mathbf{x}}}_{uu},{\overrightarrow{\mathbf{x}}}_{uv}\rangle +c\langle {\overrightarrow{\mathbf{x}}}_{uu},{\overrightarrow{\mathbf{x}}}_{vv}\rangle +d{\parallel {\overrightarrow{\mathbf{x}}}_{uv}\parallel }^2+e\langle {\overrightarrow{\mathbf{x}}}_{uv},{\overrightarrow{\mathbf{x}}}_{vv}\rangle +f{\parallel {\overrightarrow{\mathbf{x}}}_{vv}\parallel }^2)dudv,\end{array}$$

where $a,b,c,d,e,f\in \mathbb{R}$ are constants [12]. The associated Euler-Lagrange equation for the general quadratic functionals given above is

$$\begin{array}{c}\hfill a{\overrightarrow{\mathbf{x}}}_{uuuu}+b{\overrightarrow{\mathbf{x}}}_{uuuv}+(c+d){\overrightarrow{\mathbf{x}}}_{uuvv}+e{\overrightarrow{\mathbf{x}}}_{uvvv}+f{\overrightarrow{\mathbf{x}}}_{vvvv}=0.\end{array}$$

To simplify our formulation, let us consider the more general functionals which is divergence free [21], i.e.,

$$\begin{array}{c}\hfill L_{a,b,c,d,e}(\overrightarrow{\mathbf{x}})={\displaystyle \frac{1}{2}}{\int }_{\Omega }(a\parallel {\overrightarrow{\mathbf{x}}}_{uu}\parallel +b\langle {\overrightarrow{\mathbf{x}}}_{uu},{\overrightarrow{\mathbf{x}}}_{uv}\rangle +c{\parallel {\overrightarrow{\mathbf{x}}}_{uv}\parallel }^2+d\langle {\overrightarrow{\mathbf{x}}}_{uv},{\overrightarrow{\mathbf{x}}}_{vv}\rangle +e{\parallel {\overrightarrow{\mathbf{x}}}_{vv}\parallel }^2)dudv,\end{array}$$

where $a,b,c,d,e$ are constants, whose Euler-Lagrange equation [12] is

$$\begin{array}{c}\hfill a{\overrightarrow{\mathbf{x}}}_{uuuu}+b{\overrightarrow{\mathbf{x}}}_{uuuv}+c{\overrightarrow{\mathbf{x}}}_{uuvv}+d{\overrightarrow{\mathbf{x}}}_{uvvv}+e{\overrightarrow{\mathbf{x}}}_{vvvv}=0.\end{array}$$

3. B-Spline Solutions of the Second-Order PDEs

In this section, we give the existing conditions of B-spline solutions to the general second-order Euler-Lagrange equations and develop a method to get solutions. The general second-order Euler-Lagrange equations are

$$\begin{array}{c}\hfill a{\overrightarrow{\mathbf{x}}}_{uu}+b{\overrightarrow{\mathbf{x}}}_{uv}+c{\overrightarrow{\mathbf{x}}}_{vv}=0,\end{array}$$

(2)

where $a,b$ and c are constants.

Given a topological rectangular set of control points $P_{ij}\in {\mathbb{R}}^3$, $N_{i,p}(u)$ and $N_{i,q}(u)$, $i=0,1,\cdots ,m,j=0,1,\cdots ,n$, denote the B-spline basis functions of degree p and degree q defined on the knots vector $U=\{u_0,u_1,\cdots ,u_{m+p+1}\}$ and $V=\{v_0,v_1,\cdots ,v_{n+q+1}\}$ respectively [22], which can be obtained by the Cox-de Boor recursion formula [23], taking $N_{i,p}(u)$ for example,

$$\begin{array}{c}\hfill N_{i,p}(u)={\displaystyle \frac{u-u_i}{u_{i+p}-u_i}}{N}_{i,p-1}(u)+{\displaystyle \frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}}{N}_{i+1,p-1}(u).\end{array}$$

The associated B-spline surface of degree $p\times q$ is given by

$$\begin{array}{c}\hfill \overrightarrow{\mathbf{x}}(u,v)=\sum _{i=0}^m\sum _{j=0}^n{P}_{i,j}{N}_{i,p}(u)N_{j,q}(v),(u,v)\in [u_p,u_{m+1}]\times [v_q,v_{n+1}].\end{array}$$

The derivative of a B-spline basis function of degree p is simply a linear combination of two B-spline basis functions of degree $p-1$ [24].

$$\begin{array}{c}\hfill N_{i,p}^{\prime }={\displaystyle \frac{p}{u_{i+p}-u_i}}{N}_{i,p-1}(u)-{\displaystyle \frac{p}{u_{i+p+1}-u_{i+1}}}{N}_{i+1,p-1}(u).\end{array}$$

For the sake of simplicity, it is assumed that the knot vector is uniform spaced, namely the interval $[u_0,u_{m+p+1}]$ is partitioned into m intervals of spacing $h_1$. Similarly, the interval $[v_0,v_{n+q+1}]$ is partitioned into n intervals of spacing $h_2$. Thus, starting from Equation (2), substituting the B-spline surface $\overrightarrow{\mathbf{x}}(u,v)$ into the above PDE, we can obtained an equation with the control points $P_{i,j}$ to be determined.

$$\begin{array}{cc}\hfill L_{a,b,c}(\overrightarrow{\mathbf{x}})=& \frac{1}{h_1^2}\left(a\ast p(p-1)\sum _{i=0}^{m-2}\sum _{j=0}^n{∆}^{2,0}{P}_{i,j}{N}_{i,p-2}(u)N_{j,q}(v)\right)\hfill \\ & +\frac{1}{h_1{h}_2}\left(b\ast pq\sum _{i=0}^{m-1}\sum _{j=0}^{n-1}{∆}^{1,1}{P}_{i,j}{N}_{i,p-1}(u)N_{j,q-1}(v)\right)\hfill \\ & +\frac{1}{h_2^2}\left(c\ast pq(q-1)\sum _{i=0}^m\sum _{j=0}^{n-2}{∆}^{0,2}{P}_{i,j}{N}_{i,p}(u)N_{j,q-2}(v)\right).\hfill \end{array}$$

(3)

We can rewrite Equation (3) by the Cox-de Boor recursion formula as follows

$$\begin{array}{c}\hfill L_{a,b,c}(\overrightarrow{\mathbf{x}})=\sum _{i=-0}^{m+2}\sum _{j=0}^{n+2}{Q}_{i,j}{N}_{i,p-2}(u)N_{j,q-2}(v),\end{array}$$

where $Q_{i,j}$ is a linear combination of the control points obtained by the recursion formula.

The aim of the method is to find a B-spline solution

$$\begin{array}{c}\hfill \overrightarrow{\mathbf{x}}(u,v)=\sum _{i=0}^m\sum _{j=0}^n{P}_{i,j}{N}_{i,p}(u)N_{j,q}(v),\end{array}$$

where $\overrightarrow{\mathbf{x}}(u,v)$ is the solution of PDE (2), $P_{i,j}$ are control points to be determined and $N_{i,p}(u)$ and $N_{j,q}(v)$ are the tensor-product B-spline basis functions.

Since the degree elevation of B-spline basis functions is more complicated compared to the Bernstein basis functions, the expressions of $Q_{i,j}$ is complicated in the most cases. Let us look at it from another point of view. Obviously, $\overrightarrow{\mathbf{x}}(u,v)$ is a piecewise polynomial with total degree $p+q$. Therefore, the degree of $\overrightarrow{\mathbf{x}}{(u,v)}_{uu}$, $\overrightarrow{\mathbf{x}}{(u,v)}_{uv}$ or $\overrightarrow{\mathbf{x}}{(u,v)}_{vv}$ is $p+q-2$ with respect to u and v. Collecting the terms of each degree, there are finite terms, specifically, we have the following results.

Lemma 1.

In general, there are $pq+p+q-2$ terms about u and v in Equation (2).

Proof. 

For $\overrightarrow{\mathbf{x}}{(u,v)}_{uu}$, the largest degrees of u and v are $p-2$ and q respectively. Then there are $(p-2+1)(q+1)$ terms in $\overrightarrow{\mathbf{x}}{(u,v)}_{uu}$, displayed as follows,

terms	$v^q$	$v^{q-1}$	⋯	v	1
$u^{p-2}$	$u^{p-2}{v}^q$	$u^{p-2}{v}^{q-1}$	⋯	$u^{p-2}v$	$u^{p-2}$
$u^{p-3}$	$u^{p-3}{v}^q$	$u^{p-3}{v}^{q-1}$	⋯	$u^{p-2}v$	$u^{p-3}$
⋮	⋮	⋮	⋮	⋮	⋮
u	$u{v}^q$	$u{v}^{q-1}$	⋯	$uv$	u
1	$v^q$	$v^{q-1}$	⋯	v	1

Similarly, for $\overrightarrow{\mathbf{x}}{(u,v)}_{uv}$, the largest degree of u and v are $p-1$ and $q-1$ respectively, thus there are $(p-1+1)(q-1+1)$ terms in $\overrightarrow{\mathbf{x}}{(u,v)}_{uv}$, but some of them are the same as the terms in $\overrightarrow{\mathbf{x}}{(u,v)}_{uu}$. Actually, there are only q terms are new for Equation (2), to be specific, which are $u^{p-1}{v}^{q-1},u^{p-1}{v}^{q-2},\cdots ,u^{p-1}v,u^{p-1}$.

In the same way, $\overrightarrow{\mathbf{x}}{(u,v)}_{vv}$ totally has $(p+1)(q-2+1)$ terms, in which $q-1$ terms are new for (2), specifically, $u^p{v}^{q-2},\cdots ,u^{p}v,u^p$. The above analysis suggests that there are $(p-2+1)(q+1)+q+q-1$ terms in Equation (2), and then the lemma holds. □

Theorem 1.

In general, while $n\ge p$ and $m\ge q$, PDE (2) has B-spline solutions.

Proof. 

In fact, as $u\in [u_p,u_{n+1}]$, $v\in [v_q,v_{m+1}]$, $a{\overrightarrow{\mathbf{x}}}_{uu}+b{\overrightarrow{\mathbf{x}}}_{uv}+c{\overrightarrow{\mathbf{x}}}_{vv}=0$ holds if and only if the coefficient of each term about u and v vanishes. By Lemma 1, there are $pq+p+q-2$ terms about u and v in Equation (2), which implies that we get $pq+p+q-2$ linear equations about the control points. B-spline surface has $(n+1)(m+1)$ control points to be determined, while $n\ge p$ and $m\ge q$, $(n+1)(m+1)>pq+p+q-2$. This suggests that the $pq+p+q-2$ linear equations are solvable, that is to say, the PDE has B-spline solutions. □

From the previous discussion, we also get a general method to solve the PDEs Equation (2) with B-spline solutions. Actually, as the coefficients of terms about u and v vanish, $pq+p+q-2$ linear equations of the control points are obtained. While $n\ge p$ and $m\ge q$, the control points to be determined are redundancy, this means we can impose some boundary conditions on the equation. For instance, the values of B-spline surface $\overrightarrow{\mathbf{x}}(u,v)$ on the boundary and its normal derivative $\overrightarrow{\mathbf{x}}{(u,v)}_u$ and $\overrightarrow{\mathbf{x}}{(u,v)}_v$ at the corner points are specified. Thus, we conclude the method for computing the B-spline solutions of Equation (2) can be described as Algorithm 1 and explain the method in Example 1.

Algorithm 1 Algorithm to get the B-spline solution of Euler-Lagrange Equation (2)

Input:  Euler-Lagrange Equation (2) and B-spline surfaces $\overrightarrow{\mathbf{x}}(u,v)$ with the control points to be determined; Output:  The B-spline solutions $\overrightarrow{\mathbf{x}}(u,v)$ of Equation (2). Step 1: Substituting the B-spline surface $\overrightarrow{\mathbf{x}}(u,v)$ into Equation (2), collect the terms about u and v. Step 2: Collecting the $pq+p+q-2$ coefficients of each terms of equations, give $pq+p+q-2$ linear equations of the control points $\{P_{i,j},0\le i\le n,0\le j\le m\}$. Step 3: Combining with the given boundary conditions or control points and solving the linear equations, the unknown control points are obtained. Step 4: Construction of B-spline surface $\overrightarrow{\mathbf{x}}(u,v)$ by the control points.

Example 1.

For our purposes, we require that the basis functions be not only B-splines, but also at least two-times differentiable. This implies that the degree of B-spline surface is at least three. Thus, as it is not desirable to add any more complication than is strictly necessary, cubic B-splines will be used as the basis functions throughout the discussion. In the interest of simplicity, assuming that the B-splines are uniformly spaced, i.e., the interval $[u_0,u_{m+p+1}]$ is partitioned into m intervals of spacing h. Then at each point $u_i$ is centered a piecewise cubic B-spline $N_{i,3}(u)$ as

$$N_{i,3}(u,v)=\left\{\begin{array}{cc}0,\hfill & u<u_{i-2};\hfill \\ 1/6\frac{{\left(u-u_i+2h\right)}^3}{h^3},\hfill & u\in [u_{i-2},u_{i-1}];\hfill \\ 1/6\frac{4h^3-6h{u}^2+12hu{u}_i-6h{u_i}^2-3u^3+9u^2{u}_i-9u{u_i}^2+3{u_i}^3}{h^3},\hfill & u\in [u_{i-1},u_i];\hfill \\ 1/6\frac{4h^3-6h{u}^2+12hu{u}_i-6h{u_i}^2+3u^3-9u^2{u}_i+9u{u_i}^2-3{u_i}^3}{h^3},\hfill & u\in [u_i,u_{i+1}];\hfill \\ 1/6\frac{8h^3-12h^{2}u+12h^2{u}_i+6h{u}^2-12hu{u}_i+6h{u_i}^2-u^3+3u^2{u}_i-3u{u_i}^2+{u_i}^3}{h^3},\hfill & u\in [u_{i+1},u_{i+2}];\hfill \\ 0,\hfill & u>u_{i+2}.\hfill \end{array}\right.$$

$N_{i,3}(v)$ can be given in the same way. By the Theorem 1, in order to get the B-spline solutions of Equation (2), the control points set satisfies $m\ge 3$ and $n\ge 3$. Assume $m=3,n=3$, the B-spline surface to be determined can be represented by

$$\begin{array}{c}\hfill \overrightarrow{\mathbf{x}}(u,v)=\sum _{i=0}^3\sum _{j=0}^3{P}_{i,j}{N}_{i,3}(u)N_{j,3}(v),u\in [u_3,u_4],v\in [v_3,v_4],\end{array}$$

where $P_{i,j}$ are the unknown control points.

For simplicity, assuming $a=1,b=-2,c=1$ in Equation (2) and the space of interval $h=1$. Substituting the B-spline surface $\overrightarrow{\mathbf{x}}(u,v)$ into Equation (2), and collecting the terms about u and v, we get 13 terms in the equation, specifically, $u^{3}v,u^3,u^2{v}^2,u^{2}v,u^{2}v,u{v}^3,u{v}^2,uv,u,v^3,v^2,v,1.$ As $u\in [u_3,u_4]$, the PDE (2) holds if and only if the coefficients of these terms vanish, which implies the following 13 equations about control points holds.

$$\left\{\begin{array}{cc}\hfill P_{0,0}-3P_{0,1}+3P_{0,2}& -P_{0,3}-3P_{1,0}+9P_{1,1}-9P_{1,2}+3P_{1,3}+3P_{2,0}-9P_{2,1}+9P_{2,2}-3P_{2,3}-P_{3,0}\hfill \\ & +3P_{3,1}-3P_{3,2}+P_{3,3}=0;\hfill \\ \hfill -2P_{0,0}+6P_{0,1}-6P_{0,2}& +2P_{0,3}+6P_{1,0}-18P_{1,1}+18P_{1,2}-6P_{1,3}-6P_{2,0}+18P_{2,1}-18P_{2,2}+6P_{2,3}+2P_{3,0}\hfill \\ & -6P_{3,1}+6P_{3,2}-2P_{3,3}=0;\hfill \\ \hfill P_{0,0}-3P_{0,1}+3P_{0,2}& -P_{0,3}-3P_{1,0}+9P_{1,1}-9P_{1,2}+3P_{1,3}+3P_{2,0}-9P_{2,1}+9P_{2,2}-3P_{2,3}-P_{3,0}\hfill \\ & +3P_{3,1}-3P_{3,2}+P_{3,3}=0;\hfill \\ \hfill -P_{0,0}+7P_{0,1}-11P_{0,2}& +5P_{0,3}+3P_{1,0}-21P_{1,1}+33P_{1,2}-15P_{1,3}-3P_{2,0}+21P_{2,1}-33P_{2,2}+15P_{2,3}+P_{3,0}\hfill \\ & -7P_{3,1}+11P_{3,2}-5P_{3,3}=0;\hfill \\ \hfill 2P_{0,0}-14P_{0,1}+22P_{0,2}& -10P_{0,3}-2P_{1,0}+30P_{1,1}-54P_{1,2}+26P_{1,3}-2P_{2,0}-18P_{2,1}+42P_{2,2}-22P_{2,3}+2P_{3,0}\hfill \\ & +2P_{3,1}-10P_{3,2}+6P_{3,3}=0;\hfill \\ \hfill -P_{0,0}+7P_{0,1}-11P_{0,2}& +5P_{0,3}-5P_{1,0}+3P_{1,1}+9P_{1,2}-7P_{1,3}+13P_{2,0}-27P_{2,1}+15P_{2,2}-P_{2,3}-7P_{3,0}\hfill \\ & +17P_{3,1}-13P_{3,2}+3P_{3,3}=0;\hfill \\ \hfill 4P_{1,0}-12P_{1,1}+12P_{1,2}& -4P_{1,3}-8P_{2,0}+24P_{2,1}-24P_{2,2}+8P_{2,3}+4P_{3,0}\hfill \\ & -12P_{3,1}+12P_{3,2}-4P_{3,3}=0;\hfill \\ \hfill -P_{0,0}+15P_{0,1}-39P_{0,2}& +25P_{0,3}-P_{1,0}-17P_{1,1}+73P_{1,2}-55P_{1,3}+5P_{2,0}-11P_{2,1}-29P_{2,2}+35P_{2,3}-3P_{3,0}\hfill \\ & +13P_{3,1}-5P_{3,2}-5P_{3,3}=0;\hfill \\ \hfill 9/2P_{1,0}+\frac{119P_{2,1}}{2}-\frac{175P_{2,2}}{2}& +\frac{65P_{2,3}}{2}-9/2P_{2,0}-\frac{49P_{3,1}}{2}+\frac{89P_{3,2}}{2}-\frac{39P_{3,3}}{2}-1/2P_{3,0}+1/2P_{0,0}-\frac{55P_{1,1}}{2}-15/2P_{0,1}\hfill \\ & +\frac{47P_{1,2}}{2}+\frac{39P_{0,2}}{2}-1/2P_{1,3}-\frac{25}{2}{P}_{0,3}=0;\hfill \\ \hfill 16P_{1,1}-32P_{1,2}+16P_{1,3}& -12P_{2,0}+4P_{2,1}+28P_{2,2}-20P_{2,3}+12P_{3,0}-20P_{3,1}+4P_{3,2}+4P_{3,3}=0;\hfill \\ \hfill -3/2P_{1,0}-\frac{69P_{2,1}}{2}+\frac{273P_{2,2}}{2}& -\frac{175P_{2,3}}{2}-5/2P_{2,0}-7/6P_{3,1}-\frac{277P_{3,2}}{6}+\frac{235P_{3,3}}{6}+\frac{25P_{3,0}}{6}-1/6P_{0,0}+\frac{61P_{1,1}}{2}+\frac{31P_{0,1}}{6}\hfill \\ & -\frac{137P_{1,2}}{2}-\frac{131P_{0,2}}{6}+\frac{55P_{1,3}}{2}+\frac{125P_{0,3}}{6}=0;\hfill \\ \hfill -4/3P_{1,0}-58P_{2,1}+42P_{2,2}& +\frac{22P_{2,3}}{3}+\frac{26P_{2,0}}{3}+50P_{3,1}-66P_{3,2}+\frac{58P_{3,3}}{3}-10/3P_{3,0}\hfill \\ & -4P_{1,1}+36P_{1,2}-\frac{92P_{1,3}}{3}=0;\hfill \\ \hfill 2/3P_{1,0}-14/3P_{1,1}-\frac{50P_{1,2}}{3}& +\frac{110P_{1,3}}{3}-4/3P_{2,0}+\frac{172P_{2,1}}{3}-\frac{404P_{2,2}}{3}+\frac{140P_{2,3}}{3}-10/3P_{3,0}\hfill \\ & -\frac{74P_{3,1}}{3}+\frac{322P_{3,2}}{3}-\frac{190P_{3,3}}{3}=0.\hfill \end{array}\right.$$

In fact, since the coefficients of $u{v}^3$ and $u^{3}v$ are identical, the above linear system has 12 linear equations about the control points. Since $m=3,n=3$, there are 16 control points to be determined which implies there are 4 control points are free. We impose the additional values of B-spline surface or the control points on the boundary, and then we can obtain the unknown control points unique by solving the above linear equations. Supposing

$$\begin{array}{c}\hfill P_{0,0}=[1,1,1];P_{0,1}=[1,2,-8];P_{0,2}=[1,3,8];P_{0,3}=[1,4,1],\end{array}$$

by solving the above equations, the remaining control points are given as

$$\begin{array}{c}{P}_{1,0}=[2,1,8];P_{1,1}=[2,2,-5/9];P_{1,2}=[2,3,119/9];P_{1,3}=[2,4,4/3];P_{2,0}=[3,1,-8];\hfill \\ P_{2,1}=[3,2,-121/9];P_{2,2}=[3,3,7/9];P_{2,3}=[3,4,-40/3];P_{3,0}=[4,1,1];P_{3,1}=[4,2,4/3];\hfill \\ P_{3,2}=[4,3,56/3];P_{3,3}=[4,4,5].\hfill \end{array}$$

And the corresponding B-spline solutions to the PDE (2) is displayed in the left of Figure 1, whereas the other two figures are the B-spline solutions with the different given boundary control points.

4. B-Spline Solutions of the 4th Order PDEs

In this section, we represent the solutions of the following fourth-order PDE with B-spline surfaces, which can be discussed using the proposed method.

$$\begin{array}{c}\hfill a{\overrightarrow{\mathbf{x}}}_{uuuu}+b{\overrightarrow{\mathbf{x}}}_{uuuv}+c{\overrightarrow{\mathbf{x}}}_{uuvv}+d{\overrightarrow{\mathbf{x}}}_{uvvv}+e{\overrightarrow{\mathbf{x}}}_{vvvv}=0.\end{array}$$

(4)

To get the B-spline solutions of the above PDE, we can establish the algorithm in the similar fashion with the two-order PDEs in Section 3. Substituting the B-spline surface $\overrightarrow{\mathbf{x}}(u,v)$ into the above PDE, an equation with the control points $P_{i,j}$ can be obtained. Obviously, $\overrightarrow{\mathbf{x}}(u,v)$ is a piecewise polynomial surface with total degree $p+q$. Therefore, the total degree of $\overrightarrow{\mathbf{x}}{(u,v)}_{uuuu}$, $\overrightarrow{\mathbf{x}}{(u,v)}_{uuuv}$ or $\overrightarrow{\mathbf{x}}{(u,v)}_{uuvv}$ is $p+q-4$ with respect to u and v. We give the following results analogous to the two-order PDEs.

Lemma 2.

In general, there are $pq+p+q-9$ terms about u and v in Equation (4).

Proof. 

For $\overrightarrow{\mathbf{x}}{(u,v)}_{uuuu}$, the largest degrees of u and v are $p-4$ and q respectively. Then there is $(p-4+1)(q+1)$ terms in $\overrightarrow{\mathbf{x}}{(u,v)}_{uuuu}$, displayed as follows,

terms	$v^q$	$v^{q-1}$	⋯	v	1
$u^{p-4}$	$u^{p-4}{v}^q$	$u^{p-4}{v}^{q-1}$	⋯	$u^{p-4}v$	$u^{p-4}$
$u^{p-5}$	$u^{p-5}{v}^q$	$u^{p-5}{v}^{q-1}$	⋯	$u^{p-5}v$	$u^{p-5}$
⋮	⋮	⋮	⋮	⋮	⋮
u	$u{v}^q$	$u{v}^{q-1}$	⋯	$uv$	u
1	$v^q$	$v^{q-1}$	⋯	v	1

Similarly, for $\overrightarrow{\mathbf{x}}{(u,v)}_{uuuv}$, the largest degrees of u and v is $p-3$ and $q-1$ respectively, thus there are $(p-3+1)(q-1+1)$ terms, but some of them are the same as $\overrightarrow{\mathbf{x}}{(u,v)}_{uuuu}$. Actually, there are only q terms new for Equation (4), specifically, which are $u^{p-3}{v}^{q-1},u^{p-3}{v}^{q-2},\cdots ,u^{p-3}v,u^{p-3}$.

In the same way, $\overrightarrow{\mathbf{x}}{(u,v)}_{uuvv}$ totally has $(p-3+1)q$ terms, in which only $q-1$ terms are new, to be specific, $u^{p-2}{v}^{q-1},u^{p-2}{v}^{q-2},\cdots ,u^{p-2}v,u^{p-2}$. There are $q-2$ and $q-3$ new terms in $\overrightarrow{\mathbf{x}}{(u,v)}_{uvvv}$ and $\overrightarrow{\mathbf{x}}{(u,v)}_{vvvv}$, respectively. The above analysis implies that there are $(p-4+1)(q+1)+q+q-1+q-2+q-3=pq+p+q-9$ terms in Equation (4), and then the lemma holds. □

By the proposed lemma and analysis, the following result can be given, and the proof is similar with Theorem 1, which will be omitted here.

Theorem 2.

In general, while $m\ge p$ and $n\ge q$, PDE (4) has B-spline solutions.

Thus, in similar with the method to get the B-spline solutions of the second-order PDE in Section 3, we conclude the method for computing the B-spline solutions of Equation (4) as Algorithm 2 and explain the method in Example 2.

Algorithm 2 Algorithm to get the B-spline solution of the Euler-Lagrange Equation (4)

Input:  B-spline surfaces $\overrightarrow{\mathbf{x}}(u,v)$ with the control points to be determined and the general Euler-Lagrange Equation (4); Output:  The B-spline solutions $\overrightarrow{\mathbf{x}}(u,v)$ of Equation (4). Step 1: Substituting the B-spline surface $\overrightarrow{\mathbf{x}}(u,v)$ into Equation (4), collect the terms about u and v. Step 2: Collecting the $pq+p+q-9$ coefficients of each terms of equations, give $pq+p+q-9$ linear equations of the control points $\{P_{i,j},0\le i\le n,0\le j\le m\}$. Step 3: Combining with the given boundary conditions or control points and solving the linear equations, the unknown control points are obtained. Step 4: Construction of B-spline surface $\overrightarrow{\mathbf{x}}(u,v)$ by the control points.

Example 2.

In a similar way, the fourth-order PDE implies that the B-spline solution is at least four-times differentiable. Therefore, we require the degree of B-spline surface is at least quintic. As our aim is that the solution be obtained in the form of a tensor-product B-spline surface, and then $N_{i,p}(u)$ is chosen to be quintic B-spline basis functions which span the space of piecewise quintic polynomials. Assuming that the knot vector of B-splines is uniformly spaced, and the space of the intervals is h. In the interest of simplicity, suppose $h=1$. Then at each point $u_i$ is centered a piecewise quintic B-spline $N_{i,5}(u)$ as

$$N_{i,5}(u)=\left\{\begin{array}{c}0,u<u_i-3;\hfill \\ \frac{{\left(u-u_i+3\right)}^5}{120},u\in [u_i-3,u_i-2];\hfill \\ -1/20{\left(u-u_i+2\right)}^5+\frac{{\left(u-u_i+3\right)}^5}{120},u\in [u_i-2,u_i-1];\hfill \\ 1/8{\left(u-u_i+1\right)}^5-1/20{\left(u-u_i+2\right)}^5+\frac{{\left(u-u_i+3\right)}^5}{120},u\in [u_i-1,u_i];\hfill \\ -1/6{\left(u-u_i\right)}^5+1/8{\left(u-u_i+1\right)}^5-1/20{\left(u-u_i+2\right)}^5+\frac{{\left(u-u_i+3\right)}^5}{120},u\in [u_i,u_i+1];\hfill \\ 1/8{\left(u-u_i-1\right)}^5-1/6{\left(u-u_i\right)}^5+1/8{\left(u-u_i+1\right)}^5-1/20{\left(u-u_i+2\right)}^5+\frac{{\left(u-u_i+3\right)}^5}{120},u\in [u_i+1,u_i+2];\hfill \\ -1/20{\left(u-u_i-2\right)}^5+1/8{\left(u-u_i-1\right)}^5-1/6{\left(u-u_i\right)}^5+1/8{\left(u-u_i+1\right)}^5-1/20{\left(u-u_i+2\right)}^5+\frac{{\left(u-u_i+3\right)}^5}{120},\hfill \\ \hfill u\in [u_i+2,u_i+3];\\ 0,u>u_{i+3}.\hfill \end{array}\right.$$

$N_{i,5}(v)$ can be given in the similar fashion. By the Theorem 2, in order to get the B-spline solutions, the control points satisfy $m\ge p$ and $n\ge q$. Without loss of generality, for $m=5,n=5$, the B-spline surface to be determined can be represented by

$$\begin{array}{c}\hfill \overrightarrow{\mathbf{x}}(u,v)=\sum _{i=0}^5\sum _{j=0}^5{P}_{i,j}{N}_{i,5}(u)N_{j,5}(v),u\in [u_5,u_6],\end{array}$$

where $P_{i,j}$ control points to be determined.

For simplicity, assuming $a=1,b=2,c=1$ in Equation (4) and substituting the above B-spline $\overrightarrow{\mathbf{x}}(u,v)$ into Equation (4). Collecting the terms about u and v, there are 26 terms in the equation, specifically,

$$\begin{array}{c}{u}^{5}v,u^5,u^4{v}^2,u^{4}v,u^4,u^3{v}^3,u^3{v}^2,u^{3}v,u^3,u^2{v}^4,u^2{v}^3,u^2{v}^2,\hfill \\ u^{2}v,u^2,u{v}^5,u{v}^4,u{v}^3,u{v}^2,uv,u,v^5,v^4,v^3,v^2,v,1.\hfill \end{array}$$

As $u\in [u_5,u_6]$, the PDE (4) holds if and only if the coefficients of the these terms vanish, that is to say, 26 equations of the control points are obtained. We omit the tedious and laborious details here. While $m=5$ and $n=5$, the number of control points to be determined is 36, which means there are 10 control points are freely chosen, so we can impose some additional boundary conditions for the equations. For instance, the values of control points on the boundary and the normal derivatives of B-spline $\overrightarrow{\mathbf{x}}{(u,v)}_u$ and $\overrightarrow{\mathbf{x}}{(u,v)}_v$ at the corner points are specified.

We set the given boundary control points are

$$\begin{array}{c}{P}_{0,0}=[1,1,0];P_{0,1}=[1,2,10];P_{0,2}=[1,3,10];P_{0,3}=[1,4,-10];P_{0,4}=[1,5,-10];P_{0,5}=[1,6,0]\hfill \\ P_{1,0}=[2,1,-10];P_{2,0}=[3,1,-10];P_{3,0}=[4,1,10];P_{4,0}=[5,1,10];P_{5,0}=[6,1,0].\hfill \end{array}$$

And solving the linear equations of the control points, we can get the remaining control points as

$$\begin{array}{c}{P}_{1,2}=[2,3,-1];P_{1,3}=[2,4,-27];P_{1,4}=[2,5,-106/3];P_{1,5}=[2,6,-80/3]\hfill \\ P_{2,1}=[3,2,1];P_{2,2}=[3,3,0];P_{2,3}=[3,4,-33];P_{2,4}=[3,5,-164/3];P_{2,5}=[3,6,-175/3]\hfill \\ P_{3,1}=[4,2,27];P_{3,2}=[4,3,33];P_{3,3}=[4,4,0];P_{3,4}=[4,5,-30];P_{3,5}=[4,6,-45]\hfill \\ P_{4,1}=[5,2,106/3];P_{4,2}=[5,3,164/3];P_{4,3}=[5,4,30];P_{4,4}=[5,5,0];P_{4,5}=[5,6,-20]\hfill \\ P_{5,1}=[6,2,80/3];P_{5,2}=[6,3,175/3];P_{5,3}=[6,4,45];P_{5,4}=[6,5,20];P_{5,5}=[6,6,0].\hfill \end{array}$$

Then the corresponding B-spline solution to the PDE is constructed and showed in the left of Figure 2. The other two figures show the B-spline solutions with the different given boundary control points.

5. Applications in Geometrical Modeling

As mentioned previously, the PDE surfaces arising from the functional serve as a useful tool in geometrical modeling. In this section, we discuss some existing surfaces, such as Coons surfaces and biharmonic surface, are exactly the special cases of the general PDE surfaces.

While dealing with biharmonic surfaces [25], we often study the functionals associated with the standard biharmonic operators ${∆}^2$,

$$\begin{array}{c}\hfill F_{bihar}={\displaystyle \frac{1}{2}}{\int }_{\Omega }({\overrightarrow{\mathbf{x}}}_{uu}^2+2{\overrightarrow{\mathbf{x}}}_{uv}^2+{\overrightarrow{\mathbf{x}}}_{vv}^2)dudv.\end{array}$$

As the definition of the biharmonic operator itself suggests a kind of control on the shape of the resulting surface in some sense.

Other similar functionals used in CAGD are also the special cases of the previous general functional. Among them we cite the functional related to Coons patches which can be used to show that the resulting surface minimizing the surface twist [10]

$$\begin{array}{c}\hfill F_{Coons}={\displaystyle \frac{1}{2}}{\int }_{\Omega }{\parallel {\overrightarrow{\mathbf{x}}}_{uv}\parallel }^{2}dudv.\end{array}$$

The modified biharmonic functional is [26]

$$\begin{array}{c}\hfill F_a(\overrightarrow{\mathbf{x}})={\displaystyle \frac{1}{2}}{\int }_{\Omega }({\parallel {\overrightarrow{\mathbf{x}}}_{uu}\parallel }^2+2a^2{\parallel {\overrightarrow{\mathbf{x}}}_{uv}\parallel }^2+a^4{\parallel {\overrightarrow{\mathbf{x}}}_{vv}\parallel }^2)dudv.\end{array}$$

A more general form is

$$\begin{array}{c}\hfill F_{a,b,c}(\overrightarrow{\mathbf{x}})={\displaystyle \frac{1}{2}}{\int }_{\Omega }(a^2{\parallel {\overrightarrow{\mathbf{x}}}_{uu}\parallel }^2+2a^2{b}^2\langle {\overrightarrow{\mathbf{x}}}_{uv},{\overrightarrow{\mathbf{x}}}_{uv}\rangle +(b^2+2a^2{c}^2){\parallel {\overrightarrow{\mathbf{x}}}_{uv}\parallel }^2+2b^2{c}^2\langle {\overrightarrow{\mathbf{x}}}_{uv},{\overrightarrow{\mathbf{x}}}_{uv}\rangle +c^2{\parallel {\overrightarrow{\mathbf{x}}}_{vv}\parallel }^2)dudv,\end{array}$$

where $a,b$ and c are constants. This functional is applied to minimize the resulting surface area where more weights are assigned along one particular parametric coordinate direction.

The examples of functionals displayed above are some of the typical functionals commonly used in CAGD. It is noteworthy that they are the special cases of the most general functionals,

$$\begin{array}{c}\hfill L_{a,b,c,d,e}(\overrightarrow{\mathbf{x}})={\displaystyle \frac{1}{2}}{\int }_{\Omega }(a{\parallel {\overrightarrow{\mathbf{x}}}_{uu}\parallel }^2+b\langle {\overrightarrow{\mathbf{x}}}_{uu},{\overrightarrow{\mathbf{x}}}_{uv}\rangle +c{\parallel {\overrightarrow{\mathbf{x}}}_{uv}\parallel }^2+d\langle {\overrightarrow{\mathbf{x}}}_{uv},{\overrightarrow{\mathbf{x}}}_{vv}\rangle +e{\parallel {\overrightarrow{\mathbf{x}}}_{vv}\parallel }^2)dudv,\end{array}$$

where $a,b,c,d,e\in \mathbb{R}$ are constants. The associated Euler-Lagrange equation for the general functional given above is

$$\begin{array}{c}\hfill a{\overrightarrow{\mathbf{x}}}_{uuuu}+b{\overrightarrow{\mathbf{x}}}_{uuuv}+c{\overrightarrow{\mathbf{x}}}_{uuvv}+d{\overrightarrow{\mathbf{x}}}_{uvvv}+e{\overrightarrow{\mathbf{x}}}_{vvvv}=0.\end{array}$$

Please note that above surface design techniques, such as the biharmonic surface and Coons surface, are exactly the special cases of the generalized framework with appropriate choices of the constants.

$$\begin{array}{c}\hfill F_{bihar}=L_{1,0,2,0,1},F_{Coons}=L_{0,0,1,0,0},F_a=L_{1,0,2a^2,0,a^4}.\end{array}$$

6. Conclusions

PDEs have been introduced to many areas such as computer graphics, CAGD, and animation since they can solve a variety of problems in these areas. In this paper, we consider the B-spline solutions of the general Euler-Lagrange equations. It is well known that both the harmonic and biharmonic operators associated with Laplacian and bi-Laplacian equations are widely used in many application areas. We use a linear combination of these common operators and get a wide variety of PDEs which can be used in CAGD. The existing conditions of B-spline solutions to a general PDE are given. As part of this work, we present a general method for generating B-spline solutions to the second- and fourth-order Euler-Lagrange equations.

There are several extensions of this work. It is interesting to study how the various coefficients associated with the Euler-Lagrange equations affect the shape of the resulting surfaces. We could also carry out a study of the geometric properties of surfaces generated by the PDEs which will be useful in the developing interactive tools in future.