# Spectral Collocation Method for Handling Integral and Integrodifferential Equations ofn-th Order via Certain Combinations of Shifted Legendre Polynomials

Authors: Zineb Laouar; Nouria Arar; Abdellatif Ben Makhlouf
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this study, an accurate and efficient numerical method based on spectral collocation is presented to solve integral equations and integrodifferential equations ofn-th order. The method is developed using compact combinations of shifted Legendre polynomials as a spectral basis and shifted Legendre–Gauss–Lobatto nodes as collocation points to construct the appropriate algorithm that makes simple systems easy to solve. The technique treats both types of equations: linear and nonlinear equations. The study aims to provide the relevant spectral basis by the use of compact combinations, which allows us to take advantage of shifted Legendre polynomials and to reduce the dimension of the space of approximation. The reliability of the proposed algorithms is proven via different examples of several cases and the results are discussed to confirm the effectiveness of the spectral approach.


---

## Body

## 1. Introduction

In the last years, integral equations and integrodifferential equations received a lot of interest from researchers in many areas of science. Mathematicians and physicists have given importance to this kind of equations used in fundamental problems in biology, medicine [1], chemistry, electrostatics, fluid dynamics, physics [2], economics [3], engineering [4], mechanics [5], potential theory [6], and problems of gravitation [7].Integral equations which consider an unknown function that appears under the sign of integration can be classified into two categories depending on the boundaries of integration (Fredholm and Volterra equations). When the equation also contains the derivative of the unknown function, we obtain an integrodifferential equation that has the same order of the order of derivation appearing in the equation.These types of equations, generally difficult to solve directly by classical methods, were the main subjects of numerous numerical methods that aim to approximate the solution by a reliable algorithm [8]. For this reason, mathematicians are still developing numerical methods to solve both integral and integrodifferential equations [9]. Every method has its advantages and inconveniences, so scientists are always working to improve various methods, especially with the development of computer sciences [10].When talking about numerical methods for integral and integrodifferential equations, many studies have been interested in them. Doha et al. [11] used shifted Jacobi polynomials in the spectral collocation method to solve integrodifferential equations and system of integrodifferential equations. In [12], the authors developed a collocation procedure using cubic B-splines for linear and nonlinear Fredholm and Volterra integral equations. Fathy et al. [13] presented a Legendre–Galerkin method for the linear Fredholm integrodifferential equations. Nemati in [14] used the shifted Legendre polynomials to approximate the solution of Volterra-Fredholm integral equations. Others such as Černá et al. [9] and Moghaddam et al. [15,16] used b-spline wavelets, Maleknejad et al. [17] used Haar wavelets, and Lakestani et al. [18] used multiwavelets to solve integral and integrodifferential equations numerically. Moghaddam et al. in [19] developed the fractional finite differences method, Biçer et al. [10] investigated Bernoulli polynomials, Doha et al. [20] used ultraspherical polynomials, Jalilian et al. [21] used the exponential spline method, Meng et al. [22] and El-Sayed et al. [23] used alternative (shifted) Legendre polynomials, Mandal et al. [24] used Bernstein polynomials, Machado et al. [25] reduced differential transform, Loh et al. [26] used the Laplace transform and resolvent kernel method, and Mokhtary et al. used spectral methods [27] to solve numerically integrodifferential equations.Recently, spectral methods (including all its categories: Galerkin, collocation, and tau) gained a high level of interest to solve different kinds of equations, especially to solve integral equations and integrodifferential equations [28]. Due to their simple application, spectral methods provide interesting results compared to other methods. The idea is to develop the solution as a finite sum of special basis functions [29] (generally orthogonal polynomials or a combination of orthogonal polynomials) [30] to obtain the simplest system where the solution builds the coefficients of the approximation in the chosen basis [31]. Thanks to the generous properties of orthogonal polynomials (Legendre and Chebyshev, …), spectral methods have speed convergence, which means that with small data, we obtain a great rate of convergence and gain spectral accuracy [32]. The choice of basis functions takes a very important part in the application of spectral methods [33]. When choosing the suitable basis [34], the resulted systems would be simple with a special structure of matrices, easy to invert, which reduce the cost of the method and provide efficiency and accuracy.Our motivation in this study is to elaborate a spectral approximation governing both integral and integrodifferential equations of linear and nonlinear cases. We apply the shifted Legendre–Gauss collocation method by the use of a certain combination of shifted Legendre polynomials as basis functions and the nodes of shifted Legendre–Gauss interpolation as collocation points. The initial condition of the problem is invested in the choice of the basis which leads in all cases to simple systems, which will be solved directly by Gauss elimination in the linear case and by using a Newton algorithm in the nonlinear case.The study is divided into five sections. We present first the main properties of Legendre polynomials and shifted Legendre polynomials. Then, in Section3, we start describing the proposed method by introducing the studied problem and characterizing the main steps of the method. Next, we study separately the linear and the nonlinear cases and we apply the Gauss quadrature in each case. In Section 4, several numerical examples are presented and discussed to confirm the reliability of the described method. Section 5 recaps all details as conclusion.


## 2. Preliminaries

We start by presenting basic facts of Legendre polynomials, followed by introducing important properties of shifted Legendre polynomials that helps in developing the proposed method.


### 2.1. Legendre Polynomials

We denote byPN the space of polynomials of degree less than or equal to N and LNx,x∈−1,1 the Legendre polynomial.The standard Legendre polynomials are orthogonal polynomials of degreeN defined on −1,1 resulting from the following differential Legendre equation (for more details see [29–31]):(1)1−x2LN′′x−2xLN′x+NN+1LNx=0,N≥0.The first Legendre polynomials are given by(2)L0x=1,L1x=x,L2x=3x2−12.In particular, we have(3)LN1=1,LN−x=−1NLNx,LN−1=−1N.LN′1=NN+12,LN′−1=NN+12−1N+1,LN″x=∑k=0N−2k+12NN+1−kk+1Lkx.They constitute a Hilbert basis ofL2−1,1 with the weight function wx=1 with the following norm and inner product:(4)fx2=fx,fx,fx,gx=∫−11fxgxdx.They also satisfy the relation of orthogonality relative to the weightwx=1 in −1,1:(5)∀N≠M∈ℕ,∫−11LNxLMxdx=0,∫−11LN2xdx=22N+1.They satisfy the recurrence relations:(6)1−x2LN′x=−NxLNx+NLN−1x,LN+1′x=N+1LNx+xLN′x,and the Rodrigues formula(7)LNx=−1N2NN!dNdxN1−x2N.


### 2.2. Shifted Legendre Polynomials

In order to generalize the use of these polynomials in any interval of the form0,l, the shifted Legendre polynomials are introduced by implementing the change of variable x=2t/l−1; in the special case, x=2t−1 for 0,1. We denote the obtained polynomial LN2t−1 by L˜Nt. Then, the shifted Legendre polynomials can be obtained using the following recurrence formula (for more details, see [22, 23]).(8)L˜N+1t=2N+12t−1N+1L˜Nt−NN+1L˜N−1t,N=1,2,…,and the first ones are as follows:(9)L˜0t=1,L˜1t=2t−1,L˜2t=6t2−6t+1.Note the special cases as follows:(10)L˜N0=−1N,L˜N1=1.These polynomials also satisfy the orthogonality condition with respect to the weight functionwt=1 on 0,1 with the inner product:(11)L˜Nt,L˜Mt=∫01L˜NtL˜Mtdt=12N+1δNM.


## 2.1. Legendre Polynomials

We denote byPN the space of polynomials of degree less than or equal to N and LNx,x∈−1,1 the Legendre polynomial.The standard Legendre polynomials are orthogonal polynomials of degreeN defined on −1,1 resulting from the following differential Legendre equation (for more details see [29–31]):(1)1−x2LN′′x−2xLN′x+NN+1LNx=0,N≥0.The first Legendre polynomials are given by(2)L0x=1,L1x=x,L2x=3x2−12.In particular, we have(3)LN1=1,LN−x=−1NLNx,LN−1=−1N.LN′1=NN+12,LN′−1=NN+12−1N+1,LN″x=∑k=0N−2k+12NN+1−kk+1Lkx.They constitute a Hilbert basis ofL2−1,1 with the weight function wx=1 with the following norm and inner product:(4)fx2=fx,fx,fx,gx=∫−11fxgxdx.They also satisfy the relation of orthogonality relative to the weightwx=1 in −1,1:(5)∀N≠M∈ℕ,∫−11LNxLMxdx=0,∫−11LN2xdx=22N+1.They satisfy the recurrence relations:(6)1−x2LN′x=−NxLNx+NLN−1x,LN+1′x=N+1LNx+xLN′x,and the Rodrigues formula(7)LNx=−1N2NN!dNdxN1−x2N.


## 2.2. Shifted Legendre Polynomials

In order to generalize the use of these polynomials in any interval of the form0,l, the shifted Legendre polynomials are introduced by implementing the change of variable x=2t/l−1; in the special case, x=2t−1 for 0,1. We denote the obtained polynomial LN2t−1 by L˜Nt. Then, the shifted Legendre polynomials can be obtained using the following recurrence formula (for more details, see [22, 23]).(8)L˜N+1t=2N+12t−1N+1L˜Nt−NN+1L˜N−1t,N=1,2,…,and the first ones are as follows:(9)L˜0t=1,L˜1t=2t−1,L˜2t=6t2−6t+1.Note the special cases as follows:(10)L˜N0=−1N,L˜N1=1.These polynomials also satisfy the orthogonality condition with respect to the weight functionwt=1 on 0,1 with the inner product:(11)L˜Nt,L˜Mt=∫01L˜NtL˜Mtdt=12N+1δNM.


## 3. Shifted Legendre Collocation Method

We are interested in using a shifted Legendre collocation method to solve the integrodifferential equation ofn-th order:(12)a0xux+∑i=1naixuix=fx+λ∫01Kx,tFt,utdt,joined to the following initial conditions:(13)uj0=0,j=0,…,n−1,where uix=diu/dxix denotes the i-th derivative of ux, ai,i=0,…,n are the given functions in C0,1, Kx,t is the known kernel which is continuous and square integrable function, F is the given function of ut which can be linear Ft,ut=ut or nonlinear (like: Ft,ut=uts,s≠0,1), λ is a given constant, and f∈L20,1 is a known function.The order of derivationn denotes the order of the integrodifferential equation (12), when aix=0,∀x∈0,1,∀i=1,…,n, we derive the case of integral equations.We defineSN as(14)SN=spanL˜0x,L˜1x,…,L˜Nx,where L˜kx,k=0,…,N are the shifted Legendre polynomials defined on 0,1.We setVN⊂SN, the subspace where the initial conditions (13) are verified(15)VN=v∈SNsuch thatvj0=0,forj=0,…,n−1.The spectral scheme to solve (12) is to find uN∈VN, such that for all v∈VN,(16)a0xuNx,vx+∑i=1naixuNix,vx=fx,vx+λ∫01Kx,tFt,uNtdt,vx,where u,v=∫01uxvxdx denotes the inner product in the space L20,1.When applying a spectral method, one wonders about the choice of the appropriate basis such that the obtained system is the simplest possible. Therefore, we look to use compact combinations of orthogonal polynomials as basis functions to gain more efficiency. Many studies developed several types of combinations for different equations [32, 35–38]. The choice of orthogonal polynomials as basis functions allows us to benefit from their advantages: important properties, especially as a Hilbert basis for L2. On the other hand, the choice of compact combinations of orthogonal polynomials, allows us not only to take the advantage of these orthogonal polynomials but also to reduce the dimension of the space of approximation from N+1 to N−n+1.


### 3.1. The Choice of Basis Functions

In this paper, we choose a spectral basis composed of shifted Legendre polynomials of the form(17)ϕkx=L˜kx+∑i=1nαikL˜k+ix,k=0,…,N−n.The coefficientsαik are calculated, such that ϕkx verifies the initial conditions (13), which means(18)ϕkj0=0,j=0,…,n−1,k=0,…,N−n.So,(19)L˜kj0+∑i=1nαikL˜k+ij0=0,j=0,…,n−1,k=0,…,N−n.Taking into account the important properties of shifted Legendre polynomials resulting from those of Legendre polynomials (1)–(3), (5), and (6), the basis coefficients are obtained from the following system:(20)∑i=1n−1iαik=−1,j=0,∑i=1nαik−1i∏m=0j−1k+i−mk+i+m+1=−∏m=0j−1k−mk+m+1,j=1,…,n−1.Remark 1.
Ifuj0=cj,cj≠0 for a certain j∈0,…,n−1, we move to homogeneous initial conditions by a suitable change of variable.The determinant of this system is different from zero, and we note here the following special cases.(i)
Whenn=1, we obtain α1k=1(ii)
Whenn=2, we obtain(21)α1k=2k+3k+2,α2k=k+1k+2.(iii)
Whenn=3, we obtain(22)α1k=32k+32k+5,α2k=3k+1k+3,α3k=k+12k+3k+32k+5.It is obvious thatϕkx are linearly independent and dimVN=N−n+1, so we obtain VN=spanϕkx,k=0,…,N−n and we denote our approximation as(23)uNx=∑k=0N−nukϕkx.Then, (14) becomes for j=0,…,N−n(24)a0x∑k=0N−nukϕkx,ϕjx+∑i=1naix∑k=0N−nukϕkxi,ϕjx=fx,ϕjx+λ∫01Kx,tFt,∑k=0N−nukϕktdt,ϕjx.It can be written as(25)∑k=0N−nuka0xϕkx,ϕjx+∑k=0N−nuk∑i=1naixϕkix,ϕjx=fx,ϕjx+λ∫01Kx,tFt,∑k=0N−nukϕktdt,ϕjx.We consider here two cases depending on the linearity ofFt,ut.(1)
WhenFt,ut is linear which means Ft,ut=ut,then (25) becomes, for all j=0,…,N−n,(26)∑k=0N−nuk<a0xϕkx,ϕjx>+∑k=0N−nuk∑i=1n<aixϕkix,ϕjx>=<fx,ϕjx>+λ∑k=0N−nuk<∫01Kx,tϕktdt,ϕjx>.Let us denote(27)f=f0,f1,…,fN−nT,fj=<fx,ϕjx>;U=u1,u2,…,uN−nT,uNx=∑k=0N−nukϕkx;A0=ajk00≤j,k≤N−n,ajk0=<a0xϕkx,ϕjx>;Ai,n=ajki,n0≤j,k≤N−n,ajki,n=<aixϕkix,ϕjx>;M=mjk0≤j,k≤N−n,mjk=<∫01Kx,tϕktdt,ϕjx>.So that (26) is equivalent to the following linear system:(28)A0+∑i=1nAi,n−λMU=f.(2)
WhenFt,ut is nonlinear of the form Ft,ut=uts,s≠0,1, then (25) becomes(29)∑k=0N−nuk<a0xϕkx,ϕjx>+∑k=0N−nuk∑i=1n<aixϕkix,ϕjx>=<fx,ϕjx>+λ<∫01Kx,t∑k=0N−nukϕktsdt,ϕjx>.We defineF as a vectorial function by(30)F=F0u0,…,uN−n,F1u0,…,uN−n,…,FN−nu0,…,uN−nT,where(31)Fju0,…,uN−n=∑k=0N−nuka0xϕkx,ϕjx>+∑i=1n<aixϕkix,ϕjx−<fx,ϕjx>−λ<∫01Kx,t∑k=0N−nukϕktsdt,ϕjx>.The Jacobian matrix ofF is given by(32)JFu0,…,uN−n=Jjk,0≤j,k≤N−nJjk=∂Fju0,…,uN−n∂uk,where(33)Jjk=<a0xϕkx,ϕjx>+∑i=1n<aixϕkix,ϕjx>−λ<∫01Kx,tsϕkt∑k=0N−nukϕkts−1dt,ϕjx>.


### 3.2. Gauss Quadrature

Let us denote the Legendre–Gauss–Lobatto nodes of−1,1 by xp with their respective weights wp. While moving to 0,1, we can define the shifted Legendre–Gauss–Lobatto (SLGL) nodes with their respective weights by [31](34)xp¯=12xp+1,wp¯=12wp.(i)
WhenFt,ut is linear, (26) becomes for all j=0,…,N−n,(35)∑k=0N−nuk∑p=0Nxa0xp¯ϕkxp¯ϕjxp¯wp¯+∑k=0N−nuk∑i=1n∑p=0Nxaixp¯ϕkixp¯ϕjxp¯wp¯=∑p=0Nxfxp¯ϕjxp¯wp¯+λ∑k=0N−nuk∑p=0Nx∑q=0NxKxp¯,xq¯ϕkxq¯ϕjxp¯wq¯wp¯.It leads to matrix system (28) with(36)fj=∑p=0Nxfxp¯ϕjxp¯wp¯;ajk0=∑p=0Nxa0xp¯ϕkxp¯ϕjxp¯wp¯;ajki,n=∑p=0Nxaixp¯ϕkixp¯ϕjxp¯wp¯;mjk=∑p=0Nx∑q=0NxKxp¯,xq¯ϕkxq¯ϕjxp¯wq¯wp¯.To solve (28) in both cases, we proceed by a method of Gauss elimination.(ii)
WhenFt,ut is nonlinear, then (31) becomes for all j=0,…,N−n,(37)Fju0,…,uN−n=∑k=0N−nuk∑p=0Nxa0xp¯ϕkxp¯ϕjxp¯wp¯+∑k=0N−nuk∑i=1n∑p=0Nxaixp¯ϕkixp¯ϕjxp¯wp¯−∑p=0Nxfxp¯ϕjxp¯wp¯−λ∑p=0Nx∑q=0NxKxp¯,xq¯∑k=0N−nukϕkxq¯sϕjxp¯wq¯wp¯,And the Jacobian matrix ofF is(38)JFu0,…,uN−n=Jjk,0≤j,k≤N−n,Jjk=∂Fju0,…,uN−n∂uk.Here,(39)Jjk=∑p=0Nxa0xp¯ϕkxp¯ϕjxp¯wp¯+∑i=1n∑p=0Nxaixp¯ϕkixp¯ϕjxp¯wp¯−λ∑p=0Nx∑q=0NxKxp¯,xq¯sϕkxq¯∑k=0N−nukϕkxq¯s−1ϕjxp¯wq¯wp¯.Now, to solve (31)–(33) and (37)–(39), we use the following algorithm of Newton.(i)
Initialisation: let beϵ>0, an initial vector,(40)U0=u00,u10,…,uN−n0,(41)F0=F0u00,…,uN−n0,F1u00,…,uN−n0,…,FN−nu00,…,uN−n0T.(ii)
Iteration: we solve(42)JFUk+1−Uk=−Fk.(iii)
StopUk+1−Uk<ϵ.


## 3.1. The Choice of Basis Functions

In this paper, we choose a spectral basis composed of shifted Legendre polynomials of the form(17)ϕkx=L˜kx+∑i=1nαikL˜k+ix,k=0,…,N−n.The coefficientsαik are calculated, such that ϕkx verifies the initial conditions (13), which means(18)ϕkj0=0,j=0,…,n−1,k=0,…,N−n.So,(19)L˜kj0+∑i=1nαikL˜k+ij0=0,j=0,…,n−1,k=0,…,N−n.Taking into account the important properties of shifted Legendre polynomials resulting from those of Legendre polynomials (1)–(3), (5), and (6), the basis coefficients are obtained from the following system:(20)∑i=1n−1iαik=−1,j=0,∑i=1nαik−1i∏m=0j−1k+i−mk+i+m+1=−∏m=0j−1k−mk+m+1,j=1,…,n−1.Remark 1.
Ifuj0=cj,cj≠0 for a certain j∈0,…,n−1, we move to homogeneous initial conditions by a suitable change of variable.The determinant of this system is different from zero, and we note here the following special cases.(i)
Whenn=1, we obtain α1k=1(ii)
Whenn=2, we obtain(21)α1k=2k+3k+2,α2k=k+1k+2.(iii)
Whenn=3, we obtain(22)α1k=32k+32k+5,α2k=3k+1k+3,α3k=k+12k+3k+32k+5.It is obvious thatϕkx are linearly independent and dimVN=N−n+1, so we obtain VN=spanϕkx,k=0,…,N−n and we denote our approximation as(23)uNx=∑k=0N−nukϕkx.Then, (14) becomes for j=0,…,N−n(24)a0x∑k=0N−nukϕkx,ϕjx+∑i=1naix∑k=0N−nukϕkxi,ϕjx=fx,ϕjx+λ∫01Kx,tFt,∑k=0N−nukϕktdt,ϕjx.It can be written as(25)∑k=0N−nuka0xϕkx,ϕjx+∑k=0N−nuk∑i=1naixϕkix,ϕjx=fx,ϕjx+λ∫01Kx,tFt,∑k=0N−nukϕktdt,ϕjx.We consider here two cases depending on the linearity ofFt,ut.(1)
WhenFt,ut is linear which means Ft,ut=ut,then (25) becomes, for all j=0,…,N−n,(26)∑k=0N−nuk<a0xϕkx,ϕjx>+∑k=0N−nuk∑i=1n<aixϕkix,ϕjx>=<fx,ϕjx>+λ∑k=0N−nuk<∫01Kx,tϕktdt,ϕjx>.Let us denote(27)f=f0,f1,…,fN−nT,fj=<fx,ϕjx>;U=u1,u2,…,uN−nT,uNx=∑k=0N−nukϕkx;A0=ajk00≤j,k≤N−n,ajk0=<a0xϕkx,ϕjx>;Ai,n=ajki,n0≤j,k≤N−n,ajki,n=<aixϕkix,ϕjx>;M=mjk0≤j,k≤N−n,mjk=<∫01Kx,tϕktdt,ϕjx>.So that (26) is equivalent to the following linear system:(28)A0+∑i=1nAi,n−λMU=f.(2)
WhenFt,ut is nonlinear of the form Ft,ut=uts,s≠0,1, then (25) becomes(29)∑k=0N−nuk<a0xϕkx,ϕjx>+∑k=0N−nuk∑i=1n<aixϕkix,ϕjx>=<fx,ϕjx>+λ<∫01Kx,t∑k=0N−nukϕktsdt,ϕjx>.We defineF as a vectorial function by(30)F=F0u0,…,uN−n,F1u0,…,uN−n,…,FN−nu0,…,uN−nT,where(31)Fju0,…,uN−n=∑k=0N−nuka0xϕkx,ϕjx>+∑i=1n<aixϕkix,ϕjx−<fx,ϕjx>−λ<∫01Kx,t∑k=0N−nukϕktsdt,ϕjx>.The Jacobian matrix ofF is given by(32)JFu0,…,uN−n=Jjk,0≤j,k≤N−nJjk=∂Fju0,…,uN−n∂uk,where(33)Jjk=<a0xϕkx,ϕjx>+∑i=1n<aixϕkix,ϕjx>−λ<∫01Kx,tsϕkt∑k=0N−nukϕkts−1dt,ϕjx>.


## 3.2. Gauss Quadrature

Let us denote the Legendre–Gauss–Lobatto nodes of−1,1 by xp with their respective weights wp. While moving to 0,1, we can define the shifted Legendre–Gauss–Lobatto (SLGL) nodes with their respective weights by [31](34)xp¯=12xp+1,wp¯=12wp.(i)
WhenFt,ut is linear, (26) becomes for all j=0,…,N−n,(35)∑k=0N−nuk∑p=0Nxa0xp¯ϕkxp¯ϕjxp¯wp¯+∑k=0N−nuk∑i=1n∑p=0Nxaixp¯ϕkixp¯ϕjxp¯wp¯=∑p=0Nxfxp¯ϕjxp¯wp¯+λ∑k=0N−nuk∑p=0Nx∑q=0NxKxp¯,xq¯ϕkxq¯ϕjxp¯wq¯wp¯.It leads to matrix system (28) with(36)fj=∑p=0Nxfxp¯ϕjxp¯wp¯;ajk0=∑p=0Nxa0xp¯ϕkxp¯ϕjxp¯wp¯;ajki,n=∑p=0Nxaixp¯ϕkixp¯ϕjxp¯wp¯;mjk=∑p=0Nx∑q=0NxKxp¯,xq¯ϕkxq¯ϕjxp¯wq¯wp¯.To solve (28) in both cases, we proceed by a method of Gauss elimination.(ii)
WhenFt,ut is nonlinear, then (31) becomes for all j=0,…,N−n,(37)Fju0,…,uN−n=∑k=0N−nuk∑p=0Nxa0xp¯ϕkxp¯ϕjxp¯wp¯+∑k=0N−nuk∑i=1n∑p=0Nxaixp¯ϕkixp¯ϕjxp¯wp¯−∑p=0Nxfxp¯ϕjxp¯wp¯−λ∑p=0Nx∑q=0NxKxp¯,xq¯∑k=0N−nukϕkxq¯sϕjxp¯wq¯wp¯,And the Jacobian matrix ofF is(38)JFu0,…,uN−n=Jjk,0≤j,k≤N−n,Jjk=∂Fju0,…,uN−n∂uk.Here,(39)Jjk=∑p=0Nxa0xp¯ϕkxp¯ϕjxp¯wp¯+∑i=1n∑p=0Nxaixp¯ϕkixp¯ϕjxp¯wp¯−λ∑p=0Nx∑q=0NxKxp¯,xq¯sϕkxq¯∑k=0N−nukϕkxq¯s−1ϕjxp¯wq¯wp¯.Now, to solve (31)–(33) and (37)–(39), we use the following algorithm of Newton.(i)
Initialisation: let beϵ>0, an initial vector,(40)U0=u00,u10,…,uN−n0,(41)F0=F0u00,…,uN−n0,F1u00,…,uN−n0,…,FN−nu00,…,uN−n0T.(ii)
Iteration: we solve(42)JFUk+1−Uk=−Fk.(iii)
StopUk+1−Uk<ϵ.


## 4. Numerical Examples

In this section, several examples are presented to be discussed in order to illustrate the effectiveness of the described method. To compare between the exact solutionux and the approximate solution uNx, we define the absolute error by(43)Ex=ux−uNx,the maximum absolute error by(44)Emax=maxux−uNx,x∈0,1,the square error norm by(45)E2=h∑i=0Nxuxi−uNxi21/2,where Nx is the number of interval subdivisions and h is the appropriate step, and the relative error by(46)Ere=E2ux2,ux2=h∑i=0Nxuxi21/2.To prove the exponential convergence of the method, the Log slope coefficient according to relative errorO is calculated using the following formula:(47)O=logEreN1−logEreN2logN1−logN2,where EreN1 (respectively, EreN2) denotes the relative error calculated for N=N1 (resp. N=N2) and N2=N1+2,N1≥2.Example 1.
Consider the following nonlinear integral equation:(48)ux=−x2−x322−1+2+∫01xtutdt.
The exact solution isux=2−x2.
Error values were calculated for differentN≥2, and we obtain Ex=0, for all x in 0,1.Example 2.
Consider the following nonlinear integral equation:(49)ux=x+1−e−1sinhξx−1+∫01sinhξx−1cosht−1u2tdt,where ξ=2. The exact solution is ux=x+1.Example 3.
Let us consider the following linear first-order Fredholm integrodifferential equation:(50)u′x=ex−x+xex+∫01xutdt,with the initial condition u0=0. The exact solution is ux=xex.Example 4.
Consider the following nonlinear first-order Fredholm integrodifferential equation of the form(51)xu′x−ux=−16+45x2+∫01x2+tu2tdt,with the initial condition u0=0. The exact solution is ux=x2.Example 5.
Consider the following nonlinear first-order Fredholm integrodifferential equation of the form(52)u′x+ux=12e−2−1+∫01u2tdt,with the initial condition u0=1. The exact solution is ux=e−x.Example 6.
Consider the following nonlinear second-order Fredholm integrodifferential equation of the form(53)u″x+xu′x−xux=ex−sinx+∫01sinxe−2tu2tdt,with the initial condition u0=u′0=1. The exact solution is ux=ex.


### 4.1. Results and Discussion

By applying the described method in this study, we get systems of linear and nonlinear equations that can be solved directly by Gauss elimination or by Newton’s algorithm, respectively. Several results are presented in tables and figures to confirm the accuracy of our method. All computations were carried out by MATLAB R2020a on an AMD Ryzen 5 5600X 6-Core Processor Desktop.All results are obtained forNx=50 and CPU time is given in seconds. For Newton’s method, ϵ is taken 10−10 and the maximum number of iterations is 5.Tables1–3 expose different values of εmax, ε2, and εre, allowing us to compare, for each example, the exact solution and the approximate one established via the proposed method, and this is for different increasing values of N confirming that the results are very accurate and precise for small values of N. Also, they present the coefficient O that confirms the exponential convergence of the solution according to the relative error.Table 1
Emax, E2, Ere, and O for different N and respective CPU time for Example 2.
NEmaxE2EreOCPU time21.30657⋅10−25.52993⋅10−34.46979⋅10−3/10.04241.69372⋅10−45.43398⋅10−54.39223⋅10−5−6.669124.34663.20565⋅10−68.93187⋅10−77.21954⋅10−7−10.13249.40386.95769⋅10−81.79229⋅10−81.44869⋅10−8−13.586116.59101.62236⋅10−93.99763⋅10−103.23124⋅10−10−17.042210.63123.95414⋅10−119.46293⋅10−127.64880⋅10−12−20.532356.24149.93205⋅10−132.30971⋅10−131.86691⋅10−13−24.085570.10162.55351⋅10−145.71509⋅10−154.61945⋅10−15−27.702885.40188.88178⋅10−162.00619⋅10−161.62158⋅10−16−28.4371329.3Table 2
E2 and Ere for different N and CPU time for Example 4.
NE2EreCPU time26.86325⋅10−141.49646⋅10−134.250249.58120⋅10−152.08908⋅10−149.938364.43580⋅10−179.67181⋅10−1718.02585.50628⋅10−171.20059⋅10−1629.195Table 3
Emax, E2, Ere, and O for different N and respective CPU time for Example 6.
NEmaxE2EreOCPU time23.36141⋅10−11.73884⋅10−15.59325⋅10−1—3.313844.66519⋅10−32.66201⋅10−38.56278⋅10−3−6.0299.097561.72096⋅10−59.92601⋅10−63.19284⋅10−5−13.79022.37582.86663⋅10−81.65973⋅10−85.33876⋅10−8−22.22267.270102.72220⋅10−111.57900⋅10−105.07908⋅10−10−20.861172.55121.66533⋅10−149.68711⋅10−153.11600⋅10−14−53.196338.74146.66133⋅10−162.59410⋅10−168.34430⋅10−16−23.484581.84Tables4 and 5 present the values of the exact solution and the approximation for different nodes on 0,1, and the absolute error is also calculated for each node.Table 4
Exact and approximate solution values in different nodes on0,1 when N=7 for Example 3.
xExact solutionApproximate solutionEx0.22.442805e−012.442806e−017.099942e−070.45.697298e−015.697299e−018.840336e−080.61.093271e+001.093271e+008.013137e−080.81.780432e+001.780432e+001.056572e−0712.718281e+002.718281e+001.260678e−07Table 5
Exact and approximate solution values in different nodes on0,1 when N=10 for Example 5.
xExact solutionApproximate solutionEx0.1−9.51625e−02−9.51625e−026.39766e−150.2−1.81269e−01−1.81269e−017.66053e−150.4−3.29679e−01−3.29679e−013.79141e−140.5−3.93469e−01−3.93469e−015.01265e−140.7−5.03414e−01−5.03414e−017.89368e−140.8−5.50671e−01−5.50671e−019.38138e−14Figure1 presents the curve of the absolute error Ex when N=18. The figure is homogeneous with the results obtained in Table 1. Figure 2 shows the curve of the maximum absolute error calculated for different N, while Figure 3 compares different values of absolute error calculated for different nodes on 0,1 when N=4,8,16. The error is of 10−16 order for N=16. The results confirm the reliability of the method specially when taking high values of N.Figure 1
Absolute error curve forN=18 for Example 2.Figure 2
Logarithmic maximum absolute error curve for different values ofN for Example 5.Figure 3
Logarithmic absolute error curve for different values ofN for Example 6.


## 4.1. Results and Discussion

By applying the described method in this study, we get systems of linear and nonlinear equations that can be solved directly by Gauss elimination or by Newton’s algorithm, respectively. Several results are presented in tables and figures to confirm the accuracy of our method. All computations were carried out by MATLAB R2020a on an AMD Ryzen 5 5600X 6-Core Processor Desktop.All results are obtained forNx=50 and CPU time is given in seconds. For Newton’s method, ϵ is taken 10−10 and the maximum number of iterations is 5.Tables1–3 expose different values of εmax, ε2, and εre, allowing us to compare, for each example, the exact solution and the approximate one established via the proposed method, and this is for different increasing values of N confirming that the results are very accurate and precise for small values of N. Also, they present the coefficient O that confirms the exponential convergence of the solution according to the relative error.Table 1
Emax, E2, Ere, and O for different N and respective CPU time for Example 2.
NEmaxE2EreOCPU time21.30657⋅10−25.52993⋅10−34.46979⋅10−3/10.04241.69372⋅10−45.43398⋅10−54.39223⋅10−5−6.669124.34663.20565⋅10−68.93187⋅10−77.21954⋅10−7−10.13249.40386.95769⋅10−81.79229⋅10−81.44869⋅10−8−13.586116.59101.62236⋅10−93.99763⋅10−103.23124⋅10−10−17.042210.63123.95414⋅10−119.46293⋅10−127.64880⋅10−12−20.532356.24149.93205⋅10−132.30971⋅10−131.86691⋅10−13−24.085570.10162.55351⋅10−145.71509⋅10−154.61945⋅10−15−27.702885.40188.88178⋅10−162.00619⋅10−161.62158⋅10−16−28.4371329.3Table 2
E2 and Ere for different N and CPU time for Example 4.
NE2EreCPU time26.86325⋅10−141.49646⋅10−134.250249.58120⋅10−152.08908⋅10−149.938364.43580⋅10−179.67181⋅10−1718.02585.50628⋅10−171.20059⋅10−1629.195Table 3
Emax, E2, Ere, and O for different N and respective CPU time for Example 6.
NEmaxE2EreOCPU time23.36141⋅10−11.73884⋅10−15.59325⋅10−1—3.313844.66519⋅10−32.66201⋅10−38.56278⋅10−3−6.0299.097561.72096⋅10−59.92601⋅10−63.19284⋅10−5−13.79022.37582.86663⋅10−81.65973⋅10−85.33876⋅10−8−22.22267.270102.72220⋅10−111.57900⋅10−105.07908⋅10−10−20.861172.55121.66533⋅10−149.68711⋅10−153.11600⋅10−14−53.196338.74146.66133⋅10−162.59410⋅10−168.34430⋅10−16−23.484581.84Tables4 and 5 present the values of the exact solution and the approximation for different nodes on 0,1, and the absolute error is also calculated for each node.Table 4
Exact and approximate solution values in different nodes on0,1 when N=7 for Example 3.
xExact solutionApproximate solutionEx0.22.442805e−012.442806e−017.099942e−070.45.697298e−015.697299e−018.840336e−080.61.093271e+001.093271e+008.013137e−080.81.780432e+001.780432e+001.056572e−0712.718281e+002.718281e+001.260678e−07Table 5
Exact and approximate solution values in different nodes on0,1 when N=10 for Example 5.
xExact solutionApproximate solutionEx0.1−9.51625e−02−9.51625e−026.39766e−150.2−1.81269e−01−1.81269e−017.66053e−150.4−3.29679e−01−3.29679e−013.79141e−140.5−3.93469e−01−3.93469e−015.01265e−140.7−5.03414e−01−5.03414e−017.89368e−140.8−5.50671e−01−5.50671e−019.38138e−14Figure1 presents the curve of the absolute error Ex when N=18. The figure is homogeneous with the results obtained in Table 1. Figure 2 shows the curve of the maximum absolute error calculated for different N, while Figure 3 compares different values of absolute error calculated for different nodes on 0,1 when N=4,8,16. The error is of 10−16 order for N=16. The results confirm the reliability of the method specially when taking high values of N.Figure 1
Absolute error curve forN=18 for Example 2.Figure 2
Logarithmic maximum absolute error curve for different values ofN for Example 5.Figure 3
Logarithmic absolute error curve for different values ofN for Example 6.


## 5. Conclusion

In this study, the shifted Legendre spectral method is used to find approximate solutions for both linear and nonlinear integral equations and Fredholm integrodifferential equations ofn-th order with the initial condition. Through important properties of shifted Legendre polynomials, we build a new algorithm simple to use and effective for accurate results, which has been proven through several examples. Error calculations (absolute error, maximum absolute error, square error norm, and relative error) confirm for each example the reliability of the method and the exponential convergence in terms of relative error. The use of compact combinations of orthogonal polynomials as basis functions of spectral decomposition leads to simple systems with simple matrices to invert. In addition, the collocation treatment by evaluating the obtained integrals from the inner product using Gauss quadrature on shifted Legendre–Gauss–Lobatto nodes with the respective weights enables us to get more precise calculations which reduce the cost of the method and makes it accessible and simple to use for many different kinds of equations.



---

*Source: 9043428-2022-10-08.xml*