积分方程是含有对未知函数积分运算的方程,其中未知函数是以线性形式出现的, 称为线性积分方程,否则称为非线性积分方程.在工程、力学、地球物理勘探等方面许多数学物理问题都需要通过积分方程求解.同一个问题有时既可以用微分方程的定解问题又可以用积分方程来描述, 而常微分方程与偏微分方程的定解问题也可化为等价的积分方程.积分方程这一数学工具正日益受到重视,这是因为化为积分方程可以降低维数, 减少所用节点的个数,缩短计算时间, 节省费用.

Volterra积分方程常出现在许多物理力学应用中.

例如, 在弹簧振子的阻尼振动方程与初始条件为

$ \left\{\begin{array}{ll} \varphi''+2\beta \varphi'+\omega_0^2\varphi=0,\\ \varphi(0)=\varphi_0, \varphi'(0)=0, \end{array}\right. $

可化为求解位移$\varphi(t)$所满足的积分方程

$\varphi(t)=\varphi_0\cos(\omega_0t)+\frac{2\beta \varphi_0}{\omega_0}\sin(\omega_0t)-2\beta\int_0^t\cos(\omega_0(t-s))\varphi(s)ds,$

这是第二类线性Volterra积分方程.

下面是一个与土力学有关的Volterra积分方程的例子. 设一个桩子受横向力$P$作用于桩顶. 如图 1所示,一桩子顶部受一横向力$P$的作用, $p(x)$是土的反力且未知, 在截面$x$处的弯矩$M(x)$为

$M(x)=P(e+x)-\int_0^xp(t)(x-t)dt,$

(1.1)

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

$F(x)=\int_0^xp(t)k(x, t)dt,$

(1.2)

其中$k(x, t)=x-t$为已知的核函数. 式(1.2)就是一个Volterra第一种积分方程. 求出未知函数$p(x)$ 后, 便可以用来测定土壤的一些系数^{[ 1]}.

考虑求解如下两种类型积分方程

$g(x)=\lambda\int_0^xk(x,t)\varphi(t)dt, 0\leq x\leq 1,$

(1.3)

$\varphi(x)=g(x)+\lambda\int_0^xk(x,t)\varphi(t)dt, 0\leq x\leq 1,$

(1.4)

其中方程(1.3)为Volterra第一类积分方程, 方程(1.4)为Volterra第二类积分方程. 假设$g(x)\in L_2[0,1]$, $k(x,t)\in L_2[0,1]\times L_2[0,1]$是已知的实值连续函数, $\varphi(x)\in L_2[0,1]$ 是待求的未知解函数, $\lambda$是已知的常数. 这里我们假设Volterra第一、二类积分方程(1.3), (1.4)存在解析解.

一些求解积分方程(1.3)和(1.4)的数值解法被使用, 例如Galerkin法, 配置法, Taylor级数展开法, Bernstein多项式逼近法, Adomian分解法以及求积公式法^{[ 2-7]}. 一般情况下, 这些数值计算方法将积分方程转换成可以直接求解或通过迭代方法求解的线性代数系统, 进而求得积分方程的数值近似解. 然而, 多数情况下, 所得线性系统的系数矩阵并非稀疏矩阵, 尤其当离散化所得的线性系统是高维情形时, 利用这些数值计算方法求解积分方程就会很困难而且计算量也很大.

利用蒙特卡罗方法求解Fredholm第二类积分方程在文献[8]中被提出, 而且, 蒙特卡罗方法作为一种数值计算方法, 具有受几何条件限制小, 收敛速度与问题的维数无关, 误差容易确定, 程序结构简单、灵活、易于实现的优点. 本文结合确定性求解技术和蒙特卡罗方法的优点, 推广蒙特卡罗方法在求解第一类、二类Volterra积分方程(1.3), (1.4)中的应用.

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

设$\{u_i(x)\}$, $(i=0,1,2,\cdots)$是$L_2(0,1)$中的完备函数系, 它们可以是正交的, 也可以是非正交的, 但它们一定是线性无关的. 将积分方程的未知函数$\varphi(x)$按$\{u_i(x)\}$展开, 用有限项和式作为$\varphi(x)$的近似, 即

$\varphi (x) \approx \sum\limits_{i = 0}^n {{a_i}} {u_i}(x),$

(2.1)

只要确定了系数$a_i$, 积分方程的近似解$\hat{\varphi}(x)$就可以确定.

以下我们以Volterra第二类积分方程(1.4)为例阐明对Volterra积分方程的求解方法. 将式(2.1)代入积分方程(1.4), 可得

$\sum\limits_{i=0}^na_iu_i(x)-\lambda\sum\limits_{i=0}^na_i\int_0^xk(x,t)u_i(t)dt-g(x)=r(x),$

(2.2)

这里如果对任意的$x\in[0,1]$, 误差$r(x)\equiv0$, 则通过式(2.1)就得到积分方程(1.4)的精确解, 然而, 这通常难以做到. 因此, 为得到积分方程(1.4)的近似解, 可以放宽对误差$r(x)$的要求, 根据对误差$r(x)$的不同要求可以导致确定展开系数$a_i$$(i=0,1,2,\cdots)$的不同数值计算方法.

本文采用配置法将积分方程转化成线性代数系统. 此时, 要求式(2.2)中的误差$r(x)$满足

$(r(x),v_i(x))=0 (i=0,1,2,\cdots,n),$

(2.3)

其中$v_i(x)=\delta(x-x_i)$.

这里我们选取多项式$p_i(x)$作为$[0,1]$上的完备函数系$\{u_i(x)\}$$(i=0,1,2,\cdots)$. 设积分方程(1.3)和(1.4)的未知函数$\varphi(x)$是定义在区间$[0,1]$上的连续函数.

定理 2.1 设$\varphi(x)\in C[0,1]$, 则一定存在多项式$p_n(x)\in P_n$ 使得

$\max\limits_{0\leq x\leq1}|\varphi(x)-p_n(x)|<\epsilon.$

(2.4)

证 参见文献[9].

根据定理2.1以及矩量法, 选取函数$\varphi(x)$的逼近多项式$p_n(x)$为如下形式

$P_n(\varphi(x))=\sum\limits_{i=0}^na_ix^i$

(2.5)

且满足

$\max\limits_{0\leq x\leq1}|\varphi(x)-P_n(\varphi(x))|<\epsilon.$

(2.6)

将式(2.5)代回积分方程(1.4)中得

$\sum\limits_{i=0}^na_ix^i=g(x)+\lambda\sum\limits_{i=0}^na_i\int_0^xk(x,t)t^idt,$

(2.7)

此时, 式(2.7)中的积分项已可以计算, 未知量是$a_i (i=0,1,\cdots,n)$, 为求这$n+1$个未知的系数$a_i$, 可以建立如下线性系统

$\sum\limits_{i=0}^na_ix_j^i=g(x_j)+\lambda\sum\limits_{i=0}^na_i\int_0^{x_j}k(x_j,t)t^idt,$

(2.8)

针对核函数$k(x,t)$可能会存在奇异点的问题, 取节点$x_j=j/n+\epsilon, j=0,1,\cdots,n-1$, $x_n=1-\epsilon,$ 其中$\epsilon$是任意小的正数.

线性系统(2.8)可以简写为

$AX=b,$

(2.9)

其中

$A=[x_j^i-\lambda\displaystyle\int_0^{x_j}k(x_j,t)t^idt], i,j=0,1,\cdots,n,\\ X=[a_i]^t, i=0,1,\cdots,n,\\ b=[g(x_j)]^t, j=0,1,\cdots,n.$

(2.10)

求解线性系统(2.9), 可以得到积分方程(1.4)的近似解$P_n(\widehat{\varphi}(x))=\sum\limits_{i=0}^n\hat{a}_ix^i$, 其中$\hat{a}_i$是由线性系统(2.9)求解所得. 所求近似解的误差限由以下定理给出.

引理 2.1 假设

$\|A-I\|_1=c_1<1$

(2.11)

成立,其中$I$为$n+1$阶单位阵. 则有

$\|A^{-1}\|_1\leq\frac{1}{1-c_1}.$

(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}\|\leq\|I\|+\|B\| \|(I+B)^{-1}\|,\\ \|(I+B)^{-1}\|=\|A^{-1}\|\leq\frac{1}{1-\|B\|}=\frac{1}{1-c_1}.$

定理 2.2 对于在区间$[0,1]$上连续的函数$g(x)$以及连续核$k(x,t)$, 线性系统(2.9)中的矩阵$A$是可逆的. 记$\sup\limits_{x,t\in[0,1]}|\lambda k(x,t)|=c,$ 则有

$\sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\widehat{\varphi}(x_i))|\leq\epsilon\left(1+\frac{1+c}{1-c_1}\right)$

(2.13)

成立, 其中$x_i=i/n$, $i=0,1,\cdots,n$, $\varphi(x)$是方程(1.4)的精确解. $\epsilon$是一个任意小的正数.

证

$\sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\widehat{\varphi}(x_i))|=\sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\varphi(x_i))+P_n(\varphi(x_i))-P_n(\widehat{\varphi}(x_i))|\nonumber\\ \leq \sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\varphi(x_i))|+\sup\limits_{x_i\in[0,1]}|P_n(\varphi(x_i))-P_n(\widehat{\varphi}(x_i))|.$

(2.14)

下面确定(2.14)式中两项的界限.

由式(2.6)可得

$\sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\varphi(x_i))|\leq\epsilon.$

(2.15)

积分方程(1.4)可表示为

$g(x)=P_n(\varphi(x))-\lambda\int_0^xk(x,t)P_n(\varphi(t))dt,$

(2.16)

若在式(2.16)中将$P_n(\varphi(x))$替换成$P_n(\widehat{\varphi}(x))$, 则$g(x)$将被$\widehat{g}(x)$代替, 即

$\widehat{g}(x)=P_n(\widehat{\varphi}(x))-\lambda\int_0^xk(x,t)P_n(\widehat{\varphi}(t))dt,$

(2.17)

因此有

$\sup\limits_{x_i\in[0,1]}|P_n(\varphi(x_i))-P_n(\widehat{\varphi}(x_i))|=\sup\limits_{x_i\in[0,1]}\left|\sum_{j=0}^n(a_j-\hat{a}_j)x_i^j\right| \leq\|A^{-1}\|_1\max\limits_{x_i\in[0,1]}|g(x_i)-\widehat{g}(x_i)|.$

(2.18)

接下来求$\max|g(x_i)-\widehat{g}(x_i)|$, 为此令$\widehat{g}(x)=P_n(\varphi(x))-\lambda\displaystyle\int_0^xk(x,t)P_n(\varphi(t))dt$, 由表达式(1.4)可得

$\sup\limits_{x_i\in[0,1]}|g(x_i)-\widehat{g}(x_i)|=\sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\varphi(x_i))-\int_0^{x_i}\lambda k(x_i,t)[\varphi(t)-P_n(\varphi(t))]dt|\\ \leq\sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\varphi(x_i))|+\sup\limits_{x_i\in[0,1]}\left|\int_0^{x_i}\lambda k(x_i,t)[\varphi(t)-P_n(\varphi(t))]dt\right|\\ \leq\sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\varphi(x_i))|+\int_0^1\sup\limits_{x_i\in[0,1]}\left|\lambda k(x_i,t)\right|\sup\limits_{t\in[0,1]}\left|\varphi(t)-P_n(\varphi(t))\right|dt\\ \leq\epsilon(1+c).$

将这个误差限代回式(2.18)中有

$\sup\limits_{x_i\in[0,1]}|P_n(\varphi(x_i))-P_n(\widehat{\varphi}(x_i))| \leq\|A^{-1}\|_1\epsilon(1+c).$

(2.19)

根据引理2.1以及表达式(2.14), (2.15)和(2.19)有

$\sup\limits_{x_i\in[0,1]}|\varphi(x_i)-P_n(\widehat{\varphi}(x_i))|\leq\epsilon+\epsilon\frac{1+c}{1-c_1},$

此定理得证.

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

$LY=Pb$

(2.20)

以及

$UX=Y.$

(2.21)

通过这样处理可以提高线性系统的求解速度. 对于稀疏系统(2.21)利用雅阁比松弛迭代法进行求解. 记 $S=I-DU,$ 其中$D=\mbox{diag}(\frac{\gamma}{u_{11}},\frac{\gamma}{u_{22}},\cdots, \frac{\gamma}{u_{n+1,n+1}})$, $\gamma\in(0,1]$是松弛因子, $I$是单位矩阵, 则线性方程组(2.21)可表示为如下迭代形式

$X^{(k+1)}=SX^{(k)}+Y_1, k=0,1,$

(2.22)

其中$Y_1=DY$, $X^{(0)}=Y_1$.

根据计算方法理论可得如下定理:

定理 2.3 如果矩阵$S$的谱半径$\rho(S)$满足条件

$\rho(S)<1,$

(2.23)

则迭代形式(2.22)是收敛的, 即 $\lim\limits_{k\rightarrow\infty}X^{(k)}=X.$

证 参见文献[ 10].

利用待定系数逼近法求解积分方程(1.4), 根据Weierstrass第一逼近定理, 要提高近似解的精度, 需要增加逼近多项式$p_n(x)$中的项数$n$, 因而会面临求解阶数较高的线性代数系统的困难. 而且, 在利用确定性方法求解线性系统时, 求未知解向量中一个分量时常需求出解向量的全部分量. 这在一些实际问题中无疑会导致巨大的浪费. 蒙特卡罗数值求解方法不仅可以处理维数较高的问题(因为蒙特卡罗方法的收敛速度与问题的维数无关), 而且可以只对未知函数在某个点处的值或几个点处的线性组合值进行求解, 大大减少了数值求解的计算量.

下一节将利用蒙特卡罗方法求解线性方程组的迭代形式(2.22). 蒙特卡罗方法的基本思想是将所求问题的解视为一个随机变量的数学期望, 将随机变量多次实现的平均值作为所求问题的近似解.

对于收敛的迭代格式(2.22), 有

$X=Y_1+SY_1+S^2Y_1+\cdots+S^kY_1+\cdots.$

若令$s_{ij}\in S$, $X_i$是解向量$X$的第$i$个分量, 则有

$X_i=Y_{1i}+\sum\limits_{j_1=1}^{n+1}s_{ij_1}Y_{1j_1}+\sum\limits_{j_1,j_2=1}^{n+1}s_{ij_1}s_{j_1j_2}Y_{1j_2}+\cdots+\sum\limits_{j_1,j_2,\cdots, j_k=1}^{n+1}s_{ij_1}s_{j_1j_2}\cdots s_{j_{k-1}j_{k}}Y_{1j_k}+\cdots,$

(3.1)

其中$Y_{1i}$表示向量$Y_1$的第$i$个分量, $i=1,2,\cdots,n+1$.

在式(3.1)中, 令$s_{ij}=w_{ij}p_{ij}$, 则式(3.1)可表示为

$X_i=Y_{1i}+\sum\limits_{j_1=1}^{n+1}w_{ij_1}p_{ij_1}Y_{1j_1}+\cdots+\sum\limits_{j_1,j_2,\cdots,j_k=1}^{n+1}w_{ij_1}w_{j_1j_2}\cdots w_{j_{k-1}j_{k}}p_{ij_1}p_{j_1j_2}\cdots p_{j_{k-1}j_{k}}Y_{1j_k}+\cdots.$

(3.2)

考虑如下的马尔科夫链

$T=t_0\rightarrow t_1\rightarrow t_2\rightarrow\cdots\rightarrow t_k\rightarrow\cdots,$

状态空间为$\{1,2,\cdots,n+1\}.$ 记初始概率及由$i$状态转移到$j$状态的转移概率分别为

$P(t_0=i)=p_i (p_i\geq0), P(t_k=j|t_{k-1}=i)=p_{ij} (p_{ij}\geq0), i,j=1,2,\cdots,n+1$

且$p_i$和$p_{ij}$满足

$\sum\limits_{i=1}^{n+1}p_i=1, \sum\limits_{j=1}^{n+1}p_{ij}=1.$

马尔科夫链的权函数$W_k$由如下递推公式定义

$W_0=1, W_k=W_{k-1}\frac{s_{t_{k-1}t_k}}{p_{t_{k-1}t_k}}=\frac{s_{t_0t_1}s_{t_1t_2}\cdots s_{t_{k-1}t_k}}{p_{t_0t_1}p_{t_1t_2}\cdots p_{t_{k-1}t_k}}, k=1,2,\cdots,$

(3.3)

因此 马尔科夫链的权函数$W_k$为一随机变量. 进而可以构造权函数$W_k$的随机变量函数为

$\hat{X}_i=\sum\limits_{k=0}^\infty W_k{Y_1}_{t_k}, (i=1,2,\cdots,n+1).$

(3.4)

估计量$\hat{X}_i$有如下性质:

定理 3.1 如果初始概率满足$P(t_0=i)=p_i\equiv1$, 即$P(t_0=j)\equiv0$ $(j\neq i, i,j=1,2,\cdots,n+1)$, 马尔科夫链的权函数$W_k$如(3.3)式所示, 则随机变量$\hat{X}_i$ 的数学期望为

$E(\hat{X}_i)=X_i.$

(3.5)

证

$E(\hat{X}_i)=E(W_0{Y_1}_{t_0})+E(W_1{Y_1}_{t_1})+\cdots+E(W_k{Y_1}_{t_k})+\cdots\\ =Y_{1i}+\sum\limits_{j_1=1}^{n+1}w_{ij_1}Y_{1j_1}p_{ij_1}+\cdots+\sum\limits_{j_1,j_2,\cdots,j_k=1}^{n+1}w_{ij_1}w_{j_1j_2}\cdots w_{j_{k-1}j_{k}}Y_{1j_k}p_{ij_1}p_{j_1j_2}\cdots p_{j_{k-1}j_{k}}+\cdots.$

根据表达式(3.2) 有

$E(\hat{X}_i)=X_i.$

由定理3.1可见, 估计量$\hat{X}_i$是线性代数系统(2.21)解向量$X$的第$i$个分量$X_i$的无偏估计量, 因此蒙特卡罗方法可以只计算未知函数在某个点处的函数值.

若令$X_i^{(k)}$为迭代格式(2.22)的第$k$次迭代解$X^{(k)}$的第$i$个分量, 则它的蒙特卡罗逼近解为

$\hat{X}_i^{(k)}=\sum\limits_{m=0}^k W_k{Y_1}_{t_k}, (i=1,2,\cdots,n+1),$

(3.6)

满足$E(\hat{X}_i^{(k)})=X_i^{(k)}$.

为了估计$E(\hat{X}_i^{(k)})$, 对随机变量$\hat{X}_i^{(k)}$重复实现$N$次, 利用$N$次实现的算术平均估计$E(\hat{X}_i^{(k)})$. 根据初始概率$p_i$以及状态转移概率$p_{ij}$, $N$次实现链长为$k$的马尔科夫链

$t_0^x\rightarrow t_1^x\rightarrow t_2^x\rightarrow\cdots\rightarrow t_k^x, x=1, 2, \cdots, N,$

得到$N$次实现的样本均值为

$\overline{\hat{X}_i^{(k)}}=\frac{1}{N}\sum\limits_{x=1}^N\hat{X}_{ix}^{(k)}\approx E(\hat{X}_i^{(k)}).$

如果所得样本均值的标准差$\sqrt{\mbox{Var}\left(\overline{\hat{X}_i^{(k)}}\right)}$有限, 则根据中心极限定理

$P\left(\left|\overline{\hat{X}_i^{(k)}}-E(\hat{X}_i^{(k)})\right|\leq r\sqrt{\mbox{Var}\left(\overline{\hat{X}_i^{(k)}}\right)}\right)\approx2\Phi(r)-1,$

因此 估计量$\overline{\hat{X}_i^{(k)}}$在概率意义下的精度可以由方差$\mbox{Var}\left(\overline{\hat{X}_i^{(k)}}\right)$度量. 在实际计算时, 方差$\mbox{Var}\left(\overline{\hat{X}_i^{(k)}}\right)$常通过$\overline{\hat{X}_i^{(k)}}$的样本方差进行估计.

下面根据重要抽样原理给出随机路径$T$的产生方法, 即状态转移概率$P(t_k=j|t_{k-1}=i)=p_{ij}$的选择依据. 由以上讨论易知, 计算$X_i$需计算

$\sum\limits_{j=1}^{n+1}s_{ij}Y_{1j}=E\left(\frac{s_{t_{k-1}t_k}}{p_{t_{k-1}t_k}}Y_{1t_k}\right) =\frac{1}{N}\sum\limits_{x=1}^N\frac{s_{t_{k-1}^xt_k^x}}{p_{t_{k-1}^xt_k^x}}Y_{1t_k^x}\widehat{=}\hat{\theta},$

其中路径$T=t_0\rightarrow t_1\rightarrow t_2\rightarrow\cdots\rightarrow t_k\rightarrow\cdots$根据初始概率$p_i$与状态转移概率$p_{ij}$来选择. 初始概率在定理3.1中已经给出选择. 估计量$\hat{\theta}$的方差为

$\mbox{Var}(\hat{\theta})=\frac{1}{N}\left(\sum\limits_{j=1}^{n+1}\frac{(s_{ij}Y_{1j})^2}{p_{ij}}-\left(\sum\limits_{j=1}^{n+1}s_{ij}Y_{1j}\right)^2\right).$

(3.7)

由方差表达式(3.7), 估计量$\hat{\theta}$的精度与转移概率$p_{ij}$的选择有关. 为降低因随机路径引起的误差, 可选择$p_{ij}$使其接近模拟的目标 $s_{ij}$且满足归一性$\sum\limits_{j=1}^{n+1}p_{ij}=1$, 为此选择状态转移概率$p_{ij}$形式如下

$p_{ij}=\frac{|s_{ij}|}{\sum\limits_{j=1}^{n+1}|s_{ij}|}, i=1,2,\cdots,n+1,$

这就是重要抽样蒙特卡罗方法. 重要抽样体现了显著特征优先发生转移的原则.

这一节通过一些算例验证了文中所提方法对求解第一、二类 Volterra 线性积分方程的有 效性. 使用 MATLAB R2008a 软件编程实现如下算例. 计算结果与文献 [ 5] 中的结果进行了 比较. 结果表明, 文中所提方法是有效的可选择的方法.

本文使用文献[ 5]中的误差形式来度量数值近似解的精度. 误差定义如下

$\|e_n\|=\left(\int_a^be_n^2(x)dx\right)^{1/2}\approx\left(\frac{1}{n+1}\sum\limits_{i=0}^{n}e_n^2(x_i)\right)^{1/2},$

(4.1)

其中$e_n(x_i)=\varphi(x_i)-P_n(\widehat{\varphi}(x_i)),$ $x_i=i/n,$ $i=0,1,\cdots,n$.

符号说明 TMC表示基于泰勒展式下的蒙特卡罗方法, Bernstein表示伯恩斯坦多项式法.

当$n$取不同值时, 例1至例4的计算误差$\|e_n\|$列于表1.

为直观起见, 将本文提出的随机模拟方法与文献[ 5]中的伯恩斯坦多项式近似方法的对比结果绘于图 2至图 5.

例 1 考虑如下与热传导问题有关的具有弱奇异核的第二类Volterra积分方程 (见文献[ 5])

$ \varphi(x)=g(x)+\int_0^x\frac{t^{\mu-1}}{x^\mu}\varphi(t)dt, 0\leq x\leq1. $

令$\mu=1.5$,则$g(x)$是使方程的精确解为$\varphi(x)=x^5+x^{6.5}$的连续函数.

例 2 考虑如下第二类Volterra积分方程 (见文献[ 5])

$ \varphi(x)=\cos(x)-e^x\sin(x)+\int_{0}^xe^x\varphi(t)dt, 0\leq x\leq1. $

此方程的精确解为$\varphi(x)=\cos(x)$.

例 3 考虑如下Abel积分方程 (见文献[ 5])

$ \frac{2}{105}\sqrt{x}(105-56x^2+48x^3)=\int_{0}^x\frac{1}{\sqrt{x-t}}\varphi(t)dt, 0\leq x\leq1. $

此方程的精确解为$\phi(x)=x^3-x^2+1$.

例 4 考虑有关测试问题的第一类Volterra积分方程 (见文献[ 5])

$ \sin(x)=\int_0^xe^{x-t}\varphi(t)dt, 0\leq x\leq1. $

此方程的精确解为$\varphi(x)=\cos(x)-\sin(x)$.

由表 1,图 2至图 5可以得出如下结论:

$\bullet$ 由表 1可见, $n$值越大, 由蒙特卡罗随机模拟方法所求近似解愈接近积分方程精确解, 而$n$的增大并没有 给本文提出算法的编程和计算带来困难.

$\bullet$ 由图 2和图 3可见, 对于Volterra第二类线性积分方程, 利用蒙特卡罗随机模拟方法计算的多项式形式近似解 的精度接近伯恩斯坦多项式近似解的精度.

$\bullet$ 由图 4和图 5可见,对于Volterra第一类线性积分方程, 利用蒙特卡罗随机模拟方法计算的多项式形式近似解 的精度优于伯恩斯坦多项式近似解的精度.

本文利用蒙特卡罗技术与待定系数逼近法相结合的方法求解了 Volterra 第一、二类以及 具有弱奇异核的积分方程. 该方法解决了求积公式法难于处理的具有奇异核的第一类积分方 程求解问题. 算例实现表明, 蒙特卡罗方法可以求得积分方程具有一定精度的数值近似解. 最 主要的是蒙特卡罗方法可以根据实际问题需要计算未知量在某一点处的函数值, 计算耗时少.