Effective Boundary Value Problem Solver via Bézier Functions

: In engineering disciplines, many important problems are to be formed as boundary value problems (BVPs) that have conditions that are speciﬁed at the extremes. To handle such problems, the conventional shooting method that transforms BVPs into initial value problems has been extensively used, but it does not always guarantee solving the problem due to the possible failure of ﬁnding a proper initial guess. This paper proposes a universal initial guess ﬁnder that is composed of Bézier functions. Various dimensional problems that include Lambert’s problem for several orbits around the spherically symmetric Earth are studied to validate the efﬁcacy of the proposed approach, and the results are compared.


Introduction
Boundary value problems (BVPs) are important problems in engineering fields because plenty of problems are formed as BVPs that are given by the differential equations with boundary conditions (BCs). In celestial mechanics, for example, Lambert's problem is the classical problem that is to determine an orbit by solving BVPs [1,2]. The goal of this problem is to find an orbit between two prescribed positions within a given fixed time of flight. The orbit determination, which solves Lambert's problem, is usually utilized for spacecraft intercept [3], rendezvous [4], orbit transfer [5], etc. A general solving approach for BVPs is the shooting method that finds unknown initial values that satisfy the given BCs. For orbit determination, it requires finding the initial velocity that corresponds to the prescribed position at the initial time. For applying the shooting method, the main challenge is to determine a suitable initial guess, because it may fail to find the solution depending on the starting guess selected. Developing an effective approach for selecting a proper initial guess is essential, because the initial guess affects not only the convergence rate of the solution, but also the success of obtaining the solution. Hence, this research proposes a novel approach to provide an appropriate initial guess for solving BVPs.
Bézier functions (BFs) are parametric functions that can describe complicated surfaces and curves with great flexibility [6][7][8]. These parametric functions are composed of linear combinations of polynomials and the coefficients that are called the control points, which can adjust the shape of curves. BFs have been investigated for describing trajectories in a variety of studies: airfoil design [9], path planning [10][11][12][13][14][15], collision avoidance [16][17][18][19], trajectory generation [20], etc. Besides, studies that are associated with the celestial mechanics by utilizing BFs are conducted by several researchers. Dilectis and Mortari [21] used BFs to describe orbit trajectories by only utilizing observations. Pan and Ma designed a formulation of Lambert's problem with respect to the argument of periapsis by using BFs [22]. Mortari and Clocchiatti [23] utilized non-rational BFs (NBFs) by approximating Kepler's equation for the entire mean anomaly range, and a root-finding tool by using BFs for solving Kepler's equation was proposed by Mortari and Elipe [24]. In addition, Kim b n,k (s) = n k s k (1 − s) n−k , where ( n k ) is a binomial coefficient, k = 0, 1, · · · , n, and s ∈ [0, 1] is the implicit parameter. BPs have several properties as follows: [26][27][28] 1.
BPs b n,k (s) and b n,n−k (s) are mirror images of each other about the interval mid-point s = 1/2 (i.e., b n,n−k (1 − s) = b n,k (s)); 2.
there is the recurrent relation between BPs of degree n + 1 and n as b n+1,k (s) = (1 − s)b n,k + sb n,k−1 ; 3.
they can always be written as a linear combination of polynomials of higher degree, b n−1,k (s) = n − k n b n,k (s) + k + 1 n b n,k+1 (s); and, 8. they span the full space of polynomials. This means that any n-degree polynomial can be expressed as a linear combination of n-degree BPs.
Note that BPs are the parametric basis function of both NBFs and RBFs.
A NBF of degree n is defined as [8,28] p(s) = n ∑ k=0 b n,k (s)c k , (2) where c k is the k-th control point. A set of control points (n + 1 number of points) is defined from c 0 through c n . The first and last control points (c 0 and c n ) are called endpoints that indicate the start and end points of the curve. While the endpoints are located on the curve, the intermediate control points (c 1 , . . . , c n−1 ) do not lie on the curve [29]. Instead, the intermediate control points mainly control the shape of the curve. Because the endpoints lie on each end of the curve, one or more intermediate control points are required to generate a curve. That is, at least a quadratic BF is needed to approximate a curve. The r-th derivative of the NBF of degree n with respect to s is recursively defined as [28] where c (r) An RBF of degree n is defined as [8,28] where w k is the relative weight that is associated with k-th control point. By adding the weights, the RBFs have more flexibility than the NBFs because the RBFs have more parameters to be adjusted. That is, the RBFs can closely approximate arbitrary curves by adjusting both the weights and the control points. In order to determine the derivative of the RBF, the numerator and denominator are redefined as p n (s) and p d (s), and then Equation (5) is rewritten as p(s) = p n (s)/p d (s). Subsequently, the r-th derivative of the RBF of degree n with respect to s is recursively derived as [28] where the derivatives of p n (s) and p d (s) are obtained by using Equations (3) and (4).

Shooting Method
A general form of nonlinear BVPs is given by that satisfies the BCs x(t i ) = x i and x(t f ) = x f . Here, t ∈ [t i , t f ] is the time span that is set with the initial and final time t i and t f , x(t) ∈ R m is the state vector, g ∈ R m and f ∈ R m are the nonlinear functions, and the subscripts i and f indicate the initial and final parameters, which also represent the lower and upper bounds. Note that g is to be called the residual. The BVP is generally solved by using the shooting method that finds the solution of the given problem by using the solutions to a sequence of initial value problems (IVPs), and this problem is for solving Equation (7) with x(t i ) = x i andẋ(t i ) = z k . Here, the unknown initial state z k is initially set to be z 0 . Afterwards, z k is selected in a manner to ensure that [30] lim where x(t f , z k ) denotes the state at t f when z k is applied to the IVP in Equation (7), andẋ(t i ) indicates the solution to the BVP in Equation (7) with BCs. An initial guess z 0 is utilized to find the solution x(t f , z k ) that satisfies

Initial Guess Finder
An initial guess is required to solve the BVP using the shooting method, since the BVP generally provides the boundary values as initial conditions, as mentioned in Section 2.2. That is, the trajectories of the states for the BVP should contain the boundary values at the end of the trajectories. Because the endpoints of BFs should lie on the curve similarly, one can say that the BCs correspond to the endpoints. From this concept, the given BVP can be approximated by BFs, and the initial guess is obtained from the approximated trajectory at the lower bound. This work considers two types of BFs to determine the initial guess in the initial guess finder (IGF). The first IGF uses the NBFs while the second one uses the RBFs, and they are named NBF-IGF and RBF-IGF, respectively.
For the NBF-IGF, it is necessary to transform the boundary values for the given BVP into the endpoints. That is, the initial and final time in the given problem are replaced with t B 0 ∈ R and t B n ∈ R, and the initial and final states are replaced with x B 0 ∈ R m and x B n ∈ R m in BFs. Note that the unknown intermediate control points are defined as t B i ∈ R and x B i ∈ R m for i = 1, · · · , n − 1. Thus, the total number of unknown parameters for the NBF-IGF is (m + 1) × (n − 1). Using Equation (2), t and x(t) are transformed into NBFs, as follows: x B (s) =    (1 − s) n x 1,B 0 + ( n 1 )(1 − s) n−1 sx 1,B 1 + · · · + ( n n−1 )(1 − s)s n−1 x 1,B n−1 + s n x 1,B n . . .
For the RBF-IGF, it should consider additional unknown weights (w t k ∈ R and w x k ∈ R m for k = 0, · · · , n). It is important to note that any endpoint weight values can be used because the intermediate weights are determined relatively. Thus, the weights corresponding to the endpoints can be defined as 1 for simplicity (i.e., w t 0 = w t n = w j,x 0 = w j,x n = 1 for j = 1, ..., m), and the number of unknown weights is (m + 1) × (n − 1). Hence, the total number of unknown parameters for the RBF-IGF is 2 × (m + 1) × (n − 1). Using Equation (5), t and x(t) are transformed into RBFs, as follows: (1 − s) n w m,x 0 x m,B 0 + ( n 1 )(1 − s) n−1 sw m,x 1 x m,B 1 + · · · + ( n n−1 )(1 − s)s n−1 w m,x n−1 x m,B n−1 + s n w m,x n x m,B n (1 − s) n w m,x 0 + ( n 1 )(1 − s) n−1 sw m,x 1 + · · · + ( n n−1 )(1 − s)s n−1 w m,x n−1 + s n w m, To express the corresponding nonlinear functions of the given BVP, the first and second time derivatives of x B (s) are derived aṡ where s in each variable is omitted for the simple notation. Using Equations (14) and (15), g defined in Equation (7) is transformed into BF form as where f B is the transformed nonlinear function that is composed of BFs. Here, g B and f B in the NBF-IGF are the function of s, t B 1 , . . . , t B n−1 , x B 1 , . . . , x B n−1 , and g B and f B in the RBF-IGF are the function of s, t B 1 , . . . , t B n−1 , x B 1 , . . . , x B n−1 , w t 1 , · · · , w t n−1 , w x 1 , · · · , w x n−1 .
To obtain an approximated trajectory, the transformed residual should be zero in order to satisfy Equation (7), and this is achieved by adjusting the unknown parameters. Once the unknown parameters are properly determined, the approximated trajectory can be easily found. Therefore, this becomes finding the proper unknown parameters that minimize the following cost function [31]: Note that the difference between the NBF-IGF and RBF-IGF arises in the first step that is the transformation of state variables, t B (s) and x B (s), using Equations (2) and (5). Depending on the BFs, the number of unknown variables is different. Once the state variables are transformed into each BF, the following steps in Equations (14)-(17) remain the same. It is evident that the RBF-IGF provides a more accurate initial guess, because it has parametric flexibility that can more precisely approximate a curve. However, it is computationally more expensive than the NBF-IGF. Subsequently, this problem becomes to solve Equation (17) with the unknown parameters and the following equations [25]: for i = 1, · · · , n − 1. Note that, if f B is made of products only in Equation (16), g B becomes a polynomial equation in s and Equation (17) can be solved in a closed form [25]. Once the obtained control points (and weights), which are the solution of the optimization problem, are applied to each state variable expressed as BFs, the approximated trajectory of the given BVP is found, and the initial guess for solving the given BVP is determined byẋ B (s)| s=0 . The given problem is generally solved by using the shooting method with random initial guesses because of the lack of information for the initial value, as shown in Figure 1a. If a proper initial guess is selected, then the shooting method provides the initial value for the lower bound that satisfies the given BCs. Subsequently, the numerical integration is performed by using the obtained initial value to find the solution, which is the state trajectory. If the BCs are not satisfied for the obtained state trajectory, different initial guesses should be selected. In fact, many engineering problems are sensitive to initial guesses, and this makes engineers spend a huge effort to find the proper initial guess. The IGFs that were proposed in this work provide a proper initial guess, even for the problems that are sensitive to initial guesses because the information of initial guess is based on the approximated trajectory for the given problem. Figure 1b shows the entire procedure for solving BVPs using the IGF that is highlighted with blue boxes. Initially, the state variables and equations are transformed into Bézier functions. Subsequently, the optimal control points and weights are obtained by solving the optimization problem that minimizes the cost function defined in Equation (17). Finally, using the obtained control points and weights, the initial guess is computed and applied to the shooting method. In fact, the initial value is relatively easily and efficiently found by using the initial guess that is provided by the proposed approach as compared to using random initial guesses. As shown in Figure 1a, the conventional approach iteratively selects random initial guesses to find the proper solution if the obtained states at the final time do not satisfy the BCs. On the other hand, replacing the green box shown in Figure 1a with the dashed blue box in Figure 1b eliminates such an iterative process.

Simulation Study
Numerical simulation studies are performed for several problems to validate the performance of the proposed IGF approaches. The shooting method is utilized to solve the given BVPs, and the simulation results using an initial guess obtained by each IGF are compared with the one using randomly selected initial guesses.
In the celestial mechanics, for example, the shape of the orbit trajectories is ellipse, parabola, or hyperbola, and a part of the entire elliptical orbit is generally the solution trajectory of Lambert's problem. Therefore, the second order BFs can sufficiently approximate the orbital trajectory. Quadratic BFs are used in the simulation studies since the aim of this research is to find an initial guess to efficiently solve BVPs.
For the NBF-IGF, the quadratic NBFs and their first and second derivatives are given as where j = 1, . . . , m. Similarly, the quadratic RBFs and their first and second derivatives for the RBF-IGF are given as For the numerical studies, it considers 1D and 2D nonlinear BVPs, and then 3D nonlinear BVPs, in particular the orbit determination problems, are considered. For the shooting method, the bisection method [30,32] is used for the 1D example, because the 1D problem is identical to a root-finding problem. In the 2D and 3D examples, the trust-region algorithm [33] is utilized to solve systems of nonlinear equations. In addition to this, the simplex algorithm [34] is used to find the optimal control points and weights. These algorithms are run on Octave, which is free software, and the specifications of the computer for the simulation studies are as follows: Intel i7-7500U @ 2.7 GHz clock speed CPU and 12 GB memory.

Nonlinear 1D BVP
A nonlinear 1D BVP is given by [35] where BCs are x(t i = 1) = 0 and x(t f = 2) = ln 256. At first, the given problem is solved using the shooting method with random initial guesses as the conventional approach. It considers three kinds of ranges for the initial guesses due to the lack of the initial value information and performing the Monte Carlo simulations with 1000 iterations for each range. Here, the initial guesses are randomly selected in the normal distribution (e.g., N µ, σ 2 , where µ is the mean and σ is the standard deviation). The obtained initial value for each initial guess is displayed as the same value because the number obtained is the same up to the 14-th decimal place, as shown in Table 1. Regardless of the initial guesses, it finds the solution with almost zero mean of the absolute error of x(t f ). Because the initial value to be obtained is near zero, the case for N 0, 1 2 requires the least computational burden. Furthermore, it is easy to find the solution with an arbitrary initial guess, because the initial value to be obtained is scalar, not a vector in this problem. Next, the proposed approach is applied to solve the given BVP. For this, the given equations are transformed into the quadratic NBFs and RBFs using Equations (19)-(22).

NBF-IGF
Because the endpoints are defined as t B 0 (s = 0) = 1, t B 2 (s = 1) = 2, x B 0 (s = 0) = 0, and x B 2 (s = 1) = ln 256, the number of unknown parameters is 2 (e.g., t B 1 and x B 1 ). Using Equation (16), the residual equation is transformed into where each variable is defined in Equations (19) and (20). Because the given BVP should satisfy the zero residual, the approximated trajectory can be found as solving the optimization problem that minimizes the cost that is defined in Equation (17). Once the optimal intermediate control points are determined, the approximated solution trajectory over s ∈ [0, 1] is obtained, and the initial guess is easily computed at s = 0 for Equation (14).

RBF-IGF
The endpoints are determined by the boundary values similar to NBF-IGF, and the weights that correspond to the endpoints are defined as 1 (i.e., w t 0 = w t 2 = w 1,x 0 = w 1,x 2 = 1) as mentioned in Section 2.3. Therefore, the unknown variables are the intermediate control points and weights, and the number of unknowns is 4. By substituting the transformed variables and their derivatives in Equations (21) and (22) into Equation (24), the residual is obtained. Subsequently, the residual is applied to the cost function for the optimization problem in Equation (17). After solving the optimization problem, the initial guess is calculated by substituting the obtained control points and weights into Equation (14) at s = 0.
In the optimization process, two kinds of stopping criteria are used: (i) the step size defined as the norm of the difference between the parameters at the current step and the previous step and (ii) the change in the value of the cost function during the current step and previous step. When both of the stopping criteria are satisfied, it terminates the optimization process provides the optimal control points and weights. Table 2 lists the initial guesses that are obtained by the proposed approaches and the corresponding simulation results applying the shooting method. Here, (B) represents the computational time for finding the control points and weights, and (S) means the computational time for performing the shooting method to find the proper initial value that satisfies the given BVP by using the initial guess obtained by the proposed approaches. The RBF-IGF provides the closer initial guess to the actual initial value than the NBF-IGF, because the RBFs have more flexibility (due to the additional weight parameters) than the NBFs. For this reason, a small number of iterations and evaluations is required to find the initial value that satisfies the given upper bound value with almost zero error. That is, the error is close to the machine epsilon value when applying the initial guess that is obtained from the RBF-IGF. Because the solution of the shooting method using the initial guess of the RBF-IGF is precisely computed, the obtained state at the upper bound is almost same as the given BC. It also explains that the trial trajectories of the shooting method using the initial guess obtained by the RBF-IGF are very close to the solution trajectory, as shown in Figure 2b, unlike the trial trajectories generated via the NBF-IGF shown in Figure 2a.  On the other hand, applying the initial guess that is obtained by the NBF-IGF into the shooting method is computationally more efficient than applying the other one because the NBF-IGF has only two unknown parameters to be optimized. In fact, the aim of this research is to provide an effective way to find the initial guess, not the actual initial value. For this reason, in between two IGFs, it can be said that the NBF-IGF is better to find the initial guess for solving the given BVP, because it is computationally inexpensive, and the initial guess is close enough to the initial value to be determined.
The IGF reduces a large amount of the computational time for applying the shooting method thanks to the properly selected initial guess, as shown in Tables 1 and 2. When compared the computational time and the accuracy for the absolute value of the upper bound, the one using N 0, 1 2 requires less computational time than the one that is required for the proposed approach, as shown in Figure 3, due to the fact that the initial value is close to zero.
However, it is shown that the NBF-IGF applied result provides better accuracy and less computational time overall.

Nonlinear 2D BVP
A nonlinear 2D BVP is given by [35] x where The given problem is first solved using the shooting method with random initial guesses, and the results are listed in Table 3. Table 3. Monte Carlo simulation results for 2D example using random initial guesses. Here, N m (µ, Σ σ 2 ) indicates a m-dimensional normally-distributed random vector with the mean vector µ ∈ R m and the covariance matrix Σ σ 2 ∈ R m×m . In this work, the covariance matrix is defined as Σ σ 2 = σ 2 I m×m . Because the obtained initial value for each component is the same up to the 7-th decimal place, these values are displayed as the same values. Depending on the covariance of initial guesses, the number of cases that found the solution is differently obtained. In other words, this problem requires a proper initial guess to find the solution, unlike the 1D example. For the cases that found the solutions, the mean values, such as the mean of the number of iterations, number of evaluations, computational time, and upper bound norm error, are summarized in Table 3. That is, the cases that failed to find the solution are not included. The results show that the computational burden, such as the computational time and the number of iterations and evaluations, is reduced as the covariance of the initial guess decreases. In fact, these results are reasonable, because the initial guess of N 2 (0, Σ 1 2 ) is closer to the initial value obtained than the others.

Initial Guess
Next, to apply the proposed approach, the variables and nonlinear functions for the given problem are transformed into the quadratic NBFs and RBFs while using Equations (19)-(22).

NBF-IGF
The given BCs are replaced with the endpoints as t B 0 (s = 0) = 1, t B 2 (s = 1) = 2, x B 0 (s = 0) = [2, 0.5] T , and x B 2 (s = 1) = [2.5, 1/3] T . Hence, three unknown control points are required to be determined, which are t B 1 and x B 1 = [x 1,B 1 , x 2,B 1 ] T . Using Equation (16), the residual and the nonlinear function on the right side of Equation (25) are transformed into the NBFs as where each variable is defined in Equations (19) and (20), and the cost function is derived by substituting Equation (26) into Equation (17). Subsequently, the initial guess is obtained by using the same process of the 1D example.

RBF-IGF
The weights that correspond to the endpoints are set as w t 0 = w t 2 = w 1,x 0 = w 1,x 2 = w 2,x 0 = w 2,x 2 = 1 similar to the 1D example. Therefore, six unknown variables, such as three control points and three weights (w t 1 and w x 1 = [w 1,x 1 , w 2,x 1 ] T ), are considered. The residual is derived by applying the transformed variables and their derivatives in Equations (21) and (22) into Equation (26). The unknown variables are obtained by solving the optimization problem that minimizes the cost in Equation (17) composed of the residual, and the initial guess is computed by applying the obtained control points and weights into Equation (14) at s = 0.
The results of the shooting method with the initial guesses obtained via the IGFs are listed in Table 4. The initial guess generated via the NBF-IGF finds the initial value after one trial shown in Figure 4a. On the other hand, the initial guess that is obtained by the RBF-IGF requires only one iteration to find the initial value. Because the RBFs provided the closer initial guess to the initial value, no trial trajectory is required, especially for this example, as shown in Figure 4b.
On the other hand, applying the NBF-IGF provides a huge benefit in the aspect of the computational efficiency. In addition to this, the upper bound norm error is obtained as a very small value that is close to zero. Similar to the 1D example, it can be said that the NBF-IGF is the better approach than the RBF-IGF.  Both IGFs provide suitable initial guesses and, thus, the shooting method finds the solution for the given problem with the small number of iterations, whereas the conventional approach sometimes fails to find the solutions, as shown in Tables 3 and 4. In particular, applying the random initial guess with N 2 (0, Σ 100 2 ) provides a poor success rate that finds the solution. In addition, Figure 5 shows that the accuracy of the upper bound norm error and the computational time when using random initial guesses and the initial guess that was obtained by the IGFs. It is shown that the shooting method with the NBF-IGF is superior in the aspect of considering both computational cost and accuracy simultaneously.

Nonlinear 3D BVP: Orbit Determination Applications
Lambert's problem is one of the popular nonlinear BVPs in the celestial mechanics with 3D space. The two-body equation is given by [36] that satisfies r(t i ) = r i and r(t f ) = r f . Here, r(t) ∈ R 3 is the position vector, µ is the gravitational coefficient, which is defined as µ = 398, 600 km 3 /sec 2 , ||r(t)|| is the magnitude of the position vector, and r i and r f are the initial and final position vectors, respectively.
The performance of the proposed approaches for the orbit determination problem is validated by considering three orbits that have different semi-major axes and eccentricities, and they are listed in Table 5. In the orbit determination problem, the initial velocity at the initial time is found by applying the shooting method. Because each orbit has a different initial velocity, depending on the initial position, three kinds of initial positions are considered to deal with various initial velocities: near apogee, intermediate (between the apogee and perigee), and near perigee. Note that only one initial position is considered for Hubble telescope orbit, because it is almost circular orbit. Table 6 lists the BCs for each orbit. Here, the initial time t i is set as zero, so only the final time information is given. Initially, the shooting method with random initial guesses is applied to all cases, and Table 7 lists the Monte Carlo simulation results (1000 iterations) for each orbit case. Note that the initial guesses are randomly generated as N 3 (0, Σ 10 2 ), because the speed at each position is around 1.5∼10 km/sec. It sometimes fails to find the solution and, in particular, no solution is found for the near perigee case of Molniya orbit, as shown in the table. Because the initial value at the lower bound for the near perigee of Molniya orbit is relatively greater than the other cases, it is critical to an initial guess in order to find the solution. Like the 2D example, the cases that found the solution are only considered to compute the mean values of the number of iterations, the number of evaluations, computational time, and the upper bound norm error. Once the solution is found, the upper bound norm error is achieved as less than 0.05 %, and the solution is well found. However, one can say that the random initial guesses do not guarantee to find the solution, as shown in Table 7. To apply the proposed approach into this problem, four variables, which are the time t, and three components of the position vector r(t), are transformed into the quadratic NBFs and RBFs using Equations (19)-(22).

NBF-IGF
For the given BCs that are listed in Table 6, the endpoints for each case are defined for the NBFs, like the previous examples, and this problem has four unknown control points (t B 1 and r B 1 = [r 1,B 1 , r 2,B 1 , r 3,B 1 ] T ). The residual, including the nonlinear function on the right side of Equation (27), is derived as and the variables are defined in Equations (19) and (20). After applying the residual in Equation (28) into the cost function in Equation (17) and solving the optimization problem, the initial guess is computed using the obtained control points at s = 0.

RBF-IGF
The endpoints are defined as the same way of the NBF-IGF, and the weights that are associated with the endpoints are defined as w t 0 = w t 2 = w 1,r 0 = w 1,r 2 = w 2,r 0 = w 2,r 2 = w 3,r 0 = w 3,r 2 = 1. Hence, the number of unknown parameters is eight: four control points and four weights that are defined as t B 1 , r B 1 , w t 1 , and w r 1 = [w 1,r 1 , w 2,r 1 , w 3,r 1 ] T . With these unknown parameters, the cost function in Equation (17) is derived using Equation (28). The unknown intermediate control points and weights are obtained by solving the optimization problem, and the initial guess is calculated by applying the obtained parameters into Equation (14) at s = 0.
Using the simulation parameters that are listed in Table 6, the numerical simulations for the proposed approaches are conducted, and Tables 8-14 summarize the results. Overall, the proposed approaches provide the solution with less than 0.1% upper bound norm errors. Similar to the 1D and 2D examples, applying the RBF-IGF provides a closer initial guess to the initial value, and the shooting method with the RBF-IGF requires a less number of iterations and evaluations to find the solution. Because the initial guess is very close to the initial value, the trial trajectories are also similar to the solution trajectory, as shown in Figures 6-8.      However, the NBF-IGF also provides an initial guess that is close enough to the initial value to be obtained and it requires less computational time to obtain the initial guess and find the solution with the shooting method. Furthermore, the obtained upper bound norm error using the NBF-IGF is the same as the one using the RBF-IGF. Therefore, one can be said that that the NBF-IGF has better performance than the RBF-IGF in between both of the IGFs proposed.  Although the shooting method with random initial guesses provides the solution over 65 % during 1000 times simulations, except for the near perigee case of Molniya orbit, this does not guarantee finding the solution. On the other hand, the proposed approaches find the solution for all cases, including the near perigee case of Molniya orbit. In addition, applying the IGF, especially the NBF-IGF, guarantees the least computational burden to find the solution than the others, as shown in Figure 9.
Here, "Conventional" indicates the shooting method with random initial guesses, and the numbers shown in Figure 9 represent the computational time. Therefore, it is concluded that the proposed IGFs are able to provide suitable initial guesses that guarantee finding the solution of the orbit determination problem, including the case that is sensitive to the initial guess. Because the NBF-IGF not only provides an initial guess close to the initial value, but is also computationally more efficient than the RBF-IGF, the NBF-IGF is more beneficial than the RBF-IGF to determine the initial guess.

Conclusions
This work proposes efficient initial guess finders (IGFs) that provide proper initial guesses as approximating BVPs by using the non-rational and rational Bézier functions (BFs). Each IGF is named NBF-IGF and RBF-IGF, respectively. The IGFs provided proper initial guesses at the lower bound, and this value is applied to the shooting method. To validate the efficacy of the proposed approach, various dimensional problems from onedimensional (1D) to 3D, including Lambert's problem, are studied and compared with the conventional shooting method. For the 1D and 2D problems, the IGFs are as good as the conventional approach that uses random initial guesses, because the problems are relatively easy to solve. On the other hand, for the 3D Lambert's problem, the proposed IGFs outperform the conventional approach in terms of the computational burden and the ability to find the solution. In between two IGFs, although the RBF-IGF provides a more accurate initial guess, the NBF-IGF provides a close enough initial guess and it requires less computational burden. Therefore, the NBF-IGF is a better option than the RBF-IGF for providing a proper initial guess to effectively solve the BVPs. Funding: This study was supported by research fund from Chosun University, 2020.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.