Exponentially convergent symbolic algorithm of the functional-discrete method for the fourth order Sturm-Liouville problems with polynomial coefficients

Exponentially convergent symbolic algorithm of the functional-discrete method for the fourth order Sturm-Liouville problems with polynomial coefficients

Volodymyr Makarov11footnotemark: 1 Nataliia Romaniuk Department of Numerical Mathematics, Institute of Mathematics of National Academy of Sciences of Ukraine, 3 Tereshchenkivs’ka Str., 01004 Kyiv-4, Ukraine
Abstract

A new symbolic algorithmic implementation of the functional-discrete (FD-) method is developed and justified for the Sturm-Liouville problem on a finite interval in the Hilbert space for the fourth order ordinary differential equation with polynomial coefficients. The sufficient conditions of an exponential convergence rate of the proposed approach are received. The obtained estimates of the absolute errors of FD-method significantly improve the accuracy of the estimates obtained earlier. The proposed algorithm of our method is developed when the potential functions are approximated by the identically zero, and for this case, FD-method can be considered as one of the variants of the homotopy method. Our algorithm is symbolic and operates with the decomposition coefficients of the eigenfunction corrections in some basis. The number of summands in these decompositions depends on the degree of the potential coefficients and the correction number. Our method uses only the algebraic operations and basic operations on column vectors and matrices and does not need solutions of any boundary value problems and computations of any integrals, unlike the previous variants of FD-method. The corrections to eigenpairs are computed exactly as analytical expressions, and we have no rounding errors. The numerical examples illustrate the theoretical results. The numerical results obtained with the FD-method are compared with the numerical test results obtained with other existing numerical techniques.

keywords:
Fourth order Sturm-Liouville problems, Eigenvalue problems, Polynomial coefficients, Functional-discrete method, Symbolic algorithm, Exponential convergence rate
Msc:
[2010] 65L15, 65L20, 65L70, 34B09, 34B24, 34L16, 35G15
journal: Journal of Computational and Applied Mathematics, ISSN: 0377-0427, myfootnote1myfootnote1footnotetext: Email address: makarov@imath.kiev.ua (Volodymyr Makarov)

1 Introduction

Numerical solving the mathematical models of various physical phenomena, in particular, which described by spectral problems, can be reduced to the computationally expensive problems. For example, they are the following: solving stiff ordinary differential equations, solving ill-conditioned systems of linear algebraic equations, computing of large sums and recurrence formulas, execution of time-consuming and complex computations on large data arrays. For the correct solution of such problems, it is necessary to use a high-precision computation (see Bailey2005 (); BaileyBorwein2015 ()), the methods organizing of which involve symbolic computations (see Bailey2016 ()). That is why the development of more efficient and less computationally costly symbolic algorithms for the high-precision analytical methods is an important and vital task today.

There is a great number of numerical methods for Sturm-Liouville problems for the second- and higher-order ordinary differential equations. The analytical methods based on perturbation and homotopy ideas 115_GMR_Arm1979 (); 118_GMR_AllgGeorg1990 (), such as functional-discrete methods from MakarovKlymenko2007 (); GMR-GavrMakPop2010 (); MakRom2014Dir (); MakRom2014DirNeum (); Conf_MakRom2015 (); Conf_Rom2015 (); MakRom2017 (); DemGavrMak2016 (); GavrMakRom2017 (); GavrMakRom20152017 () (compare with Adomian decomposition method 65_Rach2012 (); MRL-Adomian1993 ()), have been widely used for solving the eigenvalue problems. Using these analytical methods the solutions can be found as fast convergent functional series, the properties of the solution of the original problem can be investigated with the approximation solution, and, to top it all, these approaches allow in a natural way the use of computer algebra systems for developing and implementation symbolic-numerical algorithms.

Suggested by V. Makarov Mak1991 () in 1991 and developed by V. Makarov and his disciples, the FD-method enables us to overcome the many disadvantages of the discrete methods such as follows: the accuracy degradation with the increasing of the eigenvalue index; usage of the mesh generated at the start of the numerical process; saturation of accuracy; the number of reliable numerical eigenvalues is limited and depends on a mesh step (see Pryce1993 (); Zhang2015 ()). The main advantages of FD-method are his features, which differ from many other methods: 1) the approach can be applied to operator equations in general form; 2) the approach can be applied also to eigenvalue problems with multiple eigenvalues (see, for example, GavrMakRom20152017 ()); 3) all eigenpairs can be computed in parallel; 4) the convergence rate increases as the index of the eigenpair increases; 5) it was proved that in many cases the FD-method converges exponentially or super-exponentially.

Presented modification of the traditional algorithm of the functional-discrete (FD-) method was proposed in MakarovKlymenko2007 (). The general idea of the symbolic algorithms for FD-method is the representation of the eigenfunction corrections in some basis. Then we obtain recurrence relations for the decomposition coefficients. Finally, our algorithm operates with these decomposition coefficients and, in certain cases, uses only the algebraic operations and does not need solutions of any boundary value problems and computations of any integrals. And more than this, the corrections to eigenpairs are computed exactly as analytical expressions, and we have no rounding errors. But if computational difficulties or memory overflow arise then we can avoid combinatorial explosion. In proposed approach instead of using rational arithmetic we can easily transit to floating–point arithmetic which ”represents an alternative idea: round the computation at every step, not just at the end” (see Trefethen2015 ()).

Briefly described general idea was used for developing and justification of new symbolic algorithms for FD-method for Sturm-Liouville problems on a finite interval for the Schrödinger equation with a polynomial potential in MakRom2014Dir (); MakRom2014DirNeum (); MakRom2017 (), and for several linear differential operators with fractional derivatives in DemGavrMak2016 (); GavrMakRom2017 (). In this article, we use the described idea for the fourth order Sturm-Liouville problem with polynomial coefficients in order for the traditional method from GMR-GavrMakPop2010 (); GavrMakRom20152017 () to be modified and for a new symbolic algorithm of the FD-method to be developed. Note that some of the results of this article were announced in Conf_MakRom2015 (); Conf_Rom2015 (), but unlike the symbolic algorithms from ones, the approach presented below produces explicit recursive formulas for the decomposition coefficients of the representation for the eigenfunctions corrections.

The article is organized as follows. Section 2 deals with the problem statement. Section 3 contains the traditional algorithm of the simplest variant of the FD-method for the case when the potential coefficients are approximated by the identically zero, namely, FD-method is the purely analytical method. In Section 4 a new structural representation of the eigenfunctions corrections is obtained. This representation is used in Section 6 to develop a new symbolic algorithm of the FD-method. The sufficient conditions of an exponential convergence rate of the proposed approach are received in Section 5, and the obtained absolute errors estimates of the FD-method significantly improve the accuracy of the estimates obtained earlier in GMR-GavrMakPop2010 (). The formulas for the proposed new symbolic algorithmic implementation of our method are given in Section 6 and the formulas from appendices A, B, C are directly used in the calculations. The numerical algorithm is given in Section 7. Section 8 illustrates the theoretical results by a numerical examples. The numerical results obtained with the FD-method are compared with the numerical test results obtained with other existing numerical techniques 58_AttiliLesnic2006 (); 25_SyamSiyyam2009 (); HPM2010Ex1 (); HAM2011Ex1 (); 59_Chanane2010 (); RATTANA2013144 () in Examples 1 and 2. A review of the obtained results with implementation features and advantages of our developed numerical method are given in the final Section 9.

2 Problem statement

Let us consider the regular Sturm-Liouville problem in a Hilbert space for the fourth order ordinary differential equation

 uIV(x)+q2(x)u′′(x)+q1(x)u′(x)+(q0(x)−λ)u(x)=0, (1)
 x∈(0,X),0

with the boundary conditions

 u(0)=u′′(0)=u(X)=u′′(X)=0, (2)

where is the real constant. The real-valued polynomial coefficients are

 q0(x)=r0∑l=0Alxl,q1(x)=r1∑l=0Blxl,q2(x)=r2∑l=0Clxl,ri≥1,i=0,1,2, (3)

where constants , are the positive integer, and , .

In GMR-GavrMakPop2010 () it was shown that the fourth order ordinary differential equation with all derivatives of the eigenfunction could be reduced to the form (1) using the variable transformation. That is why we consider the eigenvalue problem with equation (1).

3 Traditional algorithm of the simplest variant of the FD-method

In a current article, the simplest variant of the FD-method is applied to the Sturm-Liouville problem (1)–(3). It means that we consider the simplest case of the approximation of potential coefficients (3) by the identically zero

 ¯qs(x)≡0,s=0,1,2.

This section contains the traditional algorithm of the FD-method which in this case is the purely analytical method. In this case, FD-method can be considered as one of the variants of the homotopy method (see 118_GMR_AllgGeorg1990 (); 115_GMR_Arm1979 ()), and his idea is closely related to the ideas of the Adomian decomposition method MRL-Adomian1993 (); 65_Rach2012 (). Below the symbolic algorithm for the simplest variant of the FD-method is developed.

Note that in the case when the simplest variant of the FD-method (with , ) is divergent for the smallest eigenvalues of problem (1)–(3), the general scheme of the FD-method (usually) with piecewise-constant approximations to potential coefficients is used. Developing and justification of a symbolic algorithmic implementation of the general scheme of the FD-method for the problem (1)–(3) are slated for the near future.

The exact solution of the eigenvalue problem (1)–(3) is then represented by the series

 un(x)=∞∑j=0u(j)n(x),λn=∞∑j=0λ(j)n, (4)

provided that these series converge. The sufficient conditions for the convergence of the series (4) will be presented later in Section 5. The approximate solution to the problem (1)–(3) is represented by a pair of corresponding truncated series, namely,

 mun(x)=m∑j=0u(j)n(x),mλn=m∑j=0λ(j)n, (5)

which is called an approximation of rank . The summands of series (5) (the corrections to eigenpairs , at the –th step of the FD-method) are the solutions of the following recursive sequence of problems (see GMR-GavrMakPop2010 ())

 d4dx4u(j+1)n(x)−λ(0)nu(j+1)n(x)=F(j+1)n(x),x∈(0,X), (6)
 u(j+1)n(0)=d2u(j+1)n(0)dx=u(j+1)n(X)=d2u(j+1)n(X)dx=0, (7)
 j=0,1,...,m−1.

where

 F(j+1)n(x)=j∑p=0λ(j+1−p)nu(p)n(x)−q2(x)d2dx2u(j)n(x)−q1(x)ddxu(j)n(x)−q0(x)u(j)n(x). (8)

Using the solvability condition

 (F(j+1)n(x),u(0)n(x))L2(0,X)=∫X0F(j+1)n(x)u(0)n(x)dx=0 (9)

of problems (6)–(8) for a fixed , we obtain the following formula for the eigenvalue corrections:

 λ(j+1)n=∫X0(q2(x)d2dx2u(j)n(x)+q1(x)ddxu(j)n(x)+q0(x)u(j)n(x))u(0)n(x)dx. (10)

Solutions to the problems (6)–(8) satisfies the orthogonality condition

 (u(j+1)n(x),u(0)n(x))L2(0,X)=∫X0u(j+1)n(x)u(0)n(x)dx=0, (11)
 j=0,1,...,m−1.

The initial approximation , is the solution of the so-called base problem, that is,

 d4u(0)n(x)dx4−λ(0)nu(0)n(x)=0,x∈(0,X), (12)
 u(0)n(0)=d2u(0)n(0)dx=u(0)n(X)=d2u(0)n(X)dx=0. (13)

The solution of (12)–(13) is the following

 u(0)n(x)=√2Xsin(nπXx),λ(0)n=(nπ)4X4. (14)

4 Representation of the corrections to eigenfunctions u(j)n(x)

Let us introduce the generalized Green’s function for a linear differential operator, corresponding to the problems by (6)–(14), in the following form

 gn(x,ξ)=2X3π4∞∑p=1,p≠nsin(pπXx)sin(pπXξ)p4−n4,x,ξ∈[0,X]. (15)

The function (15) can be expressed as

 gn(x,ξ)=X32π3n3[−(xX−H(x−ξX))cos(nπXx)sin(nπXξ)−(ξX−H(ξ−xX))sin(nπXx)cos(nπXξ)++1sinh(πn)sinh(πn(xX−H(x−ξX)))×sinh(πn(ξX−H(ξ−xX)))+32πnsin(nπXx)sin(nπXξ)], (16)
 x,ξ∈[0,X],

where is the Heaviside function, and .

Lemma 1

The generalized Green’s function (15), (16) has the following properties:

 gn(x,ξ)=gn(ξ,x),gn(x,ξ)=gn(X−x,X−ξ),
 ∫X0gn(x,ξ)sin(πnXx)dx=0.

For a fixed j, the solution of a problem (6)–(9), (3), satisfyig the orthogonality condition (11), can be expressed by a formula

 u(j+1)n(x)=∫X0gn(x,ξ)F(j+1)n(ξ)dξ. (17)

Considering the properties of the problems (6)–(9), (3), of the integral representation (17), of the generalized Green’s function (15), (16) (see Lemma 1), and using the solution of the base problem (14), we prove the assertion by a method of complete induction.

Lemma 2

The solution of problem (6)–(9), (3) can be represented by

 u(j+1)n(x)=M(j+1)∑p=0xp[b(j+1)n,pcos(πnXx)+a(j+1)n,psin(πnXx)]+M(j)∑p=0xp[d(j+1)n,pcosh(πnXx)+c(j+1)n,psinh(πnXx)],j=0,1,...,n=1,2,... (18)

where , , .

In Lemma 2 the coefficients

 b(j+1)p,a(j+1)p(p=0,1,...,M(j+1)),c(j+1)s,d(j+1)s(s=0,1,...,M(j)) (19)

are the decomposition coefficients of the eigenfunction corrections in the basis

 xpcos(πnXx),xpsin(πnXx)(p=0,1,...,M(j+1)),
 xscosh(πnXx),xssinh(πnXx)(s=0,1,...,M(j))

on interval . Unlike (17), the representation (18) is used below to develop a new symbolic algorithmic implementation of the FD-method, numerical algorithm for which is fully given in Section 7. Exact explicit recursive formulas for these coefficients are found in Section 6 (see (57), (58), (59) and (60)).

5 Convergence of the FD-method

To investigate the convergence of the FD-method, using the integration by parts as well as the representation for the generalized Green’s function by a series (15), we take the expressions (10) and (17) into formulas

 (20)

and

 u(j+1)n(x)=2X3π4∞∑p=1,p≠n1p4−n4sin(pπXx)×∫X0{[(pπX)2q2(ξ)sin(pπXξ)+pπX[−2q′2(ξ)+q1(ξ)]cos(pπXξ)+[−q′′2(ξ)+q′1(ξ)−q0(ξ)]sin(pπXξ)]u(j)n(ξ)+j∑s=0λ(j+1−s)nu(s)n(ξ)sin(pπXξ)}dξ. (21)

Note that the obtained in this Section absolute errors estimates of the FD-method significantly improve the accuracy of the estimates obtained earlier in GMR-GavrMakPop2010 ().

For the eigenvalue corrections one can deduce from (20) the next estimate

 ∣∣λ(j+1)n∣∣≤ω√2X[(nπX)2+nπX+1]∥∥u(j)n∥∥, (22)

where

 ω=max{∥q2∥∞,∥∥2q′2−q1∥∥∞,∥∥q′′2−q′1+q0∥∥∞},∥v∥∞=maxx∈[0,X]|v(x)|,
 ∥v∥=√(v(x),v(x))L2(0,X)=(∫X0[v(x)]2dx)1/122.

It is easy to establish that the following inequalities are correct:

  ⎷∞∑p=1,p≠np4−2k(p4−n4)2≤n1−k2n2−2n+1,k=0,1,2,n=1,2,...

which are used to obtain from (22) the estimate for the eigenfunction corrections (21):

 ∥∥u(j+1)n∥∥≤X2π212n2−2n+1[[n∥q2∥∞+Xπ∥∥−2q′2+q1∥∥∞+X2π2∥∥−q′′2+q′1−q0∥∥∞n]∥∥u(j)n∥∥+√2Xω[n+Xπ+X2nπ2]×j∑s=0∥∥u(j−s)n∥∥∥∥u(s)n∥∥]≤X2π2ω2n2−2n+1[n+Xπ+X2nπ2]×[∥∥u(j)n∥∥+√2Xj∑s=0∥∥u(j−s)n∥∥∥∥u(s)n∥∥]≤Mnj∑s=0∥∥u(j−s)n∥∥∥∥u(s)n∥∥, (23)

where

 Mn=X2π2ω2n2−2n+1[n+Xπ+X2nπ2]max{1,√2X}. (24)

Substituting in (23)

 Uj=M−jn∥∥u(j)n∥∥,U0=∥∥u(0)n∥∥=1 (25)

and replacing the new variables by the majorant variables subject to and , we come to the majorant equation

 ¯¯¯¯Uj+1=j∑s=0¯¯¯¯Uj−s¯¯¯¯Us, (26)

which is nonlinear recurrence relation, as well as the so-called convolution-type equation. The solution of the equation (26) is (see, e.g., (Vilenkin1969, , p. 159-161,210), ReinNievDeo1977 ())

 ¯¯¯¯Uj+1=(2j+2)!(j+1)!(j+2)!=4j+12(2j+1)!!(2j+4)!!. (27)

Returning to the old variables (see substitution of variables (25)), we obtain from (27) the following estimate for the solution of (23):

 ∥∥u(j+1)n∥∥≤(4Mn)j+12(2j+1)!!(2j+4)!!≤(4Mn)j+1(j+2)√π(j+1), (28)

and then, from (22), the next estimate for the eigenvalue corrections

 ∣∣λ(j+1)n∣∣≤ω√2X[(nπX)2+nπX+1](4Mn)j2(2j−1)!!(2j+2)!!≤ω√2X[(nπX)2+nπX+1](4Mn)j(j+1)√πj. (29)

The last part of inequalities (28) and (29) was obtained using the reflections like those from the proof of the Wallis formula (see, e.g., (Fichtenholz1968, , p. 344)). From estimates (28) and (29) follows the next theorem which contains the sufficient conditions of an exponential convergence rate of the FD-method and estimates of its absolute errors.

Theorem 1

Let , and the following condition holds true:

 rn=4Mn<1,n=1,2,.... (30)

Then the FD-method for the Sturm–Liouville problem (1)–(3) converges exponentially and the following estimates of the absolute errors are valid:

 (31)
 ∥∥un−mun∥∥≤2(rn)m+11−rn(2m+1)!!(2m+4)!!≤(rn)m+1(m+2)√π(m+1). (32)

6 Basic formulas for the symbolic algorithmic implementation of the FD-method

Further, in this section, we describe and develop a new symbolic algorithm of the FD-method for the problem (1)–(3). Note that in this section some formulas of the symbolic algorithm are given in shorthand notations as a matter of convenience in operation. The formulas from Appendices A, B, C are directly used in the calculations. The computer algebra system Maple was used for a software implementation. Unlike the symbolic algorithms from Conf_MakRom2015 (); Conf_Rom2015 (), the approach presented here produces explicit recursive formulas for the coefficients in (18) at the –th step of the FD-method.

Let us substitute (18) into (8), (3) and group together the summands as follows

 F(j+1)n(x)=F(j+1)n,cos(x)cos(πnXx)+F(j+1)n,sin(x)sin(πnXx)+F(j+1)n,cosh(x)cosh(πnXx)+F(j+1)n,sinh(x)sinh(πnXx). (33)

Then we change the order of summation in analytical expressions for , , and from (33), and we group together the summands as follows

 (34)

We arrive at

 F(j+1)n,cos(x)=M(j)∑t=0xtj∑s=]]tr+1[[λ(j+1−s)nb(s)n,t−M(j+1)−3∑t=0xtmin(r,t)∑l=max(0,t−M(j)+2)b(j)n,t−l+2Cl(t−l+2)(t−l+1)+M(j+1)−1∑t=0xtmin(r,t)∑l=max(0,t−M(j))φ(j+1)n,cos,l,t(πnX)−M(j+1)−2∑t=0xtmin(r,t)∑l=max(0,t−M(j)+1)[ddξφ(j+1)n,sin,l,t+1(ξ)]ξ=πnX(t−l+1), (35)

where is the smallest integer greater than or equal to a real number (in the computer algebra system Maple is the function ceil(y)), and

 φ(j+1)n,cos,l,t(ξ)=b(j)n,t−l(−Al+Clξ2)−a(j)n,t−lBlξ,
 φ(j+1)n,sin,l,t(ξ)=a(j)n,t−l(−Al+Clξ2)+b(j)n,t−lBlξ.

In order to obtain the expression for we replace in the formula (35)

 b(s)n,tbya(s)n,t,b(j)n,t−l+2% bya(j)n,t−l+2,

and

 φ(j+1)n,cos,l,tbyφ(j+1)n,sin,l,t,ddξφ(j+1)n,sin,l,t+1(ξ)by−ddξφ(j+1)n,cos,l,t+1(ξ).

By a similar way the formula for may be written as

 F(j+1)n,cosh(x)=M(j−1)∑t=0xtj∑s=]]tr+1[[+1λ(j+1−s)nd(s)n,t−M(j)−3∑t=0xtmin(r,t)∑l=max(0,t−M(j−1)+2)d(j)n,t−l+2Cl(t−l+2)(t−l+1)−M(j)−1∑t=0xtmin(r,t)∑l=max(0,t−M(j−1))φ(j+1)n,cosh,l,t(πnX)−M(j)−2∑t=0xtmin(r,t)∑l=max(0,t−M(j−1)+1)[ddξφ(j+1)n,sinh,l,t+1(ξ)]ξ=πnX(t−l+1), (36)

where

 φ(j+1)n,cosh,l,t(ξ)=d(j)n,t−l(Al+Clξ2)+c(j)n,t−lBlξ,
 φ(j+1)n,sinh,l,t(ξ)=c(j)n,t−l(Al+Clξ2)+d(j)n,t−lBlξ.

In order to obtain the expression for we replace in the formula (36)

 d(s)n,tbyc(s)n,t,d(j)n,t−l+2% byc(j)n,t−l+2,

and

 φ(j+1)n,cosh,l,tbyφ(j+1)n,sinh,l,t,ddξφ(j+1)n,sinh,l,t+1(ξ)byddξφ(j+1)n,cosh,l,t+1(ξ).

To extract the coefficients of like

 f(j+1)n,cos,p,f(j+1)n,sin,p,p=0,1,...,M(j+1)−1

and

 f(j+1)n,cosh,p,f(j+1)n,sinh,p,p=0,1,...,M(j)−1

in the polynomials , , , (see Appendix A) noted above, we can use the function coeff with corresponding arguments in Maple. These coefficients are included in the main formulas of the proposed algorithm in Section 7. The expressions for these coefficients involve only the algebraic operations and are represented through the corresponding quantities computed at previous steps of FD-method.

We require that the polynomials in the front of corresponding trigonometric functions and hyperbolic trigonometric functions are equal on the both sides of equation (6), (34). This requirement leads to the two recurrence systems (37), (38) (with the initial conditions (39), (40)) and (41), (42) (with the initial conditions (43), (44)) for the unknown coefficients of representation (18).

The first system is the following:

 ((((t+4)b(j+1)n,t+4+4πnXa(j+1)n,t+3)(t+3)−6(πnX)2b(j+1)n,t+2)(t+2)−4(πnX)3a(j+1)n,t+1)(t+1)=f(j+1)n,cos,t, (37)
 ((((t+4)a(j+1)n,t+4−4πnXb(j+1)n,t+3)(t+3)−6(πnX)2a(j+1)n,t+2)(t+2)+4(πnX)3b(j+1)n,t+1)(t+1)=f(j+1)n,sin,t, (38)
 t=0,1,...,M(j+1)−4,j=0,1,2,...

with the initial conditions

 a(j+1)n,M(j+1)−s=s∑k=0(−1)[[k2]]+1f(j+1)n,χ(k),M(j+1)−s−1+k×(2k+1)(M(j+1)−1)k22+k(M(j+1)−s+[[k+δ2,s2]])(Xπn)3+k, (39)
 (40)

where , , , , . Here and below denotes the Kronecker delta, is the greatest integer less than or equal to a real number (in the computer algebra system Maple is the function floor(y)).

The second system is the following:

 ((((t+4)d(j+1)n,t+4+4πnXc(j+1)n,t+3)(t+3)+6(πnX)2d(j+1)n,t+2)(t+2)+4(πnX)3c(j+1)n,t+1)(t+1)=f(j+1)n,cosh,t, (41)
 ((((t+4)c(j+1)n,t+4+4πnXd(j+1)n,t+3)(t+3)+6(πnX)2c(j+1)n,t+2)(t+2)+4(πnX)3d(j+1)n,t+1)(t+1)=f(j+1)n,sinh,t, (42)
 t=0,1,...,M(j)−4,j=1,2,...

with the initial conditions

 (43)
 (44)

where , , , , .

Let us introduce the column vectors:

 →Z[a,b]n(p,j)=[a(j+1)n,M(j+1)−p,b(j+1)n,M(j+1)−p]T,<