Abstract
In this paper, we present and illustrate a frozen Jacobian multistep iterative method to solve systems of nonlinear equations associated with initial value problems (IVPs) and boundary value problems (BVPs). We have used Jacobi-Gauss-Lobatto collocation (J-GL-C) methods to discretize the IVPs and BVPs. Frozen Jacobian multistep iterative methods are computationally very efficient. They require only one inversion of the Jacobian in the form of LU-factorization. The LU factors can then be used repeatedly in the multistep part to solve other linear systems. The convergence order of the proposed iterative method is , where is the number of steps. The validity, accuracy, and efficiency of our proposed frozen Jacobian multistep iterative method is illustrated by solving fifteen IVPs and BVPs. It has been observed that, in all the test problems, with one exception in this paper, a single application of the proposed method is enough to obtain highly accurate numerical solutions. In addition, we present a comprehensive comparison of J-GL-C methods on a collection of test problems.
1. Introduction
Iterative methods for solving systems of nonlinear equations associated with IVPs and BVPs are important because, in general, it is hard to find a closed form solution. Generally, nonlinear IVPs and BVPs are solved in two steps. For the first step, we discretize the nonlinear problem to obtain a system of nonlinear equations. In the second step, we use some iterative method to solve the system of nonlinear equations.
Pseudospectral collocation methods offer excellent accuracy, and we use them for the discretization of IVPs and BVPs. These methods convert the targeted partial or ordinary differential equations into a system of algebraic equations. Recently, Doha et al. [1] used a J-GL-C method for the discretization of a nonlinear Schrödinger equation to obtain a system of ordinary differential equations and solved it by using an implicit Runge-Kutta method. Highly accurate results for four different kinds of nonlinear Schrödinger equations were obtained. In [2], nonlinear reaction-diffusion equations were solved by using a J-GL-C method. In [3], a J-GL-C method was also used to solve nonlinear complex generalized Zakharov systems of equations. Many authors (see, e.g., [4–7]) have used pseudospectral collocation method to solve nonlinear IVPs and BVPs. J-GL-C methods are attractive because they depend on two parameters and by changing the values of these parameters we can obtain different pseudospectral collocation methods, for example, the Legendre, Chebyshev, and Gegenbauer pseudospectral collocation methods. J-GL-C methods are based on Jacobi polynomials. The Jacobi polynomials are the eigenfunctions of the singular Sturm-Liouville problem [8],and can be generated by using the following recurrence relation: where The th derivative of the th-degree Jacobi polynomial is given by It should be noticed that Jacobi polynomials are orthogonal over the domain with the weight function . We are interested in approximating the differential operators by using J-GL-C methods. The easiest way to work with differential operators of different orders is to construct the differentiation matrix; for details on the construction of numerical differentiation matrix on Jacobi polynomials over the domain , see [9]. Suppose is the Jacobi-Gauss-Lobatto differentiation matrix for the first-order derivative operator over the domain . Then, a derivative of order can be approximated over the interval as
After having a setup for the discretization of IVPs and BVPs by using J-GL-C methods, we are looking for an efficient iterative method for the solution of the associated system of nonlinear equations. The first iterative method that comes into mind to solve a system of nonlinear equations is the Newton-Raphson (NR) method [10, 11], which can be written aswhere . The multistep frozen Jacobian version of the Newton-Raphson method (MNR) can be written as
Frozen Jacobian methods have the advantage of computing LU decomposition of the Jacobian, because this makes the solution of linear systems, having the Jacobian as their matrix, inexpensive. However, the convergence order of MNR as a function of the number of function evaluations is relatively low, and many researchers [12–23] have obtained frozen Jacobian multistep iterative methods with higher convergence orders for the same number of function evaluations. For example, HJ, FTUC, and MSF are presented in [22], [17], and [18], respectively. These methods can be describe as
In this paper, we propose another frozen Jacobian multistep iterative method for solving systems of nonlinear equations with higher convergence order than the previously proposed frozen Jacobian methods for the same number of function evaluations. Another goal of the paper is to compare different pseudospectral collocation methods generated by J-GL-C methods by choosing different values for the parameters and .
2. New Multistep Iterative Method
To construct the higher order frozen Jacobian multistep iterative method we first construct method with unknown coefficients as
We have unknowns: namely, . To develop our method we find , , and these unknowns to obtain the maximum convergence order. Once we have solved that problem we write down the multistep part as Notice that are already computed in the base method part and remain fixed in the multistep part. After finding the values of unknowns, we establish our proposed new multistep iterative method (NMIM) as In NMIM, the convergence order of the base method is nine, and we get an increment of five in the convergence order per additional function evaluation. The LU factors of are used in the multistep part to solve five lower and upper triangular systems of equations. NMIM requires only two Jacobian evaluations and the LU-factorization of one of the Jacobians.
3. Convergence Analysis
We use the method of mathematical induction to prove the convergence order of NMIM. We prove the convergence order for and, then, using induction, we will prove the convergence order for . In our analysis, we expand the system of nonlinear equations around its simple root by using Taylor’s series and use higher order Fréchet derivatives in the expansion. We use the symbolic algebra of Maple software to deal with symbolic computations. The Fréchet differentiability of the system of nonlinear equations is an essential constraint of our proof. Gateaux differentiability would not be enough, but that differentiability is neither enough to ensure that has a good local linear approximation, which is in soul of Newton-like methods. A function is said to be Fréchet differentiable at a point if there exists a linear operator such that The linear operator is called the first-order Fréchet derivative and we denote it by . The higher order Fréchet derivatives can be computed recursively using where is a vector independent from .
Theorem 1. Let be a sufficiently Fréchet differentiable function on an open convex neighborhood of with and , where denotes the Fréchet derivative of . Let and , for , where denotes the -order Fréchet derivative of . Then, for , with an initial guess in the neighborhood of , the sequence generated by NMIM converges to with at least local order of convergence nine and error equation where , , and is a 9-linear function: that is, with .
Proof. We define the error at the th step as . To complete the convergence proof, we perform detailed computations using Maple with the following results:
Now, we present the proof of convergence of NMIM via induction.
Theorem 2. The convergence order of NMIM for is .
Proof. All the computations are made under the assumption of Theorem 1. We know from Theorem 1 that the convergence order of NMIM method is nine for . By performing the computations in a similar manner to that of Theorem 1, we get the following error equation for : This proves that the convergence order of NMIM is for . Now, assume that the convergence order of NMIM is for a number of steps’ and and let us prove that the convergence order of NMIM for steps is for th step. If the convergence order of NMIM for steps is , then where is the asymptotic constant and symbol means “approximately equal to.” By using (17), we perform the following steps:This completes the proof that the convergence order of NMIM is for steps.
4. Computational Cost
The computational costs in terms of multiplications of different iterative methods are shown in Table 1, where is the number of steps, is the number of dimensions of the problem, is the ratio between the computational cost of a division and the computational cost of a multiplication, is the ratio between the computational cost of a function component evaluation and the computational cost of a multiplication, is the ratio between the computational cost of a component of the Jacobian and the computational cost of a multiplication, and is the ratio between the computational cost of a second derivative and the computational cost of a multiplication. Table 2 displays the difference of computational costs of methods with respect to NMIM when the methods have the same convergence order as NMIM with steps. Table 3 gives the conditions on and bounds on so that, for the same convergence order, the computational cost of NMIM is smaller than that of the other methods.
5. Numerical Simulations
In this section, we verify the convergence order of NMIM and solve fifteen nonlinear initial and boundary value problems to show the validity, accuracy, and efficiency of our proposed iterative method. For the discretization, we use four pseudospectral collocation methods. In all experiments but one, we apply NMIM once with several steps to reduce the absolute error in the computed numerical solution. To verify the convergence order computationally, we adopt the following definition of computational convergence order (CCO):
5.1. Computational Verification of Convergence Order
In general, it is hard to verify the convergence order by solving initial and boundary value problems. Therefore, we consider the following five systems of nonlinear equations .
The system of nonlinear equations is solved by Mathematica using 11,000 digits. Table 4 shows that CCOs agree with theoretical convergence orders.
5.2. Solving Initial and Boundary Value Problems
In this section, we solve different initial and boundary value problems by using different J-GL-C discretization methods. In all computations, we make spatial and temporal discretizations to translate the entire IVPs and BVPs into systems of nonlinear equations. The systems of nonlinear equations are then solved by using the proposed iterative method NMIM. In most of the problems, we use vector with initial and boundary conditions introduced in it as the initial guess. The inclusion of initial and boundary conditions in makes the initial guess nonsmooth but it works for most of the problems. When the nonsmooth initial guess does not work, we then make it smooth to get convergence.
5.2.1. The Lane-Emden Equation
The Lane-Emden equation is J-GL-C discretization methods translate the Lane-Emden equation into the following system of nonlinear equations:where denotes a diagonal matrix with values in the diagonal equal to the elements of its vector argument, , , and is the domain of the problem. We have solved the Lane-Emden equation for different values of ranging from two to five and . Figure 1 shows the initial guess and the computed numerical solutions for , , and 50 grid points. Table 5 gives the norm of the residue in the solution of the nonlinear systems by a single iteration of NMIM with varying .
5.2.2. Bratu Problem
Bratu problem is the boundary value problem: where is a parameter. The discretization is carried out using J-GL-C methods. The resulting system of nonlinear equations is where . We considered several values of , values and for the parameters of J-GL-C, 50 grid points, and solved the resulting system of nonlinear equations using a single application of NMIM with varying . The initial guess and computed numerical solutions are given in Figure 2. Table 6 gives the norm of the residue of the solution of the system of nonlinear equations for varying .
5.2.3. Frank-Kamenetskii Problem
Frank-Kamenetskii problem is where is a parameter. After discretization by J-GL-C methods we get the system of nonlinear equations:We considered several values of , values and for the parameters of J-GL-C, grid points, and solved the resulting system of nonlinear equations by using a single application of NMIM with varying . Figure 3 plots the initial guess and computed solutions and Table 7 gives the norm of the residue of the solution of the system of nonlinear equations.
5.2.4. 1 + 1 Nonlinear Schrödinger Equation
1 + 1 nonlinear Schrödinger equation is wherewith the initial and boundary conditions The function is a complex function and can be written as . The Schrödinger equation can be rewritten in terms of the real functions and aswith initial and boundary conditions Using J-GL-C methods, we can discretize (30) aswhere , , , , ,, is the Kronecker product, is the element-wise multiplication, is a matrix of all zeros, is the identity matrix, , , and is the J-GL-C operational matrix. The system of nonlinear equations (32) can be rewritten as where , , and is the vector incorporating the initial-boundary conditions. Four nonlinear Schrödinger equations are listed in Table 8 with their corresponding analytical solutions. We chose a domain and solved the equations using different pair of values for the parameters and of J-GL-C and different grids for the domain. As initial guess, we used , with boundary conditions set for problem 3. For problems 1, 2, and 4 NMIM does not converge with that initial guess and we made the initial guess smoother. We smoothed the initial guess by applying the iterationonce for problem 1 and twice for problems 2 and 4. Tables 9, 10, 11, and 12 give the norm of the error in the solution when the nonlinear system of equations is solved using a single application of NMIM with varying . Figures 4–7 display the solutions and the error in the components of and for , , and a grid .
5.2.5. Klein Gordon Equation
Klein Gordon equation iswherewith the initial and boundary conditions By discretization using J-GL-C methods we obtain where , , , and . The analytical solution of (35) is , where and . We set the parameters , , , and , picked up the domain , and solved the nonlinear system of equations by a single application of NMIM with initial guess , with the boundary conditions introduced. We used several pairs of the , parameters of J-GL-C. According to Table 13, the Chebyshev collocation method of first kind exhibits best accuracy when the grid is more refined. The solution and the error components for , , , and are plotted in Figure 8.
5.2.6. 2D Nonlinear Wave Equation
The 2D nonlinear wave equation iswherewith the initial and boundary conditions By discretization using J-GL-C methods we obtain where , , , , and . The solution of (39) with , and appropriate initial and boundary conditions is . We took a domain and solved the system of nonlinear equations using a single application of NMIM with initial guess with the initial and boundary conditions introduced. Table 14 gives the norm of the error as a function of for several pairs of values for and and several grids. The collocation method with which worked best was the Chebyshev collocation method of first kind. Figure 9 gives the error components for that collocation method and that grid.
5.2.7. 3D Nonlinear Wave Equation
The 3D nonlinear wave equation iswherewith the initial and boundary conditions By discretization using J-GL-C methods we obtain where , , , , , and . The solution of (43) with , and appropriate initial and boundary conditions is . We took a domain and solved the system of nonlinear equations using a single application of NMIM with initial guess and with initial and boundary conditions introduced. Table 15 gives the norm of the error as a function of for several pairs of values for , and several grids. The collocation method which worked best for the grid was the Legendre collocation method. Figure 10 gives the error components for that collocation method and that grid.
5.2.8. 3D Nonlinear Poisson Equation
The 3D nonlinear Poisson equation iswherewith the initial and boundary conditions By discretization using J-GL-C methods we obtain where , , , , and . The solution of (47) with with appropriate initial and boundary conditions is . We took a domain and solved the system of nonlinear equations using a single application of NMIM with initial guess and with initial and boundary conditions introduced. Table 16 gives the norm of the error as a function of for several pairs of values for and and several grids. The Legendre collocation method offers best accuracy over different grids. Figure 11 gives the error components for that collocation method and that grid.
5.2.9. The Nonlinear Complex Generalized Zakharov System
The nonlinear complex Zakharov system has importance in plasma physics [3]. The system includes two coupled nonlinear PDEs which can be written assubject to the initial and boundary conditions
Several numerical methods have been proposed recently for approximating the solution of (51) and (52) such as the homotopy method [24], the finite difference method [25, 26], and the variational iteration method [27]. Also, Bao et al. [28] suggested some high-accurate numerical methods for solving numerically (51) and (52). Bao and Sun [29] applied a new technique based on time-splitting discretization for approximating the solution of a variant of (51) and (52).
One can split (51) using the real and imaginary parts of , , and , aswith the initial and boundary conditions
The matrix form of the nonlinear system (53) iswhere where the constants , , , , and are given and the functions , and , are known. We use J-GL-C methods for discretizing (55) subject to the initial and boundary conditions (54) to reduce (55) to a system of nonlinear algebraic equations.
5.2.10. Complex Generalized Zakharov Equation I
The first complex generalized Zakharov problem [3] is with domain . The nonlinear problem (57) has the analytical solution We defineto check the numerical accuracy in the computed solutions. Problem (57) is discretized by using Chebyshev collocation method of first kind. We took the initial guess with initial and boundary conditions introduced and a grid . The achieved numerical accuracy is shown in Table 17. Our computed solutions are better than those in [3] because the maximum accuracy in [3] is of the order of . The solutions and error components are plotted in Figures 12, 13, and 14.
5.2.11. Complex Generalized Zakharov Equation II
The second complex generalized Zakharov equation [3] iswith the domain . That problem has the analytical solution
Problem (60) is discretized by using Chebyshev collocation method of first kind. We took the initial guess with initial and boundary conditions introduced and a grid . The achieved numerical accuracy is shown in Table 18. Our computed solutions are better than those in [3] because the ones in [3] achieved a maximum accuracy of the order of . The solutions and error components are plotted in Figures 15, 16, and 17.
5.2.12. Murray Equation
Murray equation [2] is The initial and boundary condition for Murray equation can be computed from its analytical solution The Chebyshev collocation method of first kind is used to discretize Murray equation over the grid . The achieved numerical accuracy is shown in Table 19. The computed accuracy in [2] is of the order of . The solution and error components are shown in Figure 18.
5.2.13. Nonlinear Reaction and Diffusion Equation
Nonlinear reaction and diffusion equation is The analytical solution is The Chebyshev collocation method of first kind is used to discretize nonlinear reaction and diffusion equation over the grid . The initial guess is with initial and boundary conditions introduced. The achieved numerical accuracy is shown in Table 20. The computed accuracy in [2] is of the order of . The solution and error components are shown in Figure 19.
5.2.14. Reaction-Diffusion Equation with Fischer-Kolmogorov Reaction Term and Density Dependent Diffusion
Reaction-diffusion equation with Fischer-Kolmogorov reaction term and density dependent diffusion is The analytical solution is The Chebyshev collocation method of first kind is used to discretize reaction-diffusion equation with Fischer-Kolmogorov reaction term and density dependent diffusion over the grid . The initial guess is with initial and boundary conditions introduced. The achieved numerical accuracy after 1 and 2 application of NMIM is shown in Table 21. The computed accuracy in [2] is only of the order of . The solution and error components are shown in Figure 20.
5.2.15. Two-Dimensional Nonlinear Reaction-Diffusion Equation
The two-dimensional nonlinear reaction-diffusion equation is Its analytical solution is The Chebyshev collocation method of first kind is used to discretize the two-dimensional nonlinear reaction-diffusion equation over grid . The initial guess is with initial and boundary conditions introduced. The achieved numerical accuracy is shown in Table 22. The computed accuracy in [2] is only of the order of . The error components are shown in Figure 21.
5.2.16. Ostrovsky-Hunter Equation
Ostrovsky-Hunter equation is where is a nonlinear coefficient and , are the dispersion coefficients. We assume a periodic boundary condition for this problem; that is,
The Chebyshev collocation method of first kind is used to discretize Ostrovsky-Hunter equation over a grid . The initial guess is with initial and boundary conditions introduced. The residue is shown in Table 23. Our solution for the grid is given in Table 24. We notice that our results do not agree with those presented in [4]. Our solution is plotted in Figure 22.
6. Conclusions
Frozen Jacobian multistep iterative methods are computationally efficient to solve systems of nonlinear equations. This is because the LU factors of the Jacobian can be reused to solve additional linear systems at a very small additional cost. We have developed a new frozen Jacobian multistep iterative method with increased order of convergence for a given number of function evaluations. The claimed order of convergence of NMIM has been confirmed computationally by solving several small systems of nonlinear equations. NMIM has been used to solve fifteen nonlinear IVPs and BVPs with different collocation methods. As a whole, the Chebyshev collocation method of the first kind seems to provide the best numerical accuracy in the set of fifteen IVPs and BVPs considered in the paper.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.