# A Triangle Algorithm of Padé-Type Approximant for Two-Dimensional Fredholm Integral Equations of the Second Kind

Authors: Meilan Sun; Chuanqing Gu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The function-valued Padé-type approximation (2DFPTA) is used to solve two-dimensional Fredholm integral equation of the second kind. In order to compute 2DFPTA, a triangle recursive algorithm based on Sylvester identity is proposed. The advantage of this algorithm is that, in the process of calculating 2DFPTA to avoid the calculation of the determinant, it can start from the initial value, from low to high order, and gradually proceeds. Compared with the original two methods, the numerical examples show that the algorithm is effective.


---

## Body

## 1. Introduction

Consider a two-dimensional Fredholm integral equation of the second kind (2DF-II):(1)fx,λ=gx+λ∫DKx,yfy,λdy,x,y∈D,where D is a bounded domain, x=(x1,x2),y=(y1,y2),dy=dy1dy2, K and f are given functions defined on D, and λ is a real parameter called a characteristic value. It is known that a sequence of approximations leading to an infinite power series can be yielded by successive substitution, which is an iterative procedure, from this kind of integral equation. Set f(x,λ) to be analytic and meromorphic in a neighborhood D~ of the origin λ=0. Let its Neumann series be given by(2)fx,λ=∑i=0∞cixλi,in which ci(x) are continuous on D and c0(x)=g(x), (3)cix=∫DKx,yci-1ydy,fori=1,2,….Let ci(x),cj(x)∈L2[a,b] and the inner product and the norm be defined by(4)ci,cj=∫Dcixcjxdx,(5)ci=ci,ci1/2=∫Dci2xdx1/2.Two-dimensional Fredholm integral equation of the second kind (2DF-II) is a useful tool to model many problems arising in fracture mechanics, transportation [1], 2D electromagnetic scattering [2], and computer graphics manipulations [3]. For domain D, some works treat the square case [4–10] and some are based on discrete Galerkin method [5], Monte Carlo methods [6], or piecewise approximating polynomials [7], Nytröm methods based on cubature rules obtained as the tensor product of two-variate Gaussian rules [8, 10], and meshless method with complex factors [11]. Rational approximation has extensive application in engineering, technology and calculation [12–18]. As efficient rational approximating techniques, a generalized inverse function-valued Padé approximants had been used in solving integral equation [13, 14], and a function-valued Padé-type approximation method was used to solve one and two-dimensional Fredholm integral equation of the second kind [15–18].In the paper, the function-valued Padé-type approximation (2DFPTA) is used to solve 2DF-II. In order to compute 2DFPTA, a recursive algorithm based on Sylvester identity is proposed. The remainder of this paper is organized as follows. In Section2, we give the definition of 2DFPTA [18]. In Section 3, we apply Sylvester identity to propose a recursive algorithm to compute 2DFPTA. In Section 4, we compare our algorithm with the two algorithms in [4, 18].


## 2. Function-Valued Padé-Type Approximants and Its Convergence Analysis

In this section, we will give the definition of the function-valued Padé-type approximation [15–18]. Let ϕ(l):P→C be a linear functional on the polynomial space P, and define it by(6)ϕlti=cl+ix1,x2,i=0,1,…,l∈Z,where ϕ(l)(ti)=0 for l+i<0. We can obtain from linear functional ϕ=ϕ(0) in (6) that (7)ϕ1-tλ-1=ϕ1+tλ+tλ2+⋯=∑i=0∞cix1,x2λi=fx,λ.Let(8)vnλ=b0+b1λ+⋯+bnλn,where vn∈Pn is a scalar polynomial of degree n and assume bn≠0.Define the polynomialw with function-valued coefficient by(9)wmx,λ=ϕtm-n+1vnt-λm-n+1vnλt-λ.Note that ϕ acts on P and wm(x,λ) is a function-valued polynomial in λ of degree m.Set(10)v~nλ=λnvnλ-1=∑j=0nbjλn-j.(11)w~mx,λ=λmwmx,λ-1.The polynomialvn is called the generating polynomial of the Padé-type approximation (m/n). It can be arbitrarily chosen and we have n degrees of freedom.It is known from (10) that bn≠0 implies v~n(0)≠0. In view of (8), (9), and (11), w~m(x,λ) is of the following form:(12)w~mx,λ=∑j=0nbjλn-j∑i=0m-n+jcixλi=∑i=0m∑j=0nbjc~ijxλi,in which c~ij(x)=ci-n+j(x). Using (10) and (12), we get(13)v~nλfx,λ-w~mx,λ=∑j=0nbjλn-j∑i=0∞cixλi-∑j=0nbjλn-j∑i=0m-n+jcixλi.=∑i=0∞∑j=0ncm-n+1+i+jxλm+1+i=∑i=m+1∞∑j=0nbjc~ijxλi=Oλm+1.Definition 1.
Letw~m(x,λ) and v~n(λ) be given by (12) and (10), respectively, and then the rational function(14)Rm,nx,λ=w~mx,λv~nλ=∑j=0nbjλn-j∑i=0m-n+jcixλi∑j=0nbjλn-j=∑i=0m∑j=0nbjc~ijxλi∑j=0nbjλn-jis defined as a two-dimensional function-valued Padé-type approximation (2DFPTA) with generating polynomial vn(λ) and is denoted by (m/n)f(x,λ).
According to the relations (9) and (11), we can from [19] obtain the following error formula with the linear functional form:(15)fx,λ-mnfx,λ=λm+1vn~λϕm-n+1vnt1-tλ=λm+1v~nλϕm-n+1vnt1+tλ+⋯+tn-1λn-1+tnλn1-tλ.
Note that the error formula (15) is disadvantageous to estimate error numerically, and it implies from (13) that the Padé-type approximation property is(16)fx,λ-mnfx,λ=Oλm+1,asλ→0forv~nλ≠0.
In order to compute conveniently, we choosemaxx∈D⁡∑j=0nbjc~ij(x) as the maximum absolute value of the coefficient of λm+1 on D in (13). We can also observe from (15) and (16) that if the coefficients are given, the order of the approximation is λm+1.


## 3. A Triangle Recursive Algorithm of(m/n)f(x,λ)

As the definition given in Section2, we find that the key to calculating 2DFPTA is to compute its the generating polynomial vn(λ). If the generator polynomial vn(λ) is determined, we can compute the 2DFPTA of type (m,n) according to (10) and (12). In the section, we will apply Sylvester’s identity [20] to propose a three-term recurrence formula for computing v~n(λ) and establish a complete recursive algorithm for calculating 2DFPTA.Lemma 2 (Sylvester identity [20]).
Let(17)A=a0,1a0,2⋯a0,na1,1a1,2⋯a1,n⋮⋮⋱⋮an-1,1an-1,2⋯an-1,n,B=a0,1a0,2⋯a0,n-1a1,1a1,2⋯a1,n-1⋮⋮⋱⋮an-2,1an-2,2⋯an-2,n-1,and then it holds that(18)x0xuAB=Ax0x′u′B-A′x′xnBu′′,where A′=u′Ban-1,0w and(19)x=x1,x2,…,xn,x′=x1,x2,…,xn-1,u=a0,0,a1,0,…,an-1,0T,u′=a0,0,a1,0,…,an-2,0,u′′=a0,n,a1,n,…,an-2,nT,w=an-1,1,an-1,2,…,an-1,n-1.Set b0=1, the generating polynomial vn(λ)=b0+b1λ+⋯+bnλn.Lemma 3 (see [15]).
Let the determinantAn≠0; then the generating polynomial vn(λ) exists and is unique and can be expressed as(20)vnλ=1λ→nη→nAnAn,where(21)An=cm-n+1,cm-n+2cm-n+1,cm-n+3⋯cm-n+1,cm+1cm-n+1,cm-n+3cm-n+1,cm-n+4⋯cm-n+1,cm+2⋮⋮⋱⋮cm-n+1,cm+1cm-n+1,cm+2⋯cm-n+1,cm+n,λ→n=λ,λ2,…,λn,η→n=cm-n+1,cm-n+1,cm-n+1,cm-n+2,…,cm-n+1,cmT.Definition 4.
For anyjth degree polynomial vj(λ)=b0+b1λ+⋯+bjλj, by the higher-order linear functional ϕ(l)(λj)=cl+j, we define the vector inner product with l order as follows:(22)cl,vjλl=cl,ϕlvjλ=cl,ϕlb0+b1λ+⋯+bjλj=cl,b0cl+b1cl+1+⋯+bjcl+j=b0cl,cl+b1cl,cl+1+⋯+bjcl,cl+j,where l=m-n+1,j=0,1,2,….
To discuss conveniently, we introduce the following notations:(23)Gki=λiλi+1⋯λi+kcl,λilcl,λi+1l⋯cl,λi+kl⋮⋮⋯⋮cl,λi+k-1lcl,λi+kl⋯cl,λi+2k-1l,(24)Dki=cl,λi+1lcl,λi+2l⋯cl,λi+klcl,λi+2lcl,λi+3l⋯cl,λi+k+1l⋮⋮⋯⋮cl,λi+klcl,λi+k+1l⋯cl,λi+2k-1l.Let(25)vki=GkiDki,where if i=0,k=n and(26)Gn0=1λ→nη→nAn,Dn0=An.From (20) and (25), we conclude that(27)vn0=vnλ.To computevn(0), we give the following theorem. According to Definition 4 and the jth degree polynomial vk-1(i)(λ)=b0+b1λ+⋯+bjλj, we gain the vector inner product with l+k-1 order as follows:(28)cl,vk-1iλl+k-1=cl,ϕl+k-1vk-1iλ=b0cl,cl+k-1+b1cl,cl+k+⋯+bjcl,cl+k-1+j,where l=m-n+1,k=1,2,…,n,i=j-(k-1).Theorem 5 (three-term recurrence formula).
LetDk(i)Dk-1(i)≠0 for all k=1,2,…,n,i=0,1,2,…,n. Then it holds that(29)vki=vk-1i-cl,ϕl+k-1vk-1icl,ϕl+k-1vk-1i+1vk-1i+1,where the initial value is v0(i)=λi.Proof.
Setx0=λi,x=(λi+1,λi+2,…,λi+k),u=([c1,λi]1,[c2,λi]2,…,[ck,λi]k) and A=Dk(i),|B|=Dk-1(i). Substitute them into the formula (18) and we conclude that(30)GkiDk-1i=Gk-1iDki-Gk-1i+1Dki-1.Due to Dk(i)Dk-1(i)≠0, we have (31)GkiDki=Gk-1iDk-1i-Gk-1i+1Dki-1Dk-1iDki.By (25) we get Gk-1(i+1)=Dk-1(i+1)vk-1(i+1) and substitute it into the above formula; then (32)GkiDki=Gk-1iDk-1i-Dki-1Dk-1i+1vk-1i+1Dk-1iDki,that is,(33)vki=vk-1i-Dki-1Dk-1i+1Dk-1iDkivk-1i+1.According to the higher-order linear functional and nature of the determinant, we have(34)cl,ϕl+k-1vk-1i=-1k-1Dki-1Dk-1i,and(35)cl,ϕl+k-1vk-1i+1=-1k-1DkiDk-1i+1.Dividing (34) by (35), we deduce that(36)cl,ϕl+k-1vk-1icl,ϕl+k-1vk-1i+1=Dki-1Dk-1i+1DkiDk-1i.From (33) and (36), we obtain that(37)vki=vk-1i-cl,ϕl+k-1vk-1icl,ϕl+k-1vk-1i+1vk-1i+1.Example 6.
Let us choose an example ofn=6 to illustrate the process to generate v6(λ). By the three-term recurrence relation (16), the calculation program of v6(λ) can be arranged according to the following bottom-down triangle: (38)v00↘v01→v10↘↘v02→v11→v20↘↘↘v03→v12→v21→v30↘↘↘↘v04→v13→v22→v31→v40↘↘↘↘↘v05→v14→v23→v32→v41→v50↘↘↘↘↘↘v06→v15→v24→v33→v42→v51→v60
When computingvk(i) we mainly use the left polynomial vk-1(i+1), the polynomial vk-1(i) in the top left corner in the table, and their relation (29). It is not difficult to find that, in every column, the number of polynomial is in decline when k is bigger and bigger.
It is noticed thatvk(i) can also be arranged by another pattern as follows:(39)1λλ2λ3λ4λ5λ6↓↙↓↙↓↙↓↙↓↙↓↙v10v11v12v13v14v15↓↙↓↙↓↙↓↙↓↙v20v21v22v23v24↓↙↓↙↓↙↓↙v30v31v32v33↓↙↓↙↓↙v40v41v42↓↙↓↙v50v51↓↙v60When computing vk(i) we mainly use the above polynomial vk-1(i), the polynomial vk-1(i+1) in the top-right corner in the table, and their relation (29). Similarly, the number of polynomials is in decline when k is bigger and bigger in every row. In this way, we can recursively get v6(λ)=v6(0).
From the three-term recurrence relation (29), we now build a complete recursive algorithm for calculating 2DFPTA. The advantage of this algorithm is that in the process of calculating 2DFPTA to avoid the calculation of the determinant, it can start from the initial value, from low to high order, and is gradually calculated. And its specific computational steps are as follows.Algorithm 7 (Sylvester-type triangle algorithm).
(1)
Set the initial value asv0(i)=λi,i=0,1,…,n compute (ci,cj) by (4). (2)
Fork=1,2,…,n-1 compute by (29) in turn vk(0),vk(1),…,vk(n-k). (3)
Computevn(λ)=vn(0)=vn-1(0)-(cl,ϕ(l+k-1)(vn-1(0)))/(cl,ϕ(l+k-1)(vn-1(1)))vn-1(1). (4)
v n ( λ ) is written as vn(λ)=b0+b1λ+⋯+bnλn(bn=1) and the denominator polynomial of function-valued Padé-type approximation is(40)v~nλ=λnvnλ-1=∑i=0nbiλn-i. (5)
From (12), numerator polynomial of (m/n)f(x,λ) is(41)w~mx,λ=∑j=0nbjλn-j∑i=0m-n+jcixλi=∑i=0m∑j=0nbjc~ijxλi. (6)
In view of Definition1, there exists (m/n)f(x,λ) such that(42)mnfx,λ=w~mx,λv~nλ.Example 8.
Consider the second Fredholm integral equation(43)fx1,x2,λ=x1cos⁡x2+43sin⁡x1-1+43sin⁡1x2+λ∫-11∫-11y2sin⁡x1+x2y1fy1,y2,λdy1dy2,x1,x2∈-1,1.Its kernel isK(x,t)=(y2sin⁡(x1)+x2y1).In the following we use the triangle algorithm to yield(2/2)f(x,λ) of the integral equation.Solution. The first few terms of the power series for f(x,λ) are given by (44)fx,λ=c0x+c1xλ+c2xλ2+c3xλ3+c4xλ4+⋯=x1cos⁡x2+43sin⁡x1-1+43sin⁡1x2+-2.829281750769593sin⁡x1+2.728194267422565x2λ+3.637592356563420sin⁡x1-3.408364188910562x2λ2+-4.544485585214082sin⁡x1+4.382115538190247x2λ3+5.842820717586996sin⁡x1-5.474626880638770x2λ4+⋯.Here m=n=1, l=m-n+1=1 and we set d=l.Step 1.
Set the initial value asv0(0)=1, v0(1)=λ, v0(2)=λ2, computing (45)c1,c1=18.654952992399824,c1,c2=-23.623503547631628.Step 2.
Computev1(0)andv1(1) by (29)(46)v10=v00-c1,ϕ1v00c1,ϕ1v01v01=1-c1,c1c1,c2λ=1+0.789677659572644λ.v11=v01-c1,ϕ1v01c1,ϕ1v02v02=λ-c1,c2c1,c3λ2=λ+0.788390924640188λ2.Step 3.
Obtain the generating polynomial off(x,λ) function-value Padé-type approximant by (27) (47)v20=v10-c1,ϕ2v10c1,ϕ2v11v11=1+1.1752×10-13λ-0.622574700198084λ2.Step 4.
From (10) get the denominator of f(x,λ) function-value Padé-type approximant(48)v~2λ=λ2v2λ-1=b0λ2+b1λ+b2=-1.606232954345608λ2-1.89×10-13λ+1Step 5.
Get the numerator off(x,λ) function-value Padé-type approximant by (12) (49)w~2x,λ=b0c0+b1c1+b2c2λ2+b1c0+b2c1λ+b2c0=a2λ2+a1λ+a0where (50)a2=1.495948417436477sin⁡x1-9.0779×10-16x2-1.606232954345608x1cos⁡x2a1=-2.829281750769845sin⁡x1+2.728194267422966x2-1.89×10-13x1cos⁡x2a0=1.3333333333333333sin⁡x1-2.121961313077195x2+x1cos⁡x2Step 6.
Obtainf(x,λ) function-valued Padé-type approximant(51)22fx,λ=a2λ2+a1λ+a0b0λ2+b1λ+b2.


## 4. Numerical Results

In this section, we will consider the square domain case ofD=[α,β]×[α,β]. We apply Algorithm 7 for 2DFPTA to approximate the corresponding 2DF-II [11]. Compared with the methods in [4, 18], numerical experiments show that our method improves the calculation accuracy. Furthermore, our method can estimate the corresponding characteristic values of f(x,λ) from the real roots of v~n(λ).Example 9.
Consider the same Fredholm integral equation with Example8 [4, 18].
Whenλ=1, the exact solution of (43) is (52)fx1,x2,λ=x1cos⁡x2-x2.By Algorithm 7, we get the following calculation results:(53)v~2λ=1-1.89×10-13λ-1.606232954345608λ2,22fx1,x2,1=-x2+x1cos⁡x2+5.8593×10-14sin⁡x1+5.6×10-14x2,where when λ=1 the maximum absolute error of (2/2)f(x1,x2,1) caused by Algorithm 7 is shown in Figure 1. From Figure 1, we can observe that the error is 1.05×10-13 and in the process we just only employ the first m+n=4 items of Neumann series (4), while in [4] the error was 8.7×10-5 when the first 5 items of Chebyshev series corresponding to this equation were used and in [18] the corresponding error was 5.3×10-13.
By Algorithm7, we also get that(54)v~3λ=1+0.999999999999785λ-1.606232954345953λ2-1.606232954345763λ3,v~4λ=1-7.837686719857387λ-2.094833520416517λ2+12.589150695272894λ3+0.784806330738950λ4,where their eigenvalue has estimates; see Table 1.Table 1
Estimates of the characteristic values of(n/n)f(x,λ),n=1,2,3,4 in Example 9.
( n / n ) f ( 1 / 1 ) f ( 2 / 2 ) f (3/3)f (4/4)f λ1 0.789677659572644 0.789034029809863 0.789034029809820 0.789034029809775 λ 2  -0.789034029809980 -0.789034029809980 -0.789034029809980 λ3   -0.999999999999958 -16.167681854028405 λ4    0.126589673465619Figure 1
λ = 1, error of (2/2)f(x,λ) in Example 9.In the computational process ofv2~(λ), we get two possible characteristic values λ1=0.7890340298098 and λ2=-0.789034029809980. After further computing we can deduce that v~3(λ) and v~4(λ) both include the two real roots. Therefore we consider the two numbers as characteristic values of this equation.In Table2, (2/2)f(x,λ) is from Algorithm 7 in the paper and (2/2)u(x,λ) is from the method in [18], where (55)22ux1,x2,1=-x2+x1cos⁡x2+2.94875×10-13sin⁡x1.We find that from the numerical results (2/2)f(x,λ) has better approximation effect compared with (2/2)u(x,λ) in [18].Table 2
Whenλ=1, the errors of (2/2)f(x,λ) and (2/2)u(x,λ) at corresponding points in Example 9.
x f ( x , λ ) f x , λ - 2 / 2 u x , λ f x , λ - 2 / 2 f x , λ (-0.1,0.1) -1.9950e-001 2.9438e-014 2.4951e-016 (-0.2,0.2) -3.9601e-001 5.8583e-014 4.4057e-016 (-0.3,0.3) -5.8660e-001 8.7142e-014 5.1532e-016 (-0.4,0.4) -7.6842e-001 1.1483e-013 4.1707e-016 (-0.5,0.5) -9.3879e-001 1.4137e-013 9.0830e-017 (0.6,-0.6) 1.0952e+000 1.6650e-013 5.1608e-016 (0.7,-0.7) 1.2354e+000 1.8996e-013 1.4536e-015 (0.8,-0.8) 1.3574e+000 2.1153e-013 2.7682e-015 (0.9,-0.9) 1.4594e+000 2.3098e-013 4.5028e-015 (1.0,-1.0) 1.5403e+000 2.4813e-013 6.6960e-015Example 10.
Consider the following [4, 18]:(56)fx1,x2,λ=x1e-x2+13e-1-712x1-12x2+λ∫01∫01x1y1+x2ey2fy1,y2,λdy1dy2,x1,x2∈0,1.When λ=1, the exact solution of (56) is (57)fx1,x2,λ=x1e-x2+x2.By Algorithm 7, the calculation of Example 10 is shown as follows:(58)v~2λ=1-1.333333333334607λ+0.118548104777520λ2,22fx1,x2,1=-3.801972000573112×10-15x1+1.000000000000003x2+x1e-x2,where when λ=1 the maximum absolute error of (2/2)f(x1,x2,1) caused by Algorithm 7 is shown in Figure 2. From Figure 2, we can observe that the error is 3.7×10-15 and in the process we just only employ the first m+n=4 items of Neumann series (4). While in [4] the error was 8.8×10-7 when the first 5 items of Chebyshev series corresponding to this equation were used and in [18] the corresponding error was 5.5×10-15.Figure 2
λ = 1, error of (2/2)f(x,λ) in Example 10.Fromv2~(λ) in (58), we get two possible characteristic values λ1=0.80805470231918,λ2=10.43913761554. By Algorithm 7, we also get that(59)v~3λ=1-0.428436108131838λ-1.087981528825263λ2+0.107273851063696λ3,v~4λ=1-1.225739535044527λ-0.111597656991711λ2+0.128338192937178λ3-0.010276622714483λ4.Note that both v~3(λ) and v~4(λ) include the two real roots. Therefore, we consider the two real roots as the characteristic values of (56). For estimates of the characteristic values of (n/n)f(x,λ),n=1,2,3,4 in Example 10, see Table 3.Table 3
Estimates of the characteristic values of(n/n)f(x,λ),n=1,2,3,4 in Example 10.
n / n f 1 / 1 f 2/2f 3 / 3 f 4/4f λ 1 0.815393639551284 0.808054702319182 0.808054702319182 0.808054702319181 λ2  10.439137615545988 10.439137615545846 10.439137615544714 λ 3   0.017857142857281 -0.0476584+0.4291069i λ 4    -0.0476584-0.4291069iIn Table4, (2/2)f(x,λ) is from Algorithm 7 in the paper and (2/2)u(x,λ) is from the method in [18], where (60)22ux1,x2,1=1.000000000000009x2+1.000000000000001x1e-x2.From the numerical results in Table 4, (2/2)f(x,λ) has better approximation effect compared with (2/2)u(x,λ) in [18].Table 4
Whenλ=1, the errors of (2/2)f(x,λ) and (2/2)u(x,λ) at corresponding points in Example 10.
x f ( x , λ ) f(x,λ)-2/2u(x,λ) f x , λ - 2 / 2 f x , λ (0.1,0.1) 1.9048e-001 9.9048e-016 8.0197e-017 (0.2,0.2) 3.6375e-001 1.9637e-015 1.6039e-016 (0.3,0.3) 5.2225e-001 2.9222e-015 2.4059e-016 (0.4,0.4) 6.6813e-001 3.8681e-015 3.2079e-016 (0.5,0.5) 8.0327e-001 4.8033e-015 4.0099e-016 (0.6,0.6) 9.2929e-001 5.7293e-015 4.8118e-016 (0.7,0.7) 1.0476e+000 6.6476e-015 5.6138e-016 (0.8,0.8) 1.1595e+000 7.5595e-015 6.4158e-016 (0.9,0.9) 1.2659e+000 8.4659e-015 7.2177e-016 (1.0,1.0) 1.3679e+000 9.3679e-015 8.0197e-016



---

*Source: 3194907-2018-07-25.xml*