Multistep Methods of the Hybrid Type and Their Application to Solve the Second Kind Volterra Integral Equation

Abstract

There are some classes of methods for solving integral equations of the variable boundaries. It is known that each method has its own advantages and disadvantages. By taking into account the disadvantages of known methods, here was constructed a new method free from them. For this, we have used multistep methods of advanced and hybrid types for the construction methods, with the best properties of the intersection of them. We also show some connection of the methods constructed here with the methods which are using solving of the initial-value problem for ODEs of the first order. Some of the constructed methods have been applied to solve model problems. A formula is proposed to determine the maximal values of the order of accuracy for the stable and unstable methods, constructed here. Note that to construct the new methods, here we propose to use the system of algebraic equations which allows us to construct methods with the best properties by using the minimal volume of the computational works at each step. For the construction of more exact methods, here we have proposed to use the multistep second derivative method, which has comparisons with the known methods. We have constructed some formulas to determine the maximal order of accuracy, and also determined the necessary and sufficient conditions for the convergence of the methods constructed here. One can proved by multistep methods, which are usually applied to solve the initial-value problem for ODE, demonstrating the applications of these methods to solve Volterra integro-differential equations. For the illustration of the results, we have constructed some concrete methods, and one of them has been applied to solve a model equation.

1. Introduction

It is known that many problems of the natural sciences are reduced to the solving of integral equations of variable boundaries, which are called integral equations of Volterra type. Vito Volterra (proud Italian) fundamentally investigated these equations and also reduced the mathematical models of many problems of the natural sciences to solve these integral equations. As is known, to solve Volterra integral equations is one of the basic directions in modern mathematics. For objectivity, let us note that scientists have met with the need to solve integral equations with variable boundaries before Vito Volterra (see, for example, [1,2,3,4,5,6]). Now, let us consider the following integral equation of Volterra type:

(1)$y\left(x\right)=f\left(x\right)+{\displaystyle \underset{\alpha \left(x\right)}{\overset{\beta \left(x\right)}{\int }}}K\left(x,s,y\left(s\right)\right)ds,xϵ\left[x_0,X\right].$

One of the popular Volterra integral equations can be written as:

(2)$y\left(x\right)=f\left(x\right)+{\displaystyle \underset{x_0}{\overset{x}{\int }}}K\left(x,s,y\left(s\right)\right)ds,x_0\le s\le x\le X.$

This equation is taken as known, if given the functions $f\left(x\right)$ and $K\left(x, s, z\right).$

Suppose that the given functions $f\left(x\right)$ and $K\left(x, s, z\right)$ are sufficiently smooth and Equation (2) has the unique solution $y\left(x\right)$, which is defined on the segment $\left[x_0, X\right]$. It follows that the solution $y\left(x\right)$ is also a sufficiently smooth function. For the construction of numerical methods to solve Equation (2), let us divide the segment $\left[x_0, X\right]$ to $N$ equal parts by using nodes $x_{i+1}=x_i+h \left(i=0, 1, \dots , N\right)$. Here, $0<h$ is the step size.

As is known, for the solving of Equation (2) there are numerous methods constructed by different authors. There exists one-step and multistep methods constructed for the solving of Equation (2). Let us note that some authors, for the solving of Equation (2), have proposed to use the spline function or collocation methods (see, for example, [7,8,9,10,11]).

There are many works dedicated to the solving of integral equations, which have used the quadrature methods (in [12] for the calculation of definite integrals proposed to use the new way). Note that for solving Volterra integral equations, many authors constructed methods which are different from the above noted (see, for example, [13,14,15,16]). It is known that in this case the number of calculations increases when going from the current point to the next. By taking into account this property, in [13] they have constructed a method which is released from the indicated disadvantages. By generalization of this method, here we have constructed more exact stable methods, which we have applied to solve Equation (2). For the presentation of the essence of these methods, let us consider the partial case of Equation (2) which is obtained when replacing $K\left(x, s, y\right)=\phi \left(s, y\right).$ In this case Equation (2) can be written as the following:

(3)$y\left(x\right)=f\left(x\right)+{\displaystyle \underset{x_0}{\overset{x}{\int }}}\phi \left(s,y\left(s\right)\right)ds.$

It is not hard to understand that solving this equation is equivalent to solving the following initial-value problem for ODEs of the first order:

(4)$y^{\prime }\left(x\right)=f^{\prime }\left(x\right)+\phi \left(x,y\left(x\right)\right),y\left(x_0\right)=f\left(x_0\right).$

It follows from here that the solution of the integral equation of (3) and the initial-value problem (4) can be found by one and the same method. It is not the only case, when the initial-value problem for ODEs and Volterra integral equations can be solved by one and the same methods. To show this, let us consider the following case, when the function of $K\left(x, s, y\right)$ is degenerate and can be presented in the following form:

(5)$K\left(x,s,y\right)={\displaystyle \sum _{j=1}^m}{a}_j\left(x\right)b_j\left(s,y\right).$

By taking this in Equation (2), we receive:

$y\left(x\right)=f\left(x\right)+{\displaystyle \sum _{j=1}^m}{a}_j\left(x\right)v_j\left(x\right).$

The function $v_j\left(x\right)$, $j=0, 1, \dots , m$ can be determined as the solution of the following problem:

(6)$v_j{}^{\prime }\left(x\right)=b_j\left(x,y\left(x\right)\right),v_j\left(x_0\right)=0\left(j=1,2,\dots ,m\right).$

It is clear that by solving the system of ODEs (6), one can find the values of the function $y\left(x\right)$ at the nodes (mesh points) by using the solution of the system (6). Note that in this case this system of ODEs and the integral equation of (2) can be solved by one and the same methods. Thus, there are some domains in which the integral equation of (2) and the initial-value problem for ODEs can be solved by one and the same methods. Here, it is shown that this domain can be extended and the error received in this case can be estimated (see, for example, [17,18,19]).

It is not difficult to prove that by using the Lagrange interpolation polynomial the function $K\left(x, s, y\right)$ can be presented as:

(7)$K\left(x,s,y\right)={\displaystyle \sum _{j=1}^k}{l}_j\left(x\right)b_j\left(s,y\right)+R_k\left(x\right),$

where $l_j\left(x\right)$  $\left(i=1, 2, \dots , m\right)$ are the basic Lagrange function and $R_m\left(x\right)$ is the remainder term. By comparison of the equality of (5) and (7) we receive some connection between them.

2. Construction of Multistep Methods to Solve Both Equations (2) and (4)

Let us note that the known multistep method with constant coefficients can be applied to solve the Volterra integral equation. In the result of which, one can constructed by the following method (see, for example, [17,18,19]):

(8)${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}={\displaystyle \sum _{i=0}^k}{\alpha }_i{f}_{n+i}+h{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}K\left(x_{n+j},x_{n+i},y_{n+i}\right).$

For the construction methods with the improved properties, here we have used the generalization of the multistep methods, which can be written as the following (see, for example, [20,21,22,23,24,25,26,27,28]):

(9)${\displaystyle \sum _{i=0}^m}{\alpha }_i{y}_{n+i}=h{\displaystyle \sum _{i=0}^k}{\beta }_i{y^{\prime }}_{n+i}\left(n=0,1,\dots ,N-l;l=\max \left(m,k\right)\right).$

For the value $m\ge k$ from the method of (9), it follows the known multistep methods, but for the value $m<k$ it follows advanced methods (formally), but in reality these methods do not depend on each other. Therefore, each of them is an independent object of investigation. The stable advanced method is more accurate $(p\le k+l+1,$ for $k\ge 3l,$ and $m=k-l,$ here, $p$ is the degree and $k$ is the order of finite-difference method (9)) (see, for example, [23,24]). Let us note that by Dahlquist’s laws, there exists stable methods of type (8) which have the degree $p_{max}=2\left[k/2\right]+2$. Here, we use the conceptions of the stability, degree, and order, defined by Dahlquist (see, for example, [21,22,23,24,25,26,27,28]). By the above-described way, we find that the stable methods of the advanced type are more exact than the stable multistep methods. However, unstable multistep methods are more exact than the advanced method. Namely, $p\le 2k$ for method (8) and $p\le 2k-l$ for method (9).

Advanced methods have been constructed by Kouella for the calculation of the return of Holliley’s comet. Note that some advanced concrete methods have been constructed by known scientists such as Laplace, Steklov, etc. However, all stable advanced methods constructed by different specialists had the degree $p\le 2\left[k/2\right]+2$. That is, they obeyed the law of Dahlquist. The advantages and disadvantages of advanced methods have been shown in [23,24], and for the correction of some of the disadvantages of the advanced method there have been constructed special predictor–corrector methods (see, for example, [29]). Note that if $m>k$ in this case method (9) will be explicit and the maximal value for stable explicit methods can be found by the formula $p\le m$. Thus, it is proven that stable advanced methods are more exact than the explicit and implicit methods of type (9). Now, let us investigate these methods by using their other properties.

If method (9) is applied to solve the problem (4), then we receive:

(10)${\displaystyle \sum _{i=0}^m}{\alpha }_i\left(y_{n+i}-f_{n+i}\right)=h{\displaystyle \sum _{i=0}^k}{\beta }_i{\phi }_{n+i}.$

where $f_m=f\left(x_m\right),{\phi }_m=\phi \left(x_m,y_m\right),\mathrm{and}{\alpha }_j,{\beta }_i\left(j=0,1,\dots ,m;i=0,1,\dots ,k\right)$ are the coefficients of method (10) or (9).

Let us input $m=k$. In this case, for the value ${\alpha }_k\ne 0$, we receive the implicit method if ${\beta }_k\ne$ 0. If we compare these methods, then we find that the explicit $\left({\beta }_k=0\right)$ methods can be applied to solve some problems. However, in the application of implicit methods arises some difficulties for elimination, of which here it is proposed to use the predictor–corrector methods. Now, let us consider the application of the advanced methods. For this, we have $k=m+l$ and ${\alpha }_{k-l}\ne 0$. In this case, method (10) can be written as:

(11)${\displaystyle \sum _{i=0}^m}{\alpha }_i\left(y_{n+i}-f_{n+i}\right)=h{\displaystyle \sum _{i=0}^{m-1}}{\beta }_i{\phi }_{n+i}+h{\displaystyle \sum _{i=m}^{m+l}}{\beta }_i{\phi }_{n+i}.$

In the application of this method arises some difficulties related to the calculation of the second part, which is located on the right hand side of the equality of (11).

It is evident that for the calculation of the second sum on the right hand side of equality (11), we need to define the values $y_{n+m,} y_{n+m+1},\dots , y_{n+m+l}.$ This difficulty can be solved by using some methods for the calculation of values $y_{n+m+j }\left(0\le j\le l\right)$. For this aim, one can use the predictor–corrector methods. Let us note that if method (11) is implicit, then it will take place that $l=0$. In this case, one can also use predictor–corrector methods (see [29]). For the sake of objectivity, let us note that for using method (11), one can use the same predictor–corrector method in both the cases $l=0$ and $l\ne 0$. Therefore, in the using of method (11), additional difficulties do not arise for the case $l>0$.

Let us note that method (10) can be applied to the solving of Equation (2), from the results of which one can receive the following:

(12)${\displaystyle \sum _{i=0}^m}{\alpha }_i\left(y_{n+i}-f_{n+i}\right)=h{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}K\left(x_{n+j},x_{n+i},y_{n+i}\right).$

It is clear that this method can be written as method (11), and in the application of them to solve some problems can be used in the above-described way. Therefore, let us consider the determination of the values of the coefficients ${\alpha }_i, {\beta }_i^{\left(j\right)} \left(i, j=0, 1, \dots , k\right),$ as the basic properties of the multistep methods depend on the values of their coefficients. For this aim let us suppose that by any methods we have found the values of the coefficients ${\alpha }_i, {\beta }_i^{\left(j\right)} \left(i, j=0, 1, \dots , k\right),$ by the choosing of which method (12) can have the degree of p. The conception of degree can be defined by the following way:

Definition 1.The integer p is called the degree for method (12), if the following holds:

(13)${\displaystyle \sum _{i=0}^m}{\alpha }_i\left(y\left(x_{n+i}\right)-f\left(x_{n+i}\right)\right)=h{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}K\left(x_{n+j},x_{n+i},y(x_{n+i})\right)+O\left(h^{p+1}\right),h\to 0.$

where $y\left(x_m\right)$ is the exact value of the solution of the problem (2) at the point $x_m \left(m\ge 0\right)$.

By the above-described way we have constructed methods (8) and (12) to solve the Volterra integral equation of the second kind, presented by the equation of (2). It is known that both theoretical and practical interests are stable methods with a high order of accuracy. Therefore, let us define the maximum values of the degree for the methods (8) and (12). For this, let us consider the following theorem:

Theorem 1.

If the methods (8) and (12) have the degree of p and p₁, respectively, then 1 ≤ p ≤ 2k and 1 ≤ p₁ ≤ k + m.

If methods (8) and (12) are stable, then the following holds:

$1\le p\le 2\left[k/2\right]+2,1\le p_1\le k+l+1(\mathrm{if}m=k-l\mathrm{and}k\ge 3l).$

Proof. 

It is obvious that the theorem also holds in the case $K\left(x, s, y\right)=\phi \left(s, y\right).$ In this case, methods (8) and (12) will match with method (9) (and in the case $m=k$). By Dahlquist’s rule we find that $p_1\le k+m$, and the method with the degree $p_{max}=m+k$ (and also for the case $m=k$) is unique. If these methods are stable, then there are methods with the degree:

$p\le 2\left[k/2\right]+2,p_1\le k+l+1,$

for all the values of $K$. □

Generally speaking, there is no uniqueness for the methods of type (8) and (12) from the corresponding conditions (see [20,23]).

From here we find that the local truncation error for this method can be presented as $O\left(h^{p+1}\right)$. It is not difficult to understand that the equality of (13) will also hold in the case when $K\left(x, s, y\right)=\phi \left(s, y\right).$ In this case we find that the integral equation of (2) will be same with the equation of (3) (see, for example, [1,30,31,32,33,34]). By the above-described way we have proved that the solution of the integral equation of (3) coincides with the solution of the initial-value problem for ODEs of the first order, which have been written as the problem of (4). It follows that to the solving of the Equation (3) can been applied the methods constructed for solving the initial-value problem for ODEs. Taking into account that in the method of (12) one can replace the function of $K\left(x, s, y\right)$ with the function of $\phi \left(s, y\right)$, then we find that in this case from method (12) one can receive the following:

(14)${\displaystyle \sum _{i=0}^m}{\alpha }_i{y}_{n+i}={\displaystyle \sum _{i=0}^m}{\alpha }_i{f}_{n+i}+h{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}\phi \left(x_{n+i},y_{n+i}\right).$

If in the method of (14) we use the next replacement:

(15)${\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}={\beta }_i\left(i=0,1,\dots ,k\right),$

then the receiving method will be same as method (10).

If we assume that the coefficients ${\beta }_i \left(0\le i\le k\right)$ are known, then it follows that (15) is the system of linear algebraic equations. Note that the solution of this system is not unique. Generally speaking, finding the solution of system (15) is not difficult. As seen from here, the value of the degree $p$ is independent from the coefficients of ${\beta }_i^{\left(j\right)}\left(i, j=0, 1, \dots , k\right).$ Let us prove that the value of $p$ depends on the values of the coefficients of ${\alpha }_i, {\beta }_i \left(i=0, 1, \dots , k\right).$ To illustrate this, here we propose to use the following Taylor series:

(16)$y\left(x+ih\right)=y\left(x\right)+ih{y}^{\prime }\left(x\right)+\frac{{\left(ih\right)}^2}{2!}{y}^{″}\left(x\right)+\cdots +\frac{{\left(ih\right)}^p}{p!}{y}^p\left(x\right)+O\left(h^{p+1}\right),$

(17)$y^{\prime }\left(x+ih\right)=y^{\prime }\left(x\right)+ih{y}^{″}\left(x\right)+\frac{{\left(ih\right)}^2}{2!}{y}^{‴}\left(x\right)+\cdots +\frac{{\left(ih\right)}^{p-1}}{\left(p-1\right)!}{y}^p\left(x\right)+O\left(h^p\right).$

By taking into account these series in (14), we receive the following:

(18)$\begin{array}{l}{\displaystyle \sum _{i=0}^m}{\alpha }_i(y(x+ih)-f(x+ih))-h{\displaystyle \sum _{i=0}^k}{\beta }_i\phi \left(x+ih\right)={\displaystyle \sum _{i=0}^m}{\alpha }_i((y(x)-f(x)+h({\displaystyle \sum _{i=0}^m}i{\alpha }_i-{\displaystyle \sum _{i=0}^k}{\beta }_i){(y(x)-f(x))}^{\prime }\\ +h^2({\displaystyle \sum _{i=0}^m}\frac{i^2}{2!}{\alpha }_i-{\displaystyle \sum _{i=0}^k}{\beta }_i){(y(x)-f(x))}^{″}+\cdots +h({\displaystyle \sum _{i=0}^m}\frac{i^2}{2!}{\alpha }_i-{\displaystyle \sum _{i=0}^k}\frac{i^{p-1}}{\left(p-1\right)!}{\beta }_i){(y(x)-f(x))}^p\\ +O(h^{p+1})=0\end{array}$

where $x=x_0+nh$ is a fixed point and ${\left(y\left(x\right)-f\left(x\right)\right)}^{\left(j\right)}={\phi }^j\left(x, y\right)\left(j=0, 1, \dots , p\right)$ (it follows from here the equality which is similar to Equation (4)).

Suppose that the following equalities hold:

(19)${\displaystyle \sum _{i=0}^m}{\alpha }_i=0;{\displaystyle \sum _{i=0}^k}{\beta }_i={\displaystyle \sum _{i=0}^m}i{\alpha }_i;{\displaystyle \sum _{i=0}^k}i{\beta }_i={\displaystyle \sum _{i=0}^m}\frac{i^2}{2!}{\alpha }_i,\dots ,{\displaystyle \sum _{i=0}^k}\frac{i^{p-1}}{\left(p-1\right)!}{\beta }_i={\displaystyle \sum _{i=0}^m}\frac{i^p}{p!}{\alpha }_i.$

Then, from (18), we receive the following:

(20)${\displaystyle \sum _{i=0}^m}{\alpha }_i\left(y\left(x+ih\right)-f\left(x+ih\right)\right)-h{\displaystyle \sum _{i=0}^k}{\beta }_i\phi \left(x+ih\right)=O(h^{p+1}),h\to 0,$

where $x=x_0+nh$ is a fixed point.

In this case we find that the method of (14) has the degree of p. Now, let us prove that if the asymptotic equality of (20) holds, then the system of algebraic Equation (19) will have a solution. It is not hard to understand that, if the asymptotic equality of (20) holds, then we find that the following also holds:

(21) $\begin{array}{l}{\displaystyle \sum _{i=0}^m}{\alpha }_i(y(x)-f(x))+h({\displaystyle \sum _{i=0}^m}i{\alpha }_i-{\displaystyle \sum _{i=0}^k}{\beta }_i)(y^{\prime }(x)-f^{\prime }(x))+h^2({\displaystyle \sum _{i=0}^m}\frac{i^2}{2!}{\alpha }_i-{\displaystyle \sum _{i=0}^k}i{\beta }_i)(y^{″}(x)-f^{″}(x))+\cdots \\ +h^p({\displaystyle \sum _{i=0}^m}\frac{i^p}{p!}{\alpha }_i-{\displaystyle \sum _{i=0}^k}\frac{i^{p-1}}{\left(p-1\right)!}{\beta }_i)(y^{(p)}(x)-f^{(p)}(x))=0\end{array}$

Let us consider the following notation:

$z\left(x\right)=y\left(x\right)-f\left(x\right).$

It is known that if $z\left(x\right)$ is a sufficiently smooth function, then $z\left(x\right), z^{\prime }\left(x\right), \dots , z^p\left(x\right)$ is the independent linear system, if $z^{\left(j\right)}\not\equiv 0 \left(0\le j\le p\right).$ If we take this into account in the equality of (21), then from that it follows the system of (19). It follows from here that the fulfillment of the condition (19) for the coefficients of method (14) is necessary, and is a sufficient condition for the holding of the asymptotic equality of (20). Thus, we have proved the following lemma:

Lemma 1.

In order for the method of (14) to have a degree p, the satisfaction of its coefficients by the system of algebraic Equation (19) is necessary and sufficient.

In the system of (19) there are $k+m+2$ unknowns and $p+1$ equations. Equation (19) is a system of linear algebraic equations and in the case $p+1=k+m+2$, the determinant of this system is nonzero (in this case receiving the Vandermond determinant). As was noted above, the condition ${\alpha }_m\ne 0$ must hold. By taking into account this condition, we receive $p+1=k+m+2$. It follows that $p\le k+m.$ In the case $m=k$, from here we receive Dahlquist’s rule, which can be written as $p\le 2k \left(p_{max}=2k or p_{max}=k+m\right)$. One can prove that the method with the degree $p_{max}=k+m$ is unique (see, for example, [33,34,35,36]).

Comment 1.

In the above-described method for the comparison of the advanced and multistep methods, we usually use the values of the variables m and k. Note that the advanced methods can not be received from the multistep methods. For the proving of this, let us consider the following k-step methods:

(22) ${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}=h{\displaystyle \sum _{i=0}^k}{\beta }_i{y}_{n+i}^{\prime }.$

We usually suppose that ${\alpha }_k\ne 0,$ which has a relation with finding the value $y_{n+k}$ as the solution of the finite-difference equation of (22). From here, we find that in the case ${\alpha }_k=0$, the equality (22) is transferable to the other class method. As was noted, if method (22) is stable and ${\alpha }_k\ne 0$, then $p\le 2\left[k/2\right]+2$, and there are stable methods with the degree $p_{max}=2\left[k/2\right]+2$ for all the values of the order $k$.

Note that here we used the following definition for the stability:

Definition 2.

Method (9) is called stable if the roots of the following polynomial

$\rho \left(\lambda \right)={\alpha }_m{\lambda }^m+{\alpha }_{m-1}{\lambda }^{m-1}+\cdots +{\alpha }_1\lambda +{\alpha }_0$

lie inside the unit circle, on the boundary of which there are no multiply roots.

Comment 2.

As was noted above for method (9), the condition p ≤ k + m holds. If m = k − l (l > 0 ) then we receive p ≤ 2k − l. However, if method (9) is stable, then p ≤ k + l + 1. It is not hard to understand that the linear parts of the methods, which are investigated here, are the same. Therefore, the conception of stability for them is defined in one and the same way. It follows from here that the methods (9) and (22) are independent from each other, because these methods are the independent objects of research.

For the construction of methods with a high order of accuracy or higher degrees, the hybrid method is often used (see, for example, [37,38,39,40,41,42,43]). Therefore, let us consider the following paragraph.

3. Construction of a Generalized Hybrid Method and Its Application

Let us remember some of the popular methods, which have been applied to solve the problem (4). Among of them are the Euler methods (explicit and implicit) and trapezoidal and midpoint rules. The midpoint rule differs from others in that this method uses the calculation of variables of the type $y\left(x_n+h/2\right).$ This variable can be written in a more general form as $y\left(x_n+{\nu }_{i}h\right), \left(\left|{\nu }_i\right|<1, i=0, 1, 2,\dots , k\right).$ By the generalization of the midpoint rule, one can construct the following hybrid method:

(23)${\displaystyle \sum _{i=0}^k}{\alpha }_i\left(y_{n+i}-f_{n+i}\right)=h{\displaystyle \sum _{i=0}^k}{\beta }_i{\phi }_{n+i+{\nu }_i}\left(\left|{\nu }_i\right|<1,i=0,1,2,\dots ,k\right).$

In the work of [42], they constructed a hybrid method with the degree $p_{max}=4$, which can be received from method (23) in the case $k=1$. However, from the multistep method (22) for the case $k=1$, one can receive the method with the degree $p_{max}=2$. By simple comparison, we find that the hybrid methods can be taken as the perspective. From Equation (23) one can receive the midpoint rule, which can be taken as the explicit method. It is known that this method has the degree $p=2.$ However, as was noted above from the method of (22), one can receive the explicit method with the degree $p_{max}=1$ for $k=1.$ This comparison shows that the hybrid methods have some advantages over all the known methods. Note that some hybrid methods have the extended region of stability and all the methods investigated here are linear multistep methods; therefore, the linear part of these methods has the same properties. It follows that the conception of stability and degree can be defined in the same way for all linear methods. Let us note that the values of the coefficients ${\alpha }_i \left(i=0, 1, \dots , k\right)$ can be different from the corresponding coefficients of the other methods.

Now, let us define the values of the coefficients ${\alpha }_i, {\beta }_i, {\nu }_i\left(i=0, 1, \dots , k\right)$. To this end, let us use the following Taylor series:

(24)$y^{\prime }\left(x+l_{i}h\right)=y^{\prime }\left(x\right)+l_{i}h{y}^{″}\left(x\right)+\frac{{\left(l_{i}h\right)}^2}{2!}{y}^{‴}\left(x\right)+\cdots +\frac{{\left(l_{i}h\right)}^{p-1}}{\left(p-1\right)!}{y}^p\left(x\right)+O(h^p)$

By taking into account Equations (16) and (24) in the method of Equation (23), we receive:

(25)$\begin{array}{l}{\displaystyle \sum _{i=0}^k}({\alpha }_i(y(x+ih)-f(x+ih))-h{\beta }_i\phi (x+(i+v_i)h,y(x+(i+v_i)h))\\ ={\displaystyle \sum _{i=0}^k}{\alpha }_i(y(x)-f(x))+h{\displaystyle \sum _{i=0}^k}(i{\alpha }_i-{\beta }_i)(y^{\prime }(x)-f^{\prime }(x))+h^2{\displaystyle \sum _{i=0}^k}(\frac{i^2}{2!}{\alpha }_i-l_i{\beta }_i)(y^{″}(x)\\ -f^{″}(x))+\cdots +h^p{\displaystyle \sum _{i=0}^k}(\frac{i^p}{p!}{\alpha }_i-\frac{l_i^{p-1}}{\left(p-1\right)!}{\beta }_i)(y^{(p)}(x)-f^{(p)}(x))+O(h^{p+1})=0\end{array}$

where $x=x_0+nh$ is a fixed point and $l_i=i+{\nu }_i\left(i=0,1,\dots ,k\right).$

By using the discussion, which we have used in the investigation of method (9), and taking into account the comparison of asymptotic equality (18) with asymptotic equality (25), we receive the following system for finding the determined values of the coefficients ${\alpha }_i, {\beta }_i, {\nu }_i\left(i=0, 1, \dots , k\right):$

(26)${\displaystyle \sum _{i=0}^k}{\alpha }_i=0;{\displaystyle \sum _{i=0}^k}{\beta }_i={\displaystyle \sum _{i=0}^k}i{\alpha }_i;{\displaystyle \sum _{i=0}^k}(i+{\nu }_i){\beta }_i={\displaystyle \sum _{i=0}^m}\frac{i^2}{2!}{\alpha }_i,\dots ,{\displaystyle \sum _{i=0}^k}\frac{{\left(i+{\nu }_i\right)}^{p-1}}{\left(p-1\right)!}{\beta }_i={\displaystyle \sum _{i=0}^m}\frac{i^p}{p!}{\alpha }_i.$

This is a nonlinear system of algebraic equations. In this system there are $3k+3$ unknowns, but the amount of equations in this system is equal to $p+1.$ Note that this is a nonlinear system of algebraic equations, because defining the exact solution of such a system is not easy. By taking this into account, scientists proposed to use some numerical methods for solving them. For this, they used Mathcard 2015. Note that by using the approximate solution of system (26), they have constructed some methods with the degree of $p$. The application of some of them to solving model problems has shown that in reality, the received results correspond to the results received by the methods with a degree less than $p$. Therefore, finding a private solution of system (26) is very important. Let us note that these results correspond to the theoretical.

The system of (26) to remember the nonlinear system of algebraic equations is used for finding the coefficients of the Gauss method. Therefore, some solutions of this system will be also solutions of the corresponding Gauss system which is used for finding Gauss nodes and coefficients (see, for example, [43,44,45,46,47,48,49]).

To solve system (26) is more simple than the corresponding Gauss system, and usually by the solution of (26) one can construct hybrid methods which are different from Gauss methods. By taking into account these properties, here we have proposed to construct stable methods with high degrees. For this, let us consider the following method:

(27)${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}={\displaystyle \sum _{i=0}^k}{\alpha }_i{f}_{n+i}+h{\displaystyle \sum _{i=0}^k}{\beta }_i{\phi }_{n+i}+h{\displaystyle \sum _{i=0}^k}{\gamma }_i{\phi }_{n+i+{\nu }_i}.$

One can consider this method as the linear combination of methods (10) and (23). It is easy to prove that the system of algebraic equations which has been constructed for finding the coefficients of method (27) can be constructed as the linear combination of systems (26) and (19), respectively. By the generalization of the system (26) or (19), we receive the following system of algebraic equations:

(28)$\begin{array}{l}{\displaystyle \sum _{i=0}^k}{\alpha }_i=0;{\displaystyle \sum _{i=0}^k}({\beta }_i+{\gamma }_i)={\displaystyle \sum _{i=0}^k}i{\alpha }_i;{\displaystyle \sum _{i=0}^k}(i{\beta }_i+l_i{\gamma }_i)={\displaystyle \sum _{i=0}^m}\frac{i^2}{2!}{\alpha }_i,\dots ,{\displaystyle \sum _{i=0}^k}\left(\frac{{\left(i\right)}^{p-1}}{\left(p-1\right)!}{\beta }_i+\frac{l_i^{p-1}}{\left(p-1\right)!}{\gamma }_i\right)\\ ={\displaystyle \sum _{i=0}^m}\frac{i^p}{p!}{\alpha }_i,\left(l_i=i+{\nu }_i,0\le i\le k\right).\end{array}$

By taking into account the series (16), (17), and (24) in the equality of (27), one can receive the system of (28), in the case when the method will have the degree of p. Now, let us investigate the solvability of the system (28).

As was noted above, the system (28) is nonlinear, therefore to find the exact solution of that is perhaps not always possible. Hence, let us consider an investigation of the system of (28). In this system participates $4k+4$ unknowns and $p+1$ nonlinear algebraic equations. Let us consider the case $p+1\le 4k+4$. Note that without breaking the generality one can take ${\alpha }_k=1$. By taking this into account, let us investigate the system of (28) for $p\le 4k+2$. It is clear that in the case of $k=1,$ there does not arise a question on the stability of the methods of type (27). Let us in the system of (28) take $k=1$. In this case, from system (28) we receive the following (by taking into account the condition ${\alpha }_1=1$, we receive ${\alpha }_0=-1)$:

(29)${\beta }_0+{\gamma }_0+{\beta }_1+{\gamma }_1=1;{\beta }_1+l_0^j{\gamma }_0+l_1^j{\gamma }_1=1/\left(j+1\right)\left(1\le j\le 5\right).$

By using the solution of system (29) one can construct the following method:

(30)$y_{n+1}=y_n+h\left(y_{n+1}^{\prime }+y_n^{\prime }\right)/12+5h\left(y_{n+\beta }^{\prime }+y_{n+1-\beta }^{\prime }\right)/12,\beta =\left(5-\sqrt{5}\right)/10.$

This method has the degree $p=6$. Similar investigations have been given by some authors (see, for example, [32,39,40,41,42,45]). For the application of this method to solve Equation (2), it can be presented in the following form:

(31)$\begin{array}{l}{y}_{n+1}-f_{n+1}=y_n-f_n+h\left(K\left(x_{n+1}, x_{n+1}, y_{n+1}\right)+K\left(x_{n+1}, x_n, y_n\right)+2K\left(x_n, x_n, y_n\right)\right)/24\\ +5h(K\left(x_{n+1}, x_{n+\beta }, y_{n+\beta }\right)+K\left(x_{n+\beta }, x_{n+\beta }, y_{n+\beta }\right)+K\left(x_{n+1}, x_{n+1-\beta }, y_{n+1-\beta }\right)\\ ++K\left(x_{n+1-\beta }, x_{n+1-\beta }, y_{n+1-\beta }\right))/24.\end{array}$

In the construction, method (30) has used the solution of system (29) and the solution of system (15) by the addition of the following:

(32)${\displaystyle \sum _{j=1}^k}{\gamma }_i^{\left(j\right)}={\gamma }_i\left(i=0,1,\dots ,k\right).$

By this way one can construct the following multistep hybrid method:

(33)$$\begin{gathered} {\displaystyle {\displaystyle \sum }_{i=0}^k}{\alpha }_i\left(y_{n+i}-f_{n+i}\right)=h{\displaystyle {\displaystyle \sum }_{i=0}^k}{\displaystyle {\displaystyle \sum }_{j=i}^k}{\beta }_i^{\left(j\right)}K\left(x_{n+j,}{x}_{n+i},y_{n+i}\right)+ \\ +h{\displaystyle {\displaystyle \sum }_{i=0}^k}{\displaystyle {\displaystyle \sum }_{j=i}^k}{\gamma }_i^{\left(j\right)}K\left(x_{n+i+{\nu }_j,}{x}_{n+i+{\nu }_i},y_{n+i+{\nu }_i}\right),\left(\left|{\nu }_i\right|<1;i=0,1,\dots ,k\right) \end{gathered}$$

by taking into account the solutions of the systems (15) and (32).

Hence, note that method (31) can be received from the method of (33) as the partial case. The solution of systems (15) and (32) is not unique, and the order of accuracy of method (33) is independent from the solution of mentioned systems. As was proved above, the exactness of the methods of type (33) depends on the values of the coefficients ${\beta }_i, {\gamma }_i$, and ${\nu }_i \left(i=0, 1, \dots , k\right),$ which can be found as the solution of system (28).

For the construction of methods of type (33), with high accuracy, let us consider the case of $k=2$.

At first, let us define the values of ${\alpha }_i\left(i=0, 1, 2\right)$ by taking into account that the constructed method must be stable. It is clear that in this case the roots of the following polynomial:

${\lambda }^2+{\alpha }_1\lambda +{\alpha }_0=0$

must satisfy the condition of stability. Here, we have considered the following variant:

${\alpha }_1=0,{\alpha }_0=-1,{\alpha }_2=1.$

In this case, by using the solution of nonlinear system (28) and taking into account the solution of systems (15) and (32), one can construct the following method:

(34)$y_{i+2}=y_i+h\left(9y_{i+2}^{\prime }+64y_{i+1}^{\prime }+9y_i^{\prime }\right)/90+49h\left(y_{n+1+\alpha }^{\prime }+y_{n+1-\alpha }^{\prime }\right)/90,\alpha =\sqrt{21}/7.$

For the application of this method to solve nonlinear Volterra integral equations of second order, that can be modified as following:

(35)$\begin{array}{l}{y}_{i+2}=y_i+f_{i+2}-f_i+h(9K\left(x_{i+2}, x_{i+2}, y_{i+2}\right)+32K\left(x_{i+2}, x_{i+1}, y_{i+1}\right)+32K\left(x_{i+1}, x_{i+1}, y_{i+1}\right)+5K\left(x_{i+1}, x_i, y_i\right)\\ +4K\left(x_i, x_i, y_i\right))/90+49h(K\left(x_{n+2}, x_{n+1+\alpha }, y_{n+1+\alpha }\right)\\ +K\left(x_{n+1+\alpha }, x_{n+1+\alpha }, y_{n+1+\alpha }\right))+K\left(x_{n+1}, x_{n+1-\alpha }, y_{n+1-\alpha }\right)+K\left(x_{n+1-\alpha }, x_{n+1-\alpha }, y_{n+1-\alpha }\right))/180\end{array}$

It is not difficult to prove that the method of (34) can be presented in another form, which will be different from the formula (35). As was noted above, for the construction of more exact methods one can use advanced (forward-jumping) methods. However, some specialists, for the construction of more exact methods, proposed using the multistep second derivative methods with constant coefficients. For receiving some information about these methods, one can use the content of the following section.

4. On Some Properties of Multistep Second Derivative Methods with Constant Coefficients

In the last sections, we have given some information about advanced methods which have comparisons with multistep methods. Note that multistep second derivative methods can also be the advanced type. By taking into account this property, let us consider the following multistep second derivative methods, which are fundamentally investigated by some authors (see, for example, [46,47,48]):

(36)${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}=h{\displaystyle \sum _{i=0}^k}\overline{{\beta }_i}{y}_{n+i}^{\prime }+h^2{\displaystyle \sum _{i=0}^k}\overline{{\gamma }_i}{y}_{n+i}^{″}.$

This method, after application to the solving of problem (2), can be presented as the following:

(37)${\displaystyle \sum _{i=0}^k}{\alpha }_i(y_{n+i}-f_{n+i})=h{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\overline{\beta }}_i^{\left(j\right)}K\left(x_{n+j},x_{n+i},y_{n+i}\right)+h^2{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\overline{\gamma }}_i^{\left(j\right)}G\left(x_{n+j},x_{n+i},y_{n+i}\right)$

where the function of $G\left(x, x, y\right)$ is defined as: $G\left(x, x, y\right)=\frac{d}{dx}K\left(x, s, y\left(s\right)\right){|}_{s=x}.$ Depending on the used way to construct method (37), the function $G\left(x, z, y\right)$ can be defined in another form, but the received results of which will be the same with the method of (37). It follows to note that if in method (37) we input $K\left(x,s, y\right)=\phi \left(s, y\right),$ then method (36) can be received from method (37) as the partial case.

To explain the above description, it is enough to apply methods (36) and (37) to solve the following problem:

$y^{\prime }=\phi \left(x,y\right),y\left(x_0\right)=y_0,x_0\le x\le X.$

In this case we receive:

(38)${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}=h{\displaystyle \sum _{i=0}^k}\overline{{\beta }_i}{\phi }_{n+i}+h^2{\displaystyle \sum _{i=0}^k}\overline{{\gamma }_i}{g}_{n+i},$

(39)${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}=h{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\overline{\beta }}_i^{\left(j\right)}{\phi }_{n+i}+h^2{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\overline{\gamma }}_i^{\left(j\right)}{g}_{n+i},$

where the function $g\left(x,y\right)$ is defined as $g\left(x,y\right)={\phi }_x^{\prime }\left(x,y\right)+{\phi }_y^{\prime }\left(x,y\right)\phi \left(x,y\right).$

It is evident that if we take

(40)${\displaystyle \sum _{j=i}^k}{\overline{\beta }}_i^{\left(j\right)}=\overline{{\beta }_i};{\displaystyle \sum _{j=i}^k}{\overline{\gamma }}_i^{\left(j\right)}=\overline{{\gamma }_i}\left(i=0,1,2,\dots ,k\right),$

then from the formula (39) follows method (38). Therefore, to define the values of the coefficients ${\overline{\beta }}_i^{\left(j\right)}, {\overline{\gamma }}_i^{\left(j\right)} \left(i, j=0, 1, 2, \dots , k\right),$ one can use the system (40), and the system of equations, which are independent from the determination of the coefficients $\overline{{\beta }_i}, \overline{{\gamma }_i} \left(i=0, 1, 2, \dots , k\right),$ participate in the formula of (36). By taking this into account, let us consider the definitions of the values of the coefficients ${\alpha }_i, \overline{{\beta }_i}, \overline{{\gamma }_i} \left(i=0, 1, 2, \dots , k\right).$ To this end, one can use the scheme which was used in the construction of system (28). In this case, the system of algebraic equations for finding the values of the coefficients ${\alpha }_i, \overline{{\beta }_i}, \overline{{\gamma }_i}$ in one variant can be written in the following form:

(41)${\displaystyle \sum _{i=0}^k}{\alpha }_i=0;{\displaystyle \sum _{i=0}^k}\overline{{\beta }_i}=i{\alpha }_i;{\displaystyle \sum _{i=0}^k}i\overline{{\beta }_i}+{\displaystyle \sum _{i=0}^k}\overline{{\gamma }_i}={\displaystyle \sum _{i=0}^m}\frac{i^2}{2!}{\alpha }_i,\dots ,{\displaystyle \sum _{i=0}^k}\frac{i^l}{l!}\overline{{\beta }_i}+{\displaystyle \sum _{i=0}^k}\frac{i^{l-1}}{\left(l-1\right)!}\overline{{\gamma }_i}={\displaystyle \sum _{i=0}^m}\frac{i^{l+1}}{\left(l+1\right)!}{\alpha }_i,\left(l=2,3,\dots ,p\right).$

Note that the system of (28) is nonlinear, but system (41) is linear. By taking into account that the determinant of system (41) is nonzero, we find that if the amount of the unknowns and of the equations are the same, then we find that system (41) has a unique solution. Note that the amount of unknowns in system (41) is equal to $3k+3,$ but the amount of equations is equal to $p+1.$ It follows that if $p<3k+1$ then the system of (41) will have any solution, but in the case $p=3k+1,$ the corresponding solution of system (41) will be unique. It follows from here that $p_{max}=3k+1$. Note that the degree $p$ for method (36) can be defined as follows (see, for example, [46,47,48]):

Definition 3.

The integer p is called the degree for method (36) if the following holds:

(42) ${\displaystyle \sum _{i=0}^k}({\alpha }_{i}y\left(x+ih\right)-h\overline{{\beta }_i}{y}^{\prime }\left(x+ih\right)-h^2\overline{{\gamma }_i}{y}^{″}\left(x+ih\right))=O(h^{p+1}),h\to 0.$

Let us note that this and definition 1 are the same.

It is not difficult to prove that if method (36) is stable then the relationship between $k$ and $p$ (order and degree) can be presented as $p\le 2k+2$, and there exists stable methods with the degree $p=2k+2$ for all the values of $k$ (order). For the value $k=1$ we find that the maximal value for stable and unstable methods are the same; in other words, $3k+1=2k+2$ for $k=1$. Note that in this case ($k=1)$ there are not any unstable methods, so the one-step method satisfies the condition of stability. By simple comparison we find that the stable methods of type (27) are more exact than the stable methods of type (39), but application of hybrid methods of type (23) or (27) is more difficult. Note that these difficulties are related to the calculation of the values $y_{n+i+{\nu }_i} \left(i=0, 1, \dots , k\right).$ It is evident that for the calculation of these values arises the necessity to construct a special method for calculating them. Therefore, to give some advantages of these methods are difficult.

If there exist methods for calculations of the values $y_{n+i+{\nu }_i}\left(i=0, 1, \dots , k\right),$ then the methods of hybrid types will have some advantages. It is not difficult to prove that the hybrid methods constructed by using the methods of type (38) will be more exact than the known methods, but will have a more complex structure. As was noted above, the system (41) is linear; therefore, the system (41) can be solved by using the known methods. However, in the increasing the values of order k the values of the calculation in the determination of the values of the coefficients ${\alpha }_i, \overline{{\beta }_i}, \overline{{\gamma }_i} \left(i=0, 1, 2, \dots , k\right)$ also increase. For decreasing the values of the calculation works, here we have proposed a new way for the calculation of the values of the coefficients ${\alpha }_i, \overline{{\beta }_i}, \overline{{\gamma }_i} \left(i=0, 1, 2, \dots , k\right)$ (see [49]). To this end, we required analyticity from the solution of the considered problem. In the work of [48], they proposed a way by which the condition of analyticity of the solution could be simplified and replaced with the conditions that usually are used in other works (see, for example, [48]). For the method to have a degree of p one can use the following way, which is received by using the above-mentioned results:

${\alpha }_0=-{\rho }_0+{\rho }_1-{\rho }_2+\cdots +{\left(-1\right)}^{k-1}{\rho }_{k-2}+{\left(-1\right)}^k{\rho }_{k-1},$

${\alpha }_i={\displaystyle \sum _{j=i-1}^{k-1}}{\left(-1\right)}^{j-i+1}\left(j+1\right)j\left(j-1\right)\dots \left(j-i+2\right){\rho }_j/i!;i=1,2,\dots ,k.$

$\overline{{\beta }_0}={\delta }_0-{\delta }_1+{\delta }_2+\cdots +{\left(-1\right)}^{k-1}{\delta }_{k-1}+{\left(-1\right)}^k{\delta }_k,$

(43)$\overline{{\beta }_i}={\displaystyle \sum _{j=1}^k}{\left(-1\right)}^{j-1}j\left(j-1\right)\dots \left(j-i+1\right){\delta }_j/i!;i=1,2,\dots ,k.$

$\overline{{\gamma }_0}=l_0+l_1-l_2+\cdots +{\left(-1\right)}^{k-1}{l}_{k-1}+{\left(-1\right)}^k{l}_k,$

$\overline{{\gamma }_i}={\displaystyle \sum _{j=1}^k}{\left(-1\right)}^{j-1}j\left(j-1\right)\dots \left(j-i+1\right)l_j/i!;i=1,2,\dots ,k.$

The variables ${\rho }_i, {\delta }_i, l_i \left(i=0, 1, 2, \dots , k\right)$ can be defined from the following system of linear algebraic equations.

(44)${\displaystyle \sum _{i=0}^j}{c}_i{\rho }_{j-i}+{\displaystyle \sum _{i=1}^j}{\left(-1\right)}^{j-i+1}{l}_{i-1}/\left(j-i+1\right)={\delta }_j,j=0,1,\dots ,k;{\rho }_k=0,$

${\displaystyle \sum _{i=j+1}^{j+k}}{c}_i{\rho }_{j+k-i}+{\displaystyle \sum _{i=j}^j}{\left(-1\right)}^i{l}_{i+k-1}/i=0,j=0,1,\dots ,k;$

${\displaystyle \sum _{i=j+k-1}^{j+2k}}{c}_i{\rho }_{j+2k-i}+{\displaystyle \sum _{i=j+k}^{j+2k}}{\left(-1\right)}^i{l}_{j+2k-i}/i=0;j=1,2,\dots ,k;$

${\displaystyle \sum _{i=2k+2}^{3k+1}}{c}_i{\rho }_{3k+1-i}+{\displaystyle \sum _{i=2k+1}^{3k+1}}{\left(-1\right)}^i{l}_{3k+1-i}=C,$

where $C$ is the constant for the coefficient of the main leading term in the expansion of the error of method (39), but the coefficients $c_i \left(i=0, 1, \dots \right)$ are defined by the following formula:

$c_i=\frac{1}{i!}{\displaystyle \underset{0}{\overset{1}{\int }}}u\left(u-1\right)\dots \left(u-i+1\right)du;\left(i=1,2,\dots \right),c_0=1,c_1=1/26,c_2=-\frac{1}{12},\dots$

It is easy to prove that the coefficients can be calculated by the following formula:

$c_m={\displaystyle \sum _{i=1}^m}{\left(-1\right)}^{i-1}{c}_{m-i}/\left(i+1\right)(m\ge 1,c_0=1).$

Note that someone may think that finding a solution to systems (43) and (44) is more difficult than finding the solution to system (41). However, it is not. As is known, scientists have mostly constructed stable methods, considering that they are convergent. By taking this into account, we find that, basically, one can assume that the values of the quantities ${\rho }_i \left(i=0, 1, \dots , k-1\right)$ are known. In this case, we find that to solve the system of (44) is simplified. As a result of which, we obtain the actual solution of one system, which is system (43). By using this solution in the system of (42) we can compute the values of the unknowns ${\delta }_j \left(j=0, 1, \dots , k\right).$ By using these values one can find the values of the coefficients ${\alpha }_i, \overline{{\beta }_i}, \overline{{\gamma }_i} \left(i=0, 1, 2, \dots , k\right).$ Taking into account the solution of the systems (42) and (43), here we have constructed stable methods with the degree $p=8$ for $k=3$. The constructed stable method with the degree $p=8$ can be written as the following:

(45)$\begin{array}{l}{y}_{n+3}=\left(y_{n+2}+y_{n+1}+y_n\right)/3+h\left(\mathrm{10,781}{y}_{n+3}^{\prime }+\mathrm{22,707}{y}_{n+2}^{\prime }+\mathrm{16,659}{y}_{n+1}^{\prime }+4285y_n^{\prime }\right)/\mathrm{27,216}\\ -h^2\left(2099y_{n+3}^{″}-7227y_{n+2}^{″}-2853y_{n+1}^{″}-979y_n^{″}\right)/\mathrm{45,360}+3h^9{y}_n^{\left(9\right)}/\mathrm{156,800}+O\left(h^{10}\right).\end{array}$

Let us consider a comparison of the hybrid method of type (34) with method (45). These methods have the degree $p=8$ and are stable. Note that method (34) has the order of k = 2, but method of (45) has the order k = 3. It follows that for using method (45) we must know the values (y_n, y_n+1 and y_n+2) but for using method (34) we must know two values (y_n and y_n+1). For the application of method (45) it is necessary to use one explicit method as the predictor formula. However, for the application of method (34) to solve some problems, it needs to use two methods for the calculation of the values of the type $y\left(x_m\pm \nu h\right)$ ($\left|\nu \right|<1$). If in the method of (36) or (38) we input k = 2 then the degree for the stable methods will hold the condition of p ≤ 6. It follows that to give some advantages of any of these methods is difficult. Each of them has its own advantages and disadvantages. Now let, us apply method (45) to solve a Volterra integral equation. In this case we receive the following:

(46)$\begin{array}{l}{y}_{n+3}=\left(y_{n+2}+y_{n+1}+y_n\right)/3+f_{n+3}-\left(f_{n+2}+f_{n+1}+f_n\right)/3+h(\mathrm{10,781}K(x_{n+3},x_{n+3}, y_{n+3})\\ +\mathrm{11,707}K\left(x_{n+3}, x_{n+2}, y_{n+2}\right)+\mathrm{11,000}K\left(x_{n+2}, x_{n+2}, y_{n+2}\right)+8659K\left(x_{n+2}, x_{n+1}, y_{n+1}\right)\\ +8000K\left(x_{n+1}, x_{n+1}, y_{n+1}\right)+2185K\left(x_{n+2}, x_n, y_n\right)+2100K\left(x_n, x_n, y_n\right)/\mathrm{27,216}\\ -h^2(2099G\left(x_{n+3}, x_{n+3}, y_{n+3}\right)-7000G\left(x_{n+3}, x_{n+2}, y_{n+2}\right)-227G\left(x_{n+2}, x_{n+2}, y_{n+2}\right)\\ -1453G\left(x_{n+2},x_{n+1}, y_{n+1}\right)-1400G\left(x_{n+1},x_{n+1}, y_{n+1}\right)-9006\left(x_{n+1},x_n, y_n\right)\\ -79G\left(x_n,x_n, y_n\right))/\mathrm{45,360}\end{array}$

It is not easy to determine the value of $y_{n+3}$ by method (46), so, as in this case, we receive the nonlinear algebraic equation. For solving this equation, here we propose to use the predictor–corrector methods (see, for example, [29]). To this end, one can use the stable explicit methods as the predictor methods. In this case, the degree for this method will satisfy the condition $p\le 6,$ but method (46) has the degree $p=8.$ It follows that for the construction of methods with suitable accuracy, the condition $k\ge 4$ must be held. However, one can use the suitable stable explicit method (46) as the predictor and corrector method.

It is obvious that someone can propose a way to increase the accuracy of the calculated values $y_{n+3},$ which differ from the above description. Note that similar difficulties arise in the application of the quadrature method to solve the nonlinear Volterra integral equations.

Thus, we have shown that by using the properties of the investigated problem one can choose a suitable method. Lately, the specialists have predominantly used hybrid methods. Here, we have a comparison of the same numerical methods by using the conception stability, degree, and the volume of the computational works on each step. However, some authors, for the comparison of numerical methods, use the conception of the region of stability. This question can be solved by taking into account the results received when using predictor–corrector methods. Usually, the values found by the predictor–corrector methods are used for the definition of the boundaries for the step size $h>0$.

As is known, for the construction of the multistep methods are usually given the amount of mesh points. Therefore, to find some relation between the order $k$ and the degree $p$ for the investigated multistep methods is very important. It is known that one of the important questions in the investigation of multistep methods is defining the necessary conditions for their convergence. These conditions for method (36) and the methods which are received from that as partial cases, for example, method (22), have been investigated by Dahlquist (see, for example, [27]). By taking this into account, in the next section let us define the necessary condition for the convergence of method (33), which was received by the application of method (27) to solve Volterra integral equations.

5. The Conditions Imposed on Coefficients of Method (33)

Let us suppose that method (33) is convergence and prove that the following conditions are satisfied.

A. The coefficients ${\alpha }_i,{\beta }_i,{\gamma }_i,{\nu }_i\left(i=0,1,\dots ,k\right)$ are real numbers and ${\alpha }_k\ne 0$;

B. The characteristic polynomials:

$\rho \left(\lambda \right)\equiv {\displaystyle \sum _{i=0}^k}{\alpha }_i{\lambda }^i; \delta \left(\lambda \right)\equiv {\displaystyle \sum _{i=0}^k}{\beta }_i{\lambda }^i; \gamma \left(\lambda \right)\equiv {\displaystyle \sum _{i=0}^k}{\gamma }_i{\lambda }^i;$

have no common factor different from the constant;

C. The conditions $P\ge 1$ and $\delta \left(1\right)+\gamma \left(1\right)\ne 0$ hold.

The necessity of the condition ${\alpha }_k\ne 0$ is proved above. Therefore, condition A is obvious. Let us consider condition B and suppose otherwise. It follows that the polynomials have a common factor, which differs from the constant. Denote that by $\phi \left(\lambda \right)$. Then, one can write $\phi \left(\lambda \right)\not\equiv const$, and by the E, denote the shift operator. In this case, the following holds:

$E^{i}y\left(x\right)=y\left(x+ih\right)$ or ${\displaystyle {\sum }_{i=0}^k}{\alpha }_i{y}_{n+i}={\displaystyle {\sum }_{i=0}^k}{\alpha }_{i}y\left(x+ih\right)={\displaystyle {\sum }_{i=0}^k}{\alpha }_i{E}^{i}y\left(x\right),$ here $x=x_0+nh$ is a fixed point.

It is easy to see that one can write the following:

(47)$E^j{E}^{i}K\left(x_n, x_n, y_n\right)=K\left(x_{n+j}, x_{n+i}, y_{n+i}\right)=K\left(x_{n+i}+\left(j-i\right)h, x_{n+i}, y_{n+i}\right)$

If for the fixed point $x=x_n+ih$, passing the limit to respect the first argument, then we receive:

$\underset{h\to 0}{\lim }K\left(x+\left(j-i\right)h, x, y\left(x\right)\right)=K\left(x, x, y\left(x\right)\right).$

By using this in the equality of (47) we receive:

$\underset{h\to 0}{\lim }{E}^j{E}^{i}K(x_n, x_n, y_n)=E^{i}K\left(x_n, x_n, y_n\right).$

Hence, $x=x_n+ih$ is fixed.

If we use these properties in the following expression:

(48)$\underset{h\to 0}{\lim }{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}K(x_{n+j}, x_{n+i}, y_{n+i})={\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}K(x_{n+i}, x_{n+i}, y_{n+i}).$

By taking into account the systems of (15) and (32) in the equality of (48), we receive:

${\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}K(x_{n+j}, x_{n+i}, y_{n+i})={\displaystyle \sum _{i=0}^k}{\beta }_{i}K(x_{n+i}, x_{n+i}, y_{n+i})+O\left(h\right).$

Thus, we find that method (33) can be written as:

(49)${\displaystyle \sum _{i=0}^k}{\alpha }_i\left(y_{n+i}-f_{n+i}\right)-h{\displaystyle \sum _{i=0}^k}{\beta }_{i}K(x_{n+i}, x_{n+i}, y_{n+i})-h{\displaystyle \sum _{i=0}^k}{\gamma }_{i}K(x_{n+i+{\nu }_i}, x_{n+i+{\nu }_i}, y_{n+i+{\nu }_i})=O\left(h^r\right), \left(r>1\right)$

From here we receive:

(50)$\rho \left(E\right)\left(y_n-f_n\right)-h\delta \left(E\right)K\left(x_n, x_n, y_n\right)-h\gamma \left(E\right)K\left(x_n, x_n, y_n\right)=0.$

By the above-described way, we prove that the finite-difference Equations (33) and (50) are equivalents. Note that Equation (49) is homogeneous, therefore the solvability theorem for the finite-difference Equations (33) and (50) are the same. By this assumption, we find that the polynomials $\rho \left(\lambda \right), \delta \left(\lambda \right)$, and $\gamma \left(\lambda \right)$ have the common factor, which is denoted by $\phi \left(\lambda \right)$. If we use the function $\phi \left(\lambda \right)$ in the equality of (50), then we receive:

(51)$\phi \left(E\right)\left({\rho }_1\left(E\right)\left(y_n-f_n\right)-h{\delta }_1\left(E\right)K\left(x_n, x_n, y_n\right)-h{\gamma }_1\left(E\right)K\left(x_n, x_n, y_n\right)\right)=0.$

By using the condition $\phi \left(\lambda \right)\ne const,$ in the equality of (51) we receive the following:

(52)${\rho }_1\left(E\right)\left( y_n-f_n\right)-h{\delta }_1\left(E\right)K\left(x_n, x_n, y_n\right)-h{\gamma }_1\left(E\right)K\left(x_n, x_n, y_n\right)=0,$

where $\phi \left(\lambda \right){\rho }_1\left(\lambda \right)=\rho \left(\lambda \right);\phi \left(\lambda \right){\delta }_1\left(\lambda \right)=\delta \left(\lambda \right);$  $\phi \left(\lambda \right){\gamma }_1\left(\lambda \right)=\gamma \left(\lambda \right).$

By the simple comparison of Equations (50) and (52) we find that these equations are equivalents. Note that the finite-difference Equation (50) has the order of $k,$ and the order of Equation (52) satisfies the condition $k_1<k$ ($k_1$- is the order of Equation (52)). It is known that for the $k_1$- initial values the finite-difference Equation (52) has the unique solution. It is easy to prove that in this case the solution of the Equation (51) will be unique for the $k_1$- initial values ($k_1<k$ satisfies). It is known that the finite-difference equation of the order $k$ has a unique solution if it must be given $k$- initial values (solvability theorem). Obtaining a contradiction shows that the condition of B takes place. Now, let us prove the validity of the condition C.

It is evident that Equation (50) can be written as:

(53)$\rho \left(E\right)\left(y\left(x\right)-f\left(x\right)\right)=O\left(h\right),$

where $x=x_0+nh$ is fixed. If, here, we pass the limit for the $h\to 0$, then we receive:

(54)$\rho \left(1\right)=0,$

so $y\left(x\right)\not\equiv f\left(x\right).$ This condition is called the necessary condition for the convergence of method (33). By taking the equality of (50), we receive the following:

$\left(E-1\right){\rho }_1\left(E\right)\left(y_n-f_n\right)-h\left(\delta \left(E\right)+\gamma \left(E\right)\right)K\left(x_n, x_n, y_n\right)=0$

or,

(55)${\rho }_1\left(E\right)\left(y_{n+1}-y_n-f_{n+1}+f_n\right)-h\left(\delta \left(E\right)+\gamma \left(E\right)\right)K\left(x_n, x_n, y_n\right)=0.$

By changing the meaning of variable $n$ from zero to $m$ and summing the received equalities, one can find the following:

(56)${\rho }_1\left(E\right)\left(y_{m+1}-y_0-f_{m+1}+f_0\right)=\left(\delta \left(E\right)+\gamma \left(E\right)\right)h{\displaystyle \sum _{l=0}^m}K(x_l, x_l, y_l).$

If, here, we pass the limit for $h\to 0$, and take into account the equality of (48), then we receive:

(57)${\rho }_1\left(1\right)\left(y\left(x\right)-y_0-f\left(x\right)+f_0\right)=\left(\delta \left(1\right)+\gamma \left(1\right)\right){\displaystyle \underset{x_0}{\overset{x}{\int }}}F\left(x, s, y\left(s\right)\right)ds,$

where $x=x_0+mh$ is a fixed point.

By the comparison of Equation (57) with Equation (2) and taking into account $y\left(x_0\right)=f\left(x_0\right),$ we receive:

(58)${\rho }_1\left(1\right)=\delta \left(1\right)+\gamma \left(1\right).$

By using $\rho \left(1\right)=0,$ one can write:

${\rho }_1\left(\lambda \right)=\frac{\rho \left(\lambda \right)-\rho \left(1\right)}{\lambda -1}.$

Hence, we have:

${\rho }_1\left(1\right)=\underset{\lambda \to 1}{\lim }\frac{\rho \left(\lambda \right)-\rho \left(1\right)}{\lambda -1} or {\rho }_1\left(1\right)={\rho }^{\prime }\left(1\right).$

By taking this into account in the equality of (58), one can write the following:

(59)${\rho }^{\prime }\left(1\right)=\delta \left(1\right)+\gamma \left(1\right).$

By comparison of equalities (54) and (58) with the first two equations of system (28), we receive the conditions that $p\ge 1$ satisfies. Now, we prove that $\delta \left(1\right)+\gamma \left(1\right)\ne 0,$ and when supposing otherwise, input $\delta \left(1\right)+\gamma \left(1\right)=0$. In this case, from Equation (59) we receive ${\rho }^{\prime }\left(1\right)=0$. Thus, by using Equation (54) we receive $\rho \left(1\right)={\rho }^{\prime }\left(1\right)=0$. It follows from here that λ=1 is twice the root of the polynomial $\rho \left(\lambda \right).$ Now, we prove that in this case method (33) is not convergence. To this end, we use the following error of method (33):

${\epsilon }_m=y\left(x_m\right)-y_m \left(m=0, 1, 2, \dots \right).$

Let us in the equality of (33) change the approximate values $y_m$ by its exact values. Then, we receive:

(60)${\displaystyle \sum _{i=0}^k}({\alpha }_i\left(y\left(x_{n+i}\right)-f\left(x_{n+i}\right)\right)-h{\displaystyle \sum _{j=i}^k}({\beta }_i^{\left(j\right)}K\left(x_{n+j}, x_{n+i}, y\left(x_{n+i}\right)+{\gamma }_i^{\left(j\right)}K\left(x_{n+j+{\nu }_i}, x_{n+i+{\nu }_i}, y\left(x_{n+i}\right)\right)\right))=R_n.$

where $R_n$- is the reminder term.

If we subtract equality (33) from (60), then we receive

(61)${\displaystyle \sum _{i=0}^k}({\alpha }_i{\epsilon }_{n+i}-h{\displaystyle \sum _{j=i}^k}({\beta }_i^{\left(j\right)}L\left({\xi }_{n+i}\right){\epsilon }_{n+i}+{\gamma }_i^{\left(j\right)}\overline{L}\left({\xi }_{n+i}\right)){\epsilon }_{n+i+{\nu }_i})=R_n,$

where

$L_i\left({\xi }_{n+i}\right)=K_y^{\prime }\left(x_{n+j}, x_{n+i}, {\xi }_{n+i}\right); {\overline{L}}_i\left({\overline{\xi }}_{n+i}\right)=K_y^{\prime }\left(x_{n+j+{\nu }_j}, x_{n+i+{\nu }_i}, {\overline{\xi }}_{n+i}\right),$

Variable ${\xi }_{n+i}$ lies between the values of $y_{n+i}$ and $y\left(x_{n+i}\right)$, but ${\overline{\xi }}_{n+i}$ lies between the values of $y_{n+i+{\nu }_i}$ and $y\left(x_{n+i+{\nu }_i}\right)$, respectively. The Equation (61) is the nonhomogeneous finite-difference equation. The corresponding homogeneous equation has the following form:

(62)${\displaystyle \sum _{i=0}^k}{\alpha }_i{\epsilon }_{n+i}=0$

.

Let us note that one can receive Equation (62) from Equation (61) by going to the limit as the $h\to 0$.

As is known, the general solution of a homogeneous finite-difference equation with constant coefficients can be written as the following:

(63)${\epsilon }_m=c_1{\lambda }_1^m+c_2{\lambda }_2^m+\cdots +c_k{\lambda }_k^m \left({\lambda }_i\ne {\lambda }_j if i\ne j\right),$

where ${\lambda }_l\left(l=0, 1, \dots , k\right)$ are the roots of the characteristic polynomial $\rho \left(\lambda \right).$ If we use the condition ${\rho }_1\left(1\right)=\rho \left(1\right)=0$, we find that λ = 1 is twice the root. In this case the solution of Equation (62) can be presented in the following form:

(64)${\epsilon }_m=c_1+c_{2}m+c_3{\lambda }_3^m+\cdots +c_k{\lambda }_k^m.$

As follows from here, the error of the investigated method is unbounded. Hence, if in Equation (64) we pass the limit as the step size tends to zero, (i.e., $h\to 0$), then we receive ${\epsilon }_m\to \infty$. By this way, we have proved that if $\delta \left(1\right)+\gamma \left(1\right)=0,$ then the method does not converge, which contradicts the assumption. Thus, we have proved that the conditions A, B, and C take place.

Note that Dahlquist’s results can be received from the results received here. It follows that the results obtained here are the development of Dahlquist’s rule.

6. Remark

By establishing the direct relation between ODE and Volterra integral equations, we have constructed numerical methods for the solving of both Volterra integral equations and ODE. Now, we want to illustrate that one can construct simple methods for solving ODE and the Volterra integral equation and also Volterra integro-differential equations. For simplicity, we input $f\left(x\right)\equiv 0$, and in this case we have:

$y\left(x_{n+1}\right)={\displaystyle \underset{x_0}{\overset{x_{n+1}}{\int }}}K\left(x_{n+1}, s, y\left(s\right)\right)ds,$

or,

(65)$y\left(x_{n+1}\right)={\displaystyle \underset{x_0}{\overset{x_n}{\int }}}K(x_{n+1}, s, y\left(s\right)ds+{\displaystyle \underset{x_n}{\overset{x_{n+1}}{\int }}}K\left(x_{n+1}, s, y\left(s\right)\right)ds.$

The first integral can be presented as the following:

${\displaystyle \underset{x_0}{\overset{x_n}{\int }}}K(x_{n+1}, s, y\left(s\right))ds={\displaystyle \underset{x_0}{\overset{x_n}{\int }}}K(x_n, s, y\left(s\right)ds+h{\displaystyle \underset{x_0}{\overset{x_n}{\int }}}{K}_X^{\prime }\left({\xi }_n, s, y\left(s\right)\right)ds (x_n<{\xi }_n<x_n+h).$

By taking into account that the function $K\left(x, s, y\right)$ and its derivatives are bounded, then for the fixed point $x_0+nh$, we receive:

$\underset{h\to 0}{\lim }{\displaystyle \underset{x_0}{\overset{x_n}{\int }}}K\left(x_{n+1}, s, y\left(s\right)\right)ds={\displaystyle \underset{x_0}{\overset{x_n}{\int }}}K\left(x_n, s, y\left(s\right)\right)ds.$

By taking this in Equation (65) for the calculation of the value $y_{n+1}\approx y\left(x_{n+1}\right)$, we receive the following:

(66)$y_{n+1}=y_n+{\displaystyle \underset{x_n}{\overset{x_{n+1}}{\int }}}K\left(x_{n+1}, s, y\left(s\right)\right)ds.$

From here one can obtain the following methods:

$y_{n+1}=y_n+hK\left(x_{n+1}, x_{n+1}, y_{n+1}\right); y_{n+1}=y_n+hK\left(x_{n+1}, x_n, y_n\right);$

$y_{n+1}=y_n+hK\left(x_{n+1}, x_n, y_n\right); y_{n+1}=y_n+h\left(K\left(x_{n+1}, x_{n+1}, y_{n+1}\right)+K\left(x_{n+1}, x_n, y_n\right)\right)/2.$

By the generalization of these methods, we receive the following known quadrature formula:

$y_{n+1}=y_n+h{\displaystyle \sum _{i=0}^1}{\displaystyle \sum _{j=i}^1}{\beta }_i^{\left(j\right)}K(x_{n+j}, x_{n+i}, y_{n+i}).$

For the construction of methods of hybrid type it is enough to apply the midpoint rule to the calculation of integrals. In this case, we receive:

$y_{n+1}=y_n+hK\left(x_{n+1}, x_{n+\frac{1}{2}}, y_{n+\frac{1}{2}}\right).$

By using the exact solution of Equation (2), one can write the following:

(67)$y^{\prime }\left(x\right)=K\left(x, x, y\left(x\right)\right)+{\displaystyle \underset{x_0}{\overset{x}{\int }}}{K}_x^{\prime }(x, s, y\left(s\right)ds, y\left(x_0\right)=0.$

If we apply the multistep method with constant coefficients to solve the problem (67), then we receive:

(68)${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}=h{\displaystyle \sum _{i=0}^k}{\beta }_{i}K(x_{n+i}, x_{n+i}, y_{n+i})+h{\displaystyle \sum _{i=0}^k}{\beta }_i\left({\displaystyle \underset{x_0}{\overset{x_n}{{\displaystyle \int }}}}{K}_x^{\prime }\left(x_{n+i}, s, y\left(s\right)\right)ds+{\displaystyle \underset{x_n}{\overset{x_{n+i}}{{\displaystyle \int }}}}{K}_x^{\prime }\left(x_{n+i}, s, y\left(s\right)\right)ds\right).$

From this equality one can write the following:

(69)$h{\displaystyle \sum _{i=0}^k}{\beta }_i{\displaystyle \underset{x_0}{\overset{x_n}{\int }}}{K}_x^{\prime }\left(x_{n+i}, s, y\left(s\right)\right)ds={\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}-h{\displaystyle \sum _{i=0}^k}{\beta }_i(K(x_{n+i}, x_{n+i}, y_{n+i})+{\displaystyle \underset{x_n}{\overset{x_{n+i}}{\int }}}{K}_x^{\prime }\left(x_{n+i}, s, y\left(s\right)\right)ds.$

The right hand side of this equality results in the application of the following method:

(70)${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}=h{\displaystyle \sum _{i=0}^k}{\beta }_i{y}_{n+i}^{\prime },$

to solve the following initial value problem:

$y^{\prime }\left(x\right)=K\left(x, x, y\left(x\right)\right)+{\displaystyle \underset{x_n}{\overset{x}{\int }}}{K}_x^{\prime }\left(x, s, y\left(s\right)\right)ds, y\left(x_n\right)=y_n$

By taking into account that method (70) has the degree of $p$, then from the equality of (69) we find that the following holds:

$h{\displaystyle \sum _{i=0}^k}{\beta }_i{\displaystyle \underset{x_0}{\overset{x_n}{\int }}}{K}_x^{\prime }\left(x_{n+i}, s, y\left(s\right)\right)ds=O(h^{p+1}).$

Thus, we prove that if method (70) has the degree of $p$, then method (68) has also the degree of $p.$

Let us approximate the function of $K_x^{\prime }\left(x_{n+i}, s, y\left(s\right)\right)$ in the following form:

$h{K}_x^{\prime }\left(x_{n+m}, s, y\left(s\right)\right)={\displaystyle \sum _{j=0}^k}{b}_{j}K\left(x_{n+j}, s, y\left(s\right)\right).$

By using this formula, one can write the following:

$h{\displaystyle \sum _{i=0}^k}{\beta }_i{\displaystyle \underset{x_0}{\overset{x_{n+i}}{\int }}}{K}_x^{\prime }\left(x_{n+i}, s, y\left(s\right)\right)ds={\displaystyle \sum _{i=0}^k}{\beta }_i{\displaystyle \underset{x_n}{\overset{x_{n+i}}{\int }}}{\displaystyle \sum _{j=0}^k}{b}_{j}K\left(x_{n+j}, s, y\left(s\right)\right)ds.$

By using some quadrature formulas for the calculation of the definite integral, one can write:

${\displaystyle \sum _{i=0}^k}{\beta }_i{\displaystyle \underset{x_n}{\overset{x_{n+i}}{\int }}}{\displaystyle \sum _{j=i}^k}{b}_{j}K\left(x_{n+j}, s, y\left(s\right)\right)ds=h{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\tilde{\beta }}^{\left(j\right)}K(x_{n+j}, x_{n+i}, y_{n+i}).$

If we take this in the equality of (68), then we receive the following method:

${\displaystyle \sum _{i=0}^k}{\alpha }_i{y}_{n+i}=h{\displaystyle \sum _{i=0}^k}{\displaystyle \sum _{j=i}^k}{\beta }_i^{\left(j\right)}K(x_{n+j}, x_{n+i}, y_{n+i}),$

where ${\beta }_i^{\left(j\right)}={\beta }_i+{\displaystyle \sum _{j=i}^k}{\tilde{\beta }}^{\left(j\right)}.$

This method is the same as method (8).

7. Numerical Results

For the illustration of receiving theoretical results, let us consider the following model:

(71)$y\left(x\right)=1+\lambda {\displaystyle \underset{0}{\overset{x}{\int }}}y\left(s\right)ds,$

the exact solution for which can be presented as $y\left(x\right)=\exp \left(\lambda x\right),$ where $\lambda =const.$

This example well describes the behavior properties of the errors received in the application of the method which we used to solve Equation (2). Note that this integral equation has a direct relation with the following problem:

$y\prime =\lambda y, y\left(0\right)=1,$

describing the behavior of the solution in the following problem:

$y^{\prime }=f\left(x, y\right), y\left(x_0\right)=y_0,$

which was fundamentally investigated by Dahlquist. Example (71) has been solved by using the methods (31), (34), and the following:

(72)$y_{n+1}=y_n+h(y_{n+\alpha }^{\prime }+y_{n+1-\alpha }^{\prime })/2, \alpha =\frac{1}{2}-\sqrt{3}/6.$

Results are tabulated in Table 1.

In Table 1, Table 2, Table 3, Table 4 and Table 5 we have tabulated the results received by the application of methods (34) and (72) to solve example (71) for the different values of step size $h>0$ and parameter λ.

By the simple comparison of the received results one can argue that obtaining results is justified. Note that the results received for the negative $m$ ($m<0$) are better than the received results for the positive $\left(m>0\right).$ It follows from the fact that the exact values of the solution will be sufficiently small for large values of the quantity of $\left|m\right|$  $\left(m<0\right).$

According to the results in Table 1, Table 2, Table 3 and Table 4, we find that the results obtained by using method (34) are better. It is naturally because method (34) is more accurate than method (72), having the degree $p=4.$ Note that these methods have the hybrid type. Now, let us compare the results obtained by the hybrid and advanced methods. For this purpose, let us solve example (71) by the application of the following method, which is received from method (10), having the degree $p=3:$

(73)$y_{n+1}=y_n+h\left(8y_{n+1}^{\prime }+5y_n^{\prime }-y_{n+2}^{\prime }\right)/12$

Results for this method are tabulated in Table 5.

The results received by method (73) can be taken as better, which is explained by the fact that in method (73) the information about the solution at the previous and the next points is used.

8. Conclusions

Here, we have considered the comparison of some numerical methods, which have been applied to solve Volterra integral equations. To this end, we used the conception of stability and degree (order of accuracy) of the investigated methods. We also constructed a formula by which one can define the maximal value for the degree of the stable and unstable multistep methods having different forms (advanced, hybrid, etc.). We prove that the multistep second derivative methods and multistep methods of hybrid type, which have been applied to solve Volterra integral equations, are more exact than the others. These results can be taken as the development of Dahlquist’s results, which were received for the multistep second derivative, and show that hybrid methods have an application to solve Volterra integral equation of the second kind. Additionally, here we find the necessary condition for the convergence of the methods proposed to solve Volterra integral equations. For the investigation of the convergence of proposed methods, here we used the theory of finite-difference equations with constant coefficients. Therefore, multistep methods, here, are investigated in a very simple form. We prove that the initial value problem for ODE and the Volterra integral equation can be solved by one and the same methods. Here, algorithms have been constructed using a similar form to (50) and (51). Some of the received results are illustrated by the model equations. The methods proposed here are promising and we hope that they will find their followers.

Author Contributions

Formulation of the problem, V.I. and M.I.; methodology, V.I.; software, M.I.; validation, V.I. and M.I.; investigation, V.I. and M.I.; writing—original draft preparation, M.I.; writing—review and editing, V.I. and M.I.; funding acquisition, V.I. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not Applicable.

Acknowledgments

The authors wish to express their thanks to academicians Telman Aliyev and Ali Abbasov for their suggestion to investigate the computational aspects of our problem and for their frequent valuable suggestions. This work was supported by the Science Development Foundation under the President of the Republic of Azerbaijan—Grant № EİF-MQM-ETS-2020-1(35)-08/01/1-M-01, (For Vagif Ibrahimov).

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Tables

Table 1

The results received for the case $h=0.05$ and λ = 1, 5, 10.

x	$\mathit{\lambda }$ = 1		$\mathit{\lambda }$ = 5		$\mathit{\lambda }$ $=10$		$\mathit{\lambda }=15$
	Method 34	Method 72	Method 34	Method 72	Method 34	Method 72	Method 34	Method 72
0.1	5.2 × 10−12	5.0 × 10−9	2.4 × 10−8	4.1 × 10−6	1.2 × 10−6	1.0 × 10−4	1.4 × 10−5	8.4 × 10−4
0.4	3.0 × 10−11	2.5 × 10−8	4.6 × 10−7	7.4 × 10−5	1.1 × 10−4	8.4 × 10−3	5.7 × 10−3	3.0 × 10−1
0.7	7.3 × 10−11	6.0 × 10−8	3.6 × 10−6	5.8 × 10−4	3.7 × 10−3	3.0 × 10−1	9.1 × 10−1	4.8 × 101
1.0	1.4 × 10−10	1.1 × 10−7	2.3 × 10−5	3.7 × 10−3	1.0 × 10−1	8.5 × 100	1.1 × 102	6.1 × 103

Table 2

The results received for the case $h=0.05$ and λ = −1, −5, −10.

x	$\mathit{\lambda }=1$		$\mathit{\lambda }=5$		$\mathit{\lambda }=10$
	Method 34	Method 72	Method 34	Method 72	Method 34	Method 72
0.1	1.9 × 10−9	5.5 × 10−7	7.9 × 10−6	4.4 × 10−4	3.5 × 10−4	1.0 × 10−2
0.4	1.7 × 10−8	3.0 × 10−6	2.3 × 10−4	8.0 × 10−3	4.7 × 10−2	7.6 × 10−1
0.7	4.2 × 10−8	7.1 × 10−6	1.9 × 10−3	6.1 × 10−2	1.7 × 100	2.7 × 101
1.0	8.3 × 10−8	1.4 × 10−5	1.2 × 10−2	4.0 × 10−1	5.1 × 101	7.6 × 102

Table 3

The results received for the case $h=0.01$ and m = 1, 5, 10, 15.

x	$\mathit{\lambda }$ $=-1$		$\mathit{\lambda }$ $=-5$		$\mathit{\lambda }$ $=-10$
	Method 34	Method 72	Method 34	Method 72	Method 34	Method 72
0.1	1.7 × 10−9	5.0 × 10−7	9.1 × 10−5	2.9 × 10−3	5.1 × 10−4	1.0 × 10−2
0.4	8.0 × 10−9	1.5 × 10−6	2.8 × 10−5	5.6 × 10−4	3.1 × 10−5	4.4 × 10−4
0.7	1.1 × 10−8	1.9 × 10−6	2.5 × 10−6	4.9 × 10−5	6.4 × 10−7	8.0 × 10−6
1.0	1.2 × 10−8	2.0 × 10−6	1.8 × 10−7	3.4 × 10−6	1.0 × 10−8	1.2 × 10−7

Table 4

The results received for the case $h=0.01$ and λ = −1, −5, −10, −15.

x	$\mathit{\lambda }$ $=-1$		$\mathit{\lambda }$ $=-5$		$\mathit{\lambda }$ $=-10$		$\mathit{\lambda }$ $=-15$
	Method 34	Method 72	Method 34	Method 72	Method 34	Method 72	Method 34	Method 72
0.1	4.3 × 10−12	3.8 × 10−9	9.3 × 10−9	1.6 × 10−6	1.9 × 10−7	1.7 × 10−5	8.8 × 10−7	5.3 × 10−5
0.4	1.4 × 10−11	1.1 × 10−8	8.9 × 10−9	1.5 × 10−6	4.0 × 10−8	3.3 × 10−6	4.2 × 10−8	2.4 × 10−6
0.7	1.9 × 10−11	1.5 × 10−8	3.5 × 10−9	5.7 × 10−7	3.5 × 10−9	2.9 × 10−7	8.2 × 10−10	4.6 × 10−8
1.0	2.0 × 10−11	1.5 × 10−8	1.1 × 10−9	1.8 × 10−7	2.5 × 10−10	2.0 × 10−8	1.3 × 10−11	7.3 × 10−10

Table 5

The results received for the case $h=0.01$ and λ = ±1, ±5, ±10, ±15.

x	m = 1	m = 5	m = 10	m = 15	m = −1	m = −5	m = −10	m = −15
0.1	1.4 × 10−8	1.2 × 10−5	3.1 × 10−4	2.5 × 10−3	1.1 × 10−8	5.0 × 10−6	5.0 × 10−5	1.6 × 10−4
0.4	7.4 × 10−8	2.2 × 10−4	2.5 × 10−2	8.9 × 10−1	3.4 × 10−8	4.4 × 10−6	1.0 × 10−5	7.2 × 10−6
0.7	1.7 × 10−7	1.7 × 10−3	8.7 × 10−1	1.4 × 102	4.4 × 10−8	1.7 × 10−6	8.8 × 10−7	1.4 × 10−7
1.0	3.4 × 10−7	1.1 × 10−2	2.5 × 101	1.8 × 104	4.6 × 10−8	5.5 × 10−7	6.2 × 10−8	2.2 × 10−9