An Efficient Semi-Analytical Solution of a One-Dimensional Curvature Equation that Describes the Human Corneal Shape

In this paper, a numerical approach is proposed to find a semi analytical solution for a prescribed anisotropic mean curvature equation modeling the human corneal shape. The method is based on an integral operator that is constructed in terms of Green’s function coupled with the implementation of Picard’s or Mann’s fixed point iteration schemes. Using the contraction principle, it will be shown that the method is convergent for both fixed point iteration schemes. Numerical examples will be presented to demonstrate the applicability, efficiency, and high accuracy of the proposed method.


Introduction
Corsato et al. [1] studied the existence, uniqueness, and regularity of solutions of the n-dimensional prescribed anisotropic mean curvature with Dirichlet boundary conditions where a and b are positive constants, n ≥ 2, and Ω ⊂ R n is a bounded Lipschitz domain.
Okrasi ński and Płociniczak modeled the topography of the human cornea [2] by a one-dimensional form of (1), where the mean curvature operator div(∇h/ 1 + |∇h| 2 ) is replaced by its linearization div(∇h) around 0, where Ω = [0, 1], and hence obtaining the nonlinear equation It was shown in [3] that Equation (2) has a unique solution for any positive constants a and b.However, as a closed form solution for the nonlinear boundary value problem (2) cannot be obtained, researchers tried to obtain approximate solutions by different numerical approaches.With 1% error of order, an approximate solution based on hyperbolic cosine function was experimentally shown to fit the original optical data of corneal shape [2].Płociniczak et al. suggested two other approximate solutions for problem (2), one method is by linearizing the nonlinear model and the other method is by applying the perturbation method [3].In [4], he suggested, in addition to the variational iteration method, the use of Taylor series method as accessible approach to solve any non-linear two-point boundary value problems.Converting the boundary value ordinary differential equation into a parabolic partial differential equation by adding an initial value derivative (t, pseudo-time), and then using the method of lines to numerically solve the resulting partial differential equation for sufficiently large values of t, was another approximate solution attempted by Płociniczak et al. [5].An interesting approach of obtaining approximate solution was due to Coelho et al. [6], where they developed a linear iterative scheme for approximating the solution by two monotone sequences of strict lower and upper solutions starting from an explicit pair of constant lower and upper solutions.A solution using a meshless method based on radial basis functions was discussed by Griffiths et al. [7].In a recent work, Okrasi ński and Płociniczak constructed lower and upper estimates that bound the components of the exact solution [8].Other interesting mathematical models for the human corneal topography can be viewed in articles [9][10][11][12].
For nonlinear boundary value problems, in general, one would consider popular methods such as Adomian decomposition [13], variational iteration [14,15], differential transformation [16], and spline interpolation [17].Green's function based methods have been intensively employed by researchers to solution of boundary value problems that arise in science and engineering.For example, Żur [18][19][20][21][22] presented a series of work to study the vibration of thin circular plates, thin circular plates with variable thickness, and elastically supported functionally graded annular plates using Green's and quasi-Green's functions.Andrade [23] calculated the exact Green's function for arbitrary rectangular potentials.Ahyoune et al. [24] used weighted combination of 2D and 3D analytical Green's function to solve quasi-static partial element equivalent circuits.A coupled Green's function with fixed point iteration method was employed to solve Bratu problems and boundary value problems that arise in heat transfer, strong nonlinear oscillation, and electroanalytical chemistry [25][26][27][28].
In this paper, a powerful method based on Green's function and fixed point iteration schemes is employed to find an approximate solution for the nonlinear boundary value problem (2).The construction of Green's function takes into account both end points of the prescribed domain.Therefore, it is natural to expect that this method will perform better in stability and uniform convergence than methods that take into account one initial point, and consequently making the possibility of convergence deterioration grows larger as we move away from the initial point towards the terminal point.The proposed method begins by identifying the linear and nonlinear parts of Equation (2) and then use the properties of Green's function to construct the Green's function G that mimics the solution of the corresponding linear part of (2) subject to the homogenous boundary conditions.A tailored linear integral operator expressed in terms of the constructed Green's function G is defined.The integrand of the integral operator is the product of G and the standard form of Equation ( 2).An iterative procedure is then obtained by applying the well-known Picard's or Mann's fixed point schemes to the integral operator.For our problem, it was observed that Picard's iterative scheme yielded a faster convergence than Mann's.The initial iterate, h 0 , was chosen to satisfy the corresponding linear Equation (2).
The rest of the paper is organized as follows.Section 2 presents the method in details from constructing the Green's function to applying Picard's and Mann's iterative schemes.In Section 3, the convergence analysis is established.In Section 4 numerical examples are presented, where the accuracy and convergence of the proposed method will be tested.A comparison between the proposed method and other approximate methods that addressed this problem will be detailed.Concluding remarks are given in Section 5.

Description of the Method
In this section, the properties of Green's function will be used to create an integral operator that will be embedded into Picard's fixed point iteration procedure.

Construction of Green's Function
Express Equation (2) in the convenient form where It is known that Green's function G(x, s) satisfies the equation subject to the corresponding homogeneous boundary conditions: where δ(x − s) is the Dirac delta function.A particular solution to Equation ( 3) can be expressed as Let h g be a solution to the homogenous equation L[h] = 0, then The Green's function can, then, be set as The constants C 1 , C 2 , C 3 , and C 4 will be determined using the following properties of Green's function: First, G(x, s) satisfies the homogeneous boundary conditions (5), and hence and Second, the continuity of G(x, s) at x = s implies that Third, G (x, s) has a jump discontinuity at x = s, that is The unique solution to system (9)-( 12) is

Green-Picard Fixed Point Iteration
Consider the linear integral operator Within the integrand of Equation ( 15), adding and subtracting F(b, x, h p ) followed by the use of Equation ( 6) leads to When applying Picard's fixed point iteration formula to Equation ( 16), the approximate solution to boundary value problem (2) is, iteratively, obtained in the form in which the initial approximation, h 0 , is the solution of the IVP

Green-Mann Fixed Point Iteration
For Mann's fixed point iteration, we replace Equation ( 18) with where α n are constants.Notice that for the special case α n = 1, Picard's fixed point iteration is obtained.Mann's iteration provides more flexibility to achieve fast convergence since the optimal choice of α n is obtained by finding the minimum L 2 norm of the residual error given by in which the residual error is defined by Remark 1.One desirable feature of the proposed method is that it gives some freedom in choosing the linear part of a prescribed differential equation.For example, in case of Equation ( 2), if we let L , and hence the Greens' function is constructed as where Remark 2. Notice that the main part of the proposed method is the construction of the Green's function which, as was pointed out in Remark 1, depended only on the linear part of (2).Hence, the proposed method imposes no restrictions or limitations on the nonlinearity of any differential equation that may be treated numerically by the proposed method.However, the nonlinear part of the differential equation, which appears in the integrand in the approximate solutions represented by Equations ( 18) and ( 20) may cause slower convergence for some strong nonlinear equations like Bratu type or strong nonlinear elliptic boundary value problems.

Convergence Analysis
In this section, the contraction principle [29] is used to show that the proposed method is convergent for Picard's fixed point iteration scheme.Consider the scheme where {h n (x)} are particular solutions and G is the Green's function for the operator L[h] = h − ah = 0 subject to the specified BCs given by where then, the iterative sequence {h n (x)} ∞ n=1 , given by ( 25) and ( 26), where x ∈ [0, 1] and using any bounded starting function on [0, 1], converges uniformly to the exact solution, h(x), of problem (2).
Proof of Theorem 1.We use the function space C[0, 1] equipped with the maximum norm defined by Using integration by parts twice and taking into account that G, G s , and G ss are continuous on [0, 1], the iterative scheme (25) becomes Since G(x, s) satisfies the boundary conditions G s (x, 0) = G(x, 1) = 0, and h n 's satisfies the boundary conditions h n (1) = h n (0) = 0, then Equation ( 29) simplifies to and hence Equation (30) reduces to Define T G : C[0, 1] → C[0, 1] to be the right side of Equation (32), that is According to Banach-Picard fixed point theorem, to prove convergence it suffices to show that T G is a contraction mapping.We have Since 1 + (v (s)) 2 ≥ 1, then from Equation (34) and the hypothesis of the theorem, we conclude and hence As the maximum value of |g(x)| on the interval [0, 1] occurs either at the critical points or at the endpoints, then straightforward calculations show that From Equations ( 35)-(37), we have where Remark 3. The convergence for Mann's fixed point iteration can be established in a similar manner.

Results and Discussions
In this section, two examples are presented, where the proposed method is applied to solve boundary value problem (2) for different parameters a and b.Since closed form solution cannot be found, the numerical solution obtained by fourth order Runge-Kutta method will be referred to as the exact solution.The proposed method will be tested against the following four approximate methods: 1.The zero-order solution based on the hyperbolic cosine function [2]: 2. The linearization approach [3]: where 3. The perturbation method [3]: where 4. Taylor series method [4]: where in which Example 1.Consider BVP (2) with a = b = 1. Figure 1 shows a comparison between the approximate solution obtained by the proposed method (using 5 iterations) and the aforementioned approximate methods.Figure 1a shows the entire time domain while a closer look at the approximate solutions using the first half of the domain is depicted in Figure 1b.It is clear from both figures that the approximate solution obtained by the proposed method is more accurate than all other methods.
In Table 1, the relative error computed by the formula Relative Error (%) = |Approx Solution − Exact Solution| Exact Solution × 100 (48) confirms that the proposed method is in better agreement with the exact solution than all other four methods.Example 2. For two data sets consisting of 123 × 123 points representing interior and exterior surfaces of cornea, it was calculated that a in ≈ 1.72281, b in ≈ 1.59638, a ex ≈ 1.37641, and b ex ≈ 1.3084 [2].In this example, we consider BVP (2) with the typical values of corneal parameters a = 1.6 and b = 1.7.
Figure 2 shows a comparison between the approximate solution obtained by the proposed method (using 5 iterations) and the aforementioned approximate methods.Figure 2a shows the entire time domain while a closer look at the approximate solutions using the first half of the domain is depicted in Figure 2b.It is clear from both figures that the approximate solution obtained by the proposed method is more accurate than all other methods.From Table 2, where the relative error is computed for all methods, we notice that the accuracy of the proposed method surpasses all other methods.

Conclusions
In this paper, a semi-analytic method was successfully applied to solve a one-dimensional curvature equation problem that models the human corneal shape.The proposed recursive method is based on the use of Green's functions with fixed point iterative procedure.Using the contraction mapping principle, the convergence of the method was established for both Picard and Mann's fixed point iterations.To assess the stability and convergence, we presented two examples to compare between the proposed method and other methods that addressed this problem.The proposed method shows a higher rate of accuracy and remained stable on the entire prescribed domain.
Tables 1 and 2 highlight a more desirable feature of the proposed method, that is the relative error resulting from the proposed method shows a uniform convergence throughout the given time domain while convergence deteriorates for the other methods as x moves away from the initial point towards the terminal point.The reason for this is that the proposed method, which is based on Green's function, takes into consideration both end points of the interval for the boundary value problem, while other methods considered in this study take into account only the initial point.
Another advantage of the proposed method is that there is no restriction on the nonlinearity of a boundary value problem.This was obvious since the construction of the Green's function in the proposed method depends only on the linear part of the differential equation.However, the presence of strong nonlinear terms may require large number of iterations to reach an acceptable small relative error.

Theorem 1 . 1 1
Assume that f x, h, h = + (h (x)) 2 is a continuous function that satisfies a Lipschitz condition with Lipschitz constant L c .Also, assume that

Table 1 .
Comparison of relative errors obtained by different numerical schemes for problem (2) with parameters a = b = 1.

Table 2 .
Comparison of relative errors obtained by different numerical schemes for problem (2) with parameters a = 1.6, b = 1.7.