Solving Variable-Coefficient Fourth-Order ODEs with Polynomial Nonlinearity by Symmetric Homotopy Method

In this paper, the eigenfunction expansion method (EEM) is applied to find numerical solutions for variablecoefficient fourth-order ordinary differential equations (ODEs) with polynomial nonlinearity. The symmetry of the solution set for the resulting system of polynomial equations obtained from EEM of the problem is analyzed. The symmetric homotopy method is constructed to calculate all solutions of the discretization system for the problem. Due to the exploitation of symmetry, the number of computations is reduced. Numerical examples are presented to demonstrate the efficiency of the presented homotopy method.


Introduction
It is known that the nonlinear differential equations can govern many phenomena in nature. Once the nonlinearity increases in these equations, the structures of the solutions may become complicated. Due to these nonlinearities, the ODEs cannot be solved by using analytical methods. Therefore, the numerical methods can be used to solve such equations. Thus, the approximate solutions are required.
Boundary value problems (BVPs) for nonlinear higherorder ODEs have significant applications in applied mathematics. Especially, nonlinear fourth-order BVPs are commonly used in a wide variety of application fields such as physics, chemical phenomena, and engineering; see for example [1,2]. Recently, a great attention has been given to solve these problems numerically.
In this paper, we consider the following type of variablecoefficient fourth-order ODEs with polynomial nonlinearity ; , ∈ Ω ≡ 0, 1 , with boundary conditions where and are given continuous functions on the interval 0, 1 .

; ≔ ⋯
is a polynomial of with given variable coefficients , , , … , and degree . The BVP (1), usually describes the deformation of an elastic beam, which has been widely studied by many researchers [3][4][5][6]. In general, the existence and multiplicity of solutions to such kind of BVPs depend on the growth conditions of the nonlinearity term ; [7][8][9][10].
The EEM is one of the numerical techniques for finding multiple solutions of differential equations when some analytical methods fail. Hence, in this paper, we use the EEM as a discretization method, which is more accurate than other numerical methods such as finite difference and finite element methods [11].
Under the assumptions that and are in ! 0,1 , the linear fourth-order operator in BVP (1)  Let 1 & denote the 9-dimensional subspace spanned by the first 9 eigenfunctions. The approximate solution of (1) can be obtained by the linear combination of '( ) * & as follows: The approximation expressed by (3) (4) can be computed explicitly while other integral terms may not be computed explicitly. We can compute them by using the numerical quadrature. Given that 9 is large, the system of polynomial equations (5) will be a complicated system. Therefore, we can solve this system by using the numerical approximation techniques such as Fixed-point iteration method, Newton's method and Homotopy continuation method [12].
The main purpose of this paper is to construct a symmetric homotopy for solving variable-coefficient fourth-order ODEs with polynomial nonlinearity. We construct a simple fourthorder ODE as a starting system and then discretize it in eigensubspaces, where the subsystems have readily available solutions. Then, the resulting systems of polynomial equations in eigensubspaces are put together in a block-wise manner to construct the starting system for a general problem. We use the symmetry to reduce the number of computations by tracking the representative solution paths of the homotopy. For spurious solutions which may appear in the solution set of the discretized system of BVPs, we removed them by using certain filters.
The remainder of this paper is organized as follows. In section 2, we introduce a homotopy continuation method for the resulting system of polynomial equations obtained from the discretization of BVPs. In section 3, we address the analysis of symmetry group in the solution set of discretized problem. Section 4 contains the construction of symmetric homotopy for solving variable-coefficient fourth-order ODEs with polynomial nonlinearity. In section 5, we provide numerical examples to verify the efficiency of the symmetric homotopy. Finally, conclusions are summarized in section 6.

Homotopy Continuation Method for System of Polynomial Equations
It is known that Newton's method may not converge when solving a large system of polynomial equations that arise in connection with the discretization of BVPs. Sometimes, the computation of associated Jacobian matrix is rather expensive. Therefore, it is necessary to provide a substitute such as a homotopy continuation method, which has global convergence properties [13][14][15].
The basic idea of using a homotopy continuation method is to deform a simple starting system to a target system and track the zero-dimensional (all isolated) solutions of the intermediate systems. For a given system of polynomial equations, we construct an appropriate start system which can be easily solved and then deform these solutions through the smooth paths (homotopy paths) to get the desired solutions of the target system. From the numerical algebraic geometry community, there are some software packages such as PHCpack, HOM4PS-2.0, and Bertini [16][17][18]. These software packages can be used to solve the given system of polynomial equations. Therefore, the system of polynomial equations (5) can be solved by using the following total degree homotopy where K & ; , ; , … , ; & is the starting system defined as follows; moreover L ∈ ℂ a generic random number is used to avoid the singularities [19].
Recently, a few studies are dealing with computing the solutions of the discretization differential equations by using the homotopy continuation method; see, for example [11,[20][21][22] and the references therein.
Allgower et al. suggested a numerical continuation method in [20] to compute all solutions for a nonlinear second order two-point BVPs by using a finite difference method. They performed a homotopy deformation on successively refined discretization systems to obtain solutions on the finer level. When the system of polynomial equations is large, they proposed some filters for removing spurious solutions, which leads to an efficient homotopy. These filters depend on the information and properties of the solutions of the original problem, and there seems no general rule. For the symmetry of solution set for the discretized system, a special filter is selected.
Zhang et al. [11], proposed the EEM to obtain multiple solutions of semilinear elliptic partial differential equations with polynomial nonlinearity. They computed the corresponding solutions for the discretized problem on a coarse level, then utilize it as initial guesses to calculate related solutions of the discretized problem on a finer level. The extension homotopy method was constructed to find all solutions of the resulting system of polynomial equations efficiently. They used the error estimates of EEM to propose a filter strategy for removing spurious solutions. The finite element Newton method was applied to refine the computed solutions more.
In [21], Zhang et al. designed the symmetric homotopy method to find solutions of the system of polynomial equations derived from the discretizations of elliptic equations with cubic and quintic nonlinearities. They analyzed the symmetry of solution set for the system of polynomial equations arising from the eigenfunction expansion discretization of the problem. This symmetry arises from the dihedral symmetry V 7 of the unit square. They proved that this kind of homotopies could preserve the symmetry and reduce the number of computations, because only the representative paths have to be traced.
The homotopy method introduced in [22] is a bootstrapping approach, which was applied to compute multiple solutions of differential equations. This method used a homotopy continuation method based on the domain decomposition. That means to decompose the domain into subdomains, and then each subdomain is solved independently in parallel. Then, the solutions from the subdomains are combined to build solutions for the original problem. They applied this approach for solving problems consisting one and two-dimensional problems.

Symmetry Group for the Solution Set of the Discretized Problem
Throughout this section, we discuss the symmetry of discretized problem (1) by using the group actions of the dihedral group V (the point group). The dihedral group V consists of the identity and a reflection about the center of the domain; The symmetry of the discretized solution set is due to the V symmetry of the domain, and it is passed over through the eigenfunctions, the eigen-base and the expansion coefficients. We will address them in more details as follows: First, recall that the eigenfunctions of the linear fourth-order operator with boundary conditions have the following form For L ∈ V , the transformation on eigenfunctions can be obtained by; The transformation g Z on expansion coefficients induced by This leads to the following result Since ℊ Z is an orthogonal matrix, we have As a result, we find that (13) Again, this leads to the following result; Finally, the transformation on coefficients of the discretized problem with coefficients of nonlinearity ; can be as; Then, Therefore, the resulting polynomial system

Construction of the Symmetric Homotopy for Discretized Problem
In this section, we focus on the construction of the symmetric homotopy for variable-coefficient fourth-order ODEs with polynomial nonlinearity. Due to the symmetry analyzed in Section 3, we used it to construct an efficiency homotopy for our problem; and we limit our study to the nonlinearity case = 3,5.
Recall that the eigenvalues of the linear fourth-order ODE operator are 0 ≤ s ≤ s ≤ ⋯ ≤ s & , and the corresponding eigenfunctions are t , t , … , t & . Denote u s = v / t , where t = √2 ./ .0 .
Let w x y z : 1 → u s denotes the -orthogonal projection which is defined as; ? − w x y z , t@ = 0, ∀ t ∈ u s .
For L ∈ V , consider the transformation on -orthogonal projection is L ∘ ?w x y z & @ = w x y z L ∘ &

The Symmetric Homotopy for Cubic Polynomial Nonlinearity
Consider the following simple case of fourth-order BVP as the starting problem for our general problem with polynomial nonlinearity = 3 The EEM for (19) in an eigensubspace corresponding to an eigenvalue can be written as follows: find = ; t ∈ u s , such that Note that the variables ; in the discretized problem (21) are separable and each equation has 3 solutions. Therefore Thus, we find that it is easy to solve the discretized problem and then can set the resulting system of polynomial equations (21) block-wise to obtain the starting system for the general problem.
By using the symmetry, we can classify the solutions of the resulting system of polynomial equations obtained from EEM into the equivalence classes. Therefore, we only need to focus on the calculation of the representative solutions. Thus, we construct the following symmetric homotopy for our problem with cubic nonlinearity where K & is defined as follows; and, • ∈ ℂ is the generic random number for avoiding the singularities.

The Additional Symmetry for Odd Cubic Nonlinearity
It is known that the additional symmetry in the solutions of (1) can be obtained when ; is an odd function of .

That means if
& is a solution of (1), then − − & is also a solution, which is named ℤ symmetry. Therefore, we denote by V × ℤ the symmetry group of the discretized solutions and the homotopy solution paths.

The Symmetric Homotopy for Quintic Polynomial Nonlinearity
Analogues to cubic polynomial nonlinearity, the following simple case of fourth-order BVP can be the starting problem for our general problem with polynomial nonlinearity = 5 + = ~, ∈ 0, 1 0 = 1 = 0 = 1 = 0 (24) Again, the EEM for (24) in an eigensubspace corresponding to an eigenvalue can be written as follows: find = ; t ∈ u s , such that

The Additional Symmetry for Odd Quintic Nonlinearity
Similarly, the additional symmetry in the solutions of (1) can be obtained when & is a solution of (1), then ±. ±. & are also a solution besides − & , where . = √−1, which is named ℤ 7 symmetry. Therefore, we denote V × ℤ 7 to the symmetry group of the discretized solutions and the homotopy solution paths.

Numerical Results
In this section, we present some numerical examples to demonstrate the efficiency of symmetric homotopy method. As we mentioned before, there are several software packages which can be used for tracking the representative paths of our symmetric homotopy such as PHCpack, HOM4PS-2.0, and Bertini. Hence, in this paper, the numerical experiments were performed by using PHCpack which has some properties such as accepting start system defined by the user. The efficiency of our symmetric homotopy is verified by comparison with the efficiency of the total degree homotopy.

Example 1
Consider the following variable-coefficient fourth-order ODE with cubic polynomial nonlinearity Note that, such BVPs (29) may have infinitely many solutions [23], and the discretization problem has only finite solutions. The details of solution for discretized problem (29) are listed in Table 1. When 9 = 10, we observe that the number of solutions for total degree homotopy is approximately equal to 2 times the number of representative solutions, and the expected time by total degree homotopy is approximately equal to 2 times the expected time by symmetric homotopy.

Example 2
Consider the following variable-coefficient fourth-order ODE with an odd cubic polynomial nonlinearity where = m − n and = m − n 7 are chosen to be symmetric about the center of domain Ω. Table 2 shows the details of solution for discretized problem (30). When 9 = 10, the number of solutions for total degree homotopy is approximately equal to 4 times the number of representative solutions, and the expected time by total degree homotopy is approximately equal to 4 times the expected time by symmetric homotopy.
Due to the additional symmetry in the solution set of (30), we expect the efficiency of symmetric homotopy is more, comparing to solution of the problem with general cubic polynomial obtained in Table 1. As previously mentioned, the number of solutions increases with the number of eigenfunctions basis 9 . The solution set may contain spurious solutions, which means that approximation solution is not closed to original solutions for ODEs. Therefore, removing spurious solutions with certain filters will be necessary. The filtering strategy used in [11] can be a possible way to remove spurious solutions. The basic idea of this filter is based on error estimates of EEM, i.e., " -& " † N ≤ ,9 f ‡ , ‖ − & ‖ ‰ | ≤ ,9 f ‡X , where & ∈ 1 & , ≥ 0 and , is a generic positive constant independent of 9 . The approximate solutions & of discretized problems on successively finer levels satisfy the error estimates, and then can be viewed as a solution path ; 9 parameterized by discretization level 9 . Thus, the true approximate solutions should lie on a solution path. When 9 is large, the Cauchy's criterion for convergence implies that ‖ & − &X ‖ is very small. Applying Newton's method with initial guesses & = &f to the discretization system, one can expect convergence for nonspurious solutions.
The representative solutions for discretized problem (30) after filtering with 9 = 20 are displayed in Figure 1.

Example 3
Consider the following variable-coefficient fourth-order ODE with quintic polynomial nonlinearity The details of solution for discretized problem (31) are listed in Table 3.

Example 4
Consider the following variable-coefficient fourth-order ODE with an odd quintic polynomial nonlinearity   Table 4 shows the details of solution for discretized problem (32). When 9 = 7, the number of solutions for total degree homotopy is approximately equal to 2 times the number of representative solutions, and the expected time by total degree homotopy is approximately equal to 2 times the expected time by symmetric homotopy.
Similarly, because of the additional symmetry in the solution set of (32), the efficiency of symmetric homotopy is more, comparing to solution of the problem with general quintic polynomial obtained in Table 3.
The representative solutions for discretized problem (32) after filtering with 9 = 14 are displayed in Figure 2.

Conclusion
In this paper, a numerical technique based on EEM is presented to obtain the approximate solutions for variablecoefficient fourth-order ODEs with (cubic/quintic) polynomial nonlinearity. We have extensively used the symmetry with dihedral group V which leads to simple computations. The symmetric homotopy constructed for solving discretization systems of ODEs can preserve the symmetry and reduce the computational cost. Numerical results demonstrated that the symmetric homotopy method is efficient and promising.