Open access peer-reviewed chapter

Relationship between Interpolation and Differential Equations: A Class of Collocation Methods

By Francesco Aldo Costabile, Maria Italia Gualtieri and Anna Napoli

Submitted: May 15th 2016Reviewed: November 21st 2016Published: March 15th 2017

DOI: 10.5772/66995

Downloaded: 929

Abstract

In this chapter, the connection between general linear interpolation and initial, boundary and multipoint value problems is explained. First, a result of a theoretical nature is given, which highlights the relationship between the interpolation problem and the Fredholm integral equation for high-order differential problems. After observing that the given problem is equivalent to a Fredholm integral equation, this relation is used in order to determine a general procedure for the numerical solution of high-order differential problems by means of appropriate collocation methods based on the integration of the Fredholm integral equation. The classical analysis of the class of the obtained methods is carried out. Some particular cases are illustrated. Numerical examples are given in order to illustrate the efficiency of the method.

Keywords

  • boundary value problem
  • initial value problem
  • collocation methods
  • interpolation
  • Birkhoff
  • Lagrange
  • Peano
  • Fredholm

1. Introduction

The relationship between interpolation and differential equations theories has already been considered. In Ref. ([1], p. 72), Davis observed that the Peano kernel in the interpolation problem

y(a)=α,y(b)=β,a,b,α,βR,E1

is the Green’s function of the differential problem

φ(x)=f(x)φ(a)=φ(b)=0

where φ(x)=y(x)P1[y](x), being P1[y](x)the unique interpolatory polynomial for Eq. (1).

He observed that “these remarks indicate the close relationship between Peano kernels and Green’s functions, and hence between interpolation theory and the theory of linear differential equations. Unfortunately, we shall not be able to pursue this relationship” [1].

Later, Agarwal ([2], p. 2), Agarwal and Wong ([3], pp. 21, 151, 186) considered some separate boundary value problems and the related Fredholm integral equation, using only polynomial interpolation, without taking into account the related Peano kernel. They used Fredholm integral equation in order to obtain existence and uniqueness results for the solution of the considered boundary value problems.

Linear interpolation has an important role also in the numerical solution of differential problems. For example, finite difference methods (see, for instance, [46] and references therein) approximate the solution y(x)of a boundary value problem by a sequence of overlapping polynomials which interpolate y(x)in a set of grid points. This is obtained by replacing the differential equation with finite difference equations on a mesh of points that covers the range of integration. The resultant algebraic system of equations is often solved with iterative processes, such as relaxation methods.

Many authors (see [710] and references therein) used linear interpolation with spline functions for the numerical solution of boundary value problems.

Here, we take into account a more general nonlinear initial/boundary/multipoint value problems for high-order differential equations

{y(r)(x)=f(x,y(x)),xI=[a,b],r1Li[y](x)=wi,i=0,,r1,xIE2

where y(x)=(y(x),y(x),,y(q)(x)), 0q<r, yCr(I), and Liare rlinearly independent functionals on Cr(I). Moreover, we suppose that the function f:[a,b]×Rq+1Ris continuous at least in the interior of the domain of interest, and it satisfies a uniform Lipschitz condition in y, which means that there exists a nonnegative constant Λ, such that, whenever (x,y0,y1,,yq)and (x,y¯0,y¯1,,y¯q)are in the domain of f, the following inequality holds

|f(x,y0,y1,,yq)f(x,y¯0,y¯1,,y¯q)|Λk=0q|yky¯k|.E3

If Li[y]=Φ(y(a)),i=0,,r1, then (2) is an initial value problem (IVP); if Li[y]=Φ(y(a),y(b)),i=0,,r1, then (2) is a boundary value problem (BVP); if Li[y]=Φ(y(xj)),i=0,,r1, j=0,,m2, then (2) is a multipoint value problem (MVP).

In this chapter,

  1. - we assume that the conditions for the existence and uniqueness of solution of problem (2) in a certain appropriate domain of [a,b]×Rq+1are satisfied and that the solution y(x)is differentiable with continuity up to what is necessary;

  2. - we get the Fredholm integral equation related to problem (2), by polynomial interpolation and the Peano kernel of the linear interpolation problem Li[y](x)=wi,i=0,,r1. In this way, we point out the close relationship between Green’s function and Peano kernel;

  3. - then, we construct a class of spectral collocation (pseudospectral) methods which are derived by a linear interpolation process.

The reason for which we prefer collocation methods is their superior accuracy for problems whose solutions are sufficiently smooth functions. Recently, Boyd ([11], p. 8) observed that “When many decimal places of accuracy are needed, the contest between pseudospectral algorithms and finite difference and finite element methods is not an even battle but a rout: pseudospectral methods win hands-down.”

2. The Fredholm integral equation for problem (2)

We consider the general differential problem (2), and we prove that it is equivalent to a Fredholm integral equation.

Proposition 1 [1, p. 35] The linear interpolation problem

Li[P](x)=wi,wi,R,i=0,,r1,PPr1,xIE4

with Li, i=0,,r1, linearly independent functionals on Cr(I), has the unique solution

Pr1(t)=1G|01ttr1w0w1Li[xj]wr1|,G=|Li[xj]|i,j=0,,r1.E5

Proof. Since the Li, i=0,,r1are linearly independent, the result follows from the general linear interpolation theory.

Proposition 2 If y ∈ C r ( I ) and L i [ y ] ( x ) = w i , i = 0, … , r − 1 , x ∈ I , then

y(x)=Pr1[y](x)+abKrx(x,t)y(r)(t)dt,xIfixed,E6

with Li[y]=Li[Pr1], i=0,,r1, Pr1[y](x)=Pr1(x), and

Krx(x,t)=1(r1)![(xt)+r1Pr1[(xt)+r1](x)],E7

where index xmeans that Krx(x,t)is considered as a function of x.

Proof. It follows by observing that Pr1[(x)+j](t)=(t)+j,j=0,,r1and from Peano kernel Theorem [1].

Theorem 1 With the above notations and under the mentioned hypothesis, problem (2) is equivalent to the Fredholm integral equation

y(x)=Pr1[y](x)+abKrx(x,t)f(t,y(t))dt.E8

Proof. The result follows from the uniqueness of the Peano kernel and from Propositions 1 and 2.

Corollary 1 It results Li[Krx]=0,i=0,,r1.

From Theorem 1, general results on the existence and uniqueness of solution of problem (2) by standard techniques [2, 3] can be obtained. In the following, we will not linger over them, but we will outline the close relationship between interpolation and differential equations. Particularly, we will use linear interpolation in order to determine a class of collocation methods for the numerical solution of problem (2).

3. A class of Birkhoff-Lagrange collocation methods

Integral Eq. (8) allows to determine a very wide class of numerical methods for Eq. (2), which we call methods of collocation for integration.

Let {xi}i=1mbe mdistinct points in [a,b]and denoted by li(t), i=1,,m,the fundamental Lagrange polynomials on the nodes xi, that is

li(t)=ωm(t)(txi)ωm(xi),where ωm(t)=k=1m(txk).E9

Theorem 2 If the solution y(x)of Eq. (8) is in Cr+m(I), then

y(x)=Pr1[y](x)+i=1mpr,i,m(x)f(xi,y(xi))+Tr,m(y,x),E10

where

pr,i,m(x)=abKrx(x,t)li(t)dt,i=1,,m,E11

and the remainder term Tr,m(y,x)is given by:

Tr,m(y,x)=1m!abKrx(x,t)ωm(t)y(r+m)(ξx)dt,E12

being ξxa suitable point of the smallest interval containing xand all xi, i=1,,m.

Proof. From Lagrange interpolation

y(r)(x)=i=1mli(x)y(r)(xi)+R¯m(y,x)E13

where

R¯m(y,x)=1m!ωm(t)y(r+m)(ξx)E14

is the remainder term. From (2), f(x,y(x))=i=1mli(x)y(r)(xi)+R¯m(y,x). Then, from Theorem 1, inserting Eq. (13) into (8), we obtain Eq. (10).

Theorem 2 suggests to consider the implicitly defined polynomial

yr,m(x)=Pr1[yr,m](x)+i=1mpr,i,m(x)f(xi,yr,m(xi)).E15

For polynomial (15), the following theorem holds.

Theorem 3 (The main Theorem). Polynomial (15), of degree r + m − 1 , satisfies the relations

Li[yr,m](x)=wi,i=0,,r1,xI,wiRyr,m(r)(xj)=f(xj,yr,m(xj))j=1,,m,E16

that is, yr,m(x)is a collocation polynomial for Eq. (2) at nodes xj, j=1,,m.

Proof. From (15), Corollary 1 and the linearity of operators Li, we get Li[yr,m](x)=wi,i=0,,r1. By Theorems 1 and 2, we obtain y(r)(xi)=yr,m(r)(xi), and from Eq. (11), pr,i,m(r)(x)=li(x). Hence, relations (16) follow.

Remark 1 (Hermite-Birkhoff-type interpolation). Theorem 3 is equivalent to the general Hermite-Birkhoff interpolation problem [12]: given w i ∈ R , i = 0, … , r − 1 , and α j ∈ R , j = 1, … , m , determine, if there exists, the polynomial Q ( x ) ∈ P m + r − 1 such that

Li[Q]=wi,i=0,,r1Q(r)(xj)=αj,j=1,,m,xjI.E17

Remark 2 In the case of IVPs, for each method (15), we can derive the corresponding implicit Runge-Kutta method. For example, for r = 2 , let b = x 0 + h and x i = x 0 + c i h with c i ∈ [ 0,1 ] . With the change of coordinates x = x 0 + t h , t ∈ [ 0,1 ] , we can write

pr,i,m(x)=pr,i,m(x0+th)=h20t0rli(s)dsdr,li(s)=k=1kimsckcick.E18

Putting f(xi,yr,m(xi))=yr,m(xi)Ki,ai,j=pr,j(xi)=h20ci(cis)lj(s)ds,we have

Ki=f(x0+cih,y0+y0th+j=1mai,jKj)E19

and

{y1(t)yr,m(x0+th)=y0+y0th+h2i=1mpr,i,m(x0+th)Kiy1(t)yr,m(x0+th)=y0h+h2i=1mpr,i,m(x0+th)Ki.E20

Eqs. (19) and (20) are the well-known continuous Runge-Kutta method for second-order differential equations. Particularly, for t=1, we have the implicit Runge-Kutta-Nystrom method.

3.1. A-priori estimation of error

In what follows, we consider the norm

f=maxatbk=0q|f(k)(t)|,fC(q)(I).E21

Moreover, we define

Qm=i=1mpr,i,m,F(x)=abKrx(x,t)dt,H=maxatb|R¯m(y,t)|,E22

where R¯m(y,t)is defined as in (14).

Theorem 4 With the previous notations, if ΛQm<1, then

yyr,mHF1ΛQm.E23

Proof. By deriving Eqs. (10) and (15), stimes, s=0,,q, we get

y(s)(x)yr,m(s)(x)=i=1mpr,i,m(s)(x)[f(xi,y(xi))f(xi,yr,m(xi))]+sxsabKrx(x,t)R¯m(y,t)dt.E24

It follows that

|y(s)(x)yr,m(s)(x)|i=1m|pr,i,m(s)(x)|Λk=0q|y(k)(xi)yr,m(k)(xi)|+H|F(s)(x)|Λyyr,mi=1m|pr,i,m(s)(x)|+H|F(s)(x)|.E25

From this, we obtain inequality (23).

4. Algorithms and implementation

To calculate the approximate solution of problem (2) by polynomial yr,m(x)at xI, we need the values yr,m(s)(xk), k=1,,m, s=0,,q. In order to get these values, we propose the following algorithm:

- Put yk(s)=yr,m(s)(xk), k=1,,m, s=0,,qand consider the following system

yk(s)=Pr1(s)[yk](xk)+i=1mpr,i(s)(xk)f(xi,yi),E26

k=1,,m, s=0,,q, where yi=(yi,yi,,yi(q)).

System (26) can be written in the form

YAF(Y)=CE27

where

A=(A0000000Aq)m(q+1)×m(q+1)E28

with

As=(a˜1,1(s)a˜1,m(s)a˜m,1(s)a˜m,m(s))m×ma˜i,j(s)=pr,j(s)(xi),s=0,,q,E29
Y=(Y¯0,,Y¯q)m(q+1)×1T,Y¯s=(y1(s),,ym(s)),E30
F(Y)=(Fm,,Fmq)T,Fm=(f1,,fm)T,fi=f(xi,yi),E31
Bs=(Pr1(s)[y1](x1),,Pr1(s)[ym](xm)),C=(B0,,Bq)m(q+1)×1T.E32

From Eq. (27), we get

Y=AF(Y)+C,E33

or, putting G(Y)=AF(Y)+C,

Y=G(Y).E34

For the existence and uniqueness of solution of system (34), we can prove, with standard technique, the following theorem.

Theorem 5 If T=ΛA<1, system (34) has a unique solution which can be calculated by an iterative method

(Ym)ν+1=G((Ym)(ν)),ν0E35

with a fixed (Ym)0Rm(q+1)and G(Ym)=AF(Ym)+C.

Moreover, if Y is the exact solution,

(Ym)ν+1YTν1T(Ym)1(Ym)0.E36

Remark 3 If f is linear, then system (27) is a linear system which can be solved by a more suitable method.

Remark 4 System (27) can be considered as a discrete method for the numerical solution of (2).

Remark 5 Method (15) can generate the polynomial sequence

(yr,m(x))ν=Pr1[yr,m](x)+i=1mpr,i,m(x)f(xi,(yr,m(xi))ν1),(yr,m)0=Pr1[y](x)E37

which is equivalent to the discretization of Picard method for differential equations.

4.1. Numerical computation of the entries of matrix A

To calculate the elements a˜i,k(s)of the matrix A in Eq. (27), we have to compute the integrals

pr,k(s)(x)=dsdxsabKrx(x,t)li(t)dtE38

for x=xi. Integrating by parts, it remains to solve the problem of the computation of

Fi1(xj)=axjli(t)dt,Fik(xj)=axjFi,k1(t)dtk=2,,nE39
Mi1(xj)=xjbli(t)dt,Mik(xj)=xjbMi,k1(t)dtk=2,,nE40

i,j=1,m. To this aim, it suffices to compute

cxj=tkctk1ct1rm,i(t)dtdt1dtk1E41

where c=aor c=b, r0,0(t)=1,

rm,i(t)=(tx1)(txi1)(txi+1)(txm)i=1,2,,m.E42

For the computation of the integral (41), we use the recursive algorithm introduced in Ref. [13]: for each i=1,,m, let us consider the new points zj(i)=xjif j<i, and zj(i)=xj+1if ji. Moreover, let us define g0,1,c(i)(x)=xcand for s=1,,m1

gs,j,c(i)(x)=cx=tjctj1ct1(tz1(i))(tz2(i))(tzs(i))dtdt1dtj1.E43

We can easily compute g0,j,c(i)(x)=(xc)jj!.For the computation of Eq. (43), the following recurrence formula [13] holds

gs,j,c(i)(x)=(xzs(i))gs1,j,c(i)(x)jgs1,j+1,c(i)(x).E44

Thus, if Wi=k=1,kim(xixk), then

Fik(xj)=gm1,k,a(i)(xj)Wi,Mik(xj)=(1)kgm1,k,b(i)(xj)Wi.E45

Remark 6 An alternative approach for the exact computation of integrals (39) and (40) is to use a quadrature formula with a suitable degree of precision.

4.2. Outline of the method

Summarizing the proposed method consists of the following steps:

  1. determine the interpolation polynomial Pr1[y](x)satisfying the boundary conditions and compute the Peano remainder;

  2. approximate y(r)(x)by Lagrange interpolation on a set of fixed nodal point;

  3. compute the elements of matrix A (28) and solve system (26);

  4. obtain polynomial (15).

5. Some particular cases

Now we consider some special cases of problem (2), and for each case, we determine Pr1[y](x)and Krx(x,t). We also show how the proposed class of methods includes the methods presented in previous works [1224].

5.1. Initial value problems

In the case of initial value problems, in Refs. [13, 17, 25], problem

y(r)(x)=f(x,y(x))E46

has been considered, while in Ref. [23], the authors introduced the more general equation

y(r)(x)=f(x,y(x),y(x),,y(r)(x)),qr1.E47

In both cases

Pr1[y](x)=i=0r1(xa)ii!y(i)(a)E48

and

Krx(x,t)=1(r1)!(xt)+r1.E49

If {xi}i=1mare the zeros of Chebyshev polynomials of first and second kind, the explicit expression for polynomials pr,i,m(x)can be obtained [13, 17, 25] for some values of r.

Particularly, for r=1and r=2, in the case of zeros of Chebyshev polynomials of first kind, we get

p1,i,m(x)=1mk=2m1{[Tk+1(x)k+1Tk1(x)k1+2(1)k1k21]cos(2i12mkπ)}+1m[x+1+cos(2i12mπ)(x21)]E50

where Tk1(x)and Tk+1(x)are the Chebyshev polynomials of the first kind and degree k1and k+1, respectively, and

p2,i,m(x)=1m{(x+1)22+x33x23(cosπ(2i1)2m+xcosπ(2i1)m)+12k=3m1coskπ(2i1)2m[Tk+2(x)(k+1)(k+2)2Tk(x)k21+Tk2(x)(k1)(k2)12k(1)kk(k21)(k24)4(1)kk21(x+1)]}.E51

In the case of zeros of Chebyshev polynomials of second kind

p1,i,m(x)=2m+1sinπim+1k=0m1sin(k+1)πim+11k+1[Tk+1(x)+(1)k]E52

and

p2,i,m(x)=1m+1sinπim+1{sinπim+1(x+1)2+k=2m1ksinkπim+1[Tk+1(x)k+1Tk1(x)k12(x+k2k21)(1)k]}E53

In Refs. [13, 25], the authors presented the corresponding implicit Runge-Kutta methods too.

In Ref. [26], Coleman and Booth used also a polynomial interpolant of degree nfor y, but they started from an identity different to Eq. (8) and derived a collocation method for which the nodes {xi}i=1mare the zeros of Chebyshev polynomials of second kind.

5.2. Boundary value problems

5.2.1. Case r=2n

For n=1, for the exact solution y(x)of the second-order BVP

y(x)=f(x,y(x),y(x)),y(1)=y0,y(1)=y1E54

x[1,1], it is known that

y(x)=y1+y02+xy1y02+11K2x(x,t)f(x,y(x),y(x))dtE55

where

K2x(x,t)={(t+1)(x1)2tx(x+1)(t1)2x<t.E56

By applying method (15), we get [16]

y2,m(x)=y1+y02+xy1y02+i=1mpr,i,m(x)f(xi,y(xi),y(xi))E57

with pr,i,m(x)=11K2x(x,t)li(t)dt.

If xi=cosπim+1, i=1,,m, we obtain explicitly the expression of pr,i,m(x)[18]

pr,i,m(x)=1m+1sinπim+1[k=2mGk(x)ksinkπim+1+(x21)sinπim+1]E58

where

Gk(x)=Tk+1(x)k+1Tk1(x)k1+{2xk21evenkk32k21oddk.E59

The same method has been presented in Ref. [24], where also stability has been studied.

Now assume [a,b]=[0,1]and r>2. Several types of boundary conditions can be considered.

-Hermite boundary conditions [22]:

y(h)(0)=αh,y(h)(1)=βh,h=0,,n1E60

with αh,βh, h=0,,n1real constants.

In this case, Pr1is the Hermite polynomial of degree 2n1

P2n1[y](x)=i=0n1(y(i)(0)Hi1(x)+y(i)(1)Hi2(x))E61

with

Hi1(x)=xi(1x)ni!s=0ni1(n+s1n1)xsHi2(x)=xn(1x)ii!s=0ni1(n+s1n1)(1x)s.E62

The kernel is

K2nx(x,t)={i=0n1(t)2ni1(2ni1)!Hi1(x)0txi=0n1(1t)2ni1(2ni1)!Hi2(x)xt1.E63

-Lidstone boundary conditions [19]:

y(2h)(0)=αh,y(2h)(1)=βh,h=0,,n1E64

where αh,βh, h=0,,nare real constants.

In this case, Pr1is the Lidstone interpolating polynomial [3] of degree 2n1

P2n1[y](x)=k=0n1[y(2k)(0)Λk(1x)+y(2k)(1)Λk(x)]E65

where Λk(x)are the Lidstone polynomials of degree 2k+1[3], and the function K2nx(x,t)is

K2nx(x,t)={k=0n1t2n2k1(2n2k1)!Λk(1x)txk=0n1(1t)2n2k1(2n2k1)!Λk(x)xt.E66

5.2.2. Case r=2n+1

If we consider the following boundary conditions

y(0)=α0,y(2h1)(0)=αh,y(2h1)(1)=βh,h=1,,nE67

with α0,αh,βh, h=1,,nreal constants, then Pr1is the complementary Lidstone interpolating polynomial [27] of degree 2n[3, 24, 27, 28].

P2n[y](x)=y(0)+k=1n[y(2k1)(0)(vk(1)vk(1x))+y(2k1)(1)(vk(x)vk(0))],E68

where vk(x)are the complementary Lidstone polynomials of degree k [27]. The kernel is

K2n1x(x,t)={t2n(2n)!+k=1nt2n2k+1(2n2k+1)!(vk(1x)vk(1))txk=1n(1t)2n2k+1(2n2k+1)!(vk(x)vk(0))xt.E69

In Ref. [19], the proposed method applied to problem (2) with conditions (64) and (67), respectively, has been examined in detail.

5.2.3. Other special boundary conditions

If r=n1and [a,b]=[0,1], we can consider Bernoulli boundary conditions [21]

y(0)=β0,y(1)=β1,y(k)(1)y(k)(0)=βk+1,k=1,,n2E70

with βk, k=0,,n1real constants. The method has been examined in [14].

5.3. Multipoint boundary value problems

Let us now consider [15] the following conditions in I=[1,1]

y(k)(1)=αk,k=0,,s1,y(s)(xi)=ωii=1,,rs.E71

In this case

Pr1[y](x)=i=0s1(x+1)ii!αi+1(s1)!k=1rsωkpr,k(x),E72

with

pr,k(x)=1x(xt)s1lk(t)dtE73

and lk(t)are the fundamental Lagrange polynomials on the points xj,j=1,,rs. Pr1(x)is the unique polynomial of degree r1which satisfies the Birkhoff interpolation problem

Pr1(k)(1)=αk,k=0,,s1,Pr1(s)(xi)=ωi,i=1,,rs,sr1E74

with 1<x1<<xrs1. Hence, the solution of problem (2), with multipoint conditions (71), is

y(x)=Pr1[y](x)+11Krx(x,t)y(r)(t)dt,E75

with Pr1[y](x)given in Eq. (72) and

Krx(x,t)=1(r1)![(xt)+r1(r1s)si=1rspr,i,m(x)(xit)+rs1].E76

Observe that Eq. (74) is a special type of Birkhoff interpolation problem with incidence matrix E=(eij)defined by e1j=eis=1,j=0,,s1,i=2,,rs+1,eij=0otherwise and r1.

In Ref. [23], Pr1[y](x)is presented in a little different form:

Pr1[y](x)=i=0s1(x+1)ii!αi+k=1rsωkEs(x,lk(x)),E77

where Es(x,lk(x))=1x1xslk(t)dtdt.

Let us now consider the following conditions [12, 20]

y(1)=ω0,y(1)=ωr1y(xi)=ωii=1,,r2.E78

The solution to the Birkhoff interpolation problem

Pr1(1)=ω0,Pr1(1)=ωr1,Pr1(xi)=ωi,i=1,,r2E79

with 1<x1<<xr2<1is [12]

Pr1[y](x)=ωr1+ω02+ωr1ω02x+i=1r2qr,i(x)ωiE80

with

qr,i(x)=11Krx(x,t)li(t)dtE81

and

Krx(x,t)={(t+1)(x1)2tx(x+1)(t1)2x<t.E82

Hence, the solution of problem (2) is

y(x)=Pr1[y](x)+11Krx(x,t)y(r)(t)dt,E83

with Pr1[y](x)given in Eq. (80) and

Krx(x,t)=1(r1)![(xt)+r1(1t)r1(1+x)2(r1)(r2)i=1r2pr,i,m(x)(xit)+r3].E84

6. Numerical examples

In this section, we present some numerical results obtained by applying method (15), which we call CGN method, to find numerical approximations yr,m(x)to the solution of some test problems. In order to solve the nonlinear system (19), we use the so-called modified Newton method [29] (the same Jacobian matrix is used for more than one iteration) and we use algorithm (44) for the computation of the entries of the matrix, when polynomials pr,i,m(x)are not explicitly known. Since the true solutions of the analyzed problems are known, we consider the error function em(x)=|y(x)yr,m(x)|.

The maximum values of em(x)over the interval [a,b]have also been calculated by using Matlab, particularly the built-in solvers

  • ode15s, a variable-step, variable-order multistep solver based on the numerical differentiation formulas of orders 1–5;

  • ode45, a single-step solver, based on an explicit Runge-Kutta (4, 5) formula, the Dormand-Prince pair

for initial value problems, and the finite difference codes;

  • bvp4c (with an optional mesh of 200 points) that implements the three-stage Lobatto IIIa formula;

  • bvp5c that implements the four-stage Lobatto IIIa formula.

for boundary value problems.

All solvers have been used with optional parameters RelTol=AbsTol=1e−17.

Moreover, the powerful tool Chebfun [30] has been used.

Example 1 Consider the following linear ninth-order BVP [28]

{y(9)(x)=9ex+y(x)x[0,1]y(j)(0)=1jj=0,,4y(j)(1)=jej=0,,3E85

with exact solution y(x)=(1x)ex.

The unique polynomial P 8 ( x ) = P 8 [ y ] ( x ) of degree 8 satisfying the boundary conditions P 8 ( j ) ( 0 ) = 1 − j for j = 0, … ,4, and P 8 ( j ) ( 1 ) = − j   e j = 0, … ,3 is

P8(x)=112x213x318x4+(3121e2536)x5+(1321128121e)x6+(712e1932)x7+(68524212e)x8.E86

From Eq. (7), we get

K9x(x,t)=18!{70t4(x44x5+6x64x7+x8)+56t5(x3+10x520x6+15x74x8)+28t6(x220x5+45x636x7+10x8)+8t7(x+35x584x6+70x720x8)+t8(156x5+140x6120x7+35x8)0txx8+8tx728t2x6+56t3x5+70t4(4x5+6x64x7+x8)+56t5(10x520x6+15x74x8)+28t6(20x5+45x636x7+10x8)+8t7(35x584x6+70x720x8)+t8(56x5+140x6120x7+35x8)xt1.E87

Now we calculate the values of the integrals (39) by using Eq. (45), and we solve system (26). Thus, we obtain the approximate solution (15) to problem (85).

Table 1 shows the numerical results. The absolute errors are compared with those obtained in Ref. [28], where a modified decomposition method is applied for the solution of problem (85). The second and third columns of Table 1 show the error, respectively, in the method in Ref. [28] and in the CGN method, using in both cases polynomials of degree 12. The last column contains the error in the approximation by a polynomial of degree 14 using CGN method. As collocation points, equidistant nodes in [0,1]are chosen. Analogous results are obtained by using Chebyshev nodes of first and second kind, and Legendre-Gauss-Lobatto points.

xMethod in [28]CGN m=4CGN m=6
0.12.0e101.45e140.00
0.22.0e103.93e131.11e16
0.32.0e102.16e129.99e15
0.42.0e105.70e122.00e15
0.52.0e109.27e122.55e15
0.66.0e101.00e112.66e15
0.71.0e97.04e122.44e15
0.82.0e92.70e122.83e15
0.93.4e92.98e134.91e15

Table 1.

Absolute error em(x)in MDM and CGN methods for problem (85).

The maximum absolute error max{em(x)}on [0,1]has also been calculated by using Matlab ( Table 2 ) .

Chebfunbvp4cbvp5c
1.461.55e124.44e16

Table 2.

Maximum absolute error in problem (85) using Matlab built-in functions.

Example 2 Consider the fifth-order initial value problem [13]

{y(5)+(32x5+120x)y=160x3ex2x[0,1]y(0)=1,y(0)=0,y(0)=2y″′(0)=0,y(4)(0)=12E88

with solution y(x)=ex2.

Table 3 shows the absolute error in some points of the interval [0,1]for CGN method in the case, respectively, of Chebyshev nodes of first kind (Cheb I), of second kind (Cheb II) and in the case of equidistant nodes (EqPts).

xCheb ICheb IIEqPts
m=4m=6m=9
0.11.11e160.000.00
0.29.54e150.000.00
0.35.47e133.33e160.00
0.49.45e121.11e164.44e16
0.58.50e114.22e151.11e16
0.65.05e103.47e142.11e15
0.72.25e92.08e131.55e15
0.88.08e99.68e131.44e14
0.92.74e83.72e129.18e15
1.06.64e81.22e111.37e14

Table 3.

Problem (88)—example 2.

The maximum absolute errors calculated by using Matlab are displayed in Table 4 .

Chebfunode15sode45
2.11e111.35e131.33e15

Table 4.

Maximum absolute error in problem (88) using Matlab built-in functions.

Example 3 Consider now the following nonlinear problem [31]

{y(4)(x)=sinx+sin2x(y(x))2x[0,1]y(0)=0y(0)=1y(1)=sin(1)y(1)=cos(1)E89

with exact solution y(x)=sin(x).

This kind of problems models several nonlinear phenomena such as traveling waves in suspension bridges [32] or the bending of an elastic beam [33].

Suspension bridges are generally susceptible to visible oscillations, due to the forces acting on the bridge (including the force due to the cables which are considered as a spring with a one-sided restoring, the gravitation force and the external force due to the wind or other external sources). frepresents the forcing term, while yrepresents the vertical displacement when the bridge is bending.

In the case of elastic beam, frepresents the force exerted on the beam by the supports. xmeasures the position along the beam (x=0is the left-hand endpoint of the beam), yand yindicate, respectively, the height and the slope of the beam at x. ymeasures the curvature of the graph of y, and, in physical terms, it measures the bending moment of the beam at x, that is, the torque that the load places on the beam at x .

The considered boundary conditions state that the beam has both endpoints simply supported. Moreover, the derivative of the deflection function is not zero at those points, and it indicates that the beam at the wall is not horizontal.

Table 5 shows the comparison between the NMD method presented in Ref. [31] and the CGN method with m=5and m=9, respectively. The approximating polynomial of NMD method has degree 11, while the polynomial considered in CGN method for m=5has degree 8.

xNMD [31]CGN m=5CGN m=9
0.17.78e84.45e101.53e15
0.22.72e75.54e103.02e15
0.35.24e78.95e117.77e16
0.47.77e72.03e106.66e16
0.59.71e73.32e115.55e17
0.61.05e61.53e100
0.79.63e79.48e110
0.86.84e75.18e101.11e16
0.92.71e74.15e100

Table 5.

Error of NMD and CGN methods—problem (89).

The maximum absolute errors calculated by using Matlab are displayed in Table 6 .

Chebfunbvp4cbvp5c
1.67e161.22e88.88e16

Table 6.

Maximum absolute error in problem (89) using Matlab build-in functions.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Francesco Aldo Costabile, Maria Italia Gualtieri and Anna Napoli (March 15th 2017). Relationship between Interpolation and Differential Equations: A Class of Collocation Methods, Dynamical Systems - Analytical and Computational Techniques, Mahmut Reyhanoglu, IntechOpen, DOI: 10.5772/66995. Available from:

chapter statistics

929total chapter downloads

1Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Integral-Equation Formulations of Plasmonic Problems in the Visible Spectrum and Beyond

By Abdulkerim Çekinmez, Barişcan Karaosmanoğlu and Özgür Ergül

Related Book

First chapter

Appell-Gibbs Approach in Dynamics of Non-Holonomic Systems

By Jiří Náprstek and Cyril Fischer

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us