# Frozen Jacobian Multistep Iterative Method for Solving Nonlinear IVPs and BVPs.

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 1 + 1 Schrodinger 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 Schrodinger 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],

(1 - [x.sup.2])[??]" (x) + ([theta] - [phi] + ([theta] + [phi] + 2) x) [??]' (x) + n(n + [theta] + [phi] +1) [??](x) = 0, (1)

and can be generated by using the following recurrence relation:

[mathematical expression not reproducible], (2)

where

[mathematical expression not reproducible], (3)

The pth derivative of the kth-degree Jacobi polynomial [J.sup.([theta], [phi].sub.k] is given by

[mathematical expression not reproducible]. (4)

It should be noticed that Jacobi polynomials are orthogonal over the domain [-1,1] with the weight function [(1 + x).sup.[theta][(1 - x).sup.[phi]. 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 Q is the Jacobi-Gauss-Lobatto differentiation matrix for the first-order derivative operator over the domain [-1,1]. Then, a derivative of order p can be approximated over the interval [a, b] as

[d.sup.P]/[dx.sup.p] [approximately equal to] [(2/b - a Q)).sup.p)]. (5)

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

[mathematical expression not reproducible], (6)

where det(F'([Y.sub.n]) [not equal to] 0. The multistep frozen Jacobian version of the Newton-Raphson method (MNR) can be written as

[mathematical expression not reproducible]. (7)

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

[mathematical expression not reproducible]. (8)

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 [theta] and [phi].

2. New Multistep Iterative Method

To construct the higher order frozen Jacobian multistep iterative method we first construct method with unknown coefficients as

[mathematical expression not reproducible]. (9)

We have r + 3 unknowns: namely, [[alpha].sub.1,1], [[alpha].sub.2,1], [[alpha].sub.3,1i, [[alpha].sub.4,1], ... [[alpha].sub.4,q], [[alpha].sub.5,q+1], ..., [[alpha].sub.5,r]. To develop our method we find q, r, and these unknowns to obtain the maximum convergence order. Once we have solved that problem we write down the multistep part as

[mathematical expression not reproducible]. (10)

end.

Notice that [[alpha].sub.5,q+1], ..., [[alpha].sub.5,r] 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

[mathematical expression not reproducible]. (11)

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.sub.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.

3. Convergence Analysis

We use the method of mathematical induction to prove the convergence order of NMIM. We prove the convergence order for m = 4 and, then, using induction, we will prove the convergence order for m > 4. In our analysis, we expand the system of nonlinear equations around its simple root by using Taylor's series and use higher order Frechet derivatives in the expansion. We use the symbolic algebra of Maple software to deal with symbolic computations. The Frechet 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(x) is said to be Frechet differentiable at a point y if there exists a linear operator A [member of] L([R.sup.n], [R.sup.q]) such that

[mathematical expression not reproducible]. (12)

The linear operator A is called the first-order Frechet derivative and we denote it by F'(y). The higher order Frechet derivatives can be computed recursively using

F' (y) = Jacobian (F (y)), [F.sup.j] (y)[v.sup.j-1] 1 = Jacobian ([F.sup.j-1](y) [v.sup.j-1), j [greater than or equal to] 2, (13)

where v is a vector independent from y.

Theorem 1. Let F : [GAMMA] [??] [R.sup.n] [right arrow] [R.sup.n] be a sufficiently Frechet differentiable function on an open convex neighborhood [GAMMA] of y* [member of] [R.sup.n] with F'(y*) = 0 and det(F'(y*)) [not equal to] 0, where F'(y) denotes the Frechet derivative of F(y). Let [A.sub.1] = F'(y*) and [A.sub.j] = (1/j!)F'[(y*).sup.-1] [F.sup.(j)](y*),for j [greater than or equal to] 2, where [F.sup.(j)] (y) denotes the j-order Frechet derivative of F(y). Then, for m = 4, with an initial guess in the neighborhood of y*, the sequence {[y.sub.m]} generated by NMIM converges to y* with at least local order of convergence nine and error equation

[e.sub.4] = [Le.sub.0.sup.9] + O([e.sub.0.sup.10]), (14)

where [e.sub.0] = [y.sub.0] - y*, [e.sub.0.sup.p] = ([p times/[e.sub.0], [e.sub.0], ..., [e.sub.0]), and L = 184[A.sup.8.sub.2] + 8[A.sup.3.sub.2] [A.sup.3.sub.2] - 24[A.sub.3] [A.sup.6.sub.2] is a 9-linear function: that is, L [member of] ([9 times/L ([R.sup.n], [R.sup.n], [R.sup.n], ..., [R.sup.n]) with [Le.sub.0.sup.9] [member of] [R.sup.n].

Proof. We define the error at the nth step as [e.sub.n] = [y.sub.n] - y*. To complete the convergence proof, we perform detailed computations using Maple with the following results:

[mathematical expression not reproducible]. (15)

Now, we present the proof of convergence of NMIM via induction.

Theorem 2. The convergence order of NMIM for m [greater than or equal to] 4 is 5m - 11.

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 m = 4. By performing the computations in a similar manner to that of Theorem 1, we get the following error equation for m =5:

[e.sub.5] = (4[A.sup.3.sub.2] [A.sub.3] - 12[A.sub.3][A.sup.3.sub.2] + 60 [A.sup.5.sub.2]) * (184[A.sup.8.sub.2] + 8[A.sup.3.sub.2] [A.sub.3][A.sup.3.sub.2] - 24[A.sub.3][A.sup.6.sub.3]) [e.sup.14.sub.0] + O ([e.sup.15.sub.0]). (16)

This proves that the convergence order of NMIM is 5m - 11 for m = 5. Now, assume that the convergence order of NMIM is 5m - 11 for a number of steps' m >5 and and let us prove that the convergence order of NMIM for m steps is 5m- 6 for (m)th step. If the convergence order of NMIM for m steps is 5m - 11, then

[mathematical expression not reproducible]. (17)

where d is the asymptotic constant and symbol ~ means "approximately equal to." By using (17), we perform the following steps:

[mathematical expression not reproducible]. (18)

This completes the proof that the convergence order of NMIM is 5m - 6 for m + 1 steps.

4. Computational Cost

The computational costs in terms of multiplications of different iterative methods are shown in Table 1, where m is the number of steps, n is the number of dimensions of the problem, l is the ratio between the computational cost of a division and the computational cost of a multiplication, a is the ratio between the computational cost of a function component evaluation and the computational cost of a multiplication, [beta] is the ratio between the computational cost of a component of the Jacobian and the computational cost of a multiplication, and [gamma] 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 m steps. Table 3 gives the conditions on a and bounds on m 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):

[mathematical expression not reproducible]. (19)

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.

[mathematical expression not reproducible]. (20)

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 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.

5.2.1. The Lane-Emden Equation. The Lane-Emden equation is

[mathematical expression not reproducible]. (21)

J-GL-C discretization methods translate the Lane-Emden equation into the following system of nonlinear equations:

[mathematical expression not reproducible], (22)

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

5.2.2. Bratu Problem. Bratu problem is the boundary value problem:

[mathematical expression not reproducible], (23)

where [alpha] is a parameter. The discretization is carried out using J-GL-C methods. The resulting system of nonlinear equations is

F (x) = [S.sub.xx] x + [[alpha]e.sup.x] = 0,

F' (x) = [S.sub.xx] + [alpha] diag ([e.sup.x]), (24)

where [S.sub.xx] = [(2[Q.sub.x]).sup.2]. We considered several values of [alpha], values [theta] = -1/2 and [phi] = -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 m. 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 m.

5.2.3. Frank-Kamenetskii Problem. Frank-Kamenetskii problem is

[mathematical expression not reproducible], (25)

where [alpha] is a parameter. After discretization by J-GL-C methods we get the system of nonlinear equations:

[mathematical expression not reproducible]. (26)

We considered several values of [alpha], values [theta] = -1/2 and [phi] = -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 m. 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 Schrodinger Equation. 1 + 1 nonlinear Schrodinger equation is

[mathematical expression not reproducible], (27)

where

[D.sub.x] = [[a.sub.x], [b.sub.x]], [D.sub.t] = [0, [t.sub.f]], (28)

with the initial and boundary conditions

[mathematical expression not reproducible]. (29)

The function [psi](x, t) is a complex function and can be written as [psi](x, t) = [[psi].sub.1](x, t) + i[[psi].sub.2] (x, t). TheSchrodinger equation can be rewritten in terms of the real functions [[psi].sub.1] (x, t) and [[psi].sub.2] (x, t) as

[mathematical expression not reproducible], (30)

with initial and boundary conditions

[mathematical expression not reproducible]. (31)

Using J-GL-C methods, we can discretize (30) as

[mathematical expression not reproducible], (32)

where [mathematical expression not reproducible], [cross product] is the Kronecker product, [??] is the element-wise multiplication, O is a matrix of all zeros, I is the identity matrix, [S.sub.t] = (2/[t.sub.f]) [Q.sub.t], [S.sub.xx] = ((2/([b.sub.x] - [a.sub.x]))[[Q.sub.x]).sup.2], and Q is the J-GL-C operational matrix. The system of nonlinear equations (32) can be rewritten as

F ([psi]) = A[psi] + g ([psi]) - p, (33)

where [mathematical expression not reproducible], and p is the vector incorporating the initial-boundary conditions. Four nonlinear Schrodinger equations are listed in Table 8 with their corresponding analytical solutions. We chose a domain (x, t) [member of] [-1, 1] x [0, 1] and solved the equations using different pair of values for the parameters d and 0 of J-GL-C and different grids for the domain. As initial guess, we used [[psi].sub.1] = 0, [[psi].sub.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

[psi] = -[A.sup.-1] (2g ([psi]) - p) (34)

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 nonlinear system of equations is solved using a single application of NMIM with varying m. Figures 4-7 display the solutions and the error in the components of [[psi].sub.1] and [[psi].sub.2] for [theta] = -4/2, [phi] = -1/2, and a grid 15 x 30.

5.2.5. Klein Gordon Equation. Klein Gordon equation is

[u.sub.tt] (x, t) - [c.sup.2] [u.sub.xx] (x, t) + ku - [gamma]u [(x, t).sup.3] = 0, (x, t) [member of] [D.sub.x] x [D.sub.t], (35)

where

[D.sub.x] = [[a.sub.x], [b.sub.x] [D.sub.t] = [0, [t.sub.f]], (36)

with the initial and boundary conditions

[mathematical expression not reproducible]. (37)

By discretization using J-GL-C methods we obtain

F (u) = [S.sub.tt] u - [c.sup.2] [S.sub.xx] u + ku - [[gamma]u.sup.3] = 0, F' (u) = [S.sub.tt] - [c.sup.2] [S.sub.xx] + kI - 3[gamma] diag ([u.sup.2]), (38)

where [mathematical expression not reproducible]. The analytical solution of (35) is u(x, t) = [delta] sech ([kappa](v- vt)), where [kappa] = [square root of (k/([c.sup.2] - [v.sup.2])] and [delta] = [square root of (2k/[gamma])]. We set the parameters c = 1, y = 1, v = 0.5, and k = 0.5, picked up the domain (x, t) [member of] [-10, 10] x [0, 1], and solved the nonlinear system of equations by a single application of NMIM with initial guess u = 0, with the boundary conditions introduced. We used several pairs of the [theta], [phi] 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 [theta] = -1/2, [phi] = -1/2, [n.sub.x] = 120, and [n.sub.t] = 20 are plotted in Figure 8.

5.2.6. 2D Nonlinear Wave Equation. The 2D nonlinear wave equation is

[u.sub.tt] (x , y , 0 - [c.sup.2] ([u.sub.xx] (x, y,> t) + [u.sub.yy] (x, y, t)) f (u) = p (t, x, y), (x, y, t) [member of] [D.sub.x] x [D.sub.y] x [D.sub.t], (39)

where

[mathematical expression not reproducible], (40)

with the initial and boundary conditions

[mathematical expression not reproducible]. (41)

By discretization using J-GL-C methods we obtain

[mathematical expression not reproducible], (42)

where [mathematical expression not reproducible]. The solution of (39) with c = 1, f(u) = [u.sup.3] and appropriate initial and boundary conditions is [e.sup.-t] sin(x + y). We took a domain (x, y, t) [member of] [-1, 1] x [-1, 1] x [0,2] and solved the system of nonlinear equations using a single application of NMIM with initial guess u = 0 with the initial and boundary conditions introduced. Table 14 gives the norm of the error as a function of m for several pairs of values for [theta] and [phi] and several grids. The collocation method with 20 x 20 x 20 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 is

[mathematical expression not reproducible], (43)

where

[mathematical expression not reproducible], (44)

with the initial and boundary conditions

[mathematical expression not reproducible]. (45)

By discretization using J-GL-C methods we obtain

[mathematical expression not reproducible], (46)

where [mathematical expression not reproducible]. The solution of (43) with c = 0.8, f(u) = [u.sup.2] and appropriate initial and boundary conditions is [e.sup.-t] sin(x + y + z). We took a domain(x, y, z, t) [member of] [0, 1] x [0, 1]x [0 , 1] x [0 , 1] and solved the system of nonlinear equations using a single application of NMIM with initial guess u = 0 and with initial and boundary conditions introduced. Table 15 gives the norm of the error as a function of m for several pairs of values for 0, 0 and several grids. The collocation method which worked best for the grid 7 x 7 x 7 x 12 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 is

[mathematical expression not reproducible], (47)

where

[mathematical expression not reproducible], (48)

with the initial and boundary conditions

[mathematical expression not reproducible], (49)

By discretization using J-GL-C methods we obtain

[mathematical expression not reproducible], (50)

where [mathematical expression not reproducible]. The solution of (47) with f(u) = [u.sup.4] with appropriate initial and boundary conditions is sin (x + y + z). We took a domain (x, y, z) e [0, 1] x [0, 1] x [0, 1] and solved the system of nonlinear equations using a single application of NMIM with initial guess u = 0 and with initial and boundary conditions introduced. Table 16 gives the norm of the error as a function of m for several pairs of values for [theta] and [phi] 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 as

[mathematical expression not reproducible], (51)

subject to the initial and boundary conditions

[mathematical expression not reproducible], (52)

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 [psi](x, t), u(x, t), and v(x, t), as

[mathematical expression not reproducible], (53)

with the initial and boundary conditions

[mathematical expression not reproducible], (54)

The matrix form of the nonlinear system (53) is

[mathematical expression not reproducible], (55)

where

[mathematical expression not reproducible], (56)

where the constants [[delta].sub.1], [[delta].sub.2], [[delta].sub.3], [[delta].sub.4], and [c.sub.s] are given and the functions [[alpha].sub.i]((t), 1 [less than or equal to] i [less than or equal to] 6 and [[beta].sub.j](x), 1 [less than or equal to] j [less than or equal to] 4 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

[mathematical expression not reproducible], (57)

with domain [-1, 1] x [0, 3.0]. The nonlinear problem (57) has the analytical solution

[mathematical expression not reproducible]. (58)

We define

[mathematical expression not reproducible], (59)

to 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 0 with initial and boundary conditions introduced and a grid 50 x 20. 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.sup.-7]. 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] is

[mathematical expression not reproducible], (60)

with the domain [-1, 1] x [0, 1]. That problem has the analytical solution

[mathematical expression not reproducible], (61)

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 x 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.sup.-6]. The solutions and error components are plotted in Figures 15, 16, and 17.

5.2.12. Murray Equation. Murray equation [2] is

[u.sub.t] = [u.sub.xx] + [[lambda].sub.1], [uu.sub.x] + [[lambda].sub.2] u - [[lambda].sub.3] [u.sup.2], (x, t) [member of] [-1, 1] x [0, 1]. (62)

The initial and boundary condition for Murray equation can be computed from its analytical solution

u (x, t) = [[lambda].sub.2]/2[[lambda].sub.3] (1 + tanh ([[lambda].sub.2]/8[[lambda].sup.2.sub.3] (2[[lambda].sub.1] [[lambda].sub.3]x + ([[lambda].sup.2.sub.1] + 4 [[lambda].sup.2.sub.3]). (63)

The Chebyshev collocation method of first kind is used to discretize Murray equation over the grid 16 x 14. The achieved numerical accuracy is shown in Table 19. The computed accuracy in [2] is of the order of [10.sup.-8]. The solution and error components are shown in Figure 18.

5.2.13. Nonlinear Reaction and Diffusion Equation. Nonlinear reaction and diffusion equation is

[u.sub.t] = [u.sub.xx] + u (1 + u - [u.sup.2]), (x, t) [member of] [-1, 1] x [0, 1]. (64)

The analytical solution is

u (x, t) = 1/2 + [square root of (5)]/2 tanh (10/4 x + 5/4 t). (65)

The Chebyshev collocation method of first kind is used to discretize nonlinear reaction and diffusion equation over the grid 22 x 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.sup.-8]. 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

[u.sub.t] = [([u.sup.2]).sub.xx] + u (1 - u), (x, t) [member of] [-1, 1] x [0, 1]. (66)

The analytical solution is

u (x, t) = 2/1 + tanh ((x + t)/4). (67)

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 x 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.sup.-9]. 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

[u.sub.t] = 1/2 [[nabla].sup.2] u + [u.sup.2] (u- 1), (x, y, t) [member of] [-1, 1] x [-1 ,1] x [0, 20]. (68)

Its analytical solution is

u (x, y, t) = 1/1 + [e.sup.(1/[square root of (2)](x + y - (1/[square root of (2)t)]. (69)

The Chebyshev collocation method of first kind is used to discretize the two-dimensional nonlinear reaction-diffusion equation over grid 30 x 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.sup.-9]. The error components are shown in Figure 21.

5.2.16. Ostrovsky-Hunter Equation. Ostrovsky-Hunter equation is

([u.sub.t]- [[beta]u.sub.xxx] + [[[alpha]uu.sub.x]).sub.x] = [gamma]u, (70)

where [alpha] is a nonlinear coefficient and [beta] [member of] R, [gamma] > 0 are the dispersion coefficients. We assume a periodic boundary condition for this problem; that is,

[mathematical expression not reproducible]. (71)

The Chebyshev collocation method of first kind is used to discretize Ostrovsky-Hunter equation over a grid 30 x 100. The initial guess is 0 with initial and boundary conditions introduced. The residue is shown in Table 23. Our solution for the grid {0, 6[pi]/11, [pi], 16 [pi]/11} x {0.1, 0.2,0.3,0.4} 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.

https://doi.org/10.1155/2017/9407656

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] E. H. Doha, A. H. Bhrawy, M. A. Abdelkawy, and R. A. Van Gorder, "Jacobi-Gauss-Lobatto collocation method for the numerical solution of 1 +1 nonlinear Schrodinger equations," Journal of Computational Physics, vol. 261, pp. 244-255, 2014.

[2] A. H. Bhrawy, E. H. Doha, M. A. Abdelkawy, and R. A. Van Gorder, "Jacobi-Gauss-Lobatto collocation method for solving nonlinear reaction-diffusion equations subject to Dirichlet boundary conditions," Applied Mathematical Modelling. Simulation and Computation for Engineering and Environmental Systems, vol. 40, no. 3, pp. 1703-1716, 2016.

[3] A. H. Bhrawy, "An efficient Jacobi pseudospectral approximation for nonlinear complex generalized Zakharov system," Applied Mathematics and Computation, vol. 247, pp. 30-46, 2014.

[4] M. Dehghan and F. Fakhar-Izadi, "The spectral collocation method with three different bases for solving a nonlinear partial differential equation arising in modeling of nonlinear waves," Mathematical and Computer Modelling, vol. 53, no. 9-10, pp. 1865-1877, 2011.

[5] E. H. Doha, A. H. Bhrawy, and S. S. Ezz-Eldien, "Efficient Chebyshev spectral methods for solving multi-term fractional orders differential equations," Applied Mathematical Modelling. Simulation and Computation for Engineering and Environmental Systems, vol. 35, no. 12, pp. 5662-5672, 2011.

[6] E. H. Doha, A. H. Bhrawy, and R. M. Hafez, "On shifted Jacobi spectral method for high-order multi-point boundary value problems," Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 10, pp. 3802-3810, 2012.

[7] E. Tohidi and S. Lotfi Noghabi, "An efficient legendre pseudo-spectral method for solving nonlinear quasi bang-bang optimal control problems," Journal of Applied Mathematics, Statistics and Informatics, vol. 8, no. 2, pp. 73-85, 2012.

[8] G. Szego, Orthogonal Polynomials, Colloquium Publications XXIII, American Mathematical Society, 1939.

[9] J. Shen, T. Tang, and L.-L. Wang, Spectral Methods: Algorithms, Analysis and Applications, vol. 41, Springer, 2011.

[10] J. F. Traub, Iterative methods for the solution of equations, Prentice-Hall Series in Automatic Computation, Prentice-Hall, Englewood Cliffs, NJ, USA, 1964.

[11] J. M. Ortega and W. C. Rheinboldt, Iterative solution of nonlinear equations in several variables, Academic Press, New York, NY, USA, 1970.

[12] A. Cordero, M. Kansal, V. Kanwar, and J. R. Torregrosa, "A stable class of improved second-derivative free Chebyshev-Halley type methods with optimal eighth order convergence," Numerical Algorithms, vol. 72, no. 4, pp. 937-958, 2016.

[13] V. Arroyo, A. Cordero, and J. R. Torregrosa, "Approximation of artificial satellites' preliminary orbits: the efficiency challenge," Mathematical and Computer Modelling, vol. 54, no. 7-8, pp. 1802-1807, 2011.

[14] D. A. Budzko, A. Cordero, and J. R. Torregrosa, "New family of iterative methods based on the Ermakov-Kalitkin scheme for solving nonlinear systems of equations," Computational Mathematics and Mathematical Physics, vol. 55, no. 12, pp. 19471959, 2015.

[15] S. Qasim, Z. Ali, F. Ahmad, S. Serra-Capizzano, M. Z. Ullah, and A. Mahmood, "Solving systems of nonlinear equations when the nonlinearity is expensive," Computers & Mathematics with Applications, vol. 71, no. 7, pp. 1464-1478, 2016.

[16] U. Qasim, Z. Ali, F. Ahmad, S. Serra-Capizzano, M. Z. Ullah, and M. Asma, "Constructing frozen Jacobian iterative methods for solving systems of nonlinear equations, associated with ODEs and PDEs using the homotopy method," Algorithms (Basel), vol. 9, no. 1,18 pages, 2016.

[17] F. Ahmad, E. Tohidi, and J. A. Carrasco, "A parameterized multi-step Newton method for solving systems of nonlinear equations," Numerical Algorithms, vol. 71, no. 3, pp. 631-653, 2016.

[18] M. Z. Ullah, S. Serra-Capizzano, and F. Ahmad, "An efficient multi-step iterative method for computing the numerical solution of systems of nonlinear equations associated with ODEs," Applied Mathematics and Computation, vol. 250, pp. 249-259, 2015.

[19] F. Ahmad, E. Tohidi, M. Z. Ullah, and J. A. Carrasco, "Higher order multi-step Jarratt-like method for solving systems of nonlinear equations: application to PDEs and ODEs," Computers & Mathematics with Applications. An International Journal, vol. 70, no. 4, pp. 624-636, 2015.

[20] E. S. Alaidarous, M. Z. Ullah, F. Ahmad, and A. S. Al-Fhaid, "An efficient higher- order quasilinearization metho d for solving nonlinear BVPs," Journal of Applied Mathematics, vol. 2013, Article ID 259371, 11 pages, 2013.

[21] M. Z. Ullah, F. Soleymani, and A. S. Al-Fhaid, "Numerical solution of nonlinear systems by a general class of iterative methods with application to nonlinear PDEs," Numerical Algorithms, vol. 67, no. 1, pp. 223-242, 2014.

[22] H. Montazeri, F. Soleymani, S. Shateyi, and S. S. Motsa, "On a new method for computing the numerical solution of systems of nonlinear equations," Journal of Applied Mathematics, vol. 2012, Article ID 751975, 15 pages, 2012.

[23] F. Soleymani, T. Lotfi, and P Bakhtiari, "A multi-step class of iterative methods for nonlinear systems," Optimization Letters, vol. 8, no. 3, pp. 1001-1015, 2014.

[24] S. Abbasbandy, E. Babolian, and M. Ashtiani, "Numerical solution of the generalized Zakharov equation by homotopy analysis method," Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 12, pp. 4114-4121, 2009.

[25] Q. S. Chang and H. Jiang, "A conservative difference scheme for the Zakharov equations," Journal of Computational Physics, vol. 113, no. 2, pp. 309-319,1994.

[26] Q. S. Chang, B. L. Guo, and H. Jiang, "Finite difference method for generalized Zakharov equations," Mathematics of Computation, vol. 64, no. 210, pp. 537-553, 1995.

[27] M. Javidi and A. Golbabai, "Exact and numerical solitary wave solutions of generalized Zakharov equation by the variational iteration method," Chaos, Solitons and Fractals, vol. 36, no. 2, pp. 309-313, 2008.

[28] W. Bao, F. Sun, and G. W. Wei, "Numerical methods for the generalized Zakharov system," Journal of Computational Physics, vol. 190, no. 1, pp. 201-228, 2003.

[29] W. Bao and F. Sun, "Efficient and stable numerical methods for the generalized and vector Zakharov system," SIAM Journal on Scientific Computing, vol. 26, no. 3, pp. 1057-1088, 2005.

Fayyaz Ahmad, (1,2,3) Shafiq Ur Rehman, (4) Malik Zaka Ullah, (2,5) Hani Moaiteq Aljahdali, (6) Shahid Ahmad, (7) Ali Saleh Alshomrani, (5) Juan A. Carrasco, (8) Shamshad Ahmad, (9) and Sivanandam Sivasankaran (10)

(1) Departament de Fisica i Enginyeria Nuclear, Universitat Politecnica de Catalunya, Eduard Maristanu 10, 08019 Barcelona, Spain

(2) Dipartimento di Scienza eAlta Tecnologia, Universita dell'Insubria, Via Valleggio 11, 22100 Como, Italy

(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 Politecnica de Catalunya, Diagonal 647, Planta 9, 08028 Barcelona, Spain

(9) Heat and Mass Transfer Technological Center, Universitat Politecnica de Catalunya, Colom 11, 08222 Terrassa, Spain

(10) Institute of Mathematical Sciences, University of Malaya, Kuala Lumpur 50603, Malaysia

Received 30 July 2016; Revised 10 November 2016; Accepted 21 November 2016; Published 3 May 2017

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

Caption: Figure 2: Initial guess and solutions of Bratu problem for different values of [alpha].

Caption: Figure 3: Initial guess and solutions of Frank-Kamenetskii problem with different values of [alpha].

Caption: Figure 4: Solution and error components for (17).

Caption: Figure 5: Solution and error components for (30).

Caption: Figure 6: Solution and error components for (32).

Caption: Figure 7: Solution and error components for (35).

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

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

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

Caption: Figure 11: Error components for the nonlinear 3D Poisson equation.

Caption: Figure 12: Solution and error components for complex generalized Zakharov equation (57).

Caption: Figure 13: Solution and error components for complex generalized Zakharov equation (57).

Caption: Figure 14: Solution and error components for complex generalized Zakharov equation (57).

Caption: Figure 15: Solution and error components for complex generalized Zakharov equation (60).

Caption: Figure 16: Solution and error components for complex generalized Zakharov equation (60).

Caption: Figure 17: Solution and error components for complex generalized Zakharov equation (60).

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

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

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

Caption: Figure 21: Error components for the two-dimensional nonlinear reaction-diffusion equation.

Caption: Figure 22: Solution for Ostrovsky-Hunter equation.
```
TABLE 1: Computational costs of different iterative methods.

Methods                       Computational cost

MNR           1/3 [n.sup.3] ([beta] + m - 1/2 + 1/2) [n.sup.2] +
[(m[alpha] + lm - m + 1/6 - 1/2).sup.n]

HJ           1/3 [n.sup.3] (2[beta] - 3/2 + 1/2 + 3m) [n.sup.2] +
[(m[alpha] - [alpha] + 7/6 - 3/2l + 2lm).sup.n]

FTUC       1/3 [n.sup.3] (2[beta] - 7/2 + 1/2 + 3m) [n.sup.2] + [(m
[alpha] - [alpha] - [alpha] + 19/6 - 5/2; - m + 2lm).sup.n]

MSF          1/3 [n.sup.3] (2[beta] - 7/2 + 1/2 + 5m) [n.sup.2] +
[(m[alpha] + [gamma]  + 19/6 - 3/2l - 2m + 3lm).sup.n]

NMIM         1/3 [n.sup.3] (2[beta] -49/2 + 1/2 + 9m) [n.sup.2] +
[[alpha]m - 17/6 - 25/2l + 5lm).sup.n]

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 m steps.

Methods                 Difference of computational costs

MNR-NMIM     4nam - [[beta]n.sup.2] - [4n.sup.2]m - 10n[alpha] + 2ln
- 52m + [14n.sup.2] + 13n

HJ-NMIM      13/2 [n.sup.2] - 3/2 [n.sup.2]m + 3/2 n[alpha]m - 13/2
n[alpha] + 4n

FTUC-NMIM   [14n.sup.2] - [4n.sup.2]m + 2/3 n[alpha]m - 10/3 n[alpha]
25n/3 + 16/3 ln - 5/3 nm - 5/3 lmn

MSF-NMIM         -2/3 [n.sup.2], + 8/3 [n.sup.2] + 2/3 n[alpha]m
-11/3 n[alpha] + [gamma]n + 40n/3 - 10/3 nm

TABLE 3: Conditions and bounds on m so that, for the same convergence
order, the computational cost of NMIM is smaller.

Computational              Bounds on m               Conditions
cost

[C.sub.MNR] >     m < n[beta] + 10[alpha] - 2l     [alpha] < n +
[C.sub.NMIM]      - 14n -13/-14n + 4[alpha] - 5          5/4

[C.sub.HJ] >      m < 1/3 - 13n + 13[alpha] - 8/    [alpha] < n
[C.sub.NMIM]               -n + [alpha]

[C.sub.FTUC] >    m < - 42n - 16 l + 10[alpha]      [alpha] < 6n
[C.sub.NMIM]       - 25/-12 - 5l + 2[alpha] -5      + 5/2 + 5/2 l

[C.sub.MSF] >       m < - 1/ 8n - 11[alpha] +      [alpha] < n + 5
[C.sub.NMIM]      3[gamma] + 40/-2 + [alpha] -5

TABLE 4: Verification of convergence order.

[[parallel]F([x.sub.k])[parallel].sub.[infinity]]

Iter.            m = 4        m = 5         m = 6          m = 7

Problem 1
(n = 4)

1              6.56e - 8    1.09e - 11    1.82e- 15     3.04e - 19
2             6.31e - 68   5.74e - 158   8.66e - 286    2.17e - 451
3             4.44e - 608  6.90e - 2206  6.32e - 5422  6.67e - 10823

CCO                9            14            19            24

Problem 2
(n = 10)

1              7.31e - 3    6.45e - 4     5.83e - 5      5.29e - 6
2             1.54e - 21    1.38e - 47    5.42e - 84    8.20e- 131
3             1.31e- 189   5.63e - 659   1.35e- 1585   3.06e - 3126

CCO                9            14            19            24

Problem 3      7.09e - 5    2.12e - 7     6.36e - 10     1.90e- 12
(n = 100)

2             5.45e - 43    1.69e - 94   1.53e - 166    1.64e - 258
3             6.87e - 348  8.80e- 1143   7.95e - 2673  8.85e - 5184

CCO                8            12            16            20

Problem 4      1.24e - 3    8.47e - 6     8.17e - 9     9.64e - 13
(n = 5)

3             2.73e - 19    2.08e - 54   1.45e - 120    2.27e - 229
4             6.73e - 148  1.79e - 645   7.35e- 1920   5.31e - 4575
CCO              8.21          12.2          16.1          20.1

Problem 5
(n = 3)

1              3.77e - 6    8.57e - 10    1.94e - 13     4.39e- 17
2             2.90e - 67    1.71e- 151   6.21e - 269    1.46e - 419
3             3.09e - 581  2.79e - 1984  1.71e - 4737  9.99e - 9298
CCO              8.41          12.9          17.5          22.1
Theoretical        9            14            19            24
convergence
order
(5m - 11)

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

[[parallel]F([X.sub.k]))[parallel].sub.[infinity]]

P

m         2             3            4            5

10    2.69e - 17   3.70e - 14    1.04e - 12   1.03e - 11
20    1.36e - 46   7.06e - 41    9.33e - 39   1.22e - 37
30    9.33e - 76   1.68e - 67    7.56e - 65   1.97 - 63
40    6.16e - 105  4.53e - 94     6.19e-91    2.26e - 89
50    4.08e - 134  1.21e - 120   5.11e - 117  2.18e - 115
60    2.71e - 163  3.22e - 147   4.24e - 143  2.00e - 141
70    1.79e - 192  8.555e - 174  3.52e - 169  1.83e - 167
80    1.19e - 221  2.27e - 200   2.92e-195    1.69e-193
90    7.89e - 251  6.03e - 227   2.43e - 221  1.58e - 219
100   5.23e - 280   1.60e-253    2.02e - 247  1.47e - 245

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

[[parallel]F([X.sub.k)])[parallel].sub.infinity]]

[alpha]

m         1            2            3

10    7.68e - 41   4.10e - 31   1.97e - 28
20    2.20e - 93   1.01e - 70   9.65e - 64
30    7.02e - 146  2.50e - 110  4.72e - 99
40    2.23e - 198  6.17e - 150  2.30e - 134
50    7.13e - 251  1.52e - 189  1.12e - 169
60    2.27e - 303  3.75e - 229  5.52e - 205
70    7.30e - 356  9.27e - 269  2.69e - 240
80    2.30e - 408  2.28e - 308  1.32e - 275
90    7.40e - 461  5.70e - 348  6.50e - 311
100   2.30e - 513  1.40e - 387  3.20e - 346

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

[[parallel]F([X.sub.k])][parallel].sub.[infinity]]

[alpha]

m         1           1.1          1.2          1.3

10    7.09e - 41   1.89e - 39   2.63e - 38   1.80e - 37
20    2.46e - 82   2.36e - 79   6.21e - 77   4.05e - 75
30    8.58e - 124  2.95e - 119  1.46e - 115  9.09e - 113
40    2.98e - 165  3.69e - 159  3.45e - 154  2.03e - 150
50    1.03e - 206  4.61e - 199  8.15e - 193  4.57e - 188
60    3.61e - 248  5.76e - 239  1.92e - 231  1.02e - 225
70    1.25e - 289  7.19e - 279  4.53e - 270  2.29e - 263
80    4.40e - 331  9.00e - 319  1.10e - 308  5.15e - 301
90    1.50e - 372  1.10e - 358  2.50e - 347  1.20e - 338
100   5.30e - 414  1.40e - 398   .00e - 386  2.60e - 376

TABLE 8: Nonlinear Schrodinger equations.

Schroodinger equation                           Analytical Solution

i[[partial derivative].sub.t][psi] +               [psi](x, t) =
[[partial derivative].sub.xx][psi] + 2           [e.sup.i(x + t)]
[[absolute value of ([psi])].sup.2]
[psi] = 0

i[[partial derivative].sub.t] [psi] +              [psi](x, t) =
[[partial derivative].sub.xx][psi] - 2           [e.sup.i(x - 3t)]
[[absolute value of ([psi])].sup.2]
[psi] - cos [(x).sup.2] [psi] = 0

i[[partial derivative].sub.t] [psi] +            [psi](x, t) =
[[partial derivative].sub.xx][psi] - 2           [e.sup.-3it/2)]
[[absolute value of ([psi])].sup.2]                   sin(x)
[psi] - cos [(x).sup.2] [psi] = 0

i[[partial derivative].sub.t] [psi] +            [psi](x, t) =
[[partial derivative].sub.xx][psi] + 2            [e.sup.i(x +
[[absolute value of ([psi])].sup.2]              (1 - 2[delta])t)
[psi] - 2[delta][psi] = 0

TABLE 9: Norm of the error in the solution for equation (17).

[[parallel][psi] - [[psi].sub.analytical][parallel]
.sub.[infinity]]

Grid size

m                               5 x 30      10 x 30        15x30

1    Legendre polynomials     5.49e - 01   6.21e - 01   6.41e - 01
4    ([theta], [phi]) =       5.36e - 02   2.35e - 01   2.44e - 01
8    (0, 0)                   2.85e - 04   4.55e - 04   4.97e - 04
12                            2.85e - 04   5.81e - 08   6.38e - 08
16                            2.85e - 04   5.69e - 10   4.38e - 12
20                            2.85e - 04   6.13e - 11   5.95e - 13

1    ([theta], [phi]) =       6.86e - 01   Divergent    6.40e - 01
4    (1/2, 0)                 1.96e - 01                2.48e - 01
8                             1.64e - 03                4.94e - 04
12                            1.59e - 03                6.16e - 08
16                            1.59e - 03                6.02e - 12
20                            1.59e - 03                2.16e - 13

1    Chebyshev polynomials    5.24e - 01    6.25e-01    6.41e - 01
4    of second kind           5.90e - 02    2.39e-01    2.46e - 01
8    ([theta], [phi]) =       1.14e - 03   5.05e - 04   5.40e - 04
12   (1/2,1/2)                1.14e - 03   7.73e - 08   8.00e - 08
16                            1.14e - 03    1.90e-09    5.28e - 12
20                            1.14e - 03   1.90e - 09   4.99e - 13

1    Chebyshev polynomials    5.95e - 01   6.20e - 01   6.41e - 01
4    of first kind            3.73e - 02   2.37e - 01   2.45e - 01
8    ([theta], [phi]) =       1.52e - 03   4.24e - 04   4.81e - 04
12   (-1/2,-1/2)              1.52e - 03   5.22e - 08   5.81e - 08
16                            1.52e - 03   4.90e - 10   1.71e - 12
20                            1.52e - 03   4.90e - 10   9.31e - 14

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

[[parallel][psi] - [[psi].sub.analytical]
[parallel].sub.[infinity]]

Grid size

m                             5 x 30        10 x 30       15 x 30

1    Legendre               1.28e + 00    1.34e + 00    1.34e + 00
4    polynomials            4.27e - 01    9.65e - 01    9.72e - 01
8    ([theta], [phi]) =     1.25e - 03    2.15e - 03    2.18e - 03
12   (0,0)                  1.46e - 04    9.17e - 07    1.05e - 06
16                          1.46e - 04    3.32e - 10    4.61e - 10
20                          1.46e - 04    2.90e - 11    1.85e - 13

1    ([theta], [phi]) =     1.27e + 00    1.35e + 00    1.34e + 00
4    (1/2, 0)               4.34e - 01    9.46e - 01    9.72e - 01
8                           1.63e - 03    2.18e - 03    2.18e - 03
12                          1.52e - 03    3.89e - 07    8.77e - 07
16                          1.52e - 03    5.28e - 09    1.95e - 10
20                          1.52e - 03    3.73e - 09    1.34e - 13

1    Chebyshev              1.26e + 00    1.34e + 00    1.34e + 00
4    polynomials of         3.98e - 01    9.24e - 01    9.40e - 01
8    second kind            2.63e - 03    1.35e - 03    1.63e - 03
12   ([theta], [phi]) =     1.00e - 03    2.63e - 06    2.85e - 06
16   (1/2, 1/2)             1.00e - 03    2.35e - 08    1.65e - 08
20                          1.00e - 03    7.88e - 10    6.29e - 11

1    Chebyshev              1.33e + 00    1.34e + 00    1.34e + 00
4    polynomials of         3.90e - 01    9.80e - 01    9.80e - 01
8    first kind             1.03e - 03    2.15e - 03    2.19e - 03
12   ([theta], [phi]) =     1.03e - 03    6.56e - 07    7.42e - 07
16   (-1/2,-1/2)            1.03e - 03    3.08e - 10    1.57e - 10
20                          1.03e - 03    2.95e - 10    5.77e - 14

TABLE 11: Norm of the error in the solution for equation (32).

[[parallel][psi] - [[psi].sub.analytical]
[parallel].sub.[infinity]]

Grid size

m                            5 x 30      10 x 30      15 x 30

1    Legendre              9.98e - 02   8.73e - 02   9.66e - 02
4    polynomial            7.62e - 04   3.16e - 04   3.25e - 04
8    ([theta], [phi]) =    4.27e - 05   3.64e - 11   3.12e - 11
12   (0,0)                 4.27e - 05   3.16e - 11   1.89e - 13
16                         4.27e - 05   3.16e - 11   1.88e - 13
20                         4.27e - 05   3.16e - 11   1.89e - 13

1    ([theta], [phi]) =    9.71e - 02   9.33e - 02   9.67e - 02
4    (1/2, 0)              1.18e - 03   3.37e - 04   3.36e - 04
8                          7.68e - 04   3.23e - 08   2.93e - 11
12                         7.68e - 04   1.89e - 10   1.07e - 13
16                         7.68e - 04   1.87e - 10   1.06e - 13
20                         7.68e - 04   1.87e - 10   1.07e - 13

1    Chebyshev             9.47e - 02   9.17e - 02   9.65e - 02
4    polynomials of        9.17e - 04   3.14e - 04   3.38e - 04
8    second kind           6.70e - 04   9.13e - 11   8.39e - 11
12   ([theta], [phi]) =    6.70e - 04   8.92e - 11   1.19e - 13
16   (1/2,1/2)             6.70e - 04   8.92e - 11   1.19e - 13
20                         6.70e - 04   8.92e - 11   1.19e- 13

1    Chebyshev             1.01e - 01   9.01e - 02   9.61e - 02
4    polynomial of         1.31e - 03   2.99e - 04   3.26e - 04
8    first kind            8.53e - 04   4.90e - 11   2.27e - 11
12   ([theta], [phi]) =    8.53e - 04   5.59e - 11   3.28e - 14
16   (-1/2, -1/2)          8.53e - 04   5.59e- 11    3.24e - 14
20                         8.53e - 04   5.59e- 11    3.20e - 14

TABLE 12: Norm of the error in the solution for equation (35).

[[parallel][psi] - [[psi].sub.analytical]
[parallel].sub.[infinity]]

Grid size

m                            5 x 30      10 x 30      15 x 30

1    Legendre              6.44e - 01   7.84e - 01   8.01e - 01
4    polynomials           6.52e - 02   4.18e - 02   4.18e - 02
8    ([theta], [phi]) =    2.77e - 04   1.36e - 06   1.33e - 06
12   (0,0)                 2.75e - 04   1.70e - 10   1.81e- 10
16                         2.75e - 04   4.14e - 11   4.87e - 13
20                         2.75e - 04   4.14e - 11   4.67e - 13

1    ([theta], [phi]) =    6.51e - 01                8.08e - 01
4    (1/2, 0)              4.52e - 02                4.15e - 02
8                          1.50e - 03                1.25e - 06
12                         1.51e - 03   Divergent    1.17e - 10
16                         1.51e - 03                1.94e - 13
20                         1.51e - 03                1.99e - 13

1    Chebyshev             7.48e - 01   7.85e - 01   8.04e - 01
4    polynomials of        5.30e - 02   4.24e - 02   4.22e - 02
8    second kind           1.02e - 03   2.20e - 06   2.61e - 06
12   ([theta], [phi]) =    1.02e - 03   2.67e - 09   1.00e - 09
16   (1/2,1/2)             1.02e - 03   2.25e - 09   4.86e - 13
20                         1.02e - 03   2.25e - 09   1.96e- 13

1    Chebyshev             7.03e - 01   7.88e - 01   8.02e - 01
4    polynomials of        3.66e - 02   4.20e - 02   4.17e - 02
8    first kind            1.46e - 03   1.27e - 06   1.18e - 06
12   ([theta], [phi]) =    1.46e - 03   5.96e - 10   2.39e - 11
16   (-1/2, -1/2)          1.46e - 03   5.84e - 10   8.56e - 14
20                         1.46e - 03   5.84e - 10   8.49e - 14

TABLE 13: Norm of the error in the Klein Gordon equation as a function
of m for several values for [theta] and [phi] and grids.

[parallel]u - [u.sub.analytical]][parallel].sub.[infinity]]

Grid size

m                          30 x 30      60 x 30      120 x 30

1   Legendre              3.80e - 01   4.10e - 01   4.09e - 01
2   polynomials           7.37e - 02   8.39e - 02   8.34e - 02
3   ([theta], [phi]) =    1.18e - 02   8.73e - 03   8.65e - 03
4   (0, 0)                8.53e - 03   1.42e - 05   1.40e - 06
5                         8.53e - 03   1.42e - 05   4.37e - 09
6                         8.53e - 03   1.42e - 05   4.56e - 09

1   ([theta], [phi]) =    3.98e - 01   4.10e - 01   4.07e - 01
2   (1/2, 0)              7.64e - 02   8.40e - 02   8.28e - 02
3                         1.01e - 02   8.76e - 03   8.55e - 03
4                         7.67e - 03   3.55e - 05   1.38e - 06
5                         7.67e - 03   3.55e - 05   1.60e - 09
6                         7.67e - 03   3.55e - 05   1.45e - 09

1   Chebyshev             3.81e - 01   4.10e - 01   4.09e - 01
2   polynomials of        7.41e - 02   8.39e - 02   8.34e - 02
3   second Kind           1.12e - 02   8.73e - 03   8.65e - 03
4   ([theta], [phi]) =    7.85e - 03   1.29e - 05   1.41e - 06
5   (1/2, 1/2)            7.85e - 03   1.29e - 05   2.38e - 09
6                         7.85e - 03   1.29e - 05   2.23e - 09

1   Chebyshev             3.78e - 01   4.10e - 01   4.09e - 01
2   polynomials           7.33e - 02   8.38e - 02   8.34e - 02
3   of first kind         1.25e - 02   8.73e - 03   8.65e - 03
4   ([theta], [phi]) =    9.28e - 03   1.58e - 05   1.41e - 06
5   (-1/2, -1/2)          9.28e - 03   1.58e - 05   3.63e - 10
6                         9.28e - 03   1.58e - 05   3.40e - 10

TABLE 14: Norm of the error for the 2D nonlinear wave equation as a
function of m for several pair of values for [theta] and [phi] and
several grids.

m                           Grid size

5 x 5 x 20   10 x 10 x 20   20 x 20 x 20

1    Legendre               4.21e - 02    Divergent       Divergent
2    polynomials            3.85e - 03
3    ([theta], [phi]) =     2.03e - 04
9    (0,0)                  1.09e - 04

1    ([theta], [phi]) =     4.45e - 02    5.12e - 02     5.14e - 02
2    (1/2,0)                3.50e - 03    3.74e - 03     1.53e - 02
5                           7.34e - 04    5.74e - 08     2.16e - 04
7                           7.34e - 04    1.40e - 09     3.76e - 06
9                           7.34e - 04    1.40e - 09     6.56e - 08

1    Chebyshev              4.40e - 02    5.13e - 02
2    polynomials of         3.43e - 03    3.81e - 03
3    second kind            5.98e - 04    1.62e - 04      Divergent
4    ([theta], [phi]) =     4.39e - 04    8.01e - 09
5    (1/2, 1/2)             4.39e - 04    7.64e - 10

1    Chebyshev              3.88e - 02    5.02e - 02     5.13e - 02
2    polynomials of         3.91e - 03    3.85e - 03     3.80e - 03
3    first kind             5.94e - 04    1.54e - 04     1.79e - 03
4    ([theta], [phi]) =     6.60e - 04    4.54e - 09     4.28e - 05
5    (-1/2, -1/2)           6.60e - 04    6.78e - 10     7.94e - 07
9                           6.60e - 04    6.78e - 10     5.39e - 11

TABLE 15: Norm of the error for the 3D nonlinear wave equation as a
function of m for several values for [theta] and [phi] and grids.

[parallel]u - [u.sub.analytical][parallel].sub.[infinity]]

Grid size

m                           5 x 5 x      6 x 6 x      7 x 7 x
5 x 12       6 x 12       7 x 12

1    Legendre              9.95e - 01   1.00e + 00   1.14e + 00
3    polynomials           1.14e - 02   8.53e - 03   1.09e - 02
5    ([theta], [phi]) =    1.25e - 01   1.12e - 01   1.73e - 01
7    (0,0)                 5.53e - 06   4.12e - 07   1.73e - 06
8                          5.54e - 06   7.54e - 08   4.24e - 09

1    ([theta], [phi]) =    1.00e + 00   1.00e + 00
2    (1/2, 0)              1.00e - 01   8.68e - 02
4                          1.31e + 00   1.29e + 00   divergent
7                          2.24e - 05   1.99e - 06
8                          2.24e - 05   2.11e - 06

1    Chebyshev             9.98e - 01   9.99e - 01   1.32e + 00
2    polynomials of        1.02e - 01   8.06e - 02   1.08e - 01
5    second kind           1.22e - 01   1.18e - 01   2.46e - 01
7    ([theta], [phi]) =    1.81e - 05   1.84e - 06   3.72e - 05
8    (1/2, 1/2)            1.81e - 05   1.76e - 06   4.87e - 08

1    Chebyshev             9.95e - 01                1.09e + 00
2    polynomials of        9.98e - 02                1.04e - 01
5    first kind            1.31e - 01   divergent    1.51e - 01
7    ([theta], [phi]) =    1.92e - 05                4.37e - 07
8    (-1/2, -1/2)          1.92e - 05                3.27e - 08

TABLE 16: Norm of the error for the 3D Nonlinear Poisson equation as
a function of m for several values for [theta] and [phi] and grids.

[[parallel]u - [u.sub.analytical][parallel].sub.[infinity]]

Grid size

m                          8 x 8 x 8    10 x 10 x 10   12 x 12 x 12

1    Legendre              4.29e - 02    4.44e - 02     4.52e - 02
4    polynomials           5.54e - 08    6.03e - 08     6.30e - 08
5    ([theta], [phi]) =    2.85e - 11    8.77e - 12     9.27e - 12
6    (0, 0)                2.95e - 11    1.51e - 14     6.44e - 15

1    ([theta], [phi]) =    4.37e - 02    4.51e - 02     4.57e - 02
3    (1/2, 0)              6.80e - 04    7.09e - 04     7.24e - 04
5                          1.33e - 09    9.88e - 12     9.75e - 12
6                          1.33e - 09    9.39e - 13     7.66e - 15

1    Chebyshev             4.33e - 02    4.46e - 02     4.53e - 02
2    polynomials of        5.39e - 03    5.63e - 03     5.76e - 03
4    secondkind            5.58e - 08    6.10e - 08     6.35e - 08
6    (0,ty) =              9.59e - 10     6.78e-13      1.88e - 14
(1/2, 1/2)

1    Chebyshev             4.23e - 02    4.41e - 02     4.50e - 02
3    polynomials of        6.33e - 04    6.78e - 04     7.01e - 04
5    first kind            5.63e - 10    8.53e - 12     9.17e - 12
6    ([theta], [phi]) =    5.63e - 10    3.26e - 13     8.55e - 15
(-1/2, -1/2)

TABLE 17: Maximum absolute error in the solution of the system of
nonlinear equations for complex generalized Zakharov equation and
varying m.

m         E
1     8.56e - 1
2      2.72e0
3      4.16e0
4      8.92e0
5      1.12e1
6      3.40e1
7      2.04e1
8      2.83e1
9      4.05e0
10     2.34e0
11    777e - 1
12    3.25e - 1
13    1.14e - 1
14    6.20e - 2
15    2.79e - 2
16    1.85e - 2
17    8.11e - 3
18    4.44e - 3
19    1.79e - 3
20    8.86e - 4
21    3.91e - 4
22    1.66e - 4
23    5.98e - 5
24    1.97e - 5
25    5.32e - 6
26    1.20e - 6
27    4.12e - 7
28    1.36e - 7
29    6.95e - 8
30    2.95e - 8
31    9.29e - 9
32    3.42e - 9
33    1.43e - 9
34    5.57e - 10
35    2.62e - 10
36    1.26e - 10
37    8.65e - 11

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

m         E

1    9.2329e - 1
2    7.0559e - 1
3    5.726e - 1
4    6.3466e - 2
5    4.3852e - 3
6    2.3221e - 4
7    1.5501e - 5
8    4.2396e - 7
9    2.8906e - 8
10   7.7742e - 10
11   8.6019e - 11

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

m        E

1    2.63e - 01
2    1.33e - 01
3    3.85e - 02
4    1.47e - 04
5    7.43e - 07
6    2.15e - 09
7    7.53e - 12
8    5.59e - 14

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

m        E

1    7.28e - 02
2    2.12e - 02
3    4.99e - 03
4    1.09e - 06
5    9.20e - 10
6    6.53e - 13
7    2.24e - 14

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
m.

m (1 Iter.)       E

1             1.49e + 00
2             1.12e + 00
3             1.07e + 00
4             4.92e - 01
5             3.73e - 01

m (2 Iter.)       E

1             1.51e - 02
2             7.84e - 03
3             6.72e - 04
4             5.17e - 07
5             1.04e - 09
6             2.35e - 12
7             5.04e - 14

TABLE 22: Maximum absolute error in the solution of the system of
nonlinear equations for two-dimensional nonlinear reaction-diffusion
equation and varying m.

m         E

1     7.82e - 02
2     1.19e - 02
3     3.63e - 03
4     4.04e - 05
5     9.89e - 07
6     2.65e - 08
7     7.48e - 10

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

m     [[parallel]F(X)[parallel].sub.[infinity]]

1                         2
2                       2.5963
3                       2.797
4                      0.72034
5                      0.15165
6                      0.025873
7                      0.002996
8                     0.00020443
9                    8.7968e - 06
10                   1.4609e - 07
11                   1.3542e - 08
12                   8.8853e - 10
13                   8.5753e - 11

TABLE 24: Solution for Ostrovsky-Hunter equation.

X                         t

0.1        0.2        0.3        0.4

0           0.094587   0.17222    0.20233    0.21608
6[pi]/11    0.96556    0.94696     0.8702    0.73399
[pi]        -0.10438   -0.22921   -0.41929   -0.67358
16[pi]/11   -0.98195   -0.96996   -0.97341   -0.97748
```