This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

On the Uniqueness of Convex Central Configurations
in the Planar 44-Body Problem

Shanzhong Sun
Department of Mathematics
and
Academy for Multidisciplinary Studies
Capital Normal University, Beijing 100048, P. R. China
[email protected]
Zhifu Xie
School of Mathematics and Natural Science
The University of Southern Mississippi
Hattiesburg, MS 39406,USA
[email protected]
Peng You
School of Mathematics and Statistics
Hebei University of Economics and Business
Shijiazhuang Hebei 050061, P. R. China
[email protected]
Abstract

In this paper, we provide a rigorous computer-assisted proof (CAP) of the conjecture that in the planar four-body problem there exists a unique convex central configuration for any four fixed positive masses in a given order belonging to a closed domain in the mass space. The proof employs the Krawczyk operator and the implicit function theorem (IFT). Notably, we demonstrate that the implicit function theorem can be combined with interval analysis, enabling us to estimate the size of the region where the implicit function exists and extend our findings from one mass point to its neighborhood.

Dédier à Alain Chenciner avec admiration et amitié

Key Words: Central Configuration, Convex Central Configuration, Uniqueness, NN-Body Problem, Krawczyk Operator, Implicit Function Theorem.
Mathematics Subject Classification: 70F10, 70F15

1 Introduction

The planar Newtonian NN-body problem describes the motion of NN celestial bodies considered as point particles in the plane under the universal gravitational law, which is governed by a system of second order ordinary differential equations:

miqi¨=j=1,jiNGmimj(qjqi)qjqi3,i=1,2,..,N.m_{i}\ddot{q_{i}}=\sum_{j=1,j\neq i}^{N}\frac{Gm_{i}m_{j}(q_{j}-q_{i})}{\left\|q_{j}-q_{i}\right\|^{3}},\hskip 14.22636pti=1,2,..,N.

Hereafter the gravitational constant GG is set to 1 for convenience and qj2q_{j}\in\mathbb{R}^{2} is the position vector of the jj-th point mass mjm_{j} in the planar inertial barycentric frame. rij=qiqjr_{ij}=\|q_{i}-q_{j}\| is the Euclidean distance between qiq_{i} and qjq_{j}. Let

(1) C=m1q1++mNqNC=m_{1}q_{1}+\cdots+m_{N}q_{N}

be the linear moment and M=m1+m2++mNM=m_{1}+m_{2}+\cdots+m_{N} be the total mass with c=C/Mc=C/M the center of mass. Given a mass vector m=(m1,m2,,mN)m=(m_{1},m_{2},\cdots,m_{N}), a configuration q=(q1,q2,,qN)q=(q_{1},q_{2},\cdots,q_{N}) is called a central configuration for mm if there exists a positive constant λ\lambda such that

(2) j=1,jiNmimj(qjqi)rij3=λmi(qic).\sum_{j=1,j\neq i}^{N}\frac{m_{i}m_{j}(q_{j}-q_{i})}{r_{ij}^{3}}=-\lambda{m_{i}(q_{i}-c)}.

We say that two central configurations are equivalent if they can be related by rigid motions or scalings with respect to the center of mass. Later when we talk about the number of central configurations, we always mean their equivalence class.

There is another equivalent characterization of central configurations. Let UU be the negative Newtonian potential energy or the self potential given by

(3) U=1i<jNNmimjrij,U=\sum_{1\leq i<j\leq N}^{N}\frac{m_{i}m_{j}}{r_{ij}},

and let II be the moment of inertia given by

(4) I=1M1i<jNNmimjrij2.I=\frac{1}{M}\sum_{1\leq i<j\leq N}^{N}m_{i}m_{j}r_{ij}^{2}.

If the center of mass is fixed at the origin of the plane, the moment of inertia is equal to I=imi|qi|2I=\sum_{i}m_{i}|q_{i}|^{2}. Due to the fact that the potential is homogeneous of degree 1-1, it is well known that λ=UI\lambda=\frac{U}{I} by the classical Euler Theorem. Now the central configurations are the critical points of UU restricted to the “shape sphere”   I=1I=1, or equivalently they are the critical points of IU2IU^{2} without constraints.

Central configurations are significant in understanding the solution structure and dynamical behavior of the NN-body problem. One of their essential roles is to provide insights into the dynamics in the vicinity of collisions. Additionally, these configurations serve as the initial positions for homographic periodic solutions that preserve the configuration’s shape. As a result, they find practical applications, for example, in spacecraft missions. For more on the background about central configurations, please refer to [22] and the references therein.

The famous and long-standing Chazy-Wintner-Smale conjecture in the NN-body problem says that for any given masses, there is only a finite number of central configurations ([22, 32]). It has been confirmed for the planar 44-body problem by Hampton and Moeckel ([14]) for any masses and for the planar 55-body problem by Albouy and Kaloshin ([5]) for any masses except for masses belonging to a codimension 22 subvariety in the mass space defined by a polynomial system. In the general case, it is still widely open.

In the current paper we will focus on central configurations in the planar 4-body problem. Though the finiteness of central configurations is already established and reconfirmed in [5, 14], a precise finite number or better upper bound on the number of central configurations is still expected.

A related problem concerns the counting of special central configurations (e.g., configurations with interesting geometry) especially for those masses with some kind of “symmetry”   (e.g., equal masses). Among them, the convex central configurations have attracted attention in many recent papers ([1, 2, 4, 7, 9, 11, 12, 18, 26, 27, 28, 29, 30, 34], just to mention a few) which are also our main focus in the current paper. By (strict) convexity we mean that no body is contained inside or on the convex hull of the other three bodies.

The study on convex central configurations dates back to the classical work [19] where MacMillan and Bartky proved that there exists a convex central configuration in the planar 44-body problem for any given cyclic ordering of four positive masses. They also derived useful information on the admissible shapes of the 44-body convex central configurations which was generalized and refined in [17] and will be helpful for us. Xia ([33]) reconfirmed it in a simpler way by proving the stronger existence of local minima of IU2IU^{2}. He also conjectured that for any N(5)N(\geq 5) positive masses and any cyclic order of NN points on S1S^{1}, there is a convex planar central configuration of minimum type with the cyclic order.

However for the planar 4-body problem it was conjectured ([3], Problem 10) that there is only one convex central configuration for any four fixed positive masses in a given order. The conjecture is true for four equal masses and some other special settings with constraints on the shapes or masses as cited above, but there is no proof for the general case. It has resisted any attempts for a while due to the fact that the Hessian matrix of the Newton potential is hard to control.

Recently, Corbera-Cors-Roberts ([10]) classified the full set of convex central configurations in the planar 4-body problem by choosing convenient coordinates, and proved it to be three dimensional. As remarked by the authors, the injectivity of the mass map defined in their paper should be established in order to confirm the uniqueness.

We confirm the uniqueness conjecture by a rigorous computer-assisted-proof for almost all masses. All symbolic and numerical computations were carried out by using MATLAB (2017b), the open-source software SageMath 9.6 (2022), and Maple 2022. More precisely, we prove the following theorem after normalizing the largest mass m2m_{2} to be 11.

Theorem 1.1.

For the planar Newtonian 4-body problem there is a unique convex central configuration for any fixed positive masses m=(m1,1,m3,m4)m=(m_{1},1,m_{3},m_{4}) such that (m1,m3,m4)[0.6,1]3(m_{1},m_{3},m_{4})\in[0.6,1]^{3} in a given order as in Figure 2.

Remark 1.2.
  • The settings and conventions will be given in detail in §2.

  • In the statement of the theorem, the number 0.60.6 is far from being optimal. In fact, for any (m1,m3,m4)(m_{1},m_{3},m_{4})\in [δ0,1]3[\delta_{0},1]^{3} with given δ0\delta_{0} a fixed positive real number, our methods should work. We will address this point in the last section.

  • The remaining masses are those near zero which we will also comment on the last section, and it seems to be an interesting and challenging problem both analytically and numerically.

  • Note that in our setting we fix the relative ordering of the masses and vary the masses. This is equivalent to prove the theorem for a given mass in any ordering.

The main idea of the proof can be summarized as follows: We initially validate the conjecture for a set of selected mass vectors m0m_{0} using interval analysis and MacMillan-Bartky’s classical result on the existence of convex central configurations in the domain of mass. Subsequently, by determining the implicit function theorem’s domain size, we can extend the central configuration’s existence and uniqueness to a neighborhood of m0m_{0}, say a mass ball B(m0,ϵ)B(m_{0},\epsilon) for the position ball B(x0,r)B(x_{0},r). Then we verify that no solutions to the central configuration equations exist outside the previously considered position domain B(x0,r)B(x_{0},r) for the mass ball B(m0,ϵ)B(m_{0},\epsilon). Lastly we proceed to establish the uniqueness within the stated closed domain [0.6,1]3[0.6,1]^{3} by iteratively implementing the aforementioned method. We call such a mass ball B(m0,ϵ)B(m_{0},\epsilon) a uniqueness mass ball if there is a unique convex central configuration for any mass in this ball. In short, the idea of the proof is to construct enough uniqueness mass balls to cover the whole mass space. Please refer to the flow chart in Figure 1 for a visual representation of the process.

Initialize masses: m2=1,m_{2}=1, (m1,m3,m4)[0.6,1]3(m_{1},m_{3},m_{4})\in[0.6,1]^{3} Initialize position x=[x1,y1,x3,y3]𝐔x=[x_{1},y_{1},x_{3},y_{3}]\in\mathbf{U} Process1: Interval Arithmetic & Krawczyk Operator, uniquness of CC for m0m_{0} Process 2: Implicit Function Theorem, uniqueness mass ball Process 3: Confirm that CCs for B(m0,ϵ)B(m_{0},\epsilon) are in B(x0,r)B(x_{0},r) Decision: uniqueness mass balls cover the mass space Stop generate more uniqueness mass ballsm0m_{0}x0x_{0} and m0m_{0}B(m0,ϵ)B(m_{0},\epsilon) and B(x0,r)B(x_{0},r)yesno
Figure 1: Flow chart to illustrate the computation process. The goal of the computation is to generate enough uniqueness mass balls to cover the whole mass space. Processes 1, 2, and 3 correspond to sections 3, 4 and 5 respectively.

Note that Lee and Santoprete ([16]) have already used the interval Krawczyk method to study the exact counting of the central configurations in the planar 55-body problem with equal masses. Moczurad and Zgliczynski ([23]) also used this method to study the central configurations for equal masses for planar 5,6,75,6,7-body problems. These papers are quite different from our situation here because we are interested in ALL masses.

After the introduction to the uniqueness conjecture and our main result establishing the conjecture with the idea of the proof in the current section, we give in §2 the form of the equations about central configurations convenient for our purpose with useful results which add suitable constraints about the convex central configuration. We introduce the basics about the interval algorithm in §3, then based on this we prove the existence of a convex central configuration for a mass vector at the meshed grid point. §4 is about the size of the domain where the implicit function theorem is valid. We use the theorem from [15], and for completeness, we give a proof. We exclude the possibilities for the existence of convex central configurations outside of the domain generated by IFT. Finally in §5 we establish a computer-assisted-proof to our main Theorem 1.1. We close the paper with some remarks and comments in §6.

As always insisted by Alain Chenciner, understanding the reasoning behind the truth is of primary importance, and we leave a more conceptual proof for the uniqueness to future research.

2 Equations for the Convex Central Configurations in the Planar 4-Body Problem

For the planar 44-body central configurations we can use the equivalent Dziobek-Laura-Andoyer equations ( [13], p.241)

(5) fij=k=1,ki,j4mk(RikRjk)Δijk=0f_{ij}=\sum_{k=1,k\neq i,j}^{4}m_{k}(R_{ik}-R_{jk})\Delta_{ijk}=0

for 1i<j41\leq i<j\leq 4, where Rij=1/rij3R_{ij}=1/r_{ij}^{3} and Δijk=(qiqk)×(qjqk)\Delta_{ijk}=(q_{i}-q_{k})\times(q_{j}-q_{k}) with ×\times the cross product of two vectors. Thus the area Δijk\Delta_{ijk} gives the twice the signed area of the triangle Δijk\Delta{ijk}. If the bodies mi,mj,mkm_{i},m_{j},m_{k} are in counter clockwise, the area of Δijk\Delta_{ijk} is positive. Otherwise, it is negative.

From the point of view of CAP (computer-assisted proofs) in the problem of finding and counting all CCs, the dimension and the size of the search domain is fundamental. The obstacles arise for the following reasons:

  • two or more bodies might be arbitrary close to each other (collisions, singularities of the system of equations);

  • bodies might be arbitrary far from each other (unbounded domain);

  • central configurations are not isolated due to the invariance of the system under rotation, scaling or reflection which could give some degeneracy in search algorithm.

Numerically solving the central configuration equations (5) directly in a naive way is not feasible due to these issues. However, these problems can be addressed either by imposing conditions on the masses or by employing alternative methods. By demanding sufficiently large masses (e.g., 0.60.6 as in Theorem 1.1), the first issue can be avoided, and we will revisit this point in the future, as discussed in §6. By leveraging existing results, we can reduce the number of variables and limit the search domain to a bounded region in the configuration space, far from collisions. This approach allows us to isolate central configurations and eliminate continuous symmetries.

Lemma 2.1 ([21], Perpendicular Bisector Theorem, p.510).

Let q=(q1,q2,,qn)q=(q_{1},q_{2},\cdots,q_{n}) be a planar central configuration for m(𝐑+)nm\in(\mathbf{R^{+}})^{n} and let qiq_{i} and qjq_{j} be any two positions. Denote by Q1,Q2,Q3,Q_{1},Q_{2},Q_{3}, and Q4Q_{4} the four open quadrants in 𝐑2\mathbf{R}^{2} in clockwise order obtained by deleting the line passing through qiq_{i} and qjq_{j} and the perpendicular bisector of the line segment qiqj¯\overline{q_{i}q_{j}}. Then if there exists some qkQ1Q3q_{k}\in Q_{1}\cup Q_{3}, there must exist at least another qlQ2Q4q_{l}\in Q_{2}\cup Q_{4}.

Lemma 2.2 ([17], Theorem 4.1).

For any m=(m1,m2,m3,m4)(𝐑+)4m=(m_{1},m_{2},m_{3},m_{4})\in(\mathbf{R}^{+})^{4}, let q=(q1,q2,q=(q_{1},q_{2}, q3,q_{3}, q4)q_{4}) be a convex central configuration for mm as in Figure 2. Let q2q4¯\overline{q_{2}q_{4}} be the line segment between points q2q_{2} and q4q_{4}. Then q1q_{1} and q3q_{3} must satisfy the following inequalities:

(6) 36r24<max{|q1q2q4¯|,|q3q2q4¯|}<32r24,\frac{\sqrt{3}}{6}r_{24}<\max\{|q_{1}-\overline{q_{2}q_{4}}|,|q_{3}-\overline{q_{2}q_{4}}|\}<\frac{\sqrt{3}}{2}r_{24},
(7) (132)r24<min{|q1q2q4¯|,|q3q2q4¯|}.\left(1-\frac{\sqrt{3}}{2}\right)r_{24}<\min\{|q_{1}-\overline{q_{2}q_{4}}|,|q_{3}-\overline{q_{2}q_{4}}|\}.
Refer to caption
Figure 2: Admissible Convex Configuration: q1,q2,q3,q4q_{1},q_{2},q_{3},q_{4} are in counter clockwise. q1=(x1,y1)q_{1}=(x_{1},y_{1}), q2=(1,0)q_{2}=(-1,0), q3=(x3,y3)q_{3}=(x_{3},y_{3}), and q4=(1,0)q_{4}=(1,0). [x1,y1,x3,y3]𝐔.[x_{1},y_{1},x_{3},y_{3}]\in\mathbf{U}.

Due to the invariance of central configurations under rotation and scaling, we can choose the coordinates for the 44-body configuration as q1=(x1,y1)q_{1}=(x_{1},y_{1}), q2=(1,0)q_{2}=(-1,0), q3=(x3,y3)q_{3}=(x_{3},y_{3}), and q4=(1,0)q_{4}=(1,0). From the above two lemmas, it is easy to prove the following theorem whose proof is omitted.

Theorem 2.3.

For any m=(m1,m2,m3,m4)(𝐑+)4m=(m_{1},m_{2},m_{3},m_{4})\in(\mathbf{R}^{+})^{4}, let q=(q1,q2,q3,q4)q=(q_{1},q_{2},q_{3},q_{4}) be a convex central configuration for mm. Without loss of generality, we assume q1=(x1,y1)q_{1}=(x_{1},y_{1}), q2=(1,0)q_{2}=(-1,0), q3=(x3,y3)q_{3}=(x_{3},y_{3}), and q4=(1,0)q_{4}=(1,0) as given in Figure 2 after rotation, translation and scaling. Then

(8) 1<x1<1,1<x3<1,x1x30,-1<x_{1}<1,-1<x_{3}<1,x_{1}x_{3}\geq 0,
(9) 2(132)<y1<3,3<y3<2(132),2\left(1-\frac{\sqrt{3}}{2}\right)<y_{1}<{\sqrt{3}},\hskip 17.07182pt-{\sqrt{3}}<y_{3}<-2\left(1-\frac{\sqrt{3}}{2}\right),
(10) (x1+1)2+y12<2,(x11)2+y12<2,\sqrt{(x_{1}+1)^{2}+y_{1}^{2}}<2,\hskip 17.07182pt\sqrt{(x_{1}-1)^{2}+y_{1}^{2}}<2,
(11) (x3+1)2+y32<2,(x31)2+y32<2.\sqrt{(x_{3}+1)^{2}+y_{3}^{2}}<2,\hskip 17.07182pt\sqrt{(x_{3}-1)^{2}+y_{3}^{2}}<2.

Denote by 𝐔\mathbf{U} the admissible region of the position vectors which satisfy the above conditions.
If m2>m4m_{2}>m_{4}, then 0<x1<10<x_{1}<1 and 0<x3<10<x_{3}<1. If m2=m4m_{2}=m_{4}, then x1=x3=0x_{1}=x_{3}=0.

For the 44-body problem, the equation (5) gives a system of six equations. Here are the expressions of the equations:

Δ123=(x3+1)y1y3(x1+1),Δ124=2y1,\Delta_{123}=\left({x_{3}}+1\right){y_{1}}-{y_{3}}\left({x_{1}}+1\right),\Delta_{124}=2y_{1},
Δ134=(x3+1)y1+y3(x11),Δ234=2y3.\Delta_{134}=\left(-{x_{3}}+1\right){y_{1}}+{y_{3}}\left({x_{1}}-1\right),\Delta_{234}=-2y_{3}.
(12) f12=m3(1((x1x3)2+(y1y3)2)321((x3+1)2+y32)32)((x3+1)y1y3(x1+1))+2m4(1((x11)2+y12)3218)y1.\begin{array}[]{ll}f_{12}=&{m_{3}}\left(\frac{1}{\left(\left({x_{1}}-{x_{3}}\right)^{2}+\left({y_{1}}-{y_{3}}\right)^{2}\right)^{\frac{3}{2}}}-\frac{1}{\left(\left({x_{3}}+1\right)^{2}+{y_{3}}^{2}\right)^{\frac{3}{2}}}\right)\left(\left({x_{3}}+1\right){y_{1}}-{y_{3}}\left({x_{1}}+1\right)\right)+\\ &2{m_{4}}\left(\frac{1}{\left(\left({x_{1}}-1\right)^{2}+{y_{1}}^{2}\right)^{\frac{3}{2}}}-\frac{1}{8}\right){y_{1}}.\end{array}
(13) f13=m2(1((x1+1)2+y12)321((x3+1)2+y32)32)((x3+1)y1+y3(x1+1))+m4(1((x11)2+y12)321((x31)2+y32)32)((x3+1)y1+y3(x11).)\begin{array}[]{ll}f_{13}=&{m_{2}}\left(\frac{1}{\left(\left({x_{1}}+1\right)^{2}+{y_{1}}^{2}\right)^{\frac{3}{2}}}-\frac{1}{\left(\left({x_{3}}+1\right)^{2}+{y_{3}}^{2}\right)^{\frac{3}{2}}}\right)\left(-\left({x_{3}}+1\right){y_{1}}+{y_{3}}\left({x_{1}}+1\right)\right)+\\ &{m_{4}}\left(\frac{1}{\left(\left({x_{1}}-1\right)^{2}+{y_{1}}^{2}\right)^{\frac{3}{2}}}-\frac{1}{\left(\left({x_{3}}-1\right)^{2}+{y_{3}}^{2}\right)^{\frac{3}{2}}}\right)\left(\left(-{x_{3}}+1\right){y_{1}}+{y_{3}}\left({x_{1}}-1\right).\right)\end{array}
(14) f14=2m2(1((x1+1)2+y12)3218)y1+m3(1((x1x3)2+(y1y3)2)321((x31)2+y32)32)((x3+1)y1y3(x11)).\begin{array}[]{ll}f_{14}=&-2{m_{2}}\left(\frac{1}{\left(\left({x_{1}}+1\right)^{2}+{y_{1}}^{2}\right)^{\frac{3}{2}}}-\frac{1}{8}\right){y_{1}}+\\ &{m_{3}}\left(\frac{1}{\left(\left({x_{1}}-{x_{3}}\right)^{2}+\left({y_{1}}-{y_{3}}\right)^{2}\right)^{\frac{3}{2}}}-\frac{1}{\left(\left({x_{3}}-1\right)^{2}+{y3}^{2}\right)^{\frac{3}{2}}}\right)\left(-\left(-{x_{3}}+1\right){y_{1}}-{y_{3}}\left({x_{1}}-1\right)\right).\end{array}
(15) f23=m1(1((x1+1)2+y12)321((x1x3)2+(y1y3)2)32)((x3+1)y1y3(x1+1))2m4(181((x31)2+y32)32)y3.\begin{array}[]{ll}f_{23}=&{m_{1}}\left(\frac{1}{\left(\left({x_{1}}+1\right)^{2}+{y_{1}}^{2}\right)^{\frac{3}{2}}}-\frac{1}{\left(\left({x_{1}}-{x_{3}}\right)^{2}+\left({y_{1}}-{y_{3}}\right)^{2}\right)^{\frac{3}{2}}}\right)\left(\left({x_{3}}+1\right){y_{1}}-{y_{3}}\left({x_{1}}+1\right)\right)-\\ &2{m_{4}}\left(\frac{1}{8}-\frac{1}{\left(\left({x_{3}}-1\right)^{2}+{y_{3}}^{2}\right)^{\frac{3}{2}}}\right){y_{3}}.\end{array}
(16) f24=2m1(1((x1+1)2+y12)321((x11)2+y12)32)y1+2m3(1((x3+1)2+y32)321((x31)2+y32)32)y3.\begin{array}[]{ll}f_{24}=&2{m_{1}}\left(\frac{1}{\left(\left({x_{1}}+1\right)^{2}+{y_{1}}^{2}\right)^{\frac{3}{2}}}-\frac{1}{\left(\left({x_{1}}-1\right)^{2}+{y_{1}}^{2}\right)^{\frac{3}{2}}}\right){y_{1}}+\\ &2{m_{3}}\left(\frac{1}{\left(\left({x_{3}}+1\right)^{2}+{y_{3}}^{2}\right)^{\frac{3}{2}}}-\frac{1}{\left(\left({x_{3}}-1\right)^{2}+{y_{3}}^{2}\right)^{\frac{3}{2}}}\right){y_{3}}.\end{array}
(17) f34=m1(1((x1x3)2+(y1y3)2)321((x11)2+y12)32)((x3+1)y1+y3(x11))2m2(1((x3+1)2+y32)3218)y3.\begin{array}[]{ll}f_{34}=&{m_{1}}\left(\frac{1}{\left(\left({x_{1}}-{x_{3}}\right)^{2}+\left({y_{1}}-{y_{3}}\right)^{2}\right)^{\frac{3}{2}}}-\frac{1}{\left(\left({x_{1}}-1\right)^{2}+{y_{1}}^{2}\right)^{\frac{3}{2}}}\right)\left(\left(-{x_{3}}+1\right){y_{1}}+{y_{3}}\left({x_{1}}-1\right)\right)-\\ &2{m_{2}}\left(\frac{1}{\left(\left({x_{3}}+1\right)^{2}+{y_{3}}^{2}\right)^{\frac{3}{2}}}-\frac{1}{8}\right){y_{3}}.\end{array}

Note that the equations fij(x,m)=0f_{ij}(x,m)=0 are homogeneous and in fact linear with respect to mm. We can label the body with the largest mass value to be m2m_{2} and the corresponding mass is normalized to one. Therefore m2=1m_{2}=1 and 0<mi10<m_{i}\leq 1 for i=1,3,4i=1,3,4. Then the problem of convex central configurations for the planar 44-body problem is to find the zeros of fij=0f_{ij}=0 in the bounded domain 𝐔\mathbf{U} with bounded mass parameters m1,m3,m_{1},m_{3}, and m4m_{4}.

It is easy to use Matlab command “fmincon”   or “fsolve”   to numerically find the convex central configuration in the domain 𝐔\mathbf{U} for any given positive masses. But the commands “fmincon”  or “fsolve”  can not give any information whether the solution is unique even locally. We will use the interval analysis and the Krawczyk operator to produce such information for some given mass values mm.

Theorem 2.4.

For each given mass vector m0=[m10,m20,m30,m40](𝐑+)4m_{0}=[m_{10},m_{20},m_{30},m_{40}]\in(\mathbf{R}^{+})^{4} at a meshed grid point such that m20=1m_{20}=1 and 0.6mi010.6\leq m_{i0}\leq 1 for i=1,3,4i=1,3,4, there exists a unique convex central configuration in the order as given in Figure 2.

The computer-assisted proof of this theorem relies on the interval analysis which we now turn to, and the algorithm will be given at the end of the next section. The precise meaning for ”a meshed grid point” in the statement and the final proof of this Theorem 2.4 will be given in Step 2 of the proof of the main Theorem 1.1 in section 5.

Remark 2.5.

For the convenience of implementing the interval arithmetic computation, we need to cover 𝐔\mathbf{U} by intervals. Let

𝐔¯=[0,1]×[0.267,1.733]×[0,1]×[1.733,0.267]\underline{\mathbf{U}}=[0,1]\times[0.267,1.733]\times[0,1]\times[-1.733,-0.267]

and

𝐔¯=[1,1]×[0,1.733]×[1,1]×[1.733,0].\bar{\mathbf{U}}=[-1,1]\times[0,1.733]\times[-1,1]\times[-1.733,0].

Then

𝐔¯𝐔𝐔¯.\underline{\mathbf{U}}\subset\mathbf{U}\subset\bar{\mathbf{U}}.

For a given mass vector m0m_{0} in Theorem 2.4, we find a unique convex central configuration x0𝐔¯x_{0}\in\underline{\mathbf{U}} by using the Krawczyk operator. Extending the uniqueness of m0m_{0} to a neighborhood B(m0,ϵ)B(m_{0},\epsilon) in the corresponding position neighborhood B(x0,r)B(x_{0},r), we check that any x𝐔¯\B(x0,r)x\in\bar{\mathbf{U}}\backslash B(x_{0},r) is not a central configuration for any mB(m0,ϵ).m\in B(m_{0},\epsilon). As a result, for any mB(m0,ϵ)m\in B(m_{0},\epsilon), there exists a unique convex central configuration.

3 Interval Analysis and Krawczyk Operator

3.1 Notations and Concepts of Interval Arithmetic

In order to determine whether the equations of central configuration fij(x,m)=0f_{ij}(x,m)=0 have any zeros, we employ interval arithmetic to test whether the functions fij(x,m)f_{ij}(x,m) contain zeros within specific intervals. Although the concept of interval arithmetic/analysis has a long history, the development of algorithms and computer programs to implement it dates back to the 1960s. A very nice survey of interval arithmetic was given by Alefeld and Mayer [6] which we will follow as a reference to provide a brief introduction to this topic. For more information, readers are referred to [6] and [24].

Let [a]=[a¯,a¯][a]=[\underline{a},\bar{a}] and [b]=[b¯,b¯][b]=[\underline{b},\bar{b}] be real closed intervals. If a¯=a¯\underline{a}=\bar{a}, i.e., if [a][a] consists only of the number aa, then we identify the interval [a,a][a,a] with the real number aa, i.e., a=[a,a].a=[a,a]. Let \circ be one of the basic arithmetic operations {+,,,/}.\{+,-,\cdot,/\}. The corresponding operations for intervals [a][a] and [b][b] are defined by

(18) [a][b]={ab|a[a],b[b]},[a]\circ[b]=\{a\circ b|a\in[a],b\in[b]\},

where we assume 0[b]0\notin[b] in case of division.

One of the advantages of interval computations is the fact that [a][b][a]\circ[b] can be represented by using only the bounds of [a][a] and [b][b]. For example:

[a]+[b]=[a¯+b¯,a¯+b¯],[a][b]=[a¯b¯,a¯b¯],[a]+[b]=[\underline{a}+\underline{b},\bar{a}+\bar{b}],\hskip 28.45274pt[a]-[b]=[\underline{a}-\bar{b},\bar{a}-\underline{b}],
[a][b]=[min{ab¯,a¯b¯,a¯b¯,a¯b¯},max{ab¯,a¯b¯,a¯b¯,a¯b¯}].[a]\cdot[b]=[\min\{\underline{ab},\underline{a}\bar{b},\bar{a}\underline{b},\bar{a}\bar{b}\},\max\{\underline{ab},\underline{a}\bar{b},\bar{a}\underline{b},\bar{a}\bar{b}\}].

If we define 1[b]={1b|b[b]}\frac{1}{[b]}=\left\{\frac{1}{b}|b\in[b]\right\} for 0[b]0\notin[b], then [a]/[b]=[a]1[b].[a]/[b]=[a]\cdot\frac{1}{[b]}.

Standard interval functions φΨ={sin,\varphi\in\Psi=\{\sin, cos\cos, tan,\tan, arctan\arctan, exp,\hbox{exp}, ln, abs,\hbox{abs}, sqr,sqrt}\hbox{sqr},\hbox{sqrt}\} are defined via their range

(19) φ([x])={φ(x)|x[x]}.\varphi([x])=\{\varphi(x)|x\in[x]\}.

Suppose that f:Df:D\subseteq\mathbb{R}\rightarrow\mathbb{R} is defined by a mathematical expression f(x)f(x) involving a finite number of elementary operations +,,,/+,-,\cdot,/ and standard functions φΨ\varphi\in\Psi. If we replace the variable xx by an interval [x]D[x]\subseteq D and evaluate the resulting expression using the rules in (18) and (19), the output is again an interval, denoted by f([x])f([x]). This is commonly referred to as the interval arithmetic evaluation of ff over [x][x]. We assume that f([x])f([x]) exists whenever it is used in this paper, without explicitly stating so each time.

Observed that

[x][y]f([x])f([y]) and x[x]f(x)f([x]),[x]\subseteq[y]\Rightarrow f([x])\subseteq f([y])\hbox{ and }x\in[x]\Rightarrow f(x)\in f([x]),

it is easy to prove the following useful property

R(f;[x])f([x]),R(f;[x])\subseteq f([x]),

where R(f;[x])R(f;[x]) denotes the range of ff over [x][x].

The above fundamental property of interval arithmetic states that we can calculate lower and upper bounds for the range of a function over an interval using only the bounds of the interval itself. However, when we use the interval arithmetic evaluation f([x])f([x]), we may overestimate this range. The accuracy of the approximation to the range of ff over an interval [x][x] strongly depends on how the expression for f(x)f(x) is formulated. To illustrate these differences, we construct the Example 3.1 to compare the interval arithmetic evaluation f([x])f([x]) with the range R(f;[x])R(f;[x]) in different expressions, and to demonstrate that the false statement that f(x)f(x) contains zeros can arise due to the overestimation of f([x])f([x]).

Example 3.1.

Let f(x)f(x) be the rational function

f(x)=2x2x2+91.f(x)=\frac{2x^{2}}{x^{2}+9}-1.
  • Observed that f(x)=36x(x2+9)2f^{\prime}(x)=\frac{36x}{(x^{2}+9)^{2}}, f(x)f(x) is increasing in the interval [x]=[1,4][x]=[1,4]. Therefore R(f,[1,4])=[0.8,0.28]R(f,[1,4])=[-0.8,0.28]. It is easy to see that the equation f(x)=0f(x)=0 has a unique solution x=3x=3 in this interval.

  • By direct computations,

    f([1,4])=2[1,4]2[1,4]2+91=2[1,16][1,16]+91=[2,32][10,25]1=[2,32]1[10,25]1f([1,4])=\frac{2\cdot[1,4]^{2}}{[1,4]^{2}+9}-1=\frac{2\cdot[1,16]}{[1,16]+9}-1=\frac{[2,32]}{[10,25]}-1=[2,32]\cdot\frac{1}{[10,25]}-1
    =[2,32][125,110]1=[225,3210]1=[0.92,2.2][0.8,0.28],=[2,32]\cdot\left[\frac{1}{25},\frac{1}{10}\right]-1=\left[\frac{2}{25},\frac{32}{10}\right]-1=[-0.92,2.2]\supseteq[-0.8,0.28],

    which confirms that R(f;[x])f([x]).R(f;[x])\subseteq f([x]). On the other hand, for x0x\not=0 we can rewrite f(x)f(x) as

    f(x)=21+9/x21.f(x)=\frac{2}{1+9/x^{2}}-1.
    f([1,4])=21+9/[1,4]21=21+9/[1,16]1=21+[9/16,9]1f([1,4])=\frac{2}{1+9/[1,4]^{2}}-1=\frac{2}{1+9/[1,16]}-1=\frac{2}{1+[9/16,9]}-1
    =2[25/16,10]1=[2/10,32/25]1=[8/10,7/25]=[0.8,0.28],=\frac{2}{[25/16,10]}-1=[2/10,32/25]-1=[-8/10,7/25]=[-0.8,0.28],

    which implies that f([1,4])=R(f,[1,4]).f([1,4])=R(f,[1,4]). This shows that the interval arithmetic evaluation depends on its expressions.

  • Now let us consider two sub-intervals [x1]=[1,5/2][x_{1}]=[1,5/2], [x2]=[5/2,4][x_{2}]=[5/2,4]. By simple computations we have

    f([x1])=f([1,5/2])=2[1,5/2]2[1,5/2]2+91=[2,25/2][10,61/4]1=[53/61,1/4].f([x_{1}])=f([1,5/2])=\frac{2\cdot[1,5/2]^{2}}{[1,5/2]^{2}+9}-1=\frac{[2,25/2]}{[10,61/4]}-1=[-53/61,1/4].
    f([x2])=f([5/2,4])=2[5/2,4]2[5/2,4]2+91=[25/2,32][61/4,25]1=[1/2,67/61].f([x_{2}])=f([5/2,4])=\frac{2\cdot[5/2,4]^{2}}{[5/2,4]^{2}+9}-1=\frac{[25/2,32]}{[61/4,25]}-1=[-1/2,67/61].

    The equation f(x)=0f(x)=0 may have solutions in both [x1][x_{1}] and [x2][x_{2}]. But we can exclude the interval [x1][x_{1}] by splitting it into two smaller intervals [x11]=[1,2][x_{11}]=[1,2] and [x12]=[2,5/2][x_{12}]=[2,5/2].

    f([x11])=f([1,2])=[1113,15],f([x12])=f([2,5/2])=[2961,126],f([x_{11}])=f([1,2])=\left[-\frac{11}{13},-\frac{1}{5}\right],\hskip 28.45274ptf([x_{12}])=f([2,5/2])=\left[-\frac{29}{61},-\frac{1}{26}\right],

    which shows that f(x)=0f(x)=0 has no solutions in the intervals [x11][x_{11}] and [x12].[x_{12}]. Therefore there is no solutions in [x1]=[x11][x12].[x_{1}]=[x_{11}]\bigcup[x_{12}].

Note that our functions fijf_{ij} only involve basic arithmetic operations and square root. We employ interval arithmetic evaluation to determine if fijf_{ij} contains zeros in given intervals. Example 3.1 demonstrates how we can rule out the existence of zeros in [x][x], as well as an example of a false statement of zeros in [x][x] due to R(f;[x])f([x])R(f;[x])\subseteq f([x]).

Interval arithmetic has been implemented on many platforms such as INTLAB and SageMath. INTLAB is a MATLAB library for interval arithmetic routines and it is written in MATLAB. SageMath is a free open-source mathematics software systems. The implementation of our computations on interval arithmetic is on SageMath 9.6 (published in 2022). Now let us introduce the Krawczyk operator which is an extension of Newton’s method to search zeros of nonlinear system by incorporating interval computations.

3.2 Krawczyk Operator

Let F:𝐑n𝐑nF:\mathbf{R}^{n}\rightarrow\mathbf{R}^{n} be a C1C^{1}-map. The Krawczyk operator (Krawczyk 1969; Neumeier 1990; Alefeld 1994) is an interval analysis tool to establish the existence and uniqueness of zero for a system of nn nonlinear equations in the same number nn of variables:

(20) F(x)=0.F(x)=0.

The method proposed by Krawczyk for finding zeros of FF is as follows. Set

  • [x]𝐑n[x]\subset\mathbf{R}^{n} be an interval set (i.e., Cartesian product of intervals).

  • x0[x].x_{0}\in[x]. Typically x0x_{0} is chosen to be the midpoint of [x][x]. We will denote this by x0=mid([x])x_{0}=mid([x]).

  • C𝐑n×nC\in\mathbf{R}^{n\times n} be a linear isomorphism.

The Krawczyk operator is defined to be

(21) K(x0,[x],F):=x0CF(x0)+(ICdF([x]))([x]x0).K(x_{0},[x],F):=x_{0}-C\cdot F(x_{0})+(I-C\cdot dF([x]))([x]-x_{0}).

The motivation and concise derivation of Krawczyk operator can be found in [23], where the operator was used to study the existence of central configurations for planar nn-body problem with equal masses for n=5,6,7n=5,6,7. More detailed discussion is referred to [24].

Lemma 3.2.
  1. 1.

    If x[x]x^{*}\in[x] and F(x)=0F(x^{*})=0, then xK(x0,[x],F).x^{*}\in K(x_{0},[x],F).

  2. 2.

    If K(x0,[x],F)int[x],K(x_{0},[x],F)\subset int[x], then there exists in [x][x] exactly one solution of equation F(x)=0F(x)=0. This solution is non-degenerate, i.e., dF(x)dF(x) is an isomorphism.

  3. 3.

    If K(x0,[x],F)[x]=K(x_{0},[x],F)\cap[x]=\emptyset, then F(x)0F(x)\not=0 for all x[x].x\in[x].

Remark 3.3.
  • For the concrete Krawczyk operator used in our problem the parameter mm is not considered as an interval, and we employ IFT to cover the neighboring mass points. The referee suggests that treating the masses parameter as an interval in the Krawczyk operator may be more effective than the IFT-based approach. Currently, it is unclear to us whether this would indeed be the case. However, we find that the proof using IFT is conceptually more appealing.

  • In the Krawczyk operator, CC is usually taken to be the inverse of dF(x0)dF(x_{0}) which is updated in each iteration. If dF(x)dF(x^{*}) is singular, we can not use the Krawczyk operator.

  • The point 2 in the above lemma gives us the way to establish the existence of unique zero of FF in [x][x], whilst the point 3 rules out the existence of zeros in [x][x]. There is a quite common case that [x]K(x0,[x],F)[x]\subseteq K(x_{0},[x],F) when the diameter of the interval [x][x] is not small enough. In this case Krawczyk operator won’t give us some useful information. In our program, we just simply split it into two small subintervals by bisecting the longest side of the interval.

  • We use Sagemath 9.6 to implement the interval computations and set RIF= RealIntervalField( ).

The mass vectors are formulated by fixing m2=1m_{2}=1 and taking mim_{i} at δ0+(1δ0)kM1\delta_{0}+\frac{(1-\delta_{0})k}{M_{1}} for k=0,,k=0,\cdots, M1M_{1} and i=1,3,4,δ0=0.6i=1,3,4,\delta_{0}=0.6. Our program run through all the cases. We are able to prove that there is a unique convex central configuration for each given mass vector. We first use the interval arithmetic computation and the Krawczyk operator to prove Theorem 2.4.

3.3 The Algorithm to Prove Theorem 2.4

Take m0=[m01,1,m03,m04]m_{0}=[m_{01},1,m_{03},m_{04}] and 0<m0i10<m_{0i}\leq 1 at the meshed grid points at δ0+(1δ0)kM1\delta_{0}+\frac{(1-\delta_{0})k}{M_{1}}. The notation fij(x)=fij(x,m0)f_{ij}(x)=f_{ij}(x,m_{0}) is used for each fixed m0m_{0}.

The algorithm consists of four steps:

Step 1. For the convenience of interval arithmetic computation, we use the interval 𝐔¯\underline{\mathbf{U}} such that [x1,y1,x3,y3]𝐔¯[x_{1},y_{1},x_{3},y_{3}]\in\underline{\mathbf{U}} for the admissible set 𝐔\mathbf{U} , i.e.

x1[0,1],y1[0.267,1.733],x3[0,1], and y3[1.733,0.267].x_{1}\in[0,1],y_{1}\in[0.267,1.733],x_{3}\in[0,1],\hbox{ and }y_{3}\in[-1.733,-0.267].

We split the configuration space into small intervals. Each side is divided into N1N_{1} equal intervals. The typical initial interval is [x]k=[x1]k1×[y1]k2×[x3]k3×[y3]k4[x]_{k}=[x_{1}]_{k_{1}}\times[y_{1}]_{k_{2}}\times[x_{3}]_{k_{3}}\times[y_{3}]_{k_{4}} where [x1]k1=[k1N1,k1+1N1][x_{1}]_{k_{1}}=[\frac{k_{1}}{N_{1}},\frac{k_{1}+1}{N_{1}}] refers to the k1k_{1}-th subinterval for x1x_{1} and the others are similar.

Step 2. The initial interval [x]k[x]_{k} is ruled out if one of the six functions fij([x]k)f_{ij}([x]_{k}) doesn’t contain zero. The rest intervals are candidates for the existence of solutions and they are labeled as [x]k[x]_{k}, k=0,1,2,,kendk=0,1,2,\cdots,k_{end}. If no intervals are left, then there is no solution for the given mass in the given configuration space and we stop. Otherwise, go to next step.

Step 3. The Krawczyk operator is used as a part of iteration process for each interval [x]k[x]_{k} from Step 2:

Let F(x)=[f12(x),f13(x),f14(x),f23(x)].F(x)=[f_{12}(x),f_{13}(x),f_{14}(x),f_{23}(x)].

  1. 1.

    given [x]k𝐔¯𝐑4[x]_{k}\subseteq\underline{\mathbf{U}}\subset\mathbf{R}^{4}

  2. 2.

    Compute [y]=K(mid[x]k,[x]k,F)[y]=K(mid[x]_{k},[x]_{k},F)

  3. 3.

    if [y]int[x]k[y]\subset int[x]_{k} and the diameter d[y]=y¯y¯d[y]=\bar{y}-\underline{y} of [y][y] is smaller than the error tolerance, then return success and [y][y] and go to Step 4

    elseif [y]int[x]k[y]\subset int[x]_{k} and the diameter d[y]=y¯y¯d[y]=\bar{y}-\underline{y} of [y][y] is not smaller than the error tolerance, then go to 2 by replacing [x]k[x]_{k} by [y][y]

    elseif [x]k[y][x]_{k}\subset[y], then bisect the longest side of the interval [x]k[x]_{k} and go to Step 2 for the two smaller intervals

    elseif [x]k[y]=[x]_{k}\cap[y]=\emptyset, then return no solution in [x]k[x]_{k}

    elseif [x]k[y][x]_{k}\cap[y]\not=\emptyset, then set [x]k+1:=[y][x]k[x]_{k+1}:=[y]\cap[x]_{k} and go to Step 2

Step 4. If Step 3 returns success and [y][y], f24([y])f_{24}([y]) and f34([y])f_{34}([y]) are evaluated. If both contain zero, then [y][y] is a central configuration for the given mass vector m0m_{0} and it returns a unique solution [y][y].

In SageMath 9.6, there are two styles for printing intervals: ‘brackets’ style and ‘question’ style. In question style, we print the “known correct” part of the number, followed by a question mark. The question mark indicates that the preceding digit is possibly wrong by ±1.\pm 1. In the report of the computation results, we may use either of them.

If only one solution is obtained at the end of the program, it implies the existence of a unique central configuration in the given order as in Figure 2. Two examples are provided to demonstrate the numerically rigorous proof of the existence and uniqueness of the convex central configurations for a given mass.

Example 3.4.

For m0=[1,1,1,1]m_{0}=[1,1,1,1], there is a unique convex central configuration and x1=1.?e20,x_{1}=-1.?e-20, y1=1.00000000000000000000?y1=1.00000000000000000000?, x3=1.?e20x_{3}=-1.?e-20,
y3=1.00000000000000000000?.y_{3}=-1.00000000000000000000?.

In this case, we take N1=20N_{1}=20. Each side of the configuration is partitioned into 20 smaller intervals with length about 0.050.05 to 0.07330.0733. There are 204=160,00020^{4}=160,000 intervals in total to check in Step 2. Most of them are ruled out because at least one of the six functions fijf_{ij} contains no zero. Among them only 52 intervals are left for Krawczyk operator in Step 3. Then only one of these intervals returns the above unique solution and others return no solution. We also note that no bisection is needed for these intervals. So we conclude that there is a unique central configurations for this mass m0m_{0} in the given ordering. This way we recover the well known result ([2]): The square is the only convex central configuration in the planar 44-body problem with equal masses.

Example 3.5.

For m0=[0.2,1,0.3,0.4]m_{0}=[0.2,1,0.3,0.4], there is a unique convex central configuration and

x1[0.15385328707521665441880634609065,0.15385328707521665448042670172252],x_{1}\in[0.15385328707521665441880634609065,0.15385328707521665448042670172252],
y1[1.4086619698548151890849667722929,1.4086619698548151891412057380232],y_{1}\in[1.4086619698548151890849667722929,1.4086619698548151891412057380232],
x3[0.11611158428853550145374041262273,0.11611158428853550155068635233120],x_{3}\in[0.11611158428853550145374041262273,0.11611158428853550155068635233120],
y3[1.484878202646704369219144317714,1.484878202646704369133519439332].y_{3}\in[-1.484878202646704369219144317714,-1.484878202646704369133519439332].

In this case, we take N1=15N_{1}=15. Each side of the configuration is partitioned into 15 smaller intervals with length about 0.06670.0667 to 0.09770.0977. There are totally 154=50,62515^{4}=50,625 intervals to check in Step 2. Only 62 intervals are left for Krawczyk operator in Step 3. But this time, 60 of them return no solutions. Two of them are bisected their intervals and went back to Step 2. Only one of them returns the above unique solution and others return no solution. We still have the conclusion that there is a unique central configurations for this mass m0m_{0} in the given ordering.

This completes the proof of Theorem 2.4 for the given mass m0m_{0} at the meshed grid points. To get the proof for all masses, we have to fill in the gaps among the discrete masses by the Cartesian product of intervals. Here we employ implicit function theorem to produce such domains. Note that for any central configuration x0=x0(m0)=[x10,y10,x30,y30]x_{0}=x_{0}(m_{0})=[x_{10},y_{10},x_{30},y_{30}] with the corresponding mass m0m_{0}, F(x0,m0)=0F(x_{0},m_{0})=0. By implicit function theorem, if a nondegeneracy condition on (x0,m0)(x_{0},m_{0}) is fulfilled, there exist two radius r>0r>0 and ϵ>0\epsilon>0, and a unique function x=x(m)x=x(m) in the neighborhood B(x0,r)B(x_{0},r) of x0x_{0} and B(m0,ϵ)B(m_{0},\epsilon) of m0m_{0} such that x0=x(m0)x_{0}=x(m_{0}) and F(x(m),m)0F(x(m),m)\equiv 0. Since m2=1m_{2}=1 is fixed once and for all, m0m_{0} means m0=(m10,m30,m40)m_{0}=(m_{10},m_{30},m_{40}) when m0m_{0} is considered as a variable. Hence

B(m0,ϵ)={(m1,m3,m4)|(m1m10)2+(m3m30)2+(m4m40)2<ϵ}.B(m_{0},\epsilon)=\{(m_{1},m_{3},m_{4})|\sqrt{(m_{1}-m_{10})^{2}+(m_{3}-m_{30})^{2}+(m_{4}-m_{40})^{2}}<\epsilon\}.

In the next section, we will refine the implicit function theorem to estimate the radius rr and ϵ\epsilon of the neighborhoods of the concerning point. For each such mass neighborhood, we prove the existence and uniqueness of convex central configuration.

4 Implicit Function Theorem: existence, uniqueness and size of valid domain

Before we introduce the implicit function theorem (IFT) and its applications to our case, we first briefly go over the notations and the definitions that are used in this paper following [15], especially the norm we used should be consistent. For x=(x1,x2,,xn)nx=(x_{1},x_{2},\cdots,x_{n})\in\mathbb{R}^{n}, x=x12+x22++xn2\|x\|=\sqrt{x_{1}^{2}+x_{2}^{2}+\cdots+x_{n}^{2}} is the standard Euclidean norm. For a given real number r>0r>0 and x0n,x_{0}\in\mathbb{R}^{n},

B(x0,r)={xn|xx0<r}B(x_{0},r)=\{x\in\mathbb{R}^{n}|\|x-x_{0}\|<r\}

denotes the open ball of radius rr centered at x0.x_{0}. For given positive integers n1,n2n_{1},n_{2}\in\mathbb{N}, the norm A\|A\| of a linear map A:n1n2A:\mathbb{R}^{n_{1}}\rightarrow\mathbb{R}^{n_{2}} is defined by

A:=sup{Ax|xB(0,1)n1},\|A\|:=\sup\{\|Ax\||x\in B(0,1)\subseteq\mathbb{R}^{n_{1}}\},

which is equivalent to

A=sup{Ax|xn1 and x=1} or A=sup{Axx|xn1}\|A\|=\sup\{\|Ax\||x\in\mathbb{R}^{n_{1}}\hbox{ and }\|x\|=1\}\hbox{ or }\|A\|=\sup\left\{\frac{\|Ax\|}{\|x\|}|x\in\mathbb{R}^{n_{1}}\right\}

with Ax\|Ax\| the norm of AxAx in n2\mathbb{R}^{n_{2}}. The standard Euclidean norm is easily computed in Matlab by using the command ‘norm’.

For our later purpose, let

f:n1×n2\displaystyle f:\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}} \displaystyle\rightarrow n1\displaystyle\mathbb{R}^{n_{1}}
(x,m)\displaystyle(x,m) \displaystyle\mapsto f(x,m)\displaystyle f(x,m)

be a multivariable vector-valued function. ff has the form f=(f1,f2,,fn1)f=(f_{1},f_{2},\cdots,f_{n_{1}}) with components f1,f2,,fn1:n1×n2f_{1},f_{2},\cdots,f_{n_{1}}:\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}\rightarrow\mathbb{R} multivariable real-valued functions of (x,m)(x,m). For our problem f(x,m)f(x,m) is a function of position xx with mass parameter mm. Dxf(x0,m0)D_{x}f(x_{0},m_{0}) and Dmf(x0,m0)D_{m}f(x_{0},m_{0}) are Jacobian matrices of ff with respect to xx and mm respectively.

The second derivatives of ff with respect to xx is an array of n1n_{1} Hessian matrices one for each component fjf_{j}:

Dx2f=[Hx(f1),Hx(f2),,Hx(fn1)],D_{x}^{2}f=[H_{x}(f_{1}),H_{x}(f_{2}),\cdots,H_{x}(f_{n_{1}})],

whose norm Dx2f(x0,m0)\|D_{x}^{2}f(x_{0},m_{0})\| is computed as

Dx2f(x0,m0)=sup{Dx2f(x0,m0)(v,w)|vn1,wn1 and v=w=1}.\|D_{x}^{2}f(x_{0},m_{0})\|=\sup\{\|D_{x}^{2}f(x_{0},m_{0})(v,w)\||v\in\mathbb{R}^{n_{1}},w\in\mathbb{R}^{n_{1}}\hbox{ and }\|v\|=\|w\|=1\}.

The second derivative DxDmf(x0,m0)D_{x}D_{m}f(x_{0},m_{0}) of ff with respect to variables mm and then variable xx is defined similarly.

Now it is ready to state IFT.

4.1 Implicit Function Theorem

We start with a useful lemma.

Lemma 4.1.

Given a continuous map φ:clB(y0,r)×Un1\varphi:\mathrm{cl}\mathrm{B}(y_{0},r)\times U\rightarrow\mathbb{R}^{n_{1}} with clB(y0,r)n1\mathrm{cl}\mathrm{B}(y_{0},r)\subset\mathbb{R}^{n_{1}} a closed ball and Un2U\subset\mathbb{R}^{n_{2}} an open subset and a fixed point x0Ux_{0}\in U, the function

l:U\displaystyle l:U \displaystyle\rightarrow \displaystyle\mathbb{R}
x\displaystyle x \displaystyle\mapsto l(x):=maxyclB(y0,r)φ(y,x)φ(y,x0)\displaystyle l(x):=\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x)-\varphi(y,x_{0})\|

is continuous on UU.

Proof.

By triangle inequality, for any x1Ux_{1}\in U

φ(y,x)φ(y,x0)\displaystyle\|\varphi(y,x)-\varphi(y,x_{0})\| \displaystyle\leq φ(y,x)φ(y,x1)+φ(y,x1)φ(y,x0)\displaystyle\|\varphi(y,x)-\varphi(y,x_{1})\|+\|\varphi(y,x_{1})-\varphi(y,x_{0})\|
\displaystyle\leq maxyclB(y0,r)φ(y,x)φ(y,x1)+maxyclB(y0,r)φ(y,x1)φ(y,x0)\displaystyle\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x)-\varphi(y,x_{1})\|+\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x_{1})-\varphi(y,x_{0})\|
\displaystyle\leq maxyclB(y0,r)φ(y,x)φ(y,x1)+l(x1).\displaystyle\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x)-\varphi(y,x_{1})\|+l(x_{1}).

So

l(x)=maxyclB(y0,r)φ(y,x)φ(y,x0)maxyclB(y0,r)φ(y,x)φ(y,x1)+l(x1),l(x)=\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x)-\varphi(y,x_{0})\|\leq\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x)-\varphi(y,x_{1})\|+l(x_{1}),

that is

l(x)l(x1)maxyclB(y0,r)φ(y,x)φ(y,x1).l(x)-l(x_{1})\leq\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x)-\varphi(y,x_{1})\|.

Similarly, we have

l(x1)l(x)maxyclB(y0,r)φ(y,x1)φ(y,x).l(x_{1})-l(x)\leq\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x_{1})-\varphi(y,x)\|.

So x1U\forall x_{1}\in U, we have

|l(x)l(x1)|maxyclB(y0,r)φ(y,x)φ(y,x1).|l(x)-l(x_{1})|\leq\max_{y\in\mathrm{cl}\mathrm{B}(y_{0},r)}\|\varphi(y,x)-\varphi(y,x_{1})\|.

Take a closed ball clB(x1,rx1)U\mathrm{cl}\mathrm{B}(x_{1},r_{x_{1}})\subset U, then φ\varphi is uniformly continuous on clB(y0,r)×clB(x1,rx1)\mathrm{cl}\mathrm{B}(y_{0},r)\times\mathrm{cl}\mathrm{B}(x_{1},r_{x_{1}}). Thus, ϵ>0\forall\epsilon>0, δ>0\exists\delta>0, if (y3,x3)(y2,x2)<δ\|(y_{3},x_{3})-(y_{2},x_{2})\|<\delta, then φ(y3,x3)φ(y2,x2)<ϵ\|\varphi(y_{3},x_{3})-\varphi(y_{2},x_{2})\|<\epsilon. So for xUx\in U such that xx1<δ\|x-x_{1}\|<\delta, we have (y,x)(y,x1)<δ\|(y,x)-(y,x_{1})\|<\delta, φ(y,x)φ(y,x1)<ϵ\|\varphi(y,x)-\varphi(y,x_{1})\|<\epsilon for any yclB(y0,r)y\in\mathrm{cl}\mathrm{B}(y_{0},r) which means that

|l(x)l(x1)|<ϵ.|l(x)-l(x_{1})|<\epsilon.

The proof is complete. ∎

Now we use the lemma to prove a version of implicit function theorem used on the estimation of the valid domain.

Theorem 4.2.

Let Yn1Y\subset\mathbb{R}^{n_{1}}, Xn2X\subset\mathbb{R}^{n_{2}} be nonempty open sets. Let Y×X(y,x)f(y,x)n1Y\times X\owns(y,x)\mapsto f(y,x)\in\mathbb{R}^{n_{1}} be a 𝒞1\mathscr{C}^{1} map. For (y0,x0)Y×X(y_{0},x_{0})\in Y\times X, suppose Dyf(y0,x0)D_{y}f(y_{0},x_{0}) is invertible. Then for any neighborhood 𝒪(y0)\mathscr{O}(y_{0}) of y0y_{0}, there is a neighborhood 𝒪(x0)\mathscr{O}(x_{0}) of x0x_{0} and a unique continuous map g:𝒪(x0)𝒪(y0)g:\mathscr{O}(x_{0})\rightarrow\mathscr{O}(y_{0}) such that

1. g(x0)=y0g(x_{0})=y_{0},

2. f(g(x),x)=f(y0,x0)f(g(x),x)=f(y_{0},x_{0}) for all x𝒪(x0)x\in\mathscr{O}(x_{0}),

3. gg is continuous at point x0x_{0}.

Proof.

The proof is based on Banach Fixed-Point Theorem.

Let T:=Dyf(y0,x0)T:=D_{y}f(y_{0},x_{0}) and ω0=f(y0,x0)\omega_{0}=f(y_{0},x_{0}). Define

φ:Y×X\displaystyle\varphi:Y\times X \displaystyle\rightarrow n1\displaystyle\mathbb{R}^{n_{1}}
(y,x)\displaystyle(y,x) \displaystyle\mapsto φ(y,x):=yT1(f(y,x)ω0).\displaystyle\varphi(y,x):=y-T^{-1}(f(y,x)-\omega_{0}).

Then φ(y,x)=y\varphi(y,x)=y iff f(y,x)=ω0f(y,x)=\omega_{0}. We will show that φ(,x)\varphi(\cdot,x) has a fixed point. From the definition of φ\varphi, we have

Dyφ(y,x)=IT1Dyf(y,x).D_{y}\varphi(y,x)=I-T^{-1}D_{y}f(y,x).

Since Dyφ(y,x)D_{y}\varphi(y,x) is continuous and Dyφ(y0,x0)=0D_{y}\varphi(y_{0},x_{0})=0, there exists r>0r>0 such that Dyφ(y,x)12\|D_{y}\varphi(y,x)\|\leq\frac{1}{2} for any yclB(y0,r)y\in\text{cl}\text{B}(y_{0},r) and xB(x0,ϵ1)x\in\mathrm{B}(x_{0},\epsilon_{1}). By the Mean Value Inequality,

φ(y,x0)y0\displaystyle\|\varphi(y,x_{0})-y_{0}\| =φ(y,x0)φ(y0,x0)\displaystyle=\|\varphi(y,x_{0})-\varphi(y_{0},x_{0})\|
sup0t1Dyφ((1t)y0+ty,x0)yy0\displaystyle\leq\sup_{0\leq t\leq 1}\|D_{y}\varphi((1-t)y_{0}+ty,x_{0})\|\|y-y_{0}\|
12yy0r2,yclB(y0,r).\displaystyle\leq\frac{1}{2}\|y-y_{0}\|\leq\frac{r}{2},\forall y\in\mathrm{cl}\mathrm{B}(y_{0},r).

For each xx, let l(x):=maxyclB(y0,r)φ(y,x)φ(y,x0)l(x):=\max_{y\in\text{cl}\text{B}(y_{0},r)}\|\varphi(y,x)-\varphi(y,x_{0})\| as in the previous Lemma 4.1. Since l(x)l(x) is continuous by the lemma and l(x0)=0l(x_{0})=0, there exists ϵ(0,ϵ1)\epsilon\in(0,\epsilon_{1}) such that

l(x)<r2,xB(x0,ϵ).l(x)<\frac{r}{2},\,\,\forall x\in\text{B}(x_{0},\epsilon).

For each xB(x0,ϵ)x\in\text{B}(x_{0},\epsilon), we have

φ(y,x)y0φ(y,x)φ(y,x0)+φ(y,x0)y0<r2+r2=r,yclB(y0,r).\|\varphi(y,x)-y_{0}\|\leq\|\varphi(y,x)-\varphi(y,x_{0})\|+\|\varphi(y,x_{0})-y_{0}\|<\frac{r}{2}+\frac{r}{2}=r,\,\,\forall y\in\text{cl}\text{B}(y_{0},r).

Thus, for each xB(x0,ϵ)x\in\text{B}(x_{0},\epsilon), φ(,x)\varphi(\cdot,x) maps

clB(y0,r)yφ(y,x)clB(y0,r).\text{cl}\text{B}(y_{0},r)\owns y\mapsto\varphi(y,x)\in\text{cl}\text{B}(y_{0},r).

For y1y_{1}, y2clB(y0,r)y_{2}\in\mathrm{cl}\mathrm{B}(y_{0},r) and xB(x0,ϵ)x\in\mathrm{B}(x_{0},\epsilon),

φ(y1,x)φ(y2,x)max0t1Dyφ((1t)y1+ty2,x)y1y212y1y2.\|\varphi(y_{1},x)-\varphi(y_{2},x)\|\leq\max_{0\leq t\leq 1}\|D_{y}\varphi((1-t)y_{1}+ty_{2},x)\|\|y_{1}-y_{2}\|\leq\frac{1}{2}\|y_{1}-y_{2}\|.

Hence by Banach Fixed-Point Theorem, for each xB(x0,ϵ)x\in\text{B}(x_{0},\epsilon), there exists a g(x)clB(y0,r)g(x)\in\text{cl}\text{B}(y_{0},r) such that φ(g(x),x)=g(x)\varphi(g(x),x)=g(x) or equivalently f(g(x),x)=ω0f(g(x),x)=\omega_{0}. ∎

Based on the proof of the implicit function theorem, one can get the following estimate about the size of the valid domain.

Theorem 4.3.

([15], Proposition 3.12) Let Yn1Y\subset\mathbb{R}^{n_{1}} and Xn2X\subset\mathbb{R}^{n_{2}} be two nonempty open sets. Let f:Y×Xn1f:Y\times X\rightarrow\mathbb{R}^{n_{1}} be a 𝒞2\mathscr{C}^{2} map. For some (y0,x0)Y×X(y_{0},x_{0})\in Y\times X, Dyf(y0,x0)D_{y}f(y_{0},x_{0}) is invertible. Define

Lx:=Dxf(y0,x0),My:=(Dyf(y0,x0))1L_{x}:=\|D_{x}f(y_{0},x_{0})\|,\ M_{y}:=\left\|(D_{y}f(y_{0},x_{0}))^{-1}\right\|

and for any given R1>0R_{1}>0, R2>0R_{2}>0, define

Kyy\displaystyle K_{yy} :=supyB(y0,R1)Dy2f(y,x0),\displaystyle:=\sup_{y\in\mathrm{B}(y_{0},R_{1})}\left\|D_{y}^{2}f(y,x_{0})\right\|,
Kxx\displaystyle K_{xx} :=sup(y,x)B((y0,x0),R2)Dx2f(y,x),\displaystyle:=\sup_{(y,x)\in\mathrm{B}((y_{0},x_{0}),R_{2})}\left\|D_{x}^{2}f(y,x)\right\|,
Kxy\displaystyle K_{xy} :=sup(y,x)B((y0,x0),R2)DxDyf(y,x).\displaystyle:=\sup_{(y,x)\in\mathrm{B}((y_{0},x_{0}),R_{2})}\left\|D_{x}D_{y}f(y,x)\right\|.

Then for all rmin{P,R2}r\leq\min\{P,R_{2}\} with P:=min{12MyKyy,R1}P:=\min\left\{\frac{1}{2M_{y}K_{yy}},R_{1}\right\}, and ϵmin{R22r2,14MyKxy}\epsilon\leq\min\left\{\sqrt{R_{2}^{2}-r^{2}},\frac{1}{4M_{y}K_{xy}}\right\} satisfying

My(Lx+Kxyr+Kxxϵ)ϵ<r2,M_{y}(L_{x}+K_{xy}r+K_{xx}\epsilon)\epsilon<\frac{r}{2},

xB(x0,ϵ)\forall x\in\mathrm{B}(x_{0},\epsilon) there exists a unique g(x)clB(y0,r)g(x)\in\mathrm{cl}\mathrm{B}(y_{0},r) such that f(g(x),x)=f(y0,x0)f(g(x),x)=f(y_{0},x_{0}).

Proof.

Using the Fundamental Theorem of Calculus, we have

Dyφ(y,x0)\displaystyle D_{y}\varphi(y,x_{0}) =01Dy2φ((1s)y0+sy,x0)(yy0)ds\displaystyle=\int_{0}^{1}D_{y}^{2}\varphi((1-s)y_{0}+sy,x_{0})(y-y_{0})\mathrm{d}s
=T101Dy2f((1s)y0+sy,x0)(yy0)ds.\displaystyle=-T^{-1}\int_{0}^{1}D_{y}^{2}f((1-s)y_{0}+sy,x_{0})(y-y_{0})\mathrm{d}s.

Thus

Dyφ(y,x0)MyKyyyy0,yB(y0,R1),\|D_{y}\varphi(y,x_{0})\|\leq M_{y}K_{yy}\|y-y_{0}\|,\ \forall y\in\mathrm{B}(y_{0},R_{1}),

by noting T=Dyf(y0,x0)T=D_{y}f(y_{0},x_{0}) as in the previous proof. From the Mean Value Inequality,

φ(y,x0)y0\displaystyle\|\varphi(y,x_{0})-y_{0}\| =φ(y,x0)φ(y0,x0)\displaystyle=\|\varphi(y,x_{0})-\varphi(y_{0},x_{0})\|
max0t1Dyφ((1t)y0+ty,x0)yy0\displaystyle\leq\max_{0\leq t\leq 1}\|D_{y}\varphi((1-t)y_{0}+ty,x_{0})\|\|y-y_{0}\|
MyKyyyy0yy0\displaystyle\leq M_{y}K_{yy}\|y-y_{0}\|\|y-y_{0}\|
MyKyyPyy0\displaystyle\leq M_{y}K_{yy}P\|y-y_{0}\|
12yy0,yB(y0,P),\displaystyle\leq\frac{1}{2}\|y-y_{0}\|,\forall y\in\mathrm{B}(y_{0},P),

where

P:=min{12MyKyy,R1}.P:=\min\left\{\frac{1}{2M_{y}K_{yy}},R_{1}\right\}.

Using again the Fundamental Theorem of Calculus we get

Dxφ(y,x)\displaystyle D_{x}\varphi(y,x) =Dxφ(y,x0)+01Dx2φ(y,(1s)x0+sx)(xx0)ds\displaystyle=D_{x}\varphi(y,x_{0})+\int_{0}^{1}D_{x}^{2}\varphi(y,(1-s)x_{0}+sx)(x-x_{0})\mathrm{d}s
=Dxφ(y0,x0)+01DyDxφ((1s)y0+sy,x0)(yy0)ds\displaystyle=D_{x}\varphi(y_{0},x_{0})+\int_{0}^{1}D_{y}D_{x}\varphi((1-s)y_{0}+sy,x_{0})(y-y_{0})\mathrm{d}s
+01Dx2φ(y,(1s)x0+sx)(xx0)ds.\displaystyle+\int_{0}^{1}D_{x}^{2}\varphi(y,(1-s)x_{0}+sx)(x-x_{0})\mathrm{d}s.

By triangle inequality and submultiplicative property of induced norms we have

Dxφ(y,x)\displaystyle\|D_{x}\varphi(y,x)\| Dxφ(y0,x0)+01DyDxφ((1s)y0+sy,x0)(yy0)ds\displaystyle\leq\|D_{x}\varphi(y_{0},x_{0})\|+\int_{0}^{1}\|D_{y}D_{x}\varphi((1-s)y_{0}+sy,x_{0})\|\|(y-y_{0})\|\mathrm{d}s
+01Dx2φ(y,(1s)x0+sx)(xx0)ds\displaystyle+\int_{0}^{1}\|D_{x}^{2}\varphi(y,(1-s)x_{0}+sx)\|\|(x-x_{0})\|\mathrm{d}s
My(Lx+Kxyyy0+Kxxxx0),(y,x)B((y0,x0),R2).\displaystyle\leq M_{y}(L_{x}+K_{xy}\|y-y_{0}\|+K_{xx}\|x-x_{0}\|),\ \forall(y,x)\in\mathrm{B}((y_{0},x_{0}),R_{2}).

From the Mean Value Inequality,

φ(y,x)φ(y,x0)\displaystyle\|\varphi(y,x)-\varphi(y,x_{0})\| max0t1Dxφ(y,(1t)x0+tx)xx0\displaystyle\leq\max_{0\leq t\leq 1}\|D_{x}\varphi(y,(1-t)x_{0}+tx)\|\|x-x_{0}\|
My(Lx+Kxyyy0+Kxxxx0)(xx0),(y,x)B((y0,x0),R2).\displaystyle\leq M_{y}(L_{x}+K_{xy}\|y-y_{0}\|+K_{xx}\|x-x_{0}\|)(\|x-x_{0}\|),\ (y,x)\in\mathrm{B}((y_{0},x_{0}),R_{2}).

For a given rmin{P,R2}r\leq\min\{P,R_{2}\}, we need the inequality

φ(y,x)φ(y,x0)r2\|\varphi(y,x)-\varphi(y,x_{0})\|\leq\frac{r}{2}

hold for all yclB(y0,r)y\in\mathrm{cl}\mathrm{B}(y_{0},r), and all xB(x0,ϵ)x\in\mathrm{B}(x_{0},\epsilon). Requring

My(Lx+Kxyr+Kxxϵ)ϵr2M_{y}(L_{x}+K_{xy}r+K_{xx}\epsilon)\epsilon\leq\frac{r}{2}

will suffice. One can solve the inequality above to get a respective

ϵmin{R22r2,14MyKxy}\epsilon\leq\min\left\{\sqrt{R_{2}^{2}-r^{2}},\frac{1}{4M_{y}K_{xy}}\right\}

for a given rmin{P,R2}r\leq\min\{P,R_{2}\}.

Dyφ(y,x)\displaystyle D_{y}\varphi(y,x) =Dyφ(y,x0)+01DxDyφ(y,(1s)x0+sx)(xx0)ds\displaystyle=D_{y}\varphi(y,x_{0})+\int_{0}^{1}D_{x}D_{y}\varphi(y,(1-s)x_{0}+sx)(x-x_{0})\mathrm{d}s
=Dyφ(y0,x0)+01Dy2φ((1s)y0+sy,x0)(yy0)ds\displaystyle=D_{y}\varphi(y_{0},x_{0})+\int_{0}^{1}D_{y}^{2}\varphi((1-s)y_{0}+sy,x_{0})(y-y_{0})\mathrm{d}s
+01DxDyφ(y,(1s)x0+sx)(xx0)ds,yclB(y0,r).\displaystyle+\int_{0}^{1}D_{x}D_{y}\varphi(y,(1-s)x_{0}+sx)(x-x_{0})\mathrm{d}s,\ \forall y\in\mathrm{cl}\mathrm{B}(y_{0},r).

Noting Dyφ(y0,x0)=0D_{y}\varphi(y_{0},x_{0})=0, we have

Dyφ(y,x)\displaystyle\|D_{y}\varphi(y,x)\| 01Dy2φ((1s)y0+sy,x0)(yy0)ds\displaystyle\leq\int_{0}^{1}\|D_{y}^{2}\varphi((1-s)y_{0}+sy,x_{0})\|\|(y-y_{0})\|\mathrm{d}s
+01DxDyφ(y,(1s)x0+sx)(xx0)ds\displaystyle+\int_{0}^{1}\|D_{x}D_{y}\varphi(y,(1-s)x_{0}+sx)\|\|(x-x_{0})\|\mathrm{d}s
MyKyyyy0+MyKxyxx0\displaystyle\leq M_{y}K_{yy}\|y-y_{0}\|+M_{y}K_{xy}\|x-x_{0}\|
34,yclB(y0,r),xB(x0,ϵ).\displaystyle\leq\frac{3}{4},\ \forall y\in\mathrm{cl}\mathrm{B}(y_{0},r),\ \forall x\in\mathrm{B}(x_{0},\epsilon).

From the Mean Value Inequality,

φ(y1,x)φ(y2,x)max0t1Dyφ((1t)y1+ty2,x)y1y234y1y2,\|\varphi(y_{1},x)-\varphi(y_{2},x)\|\leq\max_{0\leq t\leq 1}\|D_{y}\varphi((1-t)y_{1}+ty_{2},x)\|\|y_{1}-y_{2}\|\leq\frac{3}{4}\|y_{1}-y_{2}\|,

for y1y_{1}, y2clB(y0,r)y_{2}\in\mathrm{cl}\mathrm{B}(y_{0},r) and xB(x0,ϵ)x\in\mathrm{B}(x_{0},\epsilon). ∎

Combining Theorem 4.3 and our setting about the central configurations, we have immediately

Theorem 4.4.

Let F=[F1,F2,F3,F4]=[f12,f13,f14,f23]F=[F_{1},F_{2},F_{3},F_{4}]=[f_{12},f_{13},f_{14},f_{23}]. Then FF is a map from 𝐔×(+)3(x,m)\mathbf{U}\times(\mathbb{R^{+}})^{3}\owns(x,m) to 4\mathbb{R}^{4} where the admissible region 𝐔4\mathbf{U}\subset\mathbb{R}^{4} is defined in Theorem 2.3. Since m2m_{2} is fixed to be 11 as the largest value once and for all, the mass vector m(+3)m\in(\mathbb{R^{+}}^{3}) is then m=[m1,m3,m4](0,1]3m=[m_{1},m_{3},m_{4}]\in(0,1]^{3}. Let x0𝐔x_{0}\in\mathbf{U} be a central configuration for m0m_{0} such that F(x0,m0)=0F(x_{0},m_{0})=0 and DxF(x0,m0)D_{x}F(x_{0},m_{0}) is nonsingular. Define

Lm:=DmF(x0,m0),Mx:=(DxF(x0,m0))1L_{m}:=\|D_{m}F(x_{0},m_{0})\|,\hskip 14.22636ptM_{x}:=\|(D_{x}F(x_{0},m_{0}))^{-1}\|

and for given R1>0,R2>0R_{1}>0,R_{2}>0, set

Kxx:=sup{Dx2F(x,m0)|xB(x0,R1)},K_{xx}:=sup\{\|D_{x}^{2}F(x,m_{0})\||x\in B(x_{0},R_{1})\},
Kmm:=sup{Dx2F(x,m)|(x,m)B((x0,m0),R2)},K_{mm}:=sup\{\|D_{x}^{2}F(x,m)\||(x,m)\in B((x_{0},m_{0}),R_{2})\},

and

Kxm:=sup{DxDmF(x,m)|(x,m)B((x0,m0),R2)},K_{xm}:=sup\{\|D_{x}D_{m}F(x,m)\||(x,m)\in B((x_{0},m_{0}),R_{2})\},

Set P:=min{12MxKxx,R1}.P:=\min\left\{\frac{1}{2M_{x}K_{xx}},R_{1}\right\}. Then for all rmin{P,R2}r\leq\min\{P,R_{2}\} and ϵmin{R22r2,14MxKxm}\epsilon\leq\min\{\sqrt{R_{2}^{2}-r^{2}},\frac{1}{4M_{x}K_{xm}}\} satisfying

(22) Mx(Lm+Kxmr+Kmmϵ)ϵr/2,M_{x}(L_{m}+K_{xm}r+K_{mm}\epsilon)\epsilon\leq r/2,

F(x,m)=0F(x,m)=0 has a unique solution on B(x0,r)B(x_{0},r) for any given mB(m0,ϵ)m\in B(m_{0},\epsilon).

Remark 4.5.

In our computation program, in order to have a large existence and uniqueness neighborhood about (x0,m0)(x_{0},m_{0}), we choose r=min{P,R2}r=\min\{P,R_{2}\} and ϵ\epsilon be the maximum value satisfying the inequality (22) or we take

(23) ϵ=min{Mx(Lm+Kxmr)+(Mx(Lm+Kxmr))2+2MxKmmr2MxKmm,R22r2,14MxKxm}.\epsilon=\min\left\{\frac{-M_{x}(L_{m}+K_{xm}r)+\sqrt{(M_{x}(L_{m}+K_{xm}r))^{2}+2M_{x}K_{mm}r}}{2M_{x}K_{mm}},\sqrt{R_{2}^{2}-r^{2}},\frac{1}{4M_{x}K_{xm}}\right\}.

These two values strongly depend on R1R_{1} and R2R_{2}.

We give an example to illustrate our computation program to show that there exists a unique convex central configuration for any masses in the corresponding ϵ\epsilon-neighborhood.

Example 4.6.

Let x0x_{0} be the unique convex central configuration for mass m0m_{0} with m0=[0.2,1,m_{0}=[0.2,1, 0.3,0.3, 0.4]0.4] and x0=[x10,y10,x30,y30]x_{0}=[x_{10},y_{10},x_{30},y_{30}] as given in Example 3.5. To compute rr and ϵ\epsilon, we take

x10=0.1538532870752166,y10=1.4086619698548151,x_{10}=0.1538532870752166,y_{10}=1.4086619698548151,
x30=0.1161115842885355,y30=1.4848782026467043.x_{30}=0.1161115842885355,y_{30}=-1.4848782026467043.

Let R1=0.1R_{1}=0.1 and R2=0.2R_{2}=0.2. By direct computation, Mx=3.193848,M_{x}=3.193848, Lm=0.591292L_{m}=0.591292, Kxx=1.880567K_{xx}=1.880567, Kxm=1.9330632K_{xm}=1.9330632, Kmm=2.160289K_{mm}=2.160289. Then r=0.083247r=0.083247 and ϵ=0.016540.\epsilon=0.016540. Theorem 4.4 affirms that there exists a unique solution x=x(m)x=x(m) in xB(x0,r)x\in B(x_{0},r) for mB(m0,ϵ).m\in B(m_{0},\epsilon).

If we choose R1=0.2R_{1}=0.2 and R2=0.2R_{2}=0.2, then we have r=0.066053,r=0.066053, and ϵ=0.013632.\epsilon=0.013632. This illustrates that the values of rr and ϵ\epsilon depend on the choices of R1R_{1} and R2R_{2}. In our computations, we take R1=0.1R_{1}=0.1 and R2=0.2R_{2}=0.2 to have a larger ϵ\epsilon.

To complete the proof of the uniqueness of the convex central configuration for any mB(m0,ϵ)m\in B(m_{0},\epsilon), it is necessary to prove that any x𝐔¯\B(x0,r)x\in\bar{\mathbf{U}}\backslash B(x_{0},r) is not a solution of the central configuration equations for this mm. We have to change the balls to intervals in order to use interval arithmetic. Let

(m0,ϵ)=[m10ϵ,m10+ϵ]×[m30ϵ,m30+ϵ]×[m40ϵ,m40+ϵ]\mathcal{I}(m_{0},\epsilon)=[m_{10}-\epsilon,m_{10}+\epsilon]\times[m_{30}-\epsilon,m_{30}+\epsilon]\times[m_{40}-\epsilon,m_{40}+\epsilon]

be a mass interval centered at m0m_{0} with radius ϵ\epsilon. Let

(x0,r)=[x10r,x10+r]×[y10r,y10+r]×[x30r,x30+r]×[y30r,y30+r]\mathcal{I}(x_{0},r)=[x_{10}-r,x_{10}+r]\times[y_{10}-r,y_{10}+r]\times[x_{30}-r,x_{30}+r]\times[y_{30}-r,y_{30}+r]

be a position interval centered at x0x_{0} with radius rr. Then

(x0,r2)B(x0,r)\mathcal{I}\left(x_{0},\frac{r}{2}\right)\subseteq B(x_{0},r)

and

(m0,ϵ3)B(m0,ϵ)(m0,ϵ).\mathcal{I}\left(m_{0},\frac{\epsilon}{\sqrt{3}}\right)\subseteq B(m_{0},\epsilon)\subseteq\mathcal{I}(m_{0},\epsilon).

It is easy to prove the following lemma and its proof is omitted.

Lemma 4.7.

Suppose there exists a unique solution x=x(m)x=x(m) in xB(x0,r)x\in B(x_{0},r) for mB(m0,ϵ).m\in B(m_{0},\epsilon). If any x𝐔¯\(x0,r2)x\in\bar{\mathbf{U}}\backslash\mathcal{I}(x_{0},\frac{r}{2}) is not a solution for any m(m0,ϵ)m\in\mathcal{I}(m_{0},\epsilon), then there exists a unique solution x=x(m)x=x(m) in 𝐔\mathbf{U} for mB(m0,ϵ)m\in B(m_{0},\epsilon). If any x𝐔¯\(x0,r2)x\in\bar{\mathbf{U}}\backslash\mathcal{I}(x_{0},\frac{r}{2}) is not a solution for any m(m0,ϵ3)m\in\mathcal{I}(m_{0},\frac{\epsilon}{\sqrt{3}}), then there exists a unique solution x=x(m)x=x(m) in 𝐔\mathbf{U} for m(m0,ϵ3)m\in\mathcal{I}(m_{0},\frac{\epsilon}{\sqrt{3}}).

To show that any x𝐔¯\(x0,r2)x\in\bar{\mathbf{U}}\backslash\mathcal{I}(x_{0},\frac{r}{2}) is not a solution for any m(m0,ϵ)m\in\mathcal{I}(m_{0},\epsilon), we first divide the interval 𝐔¯\(x0,r2)\bar{\mathbf{U}}\backslash\mathcal{I}(x_{0},\frac{r}{2}) into smaller ones [x][x] and then we check that at least one of fij([x],(m0,ϵ))f_{ij}([x],\mathcal{I}(m_{0},\epsilon)) doesn’t contain zeros. We need to make sure that the union of [x][x] covers the whole region out of the position ball. This shows that the mass ball B(m0,ϵ)B(m_{0},\epsilon) is a uniqueness mass ball. In Example 4.6, (m0,ϵ)=[0.183460297345220,0.216539702654780]×\mathcal{I}(m_{0},\epsilon)=[0.183460297345220,0.216539702654780]\times
[0.283460297345220,0.316539702654780]×[0.383460297345220,0.416539702654780][0.283460297345220,0.316539702654780]\times[0.383460297345220,0.416539702654780]
and x0=[0,0.112229923080813]×\mathcal{I}_{x_{0}}=[0,0.112229923080813]\times [0.268,1.733]×[0.268,1.733]\times [0,1]×[0,1]\times [1.733,0.268][-1.733,-0.268] is one of the eight intervals in the closure of 𝐔¯\(x0,r2)\bar{\mathbf{U}}\backslash\mathcal{I}(x_{0},\frac{r}{2}). But the six functions fij(x0,(m0,ϵ))f_{ij}(\mathcal{I}_{x_{0}},\mathcal{I}(m_{0},\epsilon)) all contain zeros. We split the position interval x0\mathcal{I}_{x_{0}} into two subintervals by bisecting the longest sides and then check whether fijf_{ij} contains zeros on the subintervals. Repeat the process until it proves that at least one of fijf_{ij} doesn’t contain zeros in the subintervals. By this way, it confirms that B(m0,ϵ)B(m_{0},\epsilon) is a uniqueness mass ball in 𝐔\mathbf{U}.

Note that one may worry that there exists an mm in B(m0,ϵ)B(m_{0},\epsilon) near boundary whose solution xx is also near the boundary of B(x0,r)B(x_{0},r). Then such xx falls in 𝐔¯\(x0,r2)\bar{\mathbf{U}}\backslash\mathcal{I}(x_{0},\frac{r}{2}). If such case exists, our method can not confirm the uniqueness of B(m0,r)B(m_{0},r). But we can apply our method to a smaller mass intervals (m0,ϵN)B(m0,ϵ)\mathcal{I}(m_{0},\frac{\epsilon}{N})\subseteq B(m_{0},\epsilon) for an appropriate N3N\geq\sqrt{3} and prove B(m0,ϵN)B(m_{0},\frac{\epsilon}{N}) is a uniqueness mass ball. Fortunately, this case doesn’t occur in our program. All mass balls B(m0,ϵ)B(m_{0},\epsilon) in our program are uniqueness mass balls. To accelerate the numerical computations, we can produce the uniqueness mass interval (m0,ϵ3)\mathcal{I}(m_{0},\frac{\epsilon}{\sqrt{3}}) to cover the mass space by checking that at least one of fij([x],(m0,ϵ3))f_{ij}([x],\mathcal{I}(m_{0},\frac{\epsilon}{\sqrt{3}})) doesn’t contain zeros for all small intervals [x][x] out of the corresponding position ball.

5 Proof of the Theorem 1.1

To prove the uniqueness of convex central configurations in Theorem 1.1, we create finite many uniqueness mass balls B(m0,ϵ)B(m_{0},\epsilon) so that the union of such balls can cover the whole mass space [0.6,1]3[0.6,1]^{3}. Here we describe the procedures to produce the uniqueness mass balls.

Step 1: Divide the mass space [0.6,1]3[0.6,1]^{3} into a grid with equal side length [0.6,1][0.6,1] for M1M_{1} segments, such that mik=0.6+0.4kM1m_{ik}=0.6+\frac{0.4k}{M_{1}} for i=1,3,4i=1,3,4 and k=0,1,2,,M1k=0,1,2,...,M_{1}. The value of M1M_{1} is determined based on experimental results. We choose M1=40.M_{1}=40.

Step 2: For each m0=[m1k1,1,m3k3,m4k4]m_{0}=[m_{1{k_{1}}},1,m_{3{k_{3}}},m_{4{k_{4}}}] at a meshed grid point in Step 1, applying the algorithm (interval arithmetic computation and Krawczyk operator) in Section 3, we can prove there is a unique central configuration x0x_{0}. Now the proof of Theorem 2.4 is complete.

Step 3: For such central configuration x0x_{0} for m0m_{0}, applying the implicit function theorem in Section 4, we can find a uniqueness mass ball B(m0,ϵ)B(m_{0},\epsilon), where ϵ\epsilon is a small positive number such that the central configuration is unique for all masses in the ball B(m0,ϵ)B(m_{0},\epsilon) as guaranteed by Theorem 4.4 and Lemma 4.7.

Step 4: Check whether the union of all the uniqueness balls covers the mass space [0.6,1]3[0.6,1]^{3}. This follows from the fact that the compact mass space [0.6,1]3[0.6,1]^{3} can be covered by a finitely many uniqueness mass balls and the fact that we can refine our mesh points so that the uniqueness balls overlap ensuring that any mass point is in at least one uniqueness ball.

Therefore, by Step 2 and Step 3, there exists a unique central configuration for any mass point in the meshed space, and by Step 4, this uniqueness extends to the entire mass space [0.6,1]3[0.6,1]^{3}. Thus, Theorem 1.1 is proven.

Initially, set M1=5M_{1}=5. We find that the radius ϵ\epsilon of the uniqueness mass ball is between 0.00921190.0092119 and 0.0163930.016393 with an overall average radius about 0.012584360.01258436. If the minimum of the radius of the uniqueness mass balls is ϵ0\epsilon_{0}, we can cover the whole mass space by choosing M1>4320ϵ0M_{1}>\frac{4\sqrt{3}}{20\epsilon_{0}} since a uniqueness ball B(m0,ϵ0)B(m_{0},\epsilon_{0}) contains an interval centered at m0m_{0} with each side length 23ϵ03.\frac{2\sqrt{3}\epsilon_{0}}{3}. If ϵ0=0.009\epsilon_{0}=0.009, M1>38M_{1}>38. In our program, we use M1=40M_{1}=40 and totally 6400064000 uniqueness mass balls are computed to cover the whole mass space [0.6,1]3[0.6,1]^{3}. According to Lemma 4.7, our program confirms the uniqueness of the mass ball B(m0,ϵ)B(m_{0},\epsilon) for each mass m0m_{0} at the mesh grid in Step 1.

6 Conclusions and Remarks

By combining the interval analysis and implicit function theorem, we give a computer assisted proof about the uniqueness of convex central configuration for any given mass vector m[0.6,1]3m\in[0.6,1]^{3} in a fixed order in the planar 44-body problem. This closed domain can be enlarged and improved further without modifying our method and program. So the conjecture is hopefully true for all masses except for those near the mass vector with one or more component being zeros.

However we find that it becomes more and more time consuming when we approach the boundary of the mass space because ϵ\epsilon becomes smaller and smaller. This suggests the potential singular behavior when we take the limit to the boundary of the mass space.

For example, let’s consider the mass space ([0.1,1])3([0.1,1])^{3} and let M1=10M_{1}=10. We determine that the radius ϵ\epsilon of the uniqueness mass ball lies between 0.00014360.0001436 and 0.021250.02125, with an average radius of approximately 0.009984450.00998445. By using the minimum radius ϵ0=0.0001436\epsilon_{0}=0.0001436, we can cover the entire mass space with M1>5428M_{1}>5428 and a total number of uniqueness mass balls equal to approximately 1.6×10111.6\times 10^{11}. It will take a rather long time to finish the computation using the current program, and it seems unpractical without further refinement of our current program.

Reflecting further on our methods, we find the following mathematical problems may be interesting. (1) Given a function ff on a closed domain DkD\subset\mathbb{R}^{k}, for each point qDq\in D, we can find an open ball B(q,rq)B(q,r_{q}) with radius depending on the point qq and some property of the function ff (e.g. the valid domain for implicit function theorem in the current paper). In such a way, we get an open covering of the compact set DD, and in principle we can select a finite number of them to cover DD. Now the question is whether one can find an optimal one as in sphere packing problem. Here we are deliberately vague to allow refinements and variations. (2) Even more interesting problem would be trying to refine and optimize the size of the valid domain for the implicit function theorem with or without the conditions about the second order derivatives of the mapping. There is very few literature concerning this fundamental topic, and we leave to future research.

Yet from this perspective we also notice the following interesting result establishing the conjecture about the uniqueness of the convex central configurations for masses m1,m3,m4m_{1},m_{3},m_{4} in the following form near zeros in a given order.

Lemma 6.1.

([8], Theorem 2) Given μi>0(i=1,3,4)\mu_{i}>0(i=1,3,4), there exists ϵ0>0\epsilon_{0}>0 such that for 0<ϵ<ϵ0,0<\epsilon<\epsilon_{0}, the four masses m=[ϵμ1,m=[\epsilon\mu_{1}, 1,1, ϵμ3,\epsilon\mu_{3}, ϵμ4]\epsilon\mu_{4}] admit a unique convex central configuration for a given order.

Now we cannot merge the statement into our framework because ϵ0\epsilon_{0} is theoretical and not given a concrete numerical value. Note also there are other possibilities to investigate for masses approaching zero. It seems still a challenging problem even numerically to verify the conjecture.

The final remark is that even we can prove the conjecture with computer assistance, it is still desirable to have a conceptual/analytical proof.

Acknowledgements

We would like to thank the valuable suggestions from Professor Piotr Zgliczynski and the anonymous referee.

Shanzhong Sun is partially supported by National Key R&D Program of China (2020YFA 0713300), NSFC (No.s 11771303, 12171327, 11911530092, 11871045). Zhifu Xie is partially supported by Wright W. and Annie Rea Cross Endowment Funds at the University of Southern Mississippi.

References

  • [1] A. Albouy, Symétrie des configurations centrales de quatre corps, C. R. Acad. Sci. Paris 320 (1995), no.1 217–220.
  • [2] A. Albouy, The symmetric central configurations of four equal masses, Contemp. Math. 198 (1996), 131–135.
  • [3] A. Albouy, H. E. Cabral and A. A. Santos, Some problems on the classical nn-body problem, Celest. Mech. Dyn. Astron. 113 (2012), no.4, 369–375.
  • [4] A. Albouy, Y. Fu and S. Sun, Symmetry of planar four-body convex central configurations, Proc. R. Soc. A 464 (2008), 1255–1365.
  • [5] A. Albouy and V. Kaloshin, Finiteness of central configurations of five bodies in the plane, Ann. Math., 176 (2012), 1–54.
  • [6] G. Alefeld and G. Mayer, Interval analysis: theory and applications, Journal of Computational and Applied Mathematics 121 (2000), 421-464.
  • [7] M. Corbera, J. M. Cors, J. Llibre and E. Pérez-Chavela, Trapezoid central configurations, Appl. Math. Comput. 346 (2019), 127–142.
  • [8] M. Corbera, J. M. Cors, J. Llibre and R. Moeckel, Bifurcation of relative equilibria of the (1+3)(1+3)-body problem, SIAM J. Math. Anal., 47 (2015), 1377-1404.
  • [9] M. Corbera, J. M. Cors and G. E. Roberts, A four-body convex central configuration with perpendicular diagonals is necessarily a kite, Qual. Theory Dyn. Syst. 17 (2018), no.2, 367–374.
  • [10] M. Corbera, J. M. Cors and G. E. Roberts, Classifying four-body convex central configurations, Celest. Mech. Dyn. Astr. 131 (2019), 34.
  • [11] J. M. Cors and G. E. Roberts, Four-body co-circular central configurations, Nonlinearity 25 (2012), 343–370.
  • [12] A. C. Fernandes, J. Llibre and L. F. Mello, Convex central configurations of the four-body problem with two pairs of equal adjacent masses, Arch. Ration. Mech. Anal. 226 (2017), no.1, 303–320.
  • [13] Y. Hagihara, Celestial Mechanics, Volume I, Dynamical Principles and Transformation Theory, The MIT Press, 1970.
  • [14] M. Hampton and R. Moeckel, Finiteness of relative equilibria of the four body problem, Invent. Math., 163 (2006), 289–312.
  • [15] A. Jindal, R. Banavar and D. Chatterjee, Implicit function theorem: estimates on the size of the domain, 31 May 2022 (v2), arXiv:2205.12661v2.
  • [16] T. Lee and M. Santoprete, Central configurations of the five-body problem with equal masses, Celest. Mech. Dyn. Astr. 104 (2009), no.4, 369–381.
  • [17] Y. Long, Admissible shapes of 4-body non-collinear relative equilibria, Advanced Nonlinear Studies 1 (2003), 495-509.
  • [18] Y. Long and S. Sun, Four-body central configurations with some Equal Masses, Arch. Ration. Mech. Anal. 162 (2002) 25–44.
  • [19] W. MacMillan and W.D. Bartky, Permanent configurations in the problem of four bodies, Trans. Amer. Math. Soc. 34(1932), no.4, 838–875.
  • [20] MATLAB, version R2016b (9.1.0.441655) The MathWorks, Inc., Natick, Massachusetts, United States.
  • [21] R. Moeckel, On central configurations, Math. Z. 205 (1990), 499-517.
  • [22] R. Moeckel, Central configurations, in Central configurations, periodic orbits, and Hamiltonian systems, 105–167, Adv. Courses Math. CRM Barcelona, Birkhäuser/Springer, Basel, 2015.
  • [23] M. Moczurad and P. Zgliczynski, Central configurations in planar nn-body problem with equal masses for n=5,6,7n=5,6,7, Celest. Mech. Dyn. Astr. 131 (2019), 46.
  • [24] R. E. Moore, R. Baker Kearfott, Michael J. Cloud, Introduction to Interval Analysis, SIAM, 2009, ISBN: 978-0-898716-69-6.
  • [25] SageMath, the Sage Mathematics Software System (Version 9.6), The Sage Developers. http://www.sagemath. org (May 2022). Downloaded July 28, 2022.
  • [26] M. Santoprete, Four-body central configurations with one pair of opposite sides parallel, J. Math. Anal. Appl. 464(2018), 421–434.
  • [27] M. Santoprete, On the uniqueness of co-circular four body central configurations, Arch. Ration. Mech. Anal. 240 (2021), no. 2, 971–985.
  • [28] M. Santoprete, On the uniqueness of trapezoidal four-body central configurations, Nonlinearity 34 (2021), no. 1, 424–437.
  • [29] D. Schmidt, Central configurations and relative equilibria for the nn-Body Problem, Classical and Celestial Mechanics (Recife, 1993/1999), pp. 1–33. Princeton University Press, Princeton (2002)
  • [30] J. Shi and Z. Xie, Classification of four-body central configurations with three equal masses, J. Math. Anal. Appl., 363 (2010), 512–524.
  • [31] C. Simo, Relative equilibrium solutions in the four body problem, Celestial Mechanics 18 (1978), 165-184.
  • [32] S. Smale, Mathematical problems for the next century, Math. Intelligencer 20 (1998), no.2, 7–15.
  • [33] Z. Xia, Convex central configurations for the nn-body problem, J. Differential Equations 200 (2004) no.1 185-190.
  • [34] Z. Xie, Isosceles trapezoid central configurations of the Newtonian four-body problem, Proc. R. Soc. Edinb. Sect. A 142 (2012), no.3, 665–672.