Frozen Jacobian Multistep Iterative Method for Solving Nonlinear IVPs and BVPs

1 Departament de Fı́sica i Enginyeria Nuclear, Universitat Politècnica de Catalunya, Eduard Maristanu 10, 08019 Barcelona, Spain 2 Dipartimento di Scienza e Alta Tecnologia, Universita dell’Insubria, Via Valleggio 11, 22100 Como, Italy 3 UCERD, Islamabad, Pakistan 4 Department of Mathematics, University of Engineering and Technology, Lahore, Pakistan 5 Department of Mathematics, King Abdulaziz University, Jeddah 21589, Saudi Arabia 6 Department of Information Systems, Faculty of Computing and Information Technology (Rabigh), King Abdulaziz University, Rabigh 21911, Saudi Arabia 7 Department of Mathematics, Government College University Lahore, Lahore, Pakistan 8 Departament d’Enginyeria Electronica, Universitat Politècnica de Catalunya, Diagonal 647, Planta 9, 08028 Barcelona, Spain 9 Heat and Mass Transfer Technological Center, Universitat Politècnica de Catalunya, Colom 11, 08222 Terrassa, Spain 10Institute of Mathematical Sciences, University of Malaya, Kuala Lumpur 50603, Malaysia


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 1+1 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 1 + 1 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  ≥ 1,  (,) 0 ( The th derivative of the th-degree Jacobi polynomial  (,)  is given by  ()  (,) () = Γ ( +  +  +  + 1) 2  Γ ( +  +  + 1)  (+,+) −

(𝑥) . (4)
It should be noticed that Jacobi polynomials are orthogonal over the domain [−1, 1] with the weight function (1 + )  (1 − )  .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 [−1, 1], see [9].Suppose  is the Jacobi-Gauss-Lobatto differentiation matrix for the first-order derivative operator over the domain [−1 , 1].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 F(y) = 0 is the Newton-Raphson (NR) method [10,11], which can be written as where det(F  (y  )) ̸ = 0.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][13][14][15][16][17][18][19][20][21][22][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 Second-order Fréchet derivative = 1

Number of solutions of lower and upper triangular systems of equations
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 .

New Multistep Iterative Method
To construct the higher order frozen Jacobian multistep iterative method we first construct method with unknown coefficients as for  = 1, 3 We have  + 3 unknowns: namely,  1,1 ,  2,1 ,  3,1 ,  4,1 , . . .,  4, ,  5,+1 , . . .,  5, .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 for  = 5, Notice that  5,+1 , . . .,  5, 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 for  = 5, 8 for  = 10, 13 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 F  (y 0 ) 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.

Convergence Analysis
We use the method of mathematical induction to prove the convergence order of NMIM.We prove the convergence order for  = 4 and, then, using induction, we will prove the convergence order for  > 4. 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 F(x) has a good local linear approximation, which is in soul of Newton-like methods.A function F(⋅) is said to be Fréchet differentiable at a point y if there exists a linear operator The linear operator A is called the first-order Fréchet derivative and we denote it by F  (y).The higher order Fréchet derivatives can be computed recursively using where k is a vector independent from y.
Theorem 1.Let F : Γ ⊆ R  → R  be a sufficiently Fréchet differentiable function on an open convex neighborhood Γ of y * ∈ R  with F(y * ) = 0 and det(F  (y * )) ̸ = 0, where F  (y) denotes the Fréchet derivative of F(y).Let A 1 = F  (y * ) and , where F () (y) denotes the -order Fréchet derivative of F(y).Then, for  = 4, with an initial guess in the neighborhood of y * , the sequence {y  } generated by NMIM converges to y * with at least local order of convergence nine and error equation e 4 = Le 0 9 +  (e 0 10 ) , where e 0 = y 0 − y * , e 0  =   ⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ (e 0 , e 0 , . . ., e 0 ), and Proof.We define the error at the th step as e  = y  − y * .
Now, we present the proof of convergence of NMIM via induction.
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  = 4.By performing the computations in a similar manner to that of Theorem 1, we get the following error equation for  = 5: 2 ) e 14 0 +  (e 15 0 ) .
This proves that the convergence order of NMIM is 5 − 11 for  = 5.Now, assume that the convergence order of NMIM is 5 − 11 for a number of steps'  ≥ 5 and and let us prove that the convergence order of NMIM for  steps is 5 − 6 for ()th step.If the convergence order of NMIM for  steps is 5 − 11, then where d 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 5 − 6 for  + 1 steps.

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.

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 F(x) = 0.
The system of nonlinear equations is solved by Mathematica using 11,000 digits.Table 4 shows that CCOs agree with theoretical convergence orders.

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 0 with initial and boundary conditions introduced in it as the initial guess.The inclusion of initial and boundary conditions in 0 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.

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 diag(⋅) denotes a diagonal matrix with values in the diagonal equal to the elements of its vector argument, S  = 2/Q  , S  = (2/Q  ) 2 , and [0, ] is the domain of the problem.We have solved the Lane-Emden equation for different values of  ranging from two to five and  = 3.
Figure 1 shows the initial guess and the computed numerical solutions for  = 1/2,  = −1/2, 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 .

Bratu Problem.
Bratu problem is the boundary value problem: Table 3: Conditions and bounds on  so that, for the same convergence order, the computational cost of NMIM is smaller.

Computational cost Bounds on 𝑚 Conditions
where  is a parameter.The discretization is carried out using J-GL-C methods.The resulting system of nonlinear equations is where S  = (2Q  ) 2 .We considered several values of , values  = −1/2 and  = −1/2 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 .

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: Complexity   We considered several values of , values  = −1/2 and  = −1/2 for the parameters of J-GL-C, 50 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   norm of the residue of the solution of the system of nonlinear equations.

1 + 1 Nonlinear Schrödinger
where with the initial and boundary conditions The function (, ) is a complex function and can be written as (, ) =  1 (, ) +  2 (, ).The Schrödinger equation can be rewritten in terms of the real functions  1 (, ) and  2 (, ) as with initial and boundary conditions Using J-GL-C methods, we can discretize (30) as where , and Q is the J-GL-C operational matrix.The system of nonlinear equations (32) can be rewritten as where ], and p 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 (, ) ∈ [−1, 1] × [0, 1] 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  1 = 0,  2 = 0 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 iteration once 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 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 10

Complex Generalized Zakharov Equation II.
The second complex generalized Zakharov equation [3] is Problem ( 60) is discretized by using Chebyshev collocation method of first kind.We took the initial guess 0 with initial and boundary conditions introduced and a grid 21 × 28.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 10 −6 .The solutions and error components are plotted in Figures 15, 16, and 17. [2] is

Murray Equation. Murray equation
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 16×14.The achieved numerical accuracy is shown in Table 19.The computed accuracy in [2] is of the order of 10 −8 .The solution and error components are shown in Figure 18.

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 22 × 14.The initial guess is 0 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 10 −8 .The solution and error components are shown in Figure 19.

Reaction-Diffusion Equation with Fischer-Kolmogorov
Reaction Term and Density Dependent Diffusion.Reactiondiffusion 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 30 × 30.The initial guess is 0 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 10 −9 .The solution and error components are shown in Figure 20.
The Chebyshev collocation method of first kind is used to discretize the two-dimensional nonlinear reaction-diffusion equation over grid 30×30.The initial guess is 0 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 10 −9 .The error components are shown in Figure 21.
where  is a nonlinear coefficient and  ∈ R,  > 0 are the dispersion coefficients.We assume a periodic boundary  (71) The Chebyshev collocation method of first kind is used to discretize Ostrovsky-Hunter equation over a grid 30 × 100.The initial guess is 0 with initial and boundary conditions  24.We notice that our results do not agree with those presented in [4].Our solution is plotted in Figure 22.

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.

Figure 1 :
Figure 1: Initial guess and solutions of the Lane-Emden equation for different values of .

Figure 2 :
Figure 2: Initial guess and solutions of Bratu problem for different values of .

Figure 3 :
Figure 3: Initial guess and solutions of Frank-Kamenetskii problem with different values of .
ti a l a x is T e m p o r a l a x is 0
ti a l a x is T e m p o r a l a x is

Figure 8 :
Figure 8: Solution and error components for nonlinear Klein Gordon equation.

Figure 9 :
Figure 9: Error components for the nonlinear 2D wave equation.

Figure 10 :
Figure 10: Error components for the nonlinear 3D wave equation.

Figure 11 :
Figure 11: Error components for the nonlinear 3D Poisson equation.
p o r a l a x is S p a ti a l a x p o r a l a x is S p a ti a l a x is 0

Figure 18 :
Figure 18: Solution and error components for Murray equation.

Figure 19 :
Figure 19: Solution and error components for nonlinear reaction and diffusion equation.

Figure 20 : 2 ∇ 2
Figure 20: Solution and error components for reaction-diffusion equation with Fischer-Kolmogorov reaction term and density dependent diffusion.

Figure 21 :
Figure 21: Error components for the two-dimensional nonlinear reaction-diffusion equation.
p o r a l a x is S p a ti a l a x is u(x, t)

Table 1 :
Computational costs of different iterative methods.

Table 2 :
Difference of computational costs of different iterative methods with respect to NMIM when the methods have the same convergence order as of NMIM with  steps.

Table 4 :
Verification of convergence order.

Table 5 :
Norm of the residue of the system of nonlinear equation resulting from Lane-Emden equation with varying .

Table 6 :
Norm of the residue of the solution of the system of nonlinear equations for Bratu problem and varying .

Table 7 :
Norm of the residue of the solution of the system of nonlinear equations for Fran-Kamenetskii problem and varying .

Table 9 :
(17) of the error in the solution for equation(17).
the nonlinear system of equations is solved using a single application of NMIM with varying .Figures4-7display the solutions and the error in the components of  1 and  2 for  = −4/2,  = −1/2, and a grid 15 × 30.

Table 10 :
Norm of the error in the solution for equation (30).

Table 11 :
Norm of the error in the solution for equation (32).      −  analytical     ∞

Table 12 :
Norm of the error in the solution for equation (35).      −  analytical     ∞

Table 13 :
Norm of the error in the Klein Gordon equation as a function of  for several values for  and  and grids.     u − u analytical     ∞

Table 16 :
Norm of the error for the 3D Nonlinear Poisson equation as a function of  for several values for  and  and grids.     u − u analytical     ∞   ,   ] × [  ,   ] .

Table 17 :
Maximum absolute error in the solution of the system of nonlinear equations for complex generalized Zakharov equation and varying .  −   +   )  = ,

Table 18 :
Maximum absolute error in the solution of the system of nonlinear equations for complex generalized Zakharov equation and varying .

Table 19 :
Maximum absolute error in the solution of the system of nonlinear equations for Murray equation and varying .

Table 20 :
Maximum absolute error in the solution of the system of nonlinear equations for nonlinear reaction and diffusion equation and varying .

Table 21 :
Maximum absolute error in the solution of the system of nonlinear equations for reaction-diffusion equation with Fischer-Kolmogorov reaction term and density dependent diffusion and varying .

Table 22 :
Maximum absolute error in the solution of the system of nonlinear equations for two-dimensional nonlinear reactiondiffusion equation and varying .The residue is shown in Table23.Our solution for the grid {0, 6/11, , 16/11} × {0.1, 0.2, 0.3, 0.4} is given in Table

Table 23 :
Residue of the system of nonlinear equations for Ostrovsky-Hunter equation and varying .