PDF《数学杂志》

Vol.

36 (

2016 ) No.

2 数学杂志J. of Math. (PRC) Volterra 积分方程的蒙特卡罗数值求解方法 洪志敏, 闫在在 (内蒙古工业大学理学院, 内蒙古 呼和浩特 010051) 摘要: 本文讨论了第一类、 第二类以及具有奇异核的 Volterra 积分方程的数值解问题. 利用重 要抽样蒙特卡罗随机模拟方法获得积分方程解的近似计算结果. 通过对文献中算例的实现表明文中所 提方法扩展了 Volterra 型积分方程的数值求解方法, 关键词: 线性 Volterra 积分方程;

配置法;

马尔科夫链;

蒙特卡罗算法;

重要抽样 MR(2010) 主题分类号: 65R20;

65C05 中图分类号: O175.5 文献标识码: A 文章编号: 0255-7797(2016)02-0425-12

1 引言 积分方程是含有对未知函数积分运算的方程, 其中未知函数是以线性形式出现的, 称为 线性积分方程, 否则称为非线性积分方程. 在工程、 力学、 地球物理勘探等方面许多数学物理 问题都需要通过积分方程求解. 同一个问题有时既可以用微分方程的定解问题又可以用积分 方程来描述, 而常微分方程与偏微分方程的定解问题也可化为等价的积分方程. 积分方程这 一数学工具正日益受到重视, 这是因为化为积分方程可以降低维数, 减少所用节点的个数, 缩 短计算时间, 节省费用. Volterra 积分方程常出现在许多物理力学应用中. 例如, 在弹簧振子的阻尼振动方程与初始条件为 ? + 2β? + ω2 0? = 0, ?(0) = ?0, ? (0) = 0, 可化为求解位移 ?(t) 所满足的积分方程 ?(t) = ?0 cos(ω0t) + 2β?0 ω0 sin(ω0t) ? 2β t

0 cos(ω0(t ? s))?(s)ds, 这是第二类线性 Volterra 积分方程. 下面是一个与土力学有关的 Volterra 积分方程的例子. 设一个桩子受横向力 P 作用于 桩顶. 如图

1 所示, 一桩子顶部受一横向力 P 的作用, p(x) 是土的反力且未知, 在截面 x 处 的弯矩 M(x) 为M(x) = P(e + x) ? x

0 p(t)(x ? t)dt, (1.1) ? 收稿日期: 2013-04-26 接收日期: 2014-03-18 基金项目:国家自然科学基金资助项目(11161031);

内蒙古自治区自然科学基金研究项目(2010MS0116);

高等学校博士学科点专项科研基金联合资助课题 (20131514110005). 作者简介: 洪志敏 (1975C), 女, 蒙古族, 内蒙古赤峰, 副教授, 主要研究方向: 力学问题中的数学方法.

426 数学杂志Vol.

36 图1: 桩受力图 令F(x) = P(e + x) ? M(x), 则F(x) 是已知函数. 因为力 P 和距离 e 是已知的, 而弯矩 M(x) 可以通过预先埋设于桩表面的应变片用电测法测出应变再换成应力便可求得, 因而也 是已知函数. 因此, 式(1.1) 可写成 F(x) = x

0 p(t)k(x, t)dt, (1.2) 其中 k(x, t) = x ? t 为已知的核函数. 式(1.2) 就是一个 Volterra 第一种积分方程. 求出未知 函数 p(x) 后, 便可以用来测定土壤的一些系数 [1] . 考虑求解如下两种类型积分方程 g(x) = λ x

0 k(x, t)?(t)dt,

0 ≤ x ≤ 1, (1.3) ?(x) = g(x) + λ x

0 k(x, t)?(t)dt,

0 ≤ x ≤ 1, (1.4) 其中方程 (1.3) 为Volterra 第一类积分方程, 方程 (1.4) 为Volterra 第二类积分方程. 假设 g(x) ∈ L2[0, 1], k(x, t) ∈ L2[0, 1] * L2[0, 1] 是已知的实值连续函数, ?(x) ∈ L2[0, 1] 是待求的 未知解函数, λ 是已知的常数. 这里我们假设 Volterra 第

一、二类积分方程 (1.3), (1.4) 存在 解析解. 一些求解积分方程 (1.3) 和(1.4) 的数值解法被使用, 例如 Galerkin 法, 配置法, Taylor 级数展开法, Bernstein 多项式逼近法, Adomian 分解法以及求积公式法 [2C7] . 一般情况下, 这些数值计算方法将积分方程转换成可以直接求解或通过迭代方法求解的线性代数系统, 进 而求得积分方程的数值近似解. 然而, 多数情况下, 所得线性系统的系数矩阵并非稀疏矩阵, 尤其当离散化所得的线性系统是高维情形时, 利用这些数值计算方法求解积分方程就会很困 难而且计算量也很大. 利用蒙特卡罗方法求解 Fredholm 第二类积分方程在文献 [8] 中被提出, 而且, 蒙特卡罗 方法作为一种数值计算方法, 具有受几何条件限制小, 收敛速度与问题的维数无关, 误差容易 确定, 程序结构简单、 灵活、 易于实现的优点. 本文结合确定性求解技术和蒙特卡罗方法的优 点, 推广蒙特卡罗方法在求解第一类、 二类 Volterra 积分方程 (1.3), (1.4) 中的应用. No.

2 洪志敏等: Volterra 积分方程的蒙特卡罗数值求解方法

427 文中给出的方法是简单有效的, 可以为迭代法求解积分方程提供一个在其解范围内的初 始近似解. 通过模拟文献中的数值算例证明了提出方法的效率和精度.

2 待定系数逼近法 设{ui(x)}, (i = 0, 1, 2,是L2(0, 1) 中的完备函数系, 它们可以是正交的, 也可以是 非正交的, 但它们一定是线性无关的. 将积分方程的未知函数 ?(x) 按{ui(x)} 展开, 用有限 项和式作为 ?(x) 的近似, 即?(x) ≈ n i=0 aiui(x), (2.1) 只要确定了系数 ai, 积分方程的近似解 ? ?(x) 就可以确定. 以下我们以 Volterra 第二类积分方程 (1.4) 为例阐明对 Volterra 积分方程的求解方法. 将式 (2.1) 代入积分方程 (1.4), 可得 n i=0 aiui(x) ? λ n i=0 ai x

0 k(x, t)ui(t)dt ? g(x) = r(x), (2.2) 这里如果对任意的 x ∈ [0, 1], 误差 r(x) ≡ 0, 则通过式 (2.1) 就得到积分方程 (1.4) 的精确解, 然而, 这通常难以做到. 因此, 为得到积分方程 (1.4) 的近似解, 可以放宽对误差 r(x) 的要求, 根据对误差 r(x) 的不同要求可以导致确定展开系数 ai(i = 0, 1, 2,的不同数值计算方法. 本文采用配置法将积分方程转化成线性代数系统. 此时, 要求式 (2.2) 中的误差 r(x) 满足(r(x), vi(x)) =

0 (i = 0, 1, 2,n), (2.3) 其中 vi(x) = δ(x ? xi). 这里我们选取多项式 pi(x) 作为 [0, 1] 上的完备函数系 {ui(x)}(i = 0, 1, 2,设积分 方程 (1.3) 和(1.4) 的未知函数 ?(x) 是定义在区间 [0, 1] 上的连续函数. 定理 2.1 设?(x) ∈ C[0, 1], 则一定存在多项式 pn(x) ∈ Pn 使得 max 0≤x≤1 |?(x) ? pn(x)| <

. (2.4) 证 参见文献 [9]. 根据定理 2.1 以及矩量法, 选取函数 ?(x) 的逼近多项式 pn(x) 为如下形式 Pn(?(x)) = n i=0 aixi (2.5) 且满足 max 0≤x≤1 |?(x) ? Pn(?(x)2.6) 将式 (2.5) 代回积分方程 (1.4) 中得 n i=0 aixi = g(x) + λ n i=0 ai x

0 k(x, t)ti dt, (2.7)

428 数学杂志Vol.

36 此时, 式(2.7) 中的积分项已可以计算, 未知量是 ai (i = 0, 1,n), 为求这 n +

1 个未知的 系数 ai, 可以建立如下线性系统 n i=0 aixi j = g(xj) + λ n i=0 ai xj

0 k(xj, t)ti dt, (2.8) 针对核函数 k(x, t) 可能会存在奇异点的问题, 取节点 xj = j/n + , j = 0, 1,n ? 1, xn =

1 ? , 其中 是任意小的正数. 线性系统 (2.8) 可以简写为 AX = b, (2.9) 其中 A = [xi j ? λ xj

0 k(xj, t)ti dt], i, j = 0, 1,n, X = [ai]t , i = 0, 1,n, b = [g(xj)]t , j = 0, 1,n. (2.10) 求解线性系统 (2.9), 可以得到积分方程 (1.4) 的近似解 Pn(?(x)) = n i=0 ? aixi , 其中 ? ai 是 由线性系统 (2.9) 求解所得. 所求近似解的误差限由以下定理给出. 引理 2.1 假设 A ? I

1 = c1 <

1 (2.11) 成立, 其中 I 为n+1阶单位阵. 则有 A?1

1 ≤

1 1 ? c1 . (2.12) 证 由已知条件 (2.11) 有矩阵 A 是非奇异矩阵, 令B=A?I, 则矩阵 I + B 可逆, 即(I + B)(I + B)?1 = I, 进而有 (I + B)?1 = I ? B(I + B)?1 . 在等号两端取矩阵的

1 范数有 (I + B)?1 = I ? B(I + B)?1 ≤ I + B (I + B)?1 , (I + B)?1 = A?1 ≤

1 1 ? B =

1 1 ? c1 . 定理 2.2 对于在区间 [0, 1] 上连续的函数 g(x) 以及连续核 k(x, t), 线性系统 (2.9) 中的 矩阵 A 是可逆的. 记sup x,t∈[0,1] |λk(x, t)| = c, 则有 sup xi∈[0,1] |?(xi) ? Pn(?(xi))| ≤

1 +

1 + c

1 ? c1 (2.13) 成立, 其中 xi = i/n, i = 0, 1,n, ?(x) 是方程 (1.4) 的精确解. 是一个任意小的正数. 证sup xi∈[0,1] |?(xi) ? Pn(?(xi))| = sup xi∈[0,1] |?(xi) ? Pn(?(xi)) + Pn(?(xi)) ? Pn(?(xi))| ≤ sup xi∈[0,1] |?(xi) ? Pn(?(xi))| + sup xi∈[0,1] |Pn(?(xi)) ? Pn(?(xi))|. (2.14) No.

2 洪志敏等: Volterra 积分方程的蒙特卡罗数值求解方法

429 下面确定 (2.14) 式中两项的界限. 由式 (2.6) 可得 sup xi∈[0,1] |?(xi) ? Pn(?(xi)2.15) 积分方程 (1.4) 可表示为 g(x) = Pn(?(x)) ? λ x

0 k(x, t)Pn(?(t))dt, (2.16) 若在式 (2.16) 中将 Pn(?(x)) 替换成 Pn(?(x)), 则g(x) 将被 g(x) 代替, 即g(x) = Pn(?(x)) ? λ x

0 k(x, t)Pn(?(t))dt, (2.17) 因此有 sup xi∈[0,1] |Pn(?(xi)) ? Pn(?(xi))| = sup xi∈[0,1] n j=0 (aj ? ? aj)xj i ≤ A?1

1 max xi∈[0,1] |g(xi) ? g(xi)|. (2.18) 接下来求 max |g(xi) ? g(xi)|, 为此令 g(x) = Pn(?(x)) ? λ x

0 k(x, t)Pn(?(t))dt, 由表达 式(1.4) 可得 sup xi∈[0,1] |g(xi) ? g(xi)| = sup xi∈[0,1] |?(xi) ? Pn(?(xi)) ? xi

0 λk(xi, t)[?(t) ? Pn(?(t))]dt| ≤ sup xi∈[0,1] |?(xi) ? Pn(?(xi))| + sup xi∈[0,1] xi

0 λk(xi, t)[?(t) ? Pn(?(t))]dt ≤ sup xi∈[0,1] |?(xi) ? Pn(?(xi))| +

1 0 sup xi∈[0,1] |λk(xi, t)| sup t∈[0,1] |?(t) ? Pn(?(t))| dt ≤ (1 + c). 将这个误差限代回式 (2.18) 中有 sup xi∈[0,1] |Pn(?(xi)) ? Pn(?(xi))| ≤ A?1

1 (1 + c). (2.19) 根据引理 2.1 以及表达式 (2.14), (2.15) 和(2.19) 有sup xi∈[0,1] |?(xi) ? Pn(?(xi))| ≤ +

1 + c

1 ? c1 , 此定理得证. 一般情况下, 线性系统 (2.9) 的系数矩阵 A 不是稀疏矩阵, 为此首先采用列选主元的三 角分解法将非奇异矩阵 A 分解为单位下三角矩阵和上三角矩阵的乘积, 即PA = LU, 其中 L 为单位下三角阵, U 为上三角阵, P 为排列阵. 求解线性系统 (2.9) 等价于求解两个具有稀疏 系数矩阵的线性系统 LY = Pb (2.20)

430 数学杂志Vol.

36 以及 UX = Y. (2.21) 通过这样处理可以提高线性系统的求解速度. 对于稀疏系统 (2.21) 利用雅阁比松弛迭代法进 行求解. 记S=I?DU, 其中 D = diag( γ u11 , γ u22 γ un+1,n+1 ), γ ∈ (0, 1] 是松弛因子, I 是单 位矩阵, 则线性方程组 (2.21) 可表示为如下迭代形式 X(k+1) = SX(k) + Y1, k = 0, 1,2.22) 其中 Y1 = DY , X(0) = Y1. 根据计算方法理论可得如下定理: 定理 2.3 如果矩阵 S 的谱半径 ρ(S) 满足条件 ρ(S) <

1, (2.23) 则迭代形式 (2.22) 是收敛的, 即lim k→∞ X(k) = X. 证 参见文献 [10]. 利用待定系数逼近法求解积分方程 (1.4), 根据 Weierstrass 第一逼近定理, 要提高近似 解的精度, 需要增加逼近多项式 pn(x) 中的项数 n, 因而会面临求解阶数较高的线性代数系统 的困难. 而且, 在利用确定性方法求解线性系统时, 求未知解向量中一个分量时常需求出解向 量的全部分量. 这在一些实际问题中无疑会导致巨大的浪费. 蒙特卡罗数值求解方法不仅可 以处理维数较高的问题 (因为蒙特卡罗方法的收敛速度与问题的维数无关), 而且可以只对未 知函数在某个点处的值或几个点处的线性组合值进行求解, 大大减少了数值求解的计算量. 下一节将利用蒙特卡罗方法求解线性方程组的迭代形式 (2.22). 蒙特卡罗方法的基本思 想是将所求问题的解视为一个随机变量的数学期望, 将随机变量多次实现的平均值作为所求 问题的近似解.

3 重要抽样蒙特卡罗方法 对于收敛的迭代格式 (2.22), 有X=Y1 + SY1 + S2 Y1 Sk Y1 若令 sij ∈ S, Xi 是解向量 X 的第 i 个分量, 则有 Xi = Y1i + n+1 j1=1 sij1 Y1j1 + n+1 j1,j2=1 sij1 sj1j2 Y1j2 n+1 j1,j2,・・・ ,jk=1 sij1 sj1j2 ・ ・ ・ sjk?1jk Y1jk (3.1) 其中 Y1i 表示向量 Y1 的第 i 个分量, i = 1, 2,n + 1. 在式 (3.1) 中, 令sij = wijpij, 则式 (3.1) 可表示为 Xi = Y1i + n+1 j1=1 wij1 pij1 Y1j1 +・ ・ ・+ n+1 j1,j2,・・・ ,jk=1 wij1 wj1j2 ・ ・ ・ wjk?1jk pij1 pj1j2 ・ ・ ・ pjk?1jk Y1jk +・ ・ ・ . (3.2) No.

2 洪志敏等: Volterra 积分方程的蒙特卡罗数值求解方法

431 考虑如下的马尔科夫链 T = t0 → t1 → t2 tk 状态空间为 {1, 2,n + 1}. 记初始概率及由 i 状态转移到 j 状态的转移概率分别为 P(t0 = i) = pi (pi ≥ 0), P(tk = j|tk?1 = i) = pij (pij ≥ 0), i, j = 1, 2,n +

1 且pi 和pij 满足 n+1 i=1 pi = 1, n+1 j=1 pij = 1........