一类Fredholm积分–微分方程的重心Jacobi插值求解法
Barycentric Jacobi Interpolation Method for Solving a Kind of Fredholm Integral-Differential Equation
DOI: 10.12677/pm.2024.144142, PDF, HTML, XML, 下载: 76  浏览: 148  科研立项经费支持
作者: 刘高原, 张 益:西华师范大学数学与信息学院,四川 南充;陈 冲:西华师范大学公共数学学院,四川 南充
关键词: 重心Jacobi插值配置法Fredholm积分–微分方程误差估计Barycentric Jacobi Interpolation Collocation Method Fredholm Integral-Differential Equation Error Estimates
摘要: 本文运用重心Jacobi插值配置法求解一类Fredholm积分–微分方程。首先通过取消重心Gegenbauer插值中参数相等的条件,得到重心Gegenbauer插值的一般形式——重心Jacobi插值,并表明重心Jacobi插值等价于插值节点为移位Gauss-Jacobi节点的Jacobi插值。然后基于配置法,应用重心Jacobi插值构造一类带有初值条件的Fredholm积分–微分方程的数值算法。误差估计表明,在合适的条件下,该算法是收敛的。最后,数值算例验证算法的有效性。
Abstract: In this paper, the barycentric Jacobi interpolation collocation method is used to solve a kind of Fredholm Integral-Differential equation. Firstly, the general form of barycentric Gegenbauer interpolation, barycentric Jacobi interpolation, is obtained by canceling the condition that the parameters in barycentric Gegenbauer interpolation are equal, and it is shown that the barycentric Jacobi interpolation is equivalent to the Jacobi interpolation whose interpolation nodes are shifted Gauss-Jacobi nodes. Then based on the collocation method, the numerical algorithm for a kind of Fredholm Integral-Differential equation with initial value conditions is constructed by barycentric Jacobi interpolation. The result of error estimates show that the algorithm is convergent under suitable conditions. Finally, the effectiveness of the algorithm is verified by a numerical example.
文章引用:刘高原, 张益, 陈冲. 一类Fredholm积分–微分方程的重心Jacobi插值求解法[J]. 理论数学, 2024, 14(4): 342-354. https://doi.org/10.12677/pm.2024.144142

1. 引言

积分–微分方程是积分方程与微分方程的交叉学科,在众多类型的积分–微分方程中,Fredholm积分–微分方程 [1] 常出现在数学物理应用问题中,如热–粘弹性力学问题以及三维薄翼理论相关问题 [2] ,因此,Fredholm积分–微分方程在理论上具有重要的研究意义,因这类积分–微分方程不仅包括未知函数的微分,还包括未知函数的积分,同时还满足多个边值条件,故求其精确解较困难,需采用适当的数值方法求其数值解。

数十年来,Fredholm积分–微分方程的高效数值解法一直是计算数学领域中的研究热点,同伦分析法 [3] 、配置法 [4] 、Galerkin法 [5] 及Nystörm法 [6] 等都是求Fredholm积分–微分方程的数值解有效的方法。在这些方法中,多项式是函数逼近的常用工具,近年来,Bessel多项式 [7] 、Legendre多项式 [8] 以及Chebyshev多项式 [9] 等正交多项式已被应用于构造Fredholm积分–微分方程的数值解。文献 [10] 将Fredholm积分–微分方程的求解问题转化为积分方程的求解问题,并运用正交Gegenbauer插值配置法求解了一类Fredholm积分–微分方程。文献 [11] 表明以Gauss-Gegenbauer节点为插值节点的Gegenbauer插值在数学意义上等价于Lagrange插值。Lagrange插值可改写为重心Lagrange插值,文献 [12] 提出了以Gauss-Gegenbauer节点为插值节点的Gegenbauer插值的重心形式——重心Gegenbauer插值,并将其应用于求解积分、微分以及积分–微分方程,同时表明以重心Gegenbauer插值为核心的求积方法是高阶、稳定且有效的伪谱方法。从目前的研究状况来看,重心型插值并未被广泛应用于求解Fredholm积分–微分方程,并且关于重心Gegenbauer插值中参数取值范围的问题并未被进一步研究。本文考虑扩大重心Gegenbauer插值中参数的取值范围,并基于重心Gegenbauer插值的一般形式,运用配置法求解一阶线性Fredholm积分–微分方程 [3] [4] [5]

{ u ( 1 ) ( x ) + u ( x ) + a b K ( x , t ) u ( t ) d t = f ( x ) , x D = [ a , b ] , u ( a ) = c , c , (1)

其中, f ( x ) K ( x , t ) 分别是D和 D × D 上的已知连续函数, u ( x ) 是方程(1)的未知函数,且 u ( x ) 的一阶导数 u ( 1 ) ( x ) 满足 u ( 1 ) ( x ) = d u ( x ) / d x

本文共分为四节,其中,第一节的主要内容是证明重心Jacobi插值等价于以移位Gauss-Jacobi节点为插值节点的Jacobi插值,并阐明重心Jacobi插值配置法求解积分/微分方程的基本原理与步骤。第二节的主要内容是应用重心Jacobi插值以及Gauss-Legendre求积公式构造方程(1)的数值算法,算法的构造思路是将方程(1)转化为一个积分方程,然后应用重心Jacobi插值配置法将其离散为一个线性方程组。第三节的主要内容是对本文算法进行误差估计,并通过算法的误差估计式分析本文方法的收敛性质。第四节的主要内容是通过数值结果验证本文方法的有效性,并探究Jacobi多项式中两参数的取值对数值结果的影响。

2. 重心Jacobi插值

2.1. 重心Jacobi插值多项式

定义1 [13] 设 d 1 = b a , d 2 = a + b 为任意实数,区间 Ω = [ 1 , 1 ] 上的 j ( j ) 次Jacobi多项式为 J j ( α , β ) ( x ^ ) = 1 2 j k = 0 j C j + α k C j + β j k ( x ^ 1 ) j k ( x ^ + 1 ) k ( x ^ Ω , α , β > 1 ) ,则任意区间 D = [ a , b ] 上的j次Jacobi多项式定义为

J ˜ j ( α , β ) ( x ) = J j ( α , β ) ( 2 x d 2 d 1 ) , x D .

定义2 [14] 设 u ( x ) 是定义在区间D上的实值函数, x k ( k = 0 , 1 , , m + ) 是区间D上的节点且 u ( x k ) 已知,若存在实数 u ˜ j ( j = 0 , 1 , , m ) ,使得函数

I m u ( x ) = j = 0 m u ˜ j J ˜ j ( α , β ) ( x ) , x D , α , β > 1 , (2)

满足 I m u ( x k ) = u ( x k ) ,则称 I m u ( x ) 是实值函数 u ( x ) 在区间D上以 { x k } k = 0 m 为插值节点的Jacobi插值多项式。

在定义2中,若插值节点是 J ˜ m + 1 ( α , β ) ( x ) 的零点 { x m , j ( α , β ) } j = 0 m ,则系数 u ˜ j ( j = 0 , 1 , , m ) 满足 [14]

u ˜ j = i = 0 m ϖ i J ˜ j ( α , β ) ( x m , i ( α , β ) ) u ( x m , i ( α , β ) ) ρ j ( α , β ) , j = 0 , 1 , , m , (3)

其中, ϖ i ( i = 0 , 1 , , m ) 是Gauss-Jacobi求积权, ρ j ( α , β ) 是与 α , β 有关的正则化因子,满足 [14]

ρ j ( α , β ) = 2 ( α + β + 1 ) Γ ( j + α + 1 ) ( 2 j + α + β + 1 ) j ! Γ ( j + β + 1 ) Γ ( j + α + β + 1 ) , j = 0 , 1 , , m , (4)

其中, Γ ( ) 是定义在实数域 上的Gamma函数。将(3)式代入(2)式,令

L m , i ( α , β ) ( x ) = ϖ i j = 0 m ( ρ j ( α , β ) ) 1 J ˜ j ( α , β ) ( x m , i ( α , β ) ) J ˜ j ( α , β ) ( x ) , x D , i = 0 , 1 , , m , (5)

则(2)式可记为

I m u ( x ) = i = 0 m L m , i ( α , β ) ( x ) u ( x m , i ( α , β ) ) , x D , α , β > 1. (6)

定理1 若 I m u ( x ) 形如(6)式,且满足(4)式及(5)式,则(6)式是关于实值函数 u ( x ) 的以 { x m , j ( α , β ) } j = 0 m 为插值节点的Lagrange插值多项式。

证明 利用定义1,在(5)式中,将区间D上的Jacobi多项式转换为区间 Ω 上的Jacobi多项式,可得

L m , i ( α , β ) ( x ^ ) = ϖ i j = 0 m ( ρ j ( α , β ) ) 1 J j ( α , β ) ( x ^ m , i ( α , β ) ) J j ( α , β ) ( x ^ ) , x ^ Ω , i = 0 , 1 , , m , (7)

其中, x m , i ( α , β ) = ( b a ) 2 x ^ m , i ( α , β ) + ( a + b ) 2 。由文献 [14] 中的Christoffel-Darboux定理,结合(4)式,则(7)式可化为

L m , i ( α , β ) ( x ^ ) = ϖ i 2 ( α + β ) Γ ( m + 2 ) Γ ( m + α + β + 2 ) ( 2 m + α + β + 2 ) Γ ( m + α + 1 ) Γ ( m + β + 1 ) J m + 1 ( α , β ) ( x ^ ) J m ( α , β ) ( x ^ m , i ( α , β ) ) x ^ x ^ m , i ( α , β ) , (8)

因为 { x ^ m , i ( α , β ) } i = 0 m J m + 1 ( α , β ) ( x ^ ) 的零点,故 J m + 1 ( α , β ) ( x ^ m , i ( α , β ) ) = 0 ,于是(8)式可化为

L m , i ( α , β ) ( x ^ ) = 2 ( α + β ) Γ ( m + 2 ) Γ ( m + α + β + 2 ) ϖ i ( 2 m + α + β + 2 ) Γ ( m + α + 1 ) Γ ( m + β + 1 ) [ J m + 1 ( α , β ) ( x ^ ) J m + 1 ( α , β ) ( x ^ m , i ( α , β ) ) ] J m ( α , β ) ( x ^ m , i ( α , β ) ) x ^ x ^ m , i ( α , β ) , (9)

将节点 x ^ m , k ( α , β ) ( k = 0 , 1 , , m ) 代入(9)式中,可得

L m , i ( α , β ) ( x ^ m , k ( α , β ) ) = 2 ( α + β ) Γ ( m + 2 ) Γ ( m + α + β + 2 ) ϖ i ( 2 m + α + β + 2 ) Γ ( m + α + 1 ) Γ ( m + β + 1 ) [ J m + 1 ( α , β ) ( x ^ m , k ( α , β ) ) J m + 1 ( α , β ) ( x ^ m , i ( α , β ) ) ] J m ( α , β ) ( x ^ m , i ( α , β ) ) x ^ m , k ( α , β ) x ^ m , i ( α , β ) . (10)

(10)式中,当 i k 时,因为 J m + 1 ( α , β ) ( x ^ m , k ( α , β ) ) = 0 ,故 L m , i ( α , β ) ( x ^ m , k ( α , β ) ) = 0 ;当 i = k 时,由洛必达法则 [8] ,对任意 i = 0 , 1 , , m ,有

lim x x ^ m , i ( α , β ) [ J m + 1 ( α , β ) ( x ^ ) J m + 1 ( α , β ) ( x ^ m , i ( α , β ) ) ] x ^ x ^ m , i ( α , β ) = d d x ^ J m + 1 ( α , β ) ( x ^ m , i ( α , β ) ) ,

于是

L m , i ( α , β ) ( x ^ m , i ( α , β ) ) = 2 ( α + β ) Γ ( m + 2 ) Γ ( m + α + β + 2 ) ϖ i ( 2 m + α + β + 2 ) Γ ( m + α + 1 ) Γ ( m + β + 1 ) J m ( α , β ) ( x ^ m , i ( α , β ) ) d d x ^ J m + 1 ( α , β ) ( x ^ m , i ( α , β ) ) . (11)

由Gamma函数的基本性质,可得(11)式中

Γ ( m + 2 ) = ( m + 1 ) ! , Γ ( m + α + 2 ) Γ ( m + α + 1 ) = m + α + 1 ,

将该式以及Gauss-Jacobi求积权 { ϖ j } j = 0 m 的显式表达式 [15] 代入(11)式,可得

L m , i ( α , β ) ( x ^ m , i ( α , β ) ) = 2 ( m + α + 1 ) ( m + β + 1 ) J m ( α , β ) ( x ^ m , i ( α , β ) ) ( 2 m + α + β + 2 ) ( 1 ( x ^ m , i ( α , β ) ) 2 ) d d x ^ J m + 1 ( α , β ) ( x ^ m , i ( α , β ) ) , i = 0 , 1 , , m . (12)

在Jacobi多项式三项迭代式 [15] 中令 x ^ = x ^ m , i ( α , β ) ,可得

( 1 ( x ^ m , i ( α , β ) ) 2 ) d J m + 1 ( α , β ) ( x ^ m , i ( α , β ) ) d x ^ = 2 ( m + α + 1 ) ( m + β + 1 ) ( 2 m + α + β + 2 ) J m ( α , β ) ( x ^ m , i ( α , β ) ) , i = 0 , 1 , , m , (13)

将(13)式代入(12)式,可得 L m , i ( α , β ) ( x ^ m , i ( α , β ) ) = 1

综上, L m , i ( α , β ) ( x ^ m , k ( α , β ) ) = δ i k ( i , k = 0 , 1 , , m ) ,其中 δ i k 是Kronecker-Delta函数,满足当 i k 时, δ i k = 1 ;当 i = k 时, δ i k = 0 。结合(5)式、(7)式以及定义1可知, L m , i ( α , β ) ( x m , k ( α , β ) ) = δ i k ,故(6)式满足

I m u ( x m , i ( α , β ) ) = u ( x m , i ( α , β ) ) , i = 0 , 1 , , m , α , β > 1. (14)

(14)式表明,(6)式是以 { x m , j ( α , β ) } j = 0 m 为插值节点的Lagrange插值多项式。□

{ L m , i , B ( α , β ) ( x ) } i = 0 m 表示Lagrange插值基函数 { L m , i ( α , β ) ( x ) } i = 0 m 的重心形式,则实值函数 u ( x ) 在D上的m次重心Jacobi插值多项式表示为

I m , B u ( x ) = i = 0 m L m , i , B ( α , β ) ( x ) u ( x m , i ( α , β ) ) , α , β > 1 , x D , (15)

其中

L m , i , B ( α , β ) ( x ) = σ m , i ( α , β ) x x m , i ( α , β ) ( j = 0 m σ m , j ( α , β ) x x m , j ( α , β ) ) 1 , i = 0 , 1 , , m , (16)

σ m , i ( α , β ) 是与定义1中 m + 1 次Jacobi多项式 J m + 1 ( α , β ) ( x ^ ) 的零点 x ^ m , i ( α , β ) 对应的重心插值权。

特别地,文献 [9] 中的重心Gegenbauer插值多项式是 α = β 时的重心Jacobi插值多项式。

定理2 [16] (16)式中,重心插值权 σ m , i ( α , β ) 的显式表达式为

σ m , i ( α , β ) = ( 1 ( x ^ m , i ( α , β ) ) 2 ) ϖ i , i = 0 , 1 , , m .

当m的取值较大时,该式中的 ( 1 ( x ^ m , i ( α , β ) ) 2 ) 会受到舍入误差的影响,对其引入正弦变换 x ^ m , i ( α , β ) = sin ( κ ) ( i = 0 , 1 , , m , κ [ π 2 , 3 π 2 ] ) ,可消除舍入误差的影响,此时

σ m , i ( α , β ) = ( 1 ) i ϖ i cos ( sin 1 ( x ^ m , i ( α , β ) ) ) , i = 0 , 1 , , m .

根据定理1以及(16)式可知,(15)式在数学意义上等价于(6)式。由于(15)式的插值节点是定义1中区间D上 m + 1 次Jacobi多项式 J ˜ m + 1 ( α , β ) ( x ) 的零点,与区间D上的Gauss-Jacobi求积节点(又称为移位Gauss-Jacobi节点)的含义相同 [13] ,所以(15)式是与区间D上以移位Gauss-Jacobi节点为插值节点的Jacobi插值多项式等价的重心Lagrange插值多项式。

2.2. 重心Jacobi插值配置法

配置法是一种全局数值方法 [14] ,它求解积分/微分方程的基本原理是运用全局插值函数逼近原方程的未知函数,并将原方程离散化。重心Jacobi插值配置法求解积分/微分方程的主要思路是首先运用满足(15)式的m次重心Jacobi插值多项式 I m , B u ( x ) 逼近原方程的未知函数,得到原方程的近似方程,然后令原方程在(15)式的插值节点 x m , j ( α , β ) ( j = 0 , 1 , , m ) 处精确成立,并应用求积公式得到近似方程的离散方程组,最后求解离散方程组,该离散方程组的解是原方程的未知函数在节点 x m , j ( α , β ) 处的近似值,将求得的数值结果代入(15)式,即可得到原方程的重心Jacobi插值配置解。

3. 算法构造

以下详细阐述重心Jacobi配置法求解Fredholm积分–微分方程(1)的算法步骤。

第1步将原积分–微分方程改写为积分方程

结合初值条件,积分–微分方程(1)可变形为

u ( x ) + a x u ( t ) d t + a x a b K ( t 0 , t ) u ( t ) d t d t 0 = a x f ( t ) d t + u ( a ) , x D . (17)

第2步构造积分方程的近似方程

以(15)式近似积分方程(17)中的 u ( x ) ,可得方程(17)的近似方程

I m , B u ( x ) + a x I m , B u ( t ) d t + a x a b K ( t 0 , t ) I m , B u ( t ) d t d t 0 = a x f ( t ) d t + u ( a ) . (18)

第3步将近似方程离散化

令方程(18)在节点 x m , j ( α , β ) ( j = 0 , 1 , , m ) 处精确成立,则

I m , B u ( x m , j ( α , β ) ) + a x m , j ( α , β ) I m , B u ( t ) d t + a x m , j ( α , β ) a b K ( t 0 , t ) I m , B u ( t ) d t d t 0 = a x m , j ( α , β ) f ( t ) d t + u ( a ) , j = 0,1, , m . (19)

第4步近似计算离散近似方程中的积分,并得到离散方程组

采用 M ( M + ) 个求积点的Gauss-Legendre求积公式计算方程(19)中的第二项与第三项,结合定理1,(19)式可转化为如下离散方程组

u ˜ ( x m , j ( α , β ) ) + i = 0 m [ p i , j , B ( 1 ) + p ^ i , j , B ( 1 ) ] u ˜ ( x m , i ( α , β ) ) = a x m , j ( α , β ) f ( t ) d t + u ( a ) , j = 0 , 1 , , m , (20)

其中, u ˜ ( x m , j ( α , β ) ) u ( x ) 在节点 x m , j ( α , β ) 的近似值,且

p i , j , B ( 1 ) = ( x m , j ( α , β ) a ) 2 q = 1 M ω M , q ( 0 , 0 ) L m , i , B ( α , β ) ( x ^ M , q ( 0 , 0 ) ; a , x m , j ( α , β ) ) ,

p ^ ( 1 ) i , j , B = d 1 2 q = 1 M [ a x m , j ( α , β ) K ( t 0 , x ^ M , q ( 0 , 0 ) ) d t 0 ] ω M , q ( 0 , 0 ) i = 0 m L m , i , B ( α , β ) ( x M , q ( 0 , 0 ) ) , i , j = 0 , 1 , , m ,

其中, { x ^ M , q ( 0 , 0 ) } q = 1 M 是Gauss-Legendre求积点, { ω M , q ( 0 , 0 ) } q = 1 M 是Gauss-Legendre求积权。

第5步将离散方程组改写为矩阵形式,得到线性方程组

U = [ u ˜ ( x m , 0 ( α , β ) ) u ˜ ( x m , m ( α , β ) ) ] T , P B = ( p i , j , B ( 1 ) + p ^ ( 1 ) i , j , B ) i , j = 0 m ,

F B = [ a x m , 0 ( α , β ) f ( t ) d t + u ( a ) a x m , m ( α , β ) f ( t ) d t + u ( a ) ] T ,

则方程(20)的矩阵形式为

U P B U = F B , (21)

(21)式是关于 U 的线性方程组。

第6步回代,得到原方程的数值解

求解线性方程组(21)并得到 u ˜ ( x m , j ( α , β ) ) ( j = 0 , 1 , , m ) ,将其代入(15)式,可得方程(1)的重心Jacobi插值配置解为

u ( x ) I m , B u ˜ ( x ) = i = 0 m L m , i , B ( α , β ) ( x ) u ˜ ( x m , i ( α , β ) ) , α , β > 1 , x D . (22)

4. 误差估计

本节在 χ ¯ = ( C ( D ) , ) (D上全体连续函数组成的完备空间)中对算法进行误差估计,其中,对任意 u ( x ) C ( D ) ,范数 定义为 u ( x ) = max x D | u ( x ) | 。另外,在本节中,记号 Ω 与D分别表示区间 [ 1 , 1 ] 与区间 [ a , b ] ( a , b )

将(22)式代入方程(17),并设余项为 E m ( x ) ,则有

I m , B u ˜ ( x ) + a x I m , B u ˜ ( t ) d t + a x a b K ( t 0 , t ) I m , B u ˜ ( t ) d t d t 0 = a x f ( t ) d t + u ( a ) + E m ( x ) , (23)

以(17)式减(23)式,可得

E m ( x ) = u ( x ) I m , B u ˜ ( x ) + a x ( u ( t ) I m , B u ˜ ( t ) ) d t + a x a b K ( t 0 , t ) ( u ( t ) I m , B u ˜ ( t ) ) d t d t 0 . (24)

引理1 [17] 设 f ( x ) C ( Ω ) , n + ,Gauss-Legendre求积权为 ω n , k ( 0 , 0 ) ( k = 0 , 1 , , n ) ,Gauss-Legendre求积节点为 x ^ n , k ( 0 , 0 ) ,若 R n [ f ] 为求积公式 1 1 f ( x ) d x = k = 0 n ω n , k ( 0 , 0 ) f ( x ^ n , k ( 0 , 0 ) ) 的余项,则

R n [ f ] = 2 2 n + 3 [ ( n + 1 ) ! ] 4 ( 2 n + 3 ) [ ( 2 n + 2 ) ! ] 3 f ( 2 n + 2 ) ( η ) , η ( 1 , 1 ) .

推论1 若被积函数 f ( x ) 在积分区间 Ω 内存在至多 2 n + 1 阶导数,则 R n [ f ] = 0

由(16)式可看出,重心Jacobi插值基函数 L m , i , B ( α , β ) ( x ) ( i = 0 , 1 , , m ) 是m次多项式,因此结合推论1可知,采用不小于 ( m + 1 ) / 2 个求积点的Gauss-Legendre求积公式计算方程(19)中的第二项与第三项时,相应的求积余项为0,此时, E m ( x ) 由插值误差主导。

定理3 设 S m = s p a n { L m , 0 , B ( α , β ) ( x ) , L m , 1 , B ( α , β ) ( x ) , , L m , m , B ( α , β ) ( x ) } ,若 I m , B 满足(15)式,且是从 C ( D ) S m 的插值投影算子,则 I m , B 是线性有界的,且满足 I m , B Λ m < ,其中

Λ m = sup u C ( D ) I m , B u u = sup x D i = 0 m | L m , i , B ( α , β ) ( x ) | .

证明 对任意 u , g C ( D ) , ε , δ , x D

I m , B ( ε u + δ g ) ( x ) = i = 0 m L m , i , B ( α , β ) ( x ) [ ε u ( x m , i ( α , β ) ) + δ g ( x m , i ( α , β ) ) ] = ε i = 0 m L m , i , B ( α , β ) ( x ) u ( x m , i ( α , β ) ) + δ i = 0 m L m , i , B ( α , β ) ( x ) g ( x m , i ( α , β ) ) = ε I m , B u ( x ) + δ I m , B g ( x ) ,

因此 I m , B 是一个线性算子。由于

| I m , B u ( x ) | = | i = 0 m L m , i , B ( α , β ) ( x ) u ( x m , i ( α , β ) ) | i = 0 m | L m , i , B ( α , β ) ( x ) | | u ( x m , i ( α , β ) ) | Λ m u ,

因此 I m , B Λ m < 。综上, I m , B 是一个线性有界算子。□

定理3中的 Λ m 是与节点 { x m , j ( α , β ) } j = 0 m 对应的Lebesgue常数。结合(12)式,将 Λ m 中的变量 x D 运用变量代换转化为 x ^ Ω ,可以发现

Λ m = sup x D i = 0 m | L m , i , B ( α , β ) ( x ) | = sup x Ω i = 0 m | L m , i , B ( α , β ) ( x ^ ) | ,

该式表明, Λ m 也是与节点 { x ^ m , j ( α , β ) } j = 0 m 对应的Lebesgue常数。此外,文献 [16] 表明, Λ m 有如下估计式成立

Λ m = O ( m max { α , β } + 1 / 2 ) , (25)

其中, O ( m max { α , β } + 1 / 2 ) 表示与 m max { α , β } + 1 / 2 同阶的量。

引理2 [18] 设存在 p ( p + ) ,使得 d p u ( x ) / d x = u ( p ) ( x ^ ) C ( Ω ) 。令 B m u ( x ^ ) u ( x ^ ) 的m次最佳逼近多项式,则对任意 p m + 1 ,有

B m u ( x ^ ) u ( x ^ ) ( π / 2 ) u ( p ) ( x ^ ) ( m p + 3 ) p .

定理4 (误差界)设方程(1)满足条件

K ( x , t ) = max x , t D | K ( x , t ) | h (h为常量),

精确解 u ( x ) C ( p ) ( D ) ( p + )

以(17)式~(22)式求解方程(1),设 u ˜ ( x m , i ( α , β ) ) ( i = 0 , 1 , , m ) 是方程组(21)的解, u ( x m , i ( α , β ) ) 为方程(1)在节点 x m , i ( α , β ) 的精确值,若 M = m + 1 2 p m ,则对任意 i = 0 , 1 , , m ,有

E m ( x ) [ 1 + ( b a ) + ( b a ) 2 h ] ( u ( x ) I m , B u ( x ) + Λ m e i ) = O ( m 2 max { α , β } + 1 p ) ,

其中, e i = u ( x m , i ( α , β ) ) u ˜ ( x m , i ( α , β ) ) α , β 是Jacobi多项式中的参数且满足 α , β > 1

证明 由(24)式以及范数的三角形不等式,可得

E m ( x ) = u ( x ) I m , B u ˜ ( x ) + a x ( u ( t ) I m , B u ˜ ( t ) ) d t + a x a b K ( t 0 , t ) ( u ( t ) I m , B u ˜ ( t ) ) d t d t 0 u ( x ) I m , B u ˜ ( x ) + a x ( u ( t ) I m , B u ˜ ( t ) ) d t + a x a b K ( t 0 , t ) ( u ( t ) I m , B u ˜ ( t ) ) d t d t 0 ,

由于 x D ,故

a x ( u ( t ) I m , B u ˜ ( t ) ) d t ( b a ) u ( x ) I m , B u ˜ ( x ) ,

a x a b K ( t 0 , t ) ( u ( t ) I m , B u ˜ ( t ) ) d t d t 0 [ ( b a ) 2 h ] u ( x ) I m , B u ˜ ( x ) ,

从而

E m ( x ) [ 1 + ( b a ) + ( b a ) 2 h ] u ( x ) I m , B u ˜ ( x ) [ 1 + ( b a ) + ( b a ) 2 h ] ( u ( x ) I m , B u ( x ) + I m , B u ( x ) I m , B u ˜ ( x ) ) .

对任意 x ^ Ω ,有

| u ( x ^ ) I m , B u ( x ^ ) | | u ( x ^ ) I m , B ( B m u ( x ^ ) ) | + | I m , B ( B m u ( x ^ ) ) I m , B u ( x ^ ) | = | u ( x ^ ) B m u ( x ^ ) | + | I m , B ( B m u u ) ( x ^ ) | ,

结合定理3与引理2,可得

| u ( x ^ ) I m , B u ( x ^ ) | ( 1 + Λ m ) u ( x ^ ) B m u ( x ^ ) ( 1 + Λ m ) ( π / 2 ) u ( p ) ( x ^ ) ( m p + 3 ) p .

对该式运用变量代换 x = b a 2 x ^ + a + b 2 ,结合(25)式,可得 u ( x ) 的重心Jacobi插值误差估计式为

u ( x ) I m , B u ( x ) ( b a ) ( 1 + Λ m ) 2 ( π / 2 ) u ( p ) ( x ) ( m p + 3 ) p = O ( m max { α , β } + 1 / 2 p ) . (26)

由(22)式,对任意 i = 0 , 1 , , m ,有

I m , B u ( x ) I m , B u ˜ ( x ) = i = 0 m L m , i , B ( α , β ) ( x ) [ u ( x m , i ( α , β ) ) u ˜ ( x m , i ( α , β ) ) ] Λ m e i .

因为当采用 M = ( m + 1 ) / 2 个求积点的Gauss-Legendre求积公式近似积分时,相应的求积余项为0,所以近似方程与原方程之间的误差主要是未知函数的重心Jacobi插值误差,从而结合(26)式可知 e i 满足

e i u ( x ) I m , B u ( x ) O ( m max { α , β } + 1 / 2 p ) ,

于是

I m , B u ( x ) I m , B u ˜ ( x ) Λ m O ( m max { α , β } + 1 / 2 p ) = O ( m 2 max { α , β } + 1 p ) ,

因此

E m ( x ) [ 1 + d 1 + ( d 1 ) 2 h ] [ O ( m max { α , β } + 1 / 2 p ) + O ( m 2 max { α , β } + 1 p ) ] = O ( m 2 max { α , β } + 1 p ) .

根据定理4,可得到以下推论

推论2 (收敛性)若 p > 2 max { α , β } + 1 ( p + , α , β > 1 ) ,则当 m 时, E m ( x ) 呈现指数收敛的特点:按m的 2 max { α , β } + 1 p 次幂收敛到0;此外,若 u ( x ) 在D上充分光滑,且 1 < α , β ,则当 m 时, E m ( x ) 0

证明 由定理4,可知 E m ( x ) O ( m 2 max { α , β } + 1 p ) ,因此当 p > 2 max { α , β } + 1 时,有

2 max { α , β } + 1 p < 0 , (27)

因此当 m 时, m 2 max { α , β } + 1 p 0 ,从而 E m ( x ) 0 ,且 E m ( x ) 收敛到零的速率与m的 1 p + 2 max { α , β } 次幂收敛到零的速率相同,故此时 E m ( x ) 呈现指数收敛的特点。

u ( x ) 在D上充分光滑,则 p ,若 1 < α , β ,则 p α , β ,从而仍有(27)式成立,故此时当 m 时,总有 E m ( x ) 0 成立。□

推论3 当 m 时,若 E m ( x ) 0 ,则对任意 α , β > 1 E m ( x ) 收敛到0的速率随 α , β 中最大值的增大而减小。

证明 由推论2可知,当 E m ( x ) 0 时,有 2 max { α , β } + 1 p < 0 成立,此时 E m ( x ) 的收敛速率与 m 2 max { α , β } + 1 p 的收敛速率相同。固定p的取值不变,当 max { α , β } 逐渐增大时, 2 max { α , β } + 1 p 逐渐增大, m 2 max { α , β } + 1 p 的收敛速率逐渐减小,从而 E m ( x ) 的收敛速率逐渐减小。□

5. 数值算例

本节所给算例已知精确解,实验操作系统为MATLAB R2018a,机器精度 e p s = 2.220 × 10 16 。数值实验中使用double类型的变量进行计算,并将得到的数值结果保留两位小数。设本文方法求得的数值解为 I m , B u ˜ ( x ) 且满足(22)式,其中m表示(22)式中数值解 I m , B u ˜ ( x ) 的次数,并将算例中的精确解与数值解在计算节点 { x i } i = 0 N D ( N + ) 的最大绝对误差可记为

ψ m = max 0 i N | u ( x i ) I m , B u ˜ ( x i ) | , x i D , i = 0 , 1 , , N . (28)

在算例中,通过取不同的 α , β 并根据(28)式计算 ψ m 以体现Jacobi多项式中参数 α , β 的取值对最大绝对误差的影响。

例 [3] [19] [20] 考虑一阶线性Fredholm积分–微分方程初值问题

{ u ( 1 ) ( x ) = u ( x ) + 0 1 x ( 1 + t ) ( ln 2 ) 2 u ( t ) d t + f ( x ) , x [ 0 , 1 ] , u ( 0 ) = 0 , (29)

其中, f ( x ) = x / 3 + 1 / ( 1 + x ) ln ( 1 + x ) ,方程(29)的精确解为 u ( x ) = ln ( 1 + x )

运用(17)式~(20)式,将方程(29)离散为线性方程组,求解离散线性方程组并得到数值解,以(28)式计算精确解与数值解在计算节点的绝对误差与最大绝对误差。图1是本文方法在 m = 15 , α = 0.5 , β = 0.6 时求得的数值解与方程(29)的精确解在等距节点的绝对误差分布图,表1比较了精确解与数值解在等距节点的计算结果。从图1表1可看出,本文方法对方程(29)的求解有效。当计算节点为 [ 0 , 1 ] 上的等距节点时,表2列举了本文方法在 m = 2 n 1 ( n = 1 , , 5 ) α = 0.5 , β = 0.4 : 0.4 : 2 时的数值结果;在m相同的条件下,表3列举了本文方法在 β = 0.4 : 0.4 : 2 , α = 0.5 时的数值结果。此两表中数据显示,当m相对较小时,随着m的增大,本文方法的收敛速率逐渐变大,且此时,随 α , β 中最大值的不断增大,本文方法的收敛速率不断减小。总之,通过大量的数值实验可以发现,对本例,最大绝对误差随m的增大而减小,且此时本文方法的收敛速率随 α , β 中最大值的增大而减小,但当m增大到某个范围内的值后,最大绝对误差不会再随m的增大减小,而是保持相对稳定,且此时本文方法的收敛速率的随 α , β 取值的变化规律并不明显。

Figure 1. The absolute error distribution of the presented method at equidistant nodes for α = 0.5 , β = 0.6 , m = 15

图1. 当 α = 0.5 , β = 0.6 , m = 15 时,本文方法在等距节点的绝对误差分布

由(22)式可知,本文方法得到的数值解的展开式项数为 m + 1 ,与文献 [3] 方法中构造的同伦级数解的展开式的项数相同。当 α = 0.2 , β = 0.5 时,在展开式项数以及计算节点类型、个数均相同的条件下,将本文方法所求的最大绝对误差与文献 [3] 方法比较,结果见表4。从表4可看出,与文献 [3] 方法相比,本文方法求得的最大绝对误差更小。

Table 1. The calculation result of the presented method at equidistant nodes for α = − 0.5 ,   β = − 0.6 ,   m = 15

表1. 当 α = 0.5 , β = 0.6 , m = 15 时,本文方法在等距节点的计算结果

文献 [20] 在计算方程(29)的数值解时,生成一个维数为 m 3 = ( 2 m 1 + 1 ) ( m 2 + 2 ) 的线性方程组,其中 m 1 m 2 分别表示Shannon缩放函数的项数与母波展开式的项数;本文方法在计算方程(29)的数值解时,由(20)式可知,生成一个 m 4 = m + 1 维的线性代数方程组。表5列举了本文方法在 α = 0.2 , β = 0.5 时所求的数值解在Sinc节点 [19] 的最大绝对误差,并比较了本文方法与文献 [20] 方法的最大绝对误差为同一数量级时需计算的线性方程组的维数。从表5可看出,本文方法能更精确地计算方程(29)的数值解,且当这两种方法的最大绝对误差同阶时,本文方法需要计算的线性代数方程组维数更低,因此本文方法在计算时占用的内存更少,计算量更小。

Table 2. The calculation result of the presented method at different β ,   n for α = − 0.5

表2. 当 α = 0.5 时,本文方法基于不同 β , n 的数值结果

Table 3. The calculation result of the presented method at different α ,   n for β = − 0.5

表3. 当 β = 0.5 时,本文方法基于不同 α , n 的数值结果

Table 4. Numerical results of the method presented in this paper and the method presented in reference [3] for α = − 0.2 ,   β = − 0.5

表4. 当 α = 0.2 , β = 0.5 时,本文方法与文献 [3] 方法的数值结果

Table 5. Numerical results of the method presented in this paper and the method presented in reference [20] for α = 0.2 ,   β = 0.5

表5. 当 α = 0.2 , β = 0.5 时,本文方法与文献 [20] 方法的数值结果

6. 总结

本文探究了文献 [12] 中的重心Gegenbauer插值的一般形式——重心Jacobi插值,表明了它本质上是与插值节点为移位Gauss-Jacobi节点的Jacobi插值等价的重心插值,并基于重心Jacobi插值多项式,探究了一类一阶线性Fredholm型积分–微分方程初值问题的重心Jacobi插值配置法,给出了相应的误差估计式。数值算例中的数据验证了本文方法的有效性,并体现出当原方程未知函数的重心Jacobi插值多项式的次数在一定范围内变化时,本文方法求得的最大绝对误差随Jacobi多项式中参数 α , β 最大值的增大而增大。

基金项目

自然科学青年基金,项目批准号:11801456,项目名称:超奇异积分及其积分方程的数值算法和应用。

参考文献

[1] 李星. 积分方程[M]. 北京: 科学出版社, 2008: 3-22.
[2] 云天铨. 积分方程及其在力学中的作用[M]. 广州: 华南理工大学出版社, 1990: 402-417.
[3] Jaradat, H., Alsayyed, O. and Alshara, S. (2008) Numerical Solution of Linear Integro-Differential Equations. Journal of Mathematics and Statistics, 4, 250-254.
[4] Tair, B., Guebbai, H., Segni, S., et al. (2022) An Approximation Solution of Linear Fredholm Integro-Differential Equation Using Collocation and Kantorovich Methods. Journal of Applied Mathematics and Computing, 68, 3505-3525.
https://doi.org/10.1007/s12190-021-01654-2
[5] Fathy, M., El-Gamel, M. and El-Azab, M.S. (2014) Legendre-Galerkin Method for the Linear Fredholm Integro-Differential Equations. Applied Mathematics and Computation, 243, 789-800.
https://doi.org/10.1016/j.amc.2014.06.057
[6] Tair, B., Guebbai, H., Segni, S., et al. (2021) Solving Linear Fredholm Integro-Differential Equation by Nyström Method. Journal of Applied Mathematics and Computational Mechanics, 20, 53-64.
https://doi.org/10.17512/jamcm.2021.3.05
[7] Şahina, N., Yüzbaşi, Ş. and Sezer, M. (2011) A Bessel Polynomial Approach for Solving General Linear Fredholm Integro-Differential-Difference Equations. International Journal of Computer Mathematics, 89, 3093-3111.
https://doi.org/10.1080/00207160.2011.584973
[8] Yalinba, S., Sezer, M. and Sorkun, H.H. (2009) Legendre Polynomial Solutions of High-Order Linear Fredholm Integro-Differential Equations. Applied Mathematics and Computation, 210, 334-349.
https://doi.org/10.1016/j.amc.2008.12.090
[9] Maadadi, A. and Rahmoune, A. (2018) Numerical Solution of Nonlinear Fredholm Integro-Differential Equations Using Chebyshev Polynomials. International Journal of Advanced Scientific and Technical Research, 4, 85-91.
https://doi.org/10.26808/rs.st.i8v4.09
[10] Elgindy, K.T. and Smith-Miles, K.A. (2013) Solving Boundary Value Problems, Integral and Integro-Differential Equations Using Gegenbauer Integration Matrices. Journal of Computational and Applied Mathematics, 237, 307-325.
https://doi.org/10.1016/j.cam.2012.05.024
[11] Elgindy, K.T. (2016) High-Order Numerical Solution of Second-Order One-Dimensional Hyperbolic Telegraph Equation Using a Shifted Gegenbauer Pseudospectral Method. Numerical Methods for Partial Differential Equations, 32, 307-349.
https://doi.org/10.1002/num.21996
[12] Elgindy, K.T. (2017) High-Order, Stable, and Efficient Pseudo Spectral Method Using Barycentric Gegenbauer Quadratures. Applied Numerical Mathematics, 113, 1-25.
https://doi.org/10.1016/j.apnum.2016.10.014
[13] Hesthaven, J.S., Gottlieb, S. and Gottlieb, D. (2007) Spectral Methods for Time-Dependent Problems. Vol. 21 of Cambridge Monographs on the Applied and Computational Mathematics, Cambridge University Press, Cambridge, 88-89.
https://doi.org/10.1017/CBO9780511618352
[14] Shen, J., Tang, T. and Wang, L.L. (2011) Spectral Methods: Algorithms, Analysis and Applications. Springer, Berlin.
https://doi.org/10.1007/978-3-540-71041-7
[15] Mama, F. (2018) An Introduction to Orthogonal Polynomials. 2nd AIMS-Volkswagen Stiftung Workshop, Douala, 5-12 October 2018, 3-24.
[16] Wang, H., Daan, H. and Stefn, V. (2012) Explicit Barycentric Weights for Polynomial Interpolation in the Roots or Extrema of Classical Orthogonal Polynomials. Mathematics, 83, 2893-2914.
https://doi.org/10.1090/S0025-5718-2014-02821-4
[17] 李庆扬, 王能超, 易大义. 数值分析[M]. 北京: 清华大学出版, 2008: 116-123.
[18] Liang, F. and Lin, F.R. (2010) A Fast Numerical Solution Method for Two Dimensional Fredholm Integral Equations of the Second Kind Based on Piecewise Polynomial Interpolation. Applied Mathematics and Computational, 216, 3073-3088.
https://doi.org/10.1016/j.amc.2010.04.027
[19] Rashidinia, J. and Zarebnia, M. (2007) The Numerical Solution of Integro-Differential Equation by Means of the Sinc Method. Applied Mathematics and Computation, 188, 1124-1130.
https://doi.org/10.1016/j.amc.2006.10.063
[20] Maleknejad, K. and Attaryc, M. (2011) An Efficient Numerical Approximation for Linear Class of Fredhom Integro-Differential Equations Based on Cattani’s Method. Communication Nonlinear Science and Numerical Simulation, 16, 2672-2679.
https://doi.org/10.1016/j.cnsns.2010.09.037