径向基函数插值法解偏微分方程及计算渗流问题

径向基函数插值法解偏微分方程

及计算渗流问题

作者简介：魏义坤，男，1980年1月出生，师从于成都理工大学龚 灏教授，于2009年 6月毕业于成都理工大学计算数学专业，获得理学硕士学位。

摘 要

径向基函数插值法是在近十余年来发展起来的一种微分方程数值求解的无网格方法，该方法在对微分方程数值离散时不需要网格，因此不仅避免了网格生成的复杂过程，还可以显著减少传统网格方法（如有限元法、有限差分法）等中因网格畸变带来的不利影响。本文概括了径向基函数插值法对一维、二维函数插值拟合，并用径向基函数法解泊松和Hemholtz方程，将这种方法应用于计算简单渗流问题中，得到了比较满意的结果。全文共分四章。第一章是引言，主要介绍了无网格方法及径向基函数近年来的发展和研究现状。第二章是预备知识，首先介绍了径向基函数插值的基本理论和方法,然后介绍了一种特殊的径向基函数MQ函数及径向基函数插值法偏微分方程数值求解的理论和方法。第三章将这种方法应用于一维、二维的插值算例和计算经典偏微分方程算例；从数值实验得到的结果可以看出，径向基函数插值法在易用性和精度方面较传统的网格方法都有了很大的提高。第四章为渗流数学模型及定解条件。文中给出了描述渗流运动的基本方程以及基本方程的初始条件和边界条件，并针对于每一种具体的渗流模型给出了初始条件和边界条件。最后，将径向基函数插值法应用一个简单二维渗流问题。最后是结论和展望。 关键字：无网格法 径向基函数 插值法 渗流问题,

I

成都理工大学硕士学位论文

Interpolation Method of Radial Basis Function for Partial Differential Equations And calculation of seepage problem

Introduction of the author: Wei YiKun, male, was born in April, 1980, whose tutor was Professor Gong Hao. He graduated from Chengdu University of Technology in Mathematics calculation major and was granted the Master Degree in June, 2009.

Abstract

Radial basis function interpolation in the past dozen years to develop a numerical solution of differential equations of the meshless method, the method of differential equations without numerical grid dispersion, thus avoiding not only the mesh generation complex process, but also can significantly reduce the traditional grid method (such as the finite element method, finite difference method), such as mesh distortion due to the adverse impact.

This article summarizes the radial basis function interpolation method for one-dimensional, two-dimensional interpolation function fitting, and radial basis function solution before mooring and equations, this method has been applied to simple flow problems, have been relatively satisfied with the results.

The full text is divided into five chapters. The first chapter is the introduction, introduces the meshless radial basis function method and the development and research in recent years, the status quo. Chapter II is to prepare the knowledge, first introduced the radial basis function interpolation of the basic theory and methods, and then introduced a special kind of radial basis function MQ function and radial basis function interpolation of the numerical solution of partial differential equations theory and methods. Chapter III of this method is applied to one-dimensional, two-dimensional interpolation of the classic examples and calculation examples of partial differential equations; from the results of numerical experiments we can see that the radial basis function interpolation in the ease of use and accuracy than the traditional grid method has been greatly improved. Chapter IV of the mathematical model for the flow and boundary conditions. This paper describes the basic equation of flow movement and the basic equation of initial conditions and boundary conditions, and for each specific flow model of a given initial conditions and

II

A

boundary conditions. Finally, the radial basis function interpolation applied a simple two-dimensional seepage problem.

Finally, conclusions and outlook.

Keywords: Meshless method radial basis function Interpolation method Seepage problems

III

独创性声明

本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的研究成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得 成都理工大学 或其他教育机构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示谢意。

学位论文作者签名：

年 月 日

学位论文版权使用授权书

本学位论文作者完全了解 成都理工大学 有关保留、使用学位论文的规定，有权保留并向国家有关部门或机构送交论文的复印件和磁盘，允许论文被查阅和借阅。本人授权 成都理工大学 可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。

（保密的学位论文在解密后适用本授权书）

学位论文作者签名： 学位论文作者导师签名：

年 月 日

第1章 引言

第1章 引 言 1.1历史概述及研究背景

在工程和科技等领域中的许多问题,人们可以用数学模型进行描述。这些模型通常是常微分方程或偏微分方程及定解条件。但能用解析方法求出精确解的只是少数方程比较简单,且几何形状相当规则的情况。对于大多数问题,由于方程本身的复杂性,或由于求解域的几何形状比较复杂,只能采用数值方法求解。20世纪年代以来,随着电子计算机的出现,数值计算和分析方法已经成为求解科学技术问题功能强大的有利工具。

随着计算机技术的高速发展，有限元方法成为本世纪数值分析[5][6]的重要成果，并且在科学技术的各个领域得到了充分的应用 。然而，人们一直没有停止过寻找新的数值计算方法的研究。这主要是由于现已比较成熟的计算方法并非完美无缺，由于自身的缺陷，它们在解决某些问题时会遇到难于克服的困难。对于传统的计算方法[1][2]（有限元，有限差分等），在处理某些问题时会显得有些笨拙，如：

1、工业材料冲压成型过程中，模拟大变形情况下的材料流动； 2、裂纹沿任意复杂路径动态扩展的模拟；

3、高速撞击过程中产生的结构的几何畸变(如穿透、侵彻等)； 4、奇异性等特殊问题等。

基于网格的方法，在计算过程中如果网格畸变，将导致计算失效。也就是说，网格的存在妨碍了处理与原始网格线不一致的不连续性和大变形。为了处理大变形或随时间变化的不连续性，通常需要进行网格重构。然而，这样不仅计算费用昂贵，而且会使计算精度严重受损。因此这类问题需要更加有效的数值分析方法，最近几年正迅速发展的无网格数值计算方法（Meshless Method）以其特有的优点适合这类问题。无网格方法抛开网格，直接基于结点形成近似，具有重要的研究价值和应用前景。

1.2无网格方法的提出和发展

[8]

无网格方法产生于大约三十年前，但当时发展很慢直到最近十年来,

无网格方法才受到了众多专家学者的重视，并得到了迅速的发展。

对无网格方法的研究可以追溯到20世纪70年代对非规则网格有限差分法的研究，由于当时有限元法的巨大成功，这类方法没有受到重视。1977

1

成都理工大学硕士学位论文

年Lucyt 和Gingold等光滑质点流体动力学 (smoothed particle hydrodynamics, SPH)方法，并且成功地应用于天体物理领域中。等提出了归一化光滑函数算法，使其通过分片试验，可以正确模拟常应变状态，提高了SPH法的精度；Swegle,Dykah 和Chen等提出了SPH方法不稳定的起因及稳定化方案；Vignjevic等提出了克服零能模态的具体方案；Monaghan对SPH方法进行了总结；SPH方法已被应用于冲击波模拟、流体动力学、水下爆炸仿真模拟、高速碰撞等材料动态响应的数值模拟的领域。然而由于精度及稳定性问题,该方法未能得到广泛应用。

1992年，Nayrdesd等人提出了一种全新的方法“Diffuse

approximation and Diffuse Element Method”(散射元法,简称DEM),并用此方法分析了Possion方程和弹性问题.在该方法中,首次将移动最小二乘近似(Moving Least Square Method,简称MLS)引入Galerkin法中,得出具有C1连续性的近似解。

1994年,美国Northwestern University的学者Belytschko等人对方法进行了修改和发展,基于MLS近似和Galerkin原理,,利用背景风格进行积分计算,形成”The Element-free Galerkin Method” (无单元Galerkin法,简EFG)。这类方法比SPH方法计算费用高,但具有较好的协调性及稳定性。在众多学者努力下,围绕此类方法已经进行了许多的研究,并对实际遇到的的许多问题进行了分析。比如,Belytschko和Tabbara等人将方法用于动态裂纹扩展的数值模拟,克服了有限元方法在模拟裂纹扩展时需不断进行网格重新划分的缺点，在计算中可以连续地进行模拟；Krysl等使用方法对板壳进行了EFG分析:Belytschko,Organ和Hegen等EFG对方法和有限元方法的祸合使用进行了研究，并将其用于边界条件的处理：此

外,Belytschko和Krysl又将EFG用于三维撞击和流体晃动分析等人进行了方法用于相变问题的研究。

Northwestern University的另一位学者Liu W. K,根据函数积分变换的思想，提出ReproducingKernel Particle Method（再生核粒子法,简记为RKPM）并结合小波的概念,构造了Multi Scale Reproducing Kernel Particle Method，（多尺度再生核粒子法，简记为MRKPM）。小波分析最初被广泛应用于数字信号处理领域，其具有尺度伸缩平移、多分辨率等特点。MRKPM利用小波函数的多尺度分析思想,构造了一系列可现时伸缩和平移的窗函数,实现了RKPM的自适应分析,可用于对局部进行细致的数值分析。使用RKPM方法,Liu等人对大量问题进行了数值分析，如结构动力学问题、流体力学问题、橡胶材料的几何大变形问题等。

1995年, Oden和Duarte等提出了”Hp-Clouds”无网格方法，该方法

2

第1章 引言

利用移动最小二乘原理建立单位分解函数，用来进行场量模拟，然后通过变分建立离散模型。Oden严格证明了通过在Hp-Clouds空间引入单位分解函数，可以建立任意阶连续光滑的基函数。Oden, Duatte,Zienkiewicz等又将有限元形函数作为单位分解函数,提出了” Nen Clouds-Based FEM”这种方法需要借助于有限元网格,虽然破坏了”无网格”的部分特,但给解决问题带来了方便。

波兰的学者Liszka等提出了Hp-Meshless clouds Method,,,该方法与oden等人提出的Hp-Clouds方法的不同之处在于它采用配点格式,不需要Galerkin格式中用于进行积分计算的背景网格，是一种完全的无网格方法。

1996年Babuska及Melenk提出了The Partition of Unity Method(单位分解法,简称PUM)。Babusa利用现代数学单位分解的概念,对单位分解法进行了严格的数学论证,并提出了使用特征解析解进行增强的基函数。

西班牙的学者Onate和Idelsohn等在1996年提出了The Finite Point Method(有限点法,简称FPM)。该方法采用移动最小二乘原理来构造形函数，采用配点格式进行离散，是一种不需要背景网格的完全无网格方法，主要应用于流体力学领域。

后来,著名学者Attuti等也对无网格方法进行了大量研究,提出了Local boundary integral equation method (局部边界积分法，简称LBIE)和Meshless Local Petrow-Galerkin Method（简称MLPG）这两种方法都是基于移动最小二乘原理来建立对场函数的近似，并目在积分时不需要背景网格，是完全的无网格计算方法。MLPG与LBIE的不同之处在于，LBIE是建立于局部边界积分方程，具有奇异积分的计算；而在MLPG计算中不涉及奇异积分。

随着小波应用范围的扩展,许多学者将小波用于场量的近似，并产生了Wavelet-GalerkinMethod(小波伽辽金法，简称WGM)，该方法是利用小波级数对场量进行近似，通过Galerkin变分对偏微分方程进行数值离散。这种方法在处理局部化现象时，如塑性变形局部化等问题，具有较大的优势。

总之，无网格方法在近十几年来得到了迅猛的发展。随着越来越多的学者对其研究的投入，无网格方法势必会得到更大的发展。

1.3无网格方法的优势和缺点

我们就无网格方法与传统网格方法中最有代表性的有限元法做以比较,说明无网格方法现在的优势和缺点。

3

成都理工大学硕士学位论文

建立近似函数时不借助网格，基于函数逼近近似而非插值是无网格法与有限元方法的主要区别。无网格法具有以下优点：

(1) 无网格法的近似函数没有网格依赖性，减少了因网格畸变而引起的困难，尤其对于处理高速碰撞、!动态断裂!塑性流动、流固祸合等涉及大变形和需要动态高速节点位置（网格）的各类应用问题时，具有较大优势。

（2）无网格法基函数可以包含能够反映待求问题特性的函数系列，适用于分析各类具有高梯度、奇异性等特殊性质的应用问题。

（3）采用紧支函数的无网格法和传统网格方法一样具有带状稀疏系数矩阵的特点。适用于求解大型科学与工程问题。

（4）无网格法的自适应很强。在自适应分析中不需要重新划分网格，民极易实现p自适应分析，若引进小波函数还具有多尺度分析功能。

(5) 无网格法的前处理只要节点位置信息，不用网格信息,容易分析复杂三维结构。

（6）无网格计算的结果是光滑连续的，不必再进行应力光顺化等后处理。 无网格方法相对于有限元在处理大变形问题上面尽管有很大的优势，但无网格法才刚刚起步，在严格的数学论证。计算效率、边界条件处理和大量应用实例方面都不能和成熟的有限元方法相媲美，更未形成有效的通用软件。现在的无网格方法主要有以下缺点：

（1） 用MLS和RKPM建立无网格近似函数时，都涉及到矩阵求逆,计算量

较大。

（2）与有限元方法不同，无网格的近似函数大都不是多项式，因而基于伽辽金法的无网格方法（比如EFG,RKPM等），等需要在每个背景网格或节点子域中使用高阶高斯积分以保证其精度，因此伽辽金型无网格法的计算量一般远大于有限元法。

（3）基于配点法的无网格法（比如SPH,FPM）等只是在节点上形成方程，计算效率高,但是精度和稳定性都较差。

(4) 无网格近似函数大都不具有占函数的性质，施加本质边界条件的时候不如有限元方便。

1.4径向基函数的提出和发展

径向基函数的研究是从径向基函数插值开始的。Buhmann、Dyn、Lingt、Powell 、Schaback等分别写了综述性的文章及Golberg收集的参考目录

首先考虑径向基函数插值的一些不同领域的来源。最早的可能是Krige，他在1951年矿藏的沉积看成是一个各向同性的稳定的随机函数的实现，从而导

4

第1章 引言

出了广泛应用于矿藏分析的Kriging方法。在这方面的理论深入主要是Matheron实现的。1971年Hardy利用径向基函数Multi-Quadric来处理飞机外形设计曲面拟合问题，取得了非常好的效果1975年，Duchon从样条弯曲能最小的理论出发导出了多元问题的薄板样条。

径向基函数（Radial Basis Funtion,简称RBF）简称为具有形式简单、各向同性等优点，数学界对其进行了大量的研究，已被成功应用于多变量插值中。Frank比较了29种离散数据插值方法,发现Hardy提出的Mulquadric(MQ)函数和Duchcn提出的TPS插值函数精度最高。这两类函数是常用的径向基函数。Kansa将径向基函数引入配点法中，用以求解双曲、椭圆和抛物型问题。Kansa和Sharan等人都发现函数具有指数收敛性，而目计算效率极高。

在Kansa的配点法中，系数矩阵是非对称的。Fasshauer采用Hermite型方法，得到的系数矩阵在某些特殊的偏微分中是对称或反对称的。Wendland将径向基函数引入Galerkin法中，建立了相应的无网格格式。Coleman用径向基函数求解椭圆型边值问题。吴宗敏和Franke等人证明了用径向基函数进行离散数据插和求解偏微分方程的优越性，并给出了误差估计。陈文采用径向基函数提了边界点法，只需要对求解域的边界用节点进行离散，而不需要离散求解域。Hon等人采用径向基函数求解了两相流问题。将径向基函数引入配点法以求解偏微分方程具有很多优点，例如：该方法是真正的无网格方法，不需要任何网格：该方法与空间维数无关，其收敛率为O(hd+1)，其中为配点密度，d为空间维数。然而径向基函数一般是定义在全求解域的，得到的系数矩阵的条件数很大。目前己经提出许多种方法来解决这一问题，如优化调整MQ函数的形状参数以及采用紧支径向基函数。吴宗敏给出了正定紧支径向基函数的构造法则，并构造了一系列正定函数：Buhmann给出了基于TPS函数的正定紧支径向基函数：Wendland给出的正定紧支径向基函数与吴宗敏的相似，但他们具有最小的自由度。张雄等将紧支距离函数应用于配点法中，建立了相应的无网格法，用于求解固体力学问题。虽然紧支距离函数具有很多优势，但他们不能满足完备性条件，甚至不能正确描述常应的配点型无网格格式，宋康祖等提出了满足完备性的紧支距离函数，并建立了相应的配点型无网格格式，大幅度提高了计算精度,但也增加了计算量。径向基函数方法除在大地测量学，地球物理学，测绘学等诸方面有应用外,在神经网络一同、材料力学哪刘等方面也有应用。随着研究的进一步深入，径向基函数的应用也会越来越广泛。

Frank曾做了大量的各种散乱数据插值方法的实际比较，认为径向基函数插值的结果最能使人满意。

5

成都理工大学硕士学位论文

第2章 径向基函数插值法及解偏微分方程相关知识

2.1径向基函数插值

所谓径向基插值[29]就是对给定的多元散乱数据xj,fj

⎛基函数φ:R→R，利用平移构造基函数系⎧⎨φ⎜⎩⎝

+

{}mj=1

选取径向∈Rn×R，

⎫

x−xj⎞⎟⎬并寻找插值函数S(x)

⎠⎭j=1

m

形如 满足 记

S(x)=∑ajφx−xj

() (2-1)

S(xj)=fj,

j=1,2,L,m

fT=(f1,f2,L,fm)

φT(x)=(φ(x−x1),φ(x−x2),Lφ(x−xm))，

aT=(a1,a2,L,am)

A=⎡φ(x1),φ(x2)Lφ(xm)⎤

⎣⎦

如果A为非奇异，那么（2-1）可写成

S(x)=φT(x)A−1f

（2-2）

有时候为了某种目的还添加上一个多项式。即寻找

⎞ S(x)=ajφ⎛x−xbαxα （2-3） ⎜j⎟+

⎝⎠α∑∑的函数满足插值条件，这里

xαα是标准的多元记号

α=(α1,α2,L,αm)， αj是非负整数，

α1α2αm

=x1x2Lxm

α=∑αj ，x=maxαj

f

(α)

关于还按字典排列规定一个序，α∂f

=α1α2

αn

∂x1∂x2L∂xn

α<β，如果有α1<β1或者有一个

j,αj<βj当iE=那么当增加约束条件

Tα(xαj)α≤m, X=(1,Lx,L)

∑αjxαj=0 (α6

≤d

)

第2章 径向基函数插值法及解偏微分相关知识

时就可以得到

s(x)=(φ是插值问题的解。

注意到（2-2）及（2-4）的解都是关于f线形齐次的，所以S(x)可以写成 其中AλT

⎛A,x)⎜

⎜ET⎝

TE⎞⎛f⎞

（2-4） ⎟⎜⎟⎜⎟⎟00⎠⎝⎠

=φ或者

s(x)=∑λj(x)fj （2-5）

E⎞⎛λ⎞⎛φ⎞

⎟⎜⎟=⎜⎟ （2-6） ⎜⎟⎜⎟0⎟⎠⎝μ⎠⎝X⎠

m

⎛A

⎜⎜ET⎝

每个插值点的导数

'

S(x)=∑ajφ'(xj),i=1,2,Lm。

j=1

''

mj=1

S(x)=∑ajφ''(xj),i=1,2,Lm。

22

常用的径向基函数有：

（a） Kriging方法的 Gauss分布函数：φ(r)=e−cr； （b） Kriging方法的Markoff分布函数：φ(x)=e

布函数；

（c） Hardy的Multi-Quadrik函数：φ(r)=(c2+r2)β； （d） Hardy的逆Multi-Quadrik函数：φ(r)=(c2+r2)−β； （e） Duchon的薄板样条：φ(r)=r2rlogr,φ(r)=r2k+1； （f） 紧支柱正定径向基函数[4]

下面一些是Wu’s函数 具体表达式如下：

−cr

，及其他概率分

φ0,0(r)=(1+r)∈C0∩PD1

22

φ1,0(r)=(1−r)3+(1+3r+r)∈C∩PD1

0φ1,1(r)=(1−r)2+(2+r)∈C∩PD3

.

2344φ2,0(r)=(1−r)5+(1+5r+9r+5r+r)∈C∩PD1

232φ2,1(r)=(1−r)4+(4+16r+12r+3r)∈C∩PD3

.

.

20φ2,2(r)=(1−r)3+(8+9r+3r)∈C∩PD5

234560φ3.0(r)=(1−r)7+(5+35r+101r+147r+101r+35r+5r)∈C∩PD5

23454

φ3.1(r)=(1−r)6(6+36r+82r+72r+30r+5r)∈C∩PD3+

7

成都理工大学硕士学位论文

2342

φ3.2(r)=(1−r)5+(8+40r+48r+25r+5r)∈C∩PD5

230

φ3.3(r)=(1−r)4+(16+29r+20r+5r)∈C∩PD7

2.2 径向基函数插值的存在性问题

给定数据xj,fj，关于径向基函数最简单的插值公式可以表示为

{}f*(x)=∑λjφ(x−xj)

且满足插值条件

∑λjφ(x

k

−xj)=fk （2-7）

显然上述插值方程对任何数据xj,fj∈Rd⊗R，当xj两两不同时都有解的充分必要条件是：对任何的两两不同的xj，对称矩阵（(φ(

{}xk−xj))）都是非

奇异的。如果假设径向基函数是连续的，且φ(∞)=0，那么，可以利用归纳方法得到：适当取符号假设φ(0)>0。当两个数据点时，设想两个点的距离非常远，那么φ(∞)=0得到这个两阶的正定矩阵的行列式大于零。然后让这俩个点移近，由函数的连续性得到这是一个两阶的正定矩阵。一般地，如果已经有n个数据点，方程（2-10）的系数矩阵，其直至n的各阶主子式都大于0.当再增加一个数据点时，让这个新增加的点从很远的地方向所在位置移动。那么n+1个数据点的n+1阶系数矩阵行列式在第n+1个两两不同数据点的n+1阶系数矩阵的各阶主子式都大于零，从而由归纳法得到，任何形式为(φ(阵。有定理[4]。

定理2-1 函数φ:R+

xk−xj))的矩阵都必须是正定矩

→R是连续的，limr→∞φ(r)=0。那么对d元的径

向基函数插值是存在唯一解的条件是：矩阵(φ(

xk−xj))都是正定矩阵。

对d维空间的变元x,定义Φ(x)=φ(x)，并且假设这个函数是可以由Fourier积分Φ(x)=吗、，二次型

∫R

d

eΦ(t)dt表示。那么当xj两两不同，λj不全为零时

ixt

^

ixjt2

^

∑λjλkφ(xk−xj)=∫∑λe

j

Φdt>0。

而由经典的三角函数逼近知识知道，

∑λe

j

ixjt2

可以逼近几乎任何的正函数，所

以径向基函数插值对任何两两不同的数据点都存在唯一解的充分必要条件是：

8

第2章 径向基函数插值法及解偏微分相关知识

Φ(t)几乎处处非负，且至少在一个正测度上大于零。

满足这种性质的函数φ(.)称为d变元的正定函数。利用上述关于正定函数的讨论得到了一个重要的定理

定理2-2 （Bochner） 函数Φ(x)是正定函数的充分必要条件是Fourier变换Φ(t)几乎处处非负，且至少在一个正测度上大于零。

这个定理充分说明，要判断函数φ是否可以用来构造d元的正定函数Φ，需要求一个d元的积分。

一般地，函数Φd(x)=φ(

^

^

x)，x∈Rd的一个正定函数，并不能导得

^

Φd+1(x)=φ(x)，x∈R

d+1

也是一个正定函数，这个定理是由Fourier变换描

2ϕ(r)=ϕ(r)。eϕ(t)dt，其中∫

述的。如果我们采用Laplace变换表示ϕ(r)=

−rt

那么对任意空间函数φ(x)都是正定的充分必要条件是ϕ几乎处处非负，且至少在一个正测度上大于零，还可以相似不加证明得到另一个重要定理。

定理 2-3（Micchelli） 如果定义ϕ(r2)=φ(r)，那么φ(x)对任何维空间的x都是正定函数的充分必要条件是ϕ(r)是全单调函数，即ϕ∈C∞且

^

(−1)kϕ(k)(r)>0。

可见Gauss函数，逆Multi-Quadric在任意维空间都是正定函数。Multi-Quadric薄板样条不是正定函数。但是，他们的广义Fourier变换还是大于零，只是在零点有一个极点。如果极点是γ阶的，那么当满足条件二次型仍然可以用积分表示

^

ixjt2

∑λjxαj=0，∀α≤γ时，

∑λjλφ(x

k

k

−xj)=∫∑λje

Φ(t)dt>0.

这类函数称为γ阶的条件正定函数。定义矩阵

⎛(φ(⎜

A=⎜

⎜⎝其中0≤α⎞xj−xk))(xα)j⎟

xk

β0

⎟ ⎟⎠

,β≤γ。还要证明，这个矩阵当(xαj)满秩时，矩阵A也是非奇

同时得到∑λjxα=0会导致μ=0。j=0,∀(α)≤k。

异的。用反正法，如若不然，那么存在(λT,μT)，满足(λT,μT)A=0。由(xαj)满秩得到λT≠0，否则λ计算有(λT,μT)A(λT,μT)T

=λT(φ(xj−xk))λ>0，与反正假设矛盾。

对条件正定径向基函数，如果把插值方程作适当的修改，添加上一个多项

9

成都理工大学硕士学位论文

式项

满足

f*(x)=∑λjφ(x−xj)+∑bαxα+∑cjβj(x) fk=∑ajφ(xk−xj)+∑bαxα

及附加条件

的解唯一存在，容易验证这也是一个插值公式。上式中的α的选取范围与条件正定的阶数有关。确定地说如果要包含条件正定的阶数要求，需要选取所有的

∑

ajxαj=0

{α/α<γ}，并且还可以再添加上一些，那么上述方程都唯一可解。

对于条件正定函数也有两个与正定函数相对应的定理

定理2-4

^

Φ(x)作为d变元的函数是k阶正定函数的充分必要条件是其

Fourier变换Φ(t)几乎处处非负，且至少在一个正测度上大于零，且在零点是一个不超过k+d-1阶的极点。

定理 2-5如果定义ϕ(r2)=φ(r)，那么φ(x)对任何维空间的x都是k阶条件正定函数的充分必要条件是：ϕ(r)是阶的全单调函数，即ϕ∈C∞且

(−1)lϕ(l)(r)>0， ∀l>k。

一般地，如果再选择一些函数βj(x),j=1,2,L,b,那么当矩阵(xαj/βk(xj))

满秩时

满足

f*(x)=∑λjφ(x−xj)+∑bαxα+∑cjβj(x)，

fk=∑ajφ(xk−xj)+∑bαxα+∑cjβj(xk),∀k，

j

及附加条件

∑aβj

k

(xj)=0，

∀k，∑ajxαj=0，∀α。

的解也是唯一存在的。而且这个公式也是一个插值公式。并且进一步地由方程的系数矩阵非奇异性得到解的唯一性定理。容易验证，如果数据是采自一个形如

f(x)=∑bαxα+∑cjβj(x)

αfk(x)=∑bαxk+∑cjβj(xk)，那么 f*(x)=f(x)=∑bαxα+∑cjβj(x)

的函数，即

这是说如果希望插值方法有某类函数空间再生性质，那么只要将这类函数加入到

βj这类里就可以了。这是一个十分重要的性质，特别对实际问题有着非常广泛

的应用前景。

10

第2章 径向基函数插值法及解偏微分相关知识

2.3 径向基函数插值的收性问题

径向基函数插值，不仅有用事实上的一元函数描述多元函数的计算就表示简单的优点，而且还有最小方差无偏估计这个统计意义上的最优解及薄板样条弯曲能两次泛函最这个物理意义上的最优解，这两个非常明确合理的物理意义，如果适当选取自相关函数或者合适的两次泛函量度，那么在某种意义下，径向基函数插值有某种最优性。有了这样的插值方法，还希望从纯数学的角度估计插值的误差，或者讨论当数据密集化时的收剑速度。

只要数据点两两不同，散乱数据径向基函数插值问题是唯一可解的[4]。它的解可以写成

满足

f*(x)=∑ajφ(x−xj)+∑bjpj(x)，

f*(xk)=fk=∑ajφ(xk=xj)+∑bjpj(xk)

j

j

及附加条件

∑ajpk(xj)=0

Mj=1

这个函数同时有一个Lagrange表示 如果定义

f(x)=∑f(xj)uj(x).

*

R(x)=(φ(x−x1),L,φ(x−xM))T， U(x)=(μ1(x),L,μM(x))T， V(x)=(v1(x),LvQ(x))T， S(x)=(p1(x),LpQ(x))T，

A=(φ(xj−xk)) ， P=(pk(xj))jk

那么径向基插值的Lagrange基满足

⎛AP⎞⎛U(x)⎞⎛R(x)⎞

⎟⎜ ⎜T⎟=⎜⎟。 ⎜⎟⎜⎜P⎟0⎠⎝V(x)⎠⎝S(x)⎟⎠⎝

把V(x)看成是Lagrange乘子，那么径向基函数插值的Lagrange表示U(x)同时是在给定x时两次条件最优问题

minUTAU−2UTR(x)+φ(0)/U∈RM,PTU=S(x)

{}的解。这个两次条件最优问题是唯一可解的，并且它的解就是径向基函数插值的解。这里给出了径向基函数插值的另一个物理意义（Kriging意义）。如果φ是2m

11

成都理工大学硕士学位论文

次可为的，定义Φ(x)=φ(x)，那么U(x)满足下列线形方程组

⎛AP⎞⎛U(μ)(x)⎞⎛R(μ)(x)⎞

⎟⎜(μ) ⎜T⎟=⎜(μ)⎟， ⎜⎟⎜⎜P⎟0⎠⎝V(x)⎠⎝S(x)⎟⎠⎝

从而它是在给定x时两次条件最优问题

minU(μ)AU−2U(μ)R(μ)(x)+Φ(2μ)(0)

T

T

{}， （2-8）

有限制条件

U(μ)∈RM， PTU(μ)=S(μ)(x)

下的解。为了深入讨论散乱数据的径向基函数插值问题，给下列定义

定义2-1 φ(x)是q阶条件正定的函数，并假设φ∈C2m[0,∞)，那么对于0≤μ≤m及数据点x1,L,xM，函数

k(x):=minU

(μ)q

{(μ)T

AU

(μ)

−2U

(μ)T

R(μ)(x)+Φ(2μ)(0)

}，

其中U满足下述条件限制

U(μ)∈Kμ,q(x)， 称为在点x的μ阶的Kriging函数。这里集合

M⎧⎫⎪⎪M(μ)

Kμ,q(x):=⎨μ∈R/∑μjp(xj)=p,∀p∈Pq⎬ （2-9）

j=1⎪⎪⎩⎭

称之为允许向量集合。从而条件最优问题在允许向量集合中使得Kriging函数取

到最小值。事实上，Kriging函数是径向基函数插值误差的在一个Hilbert再生核空间下的范数。

径向基函数插值的误差可以由Kriging函数控制 Kriging函数有一二积分表示，利用Fourier变换得到，利用合适允许向量U∈Kμ,q(x)放入广义Fourier变换。

如果φ(.)是一个绝对可积函数，并且有非负的Fourier变换Φ，满足 φ(x

可以得到恒等式

^

)=∫R

(μ)T

n

e

ix,t

Φ(t)dt，

(μ)T

^

x∈Rn，

U

AU−2U

M

ixj,tjj=1

R(μ)(x)+Φ(0)

2

=∫R

T

n

∑μe

−e

ix,t

Φdt

^

， （2-10）

这是两次最优问题或者Kriging范数的Fourier表示。进一步有

U(μ)AU(μ)−2U(μ)R(μ)(x)+Φ(2μ)(0)

T

=∫R

n

∑u

j=0

M

2

(μ)

j

e

ixj,t

−(it)e

ix,t

Φdt

^

（2-11）

12

第2章 径向基函数插值法及解偏微分相关知识

这样就用一个积分表示了Kriging函数

在一般情形，函数φ并不一定是绝对可积的（只是局部可积）

M

g(t):=∑uje

j=1

i

−(it)μei<.x,t> （2-12）

（当uj,x及μ给定时）并不是分布函数的一个试验函数。（2-10）或（2-11）的积分可能并不存在，所以不能直接套用广义Fourier变换来解决这个问题。其中g(t)有个性质：

引理 2-1 对每个允许向量u∈Kμ,q(x)，当x及数据点xj都给定，并且

0≤μ≤q，函数g(t)在零点的行为是t

证明 令 ex

q+1

，而且在无穷远处的行为是t

μ。

=p(x)+xq+1ρ(x),， ρ(x)≤ex

为ex的q阶的Taylor展开式，而μ是允许向量μ∈Kμ,q(x)，那么

g(t)e

Mj=1

−ixj−x,t

=∑uje

j=1

M

ixj−x,t

−(it)μ

μM

=∑ujp(ixj−x,t)−(it) +∑uj(ixj−x,t)q+1ρ(ixj−x,t)

j=1

=∑uj(ixj−x,t)q+1ρ(ixj−x,t)

j=1

M

因为ρ(ixj−x,t)是一个关于y的多项式，它的μ阶导数在点x是(it)μ，所以当x，所有的xj,0≤

μ≤q及u都给定时

⎧Ο(t⎪

g(t)≤⎨

⎪⎩Ο(t

q+1

),t=0

μ),t=∞

2.4 误差估计问题

如果数据xj,fj采自于一个可以由Fourier积分

{}表示的函数，可以利用式（2-11）中的u(μ)(x)得到误差表示

^⎛m(μ)ixj,tμix,t⎞ f(x)−f(x)=∫n⎜∑uj(x)e −(it)e⎟f(t)da （2-13）R

⎝j=1⎠

这个式子与Kriging函数有一些相像的地方。如果

*(μ)

(μ)

2

f(x)=∫ne

R

ix,t^

f(t)dt

13

成都理工大学硕士学位论文

⎛^⎞2

cf:=∫nf(t)⎜φ(t)⎟dt<∞ （2-14）

R

⎝⎠

利用Cauchy-Schwarz不等式，得到

^

2−1

f*(μ)(x)−fμ(x)

≤∫

2

⎛≤⎜⎜∫Rn⎝

∑u(jμ)(x)e

j=1

m

ixj,t

−(it)μe

^

2

ix,t

⎞f(t)dt⎟

⎟⎠

^

2

R

n

∑uμ(x)e

()jj=1

m

2

ixj,t

−(it)μe

ix,t

φ(t)dt⋅∫

^

Rn

⎛^⎞

f(t)⎜φ(t)⎟dt

⎝⎠

−1

=kqμ(x)⋅c2f，

从而插值误差是点点受控于Kriging函数。

定义 2-1 称函数f: Rn→R受控于一个在Rn\\{0}满足2μ<β<2q+2及

0≤φ(t)≤ct

n

^

−n−β放入径向基函数φ，当且仅当，如果f有一个广义的Fourier变

^

换在R\\{0}表示f，满足（2-14）且 f(x)≤cx

^

−n−γ2

,x∈Rn\\{0}， 及 2γ2=β

^

−n−β定理 2-6 如果f受控于一个满足2μ<β<2q+2及0≤φ(t)≤ct径向基函数φ，那么插值误差

f

*(μ)

的

−f

(μ)

2

μ≤κq(x)⋅c2f

受控于Kriging函数.

现在知道径向基函数插值的误差可由Kriging函数控制，而径向基函数插值[25]又使得Kriging函数在允许向量空间取最小，以下的问题是如何构造向量，使得Kriging函数尽可能小。有

引理 2-2 对任何k>0,r>0及任何的点x∈Ω，存在正的常数h0,c1，使得在任何的具有局部密度

hr(x):=max

y−xmin

1≤i≤M

y−xi

≤h0

的数据分布xi∈Ω,1≤i≤M中，可以找到一个子集zβ∈X中β∈N0n是满足0≤

:={xi/1≤i≤M},其

β≤k的多元记号，成立

zβ−x≤c1⋅hr(x)⎛

Z:=⎜zβ−x

⎝

⎝

()μ⎞0

⎟满足detZ≠0其中0:=1， ⎠μ,βZ−1=B*⎛⎜δμβ⋅hr

−μ(x)⎞⎟

⎠μ,β=BC−1(N

−μ−μδμβ)⎛hδμβ⎞⎜⎟，

⎝⎠

14

第2章 径向基函数插值法及解偏微分相关知识

B

*B≤2， ∞

∞

≤c。 （2-16）

这里及以下的多元记号μ和β总满足0≤

β,μ≤k。如果多元记号发生在矩阵

中，总是假设以字典顺序排列。而在矩阵括号外的下标的第一个元代表矩阵的行。譬如B=(bμβ)μ,β表示矩阵的行以μ排，列以β排。

证明 矩阵C=(βμ)μ,β是通常的Vondermonde矩阵的多元形式。因为次数不超过K的多项式在数据点

{β/0≤β≤k}上的插值的唯一存在的，所以矩

∞,

阵C是一个非奇异的矩阵。利用矩阵C定义常数

N≥2k(2k+1)k+1Q

C−1

r

。 Nk

现在假设数据分布稠密，使得h=hr(x)≤h0在点x+Nβh的邻近存在点

c1≥(1+Nk)n，h0<

zβ∈X

满足

zβ−x−Nβh≤h （2-17） 尽快（2-15）

⎛

利用（2-17）和（2-15）来估计矩阵⎜

⎝

(β⋅N⋅h)μ⎞

⎛和⎜⎟

⎠μ,β⎝

μ−1

(z

β−x

)μ⎞

差⎟⎠μ,β的上界，利用多项式展开式，得到

⎛⎞μμ t−(t+y)≤μy⎜max(t+τy)∞⎟

∞⎜⎟0≤τ≤1⎝⎠

从而对矩阵Z-V的元素得到估计式

，

(zβ−x)−(βNh)≤μh(c1+kN)h

μμ()μ−1μ⎞

，

⎛

这里变形的Vondermonde矩阵V=⎜

⎝V=

(βNh)有形式 ⎟⎠μ,β((βNh)μ)∞

μ,β=δμβ⋅(Nh)(μ)μ,β⋅C，

可以选择N足够大，使得 V−1(Z−V)=C−1⋅⎛⎜δμ,β(Nh)⎝

−μ⎞

⎟(Z−V) ⎠μ,β∞

15

成都理工大学硕士学位论文

≤C−1max(Nh)∞

−μμhμ(c1+kN)k−1

μ−1

Q

≤C−1≤C−1

−1

∞

∞

1

Q

N11⋅k(2k+1)Q≤N2⎞

+k⎟

⎟⎠

⎛c1⋅k⎜⎜N⎝

有界。这样Neumann序列矩阵

B=ZV=VZ

(−1

)=(I−V(V−Z))∑

−1

󰀀

−1−1

⎡V−1=

⎢⎣j=0

∞

(V−Z)⎤

⎥⎦

j

存在，且它的范数小于2.从而矩阵Z−1=BV−1满足引理。 现在用引理的结果为（2-8）构造允许向量Uμ，如果B*=

[4]

(b)βμβ,μ，对

任何的r,满足0≤

r≤k,得到一个关于y的多项式

−μr

bβμh(zβ−x+y)∑β−μ

=∑bβμh

β1

r!∑zβ−x

r!r!r1+r2=r12

()r1

yr2

r11−μbβμh(zβ−x) =r!∑y∑r1!r2!βr1+r2=r

r!1r(μ)r−μr2

=

μ!(r−μ)!

y=

μ!(y

)

上式中令y=x,得到 其中

󰀀(μ)⎪μ!bβμhui:=⎨

(x)

r(μ)=μ!∑bβμh

β−μr

zβ=∑uixir， （2-18）

i=1

m󰀀

⎧⎪⎩

,xi=zβ

0,xi≠zβ⎛󰀀μ⎞

=⎜ui⎟ ，可以成为式（2-18）的允许向量，所以 ⎝⎠i

−μ由方程（2-18），U g(t)=

󰀀

󰀀

(μ)∑

M󰀀(μ)ix,t

ujejj=1

−(it)μe

ix,t

μix,t

=∑μ!bβμhe

β−μizβ,t

−(it)e

16

第2章 径向基函数插值法及解偏微分相关知识

⎛μ⎞−μizβ−x,t

=e⎜μ!bβμhe−it⎟

⎜⎟β⎝⎠

∞jjμ⎞ix,t⎛−μi

bβμhzβ−x,t−it⎟=e⎜μ!

⎜⎟j=0j!β⎝⎠

ix,t

∑

()∑∑

()=e

ix,t

⎧∞

ij⎪

⎨μ!

j=k+1j!β⎪⎩

∞

∑∑bβμh

−μ⎫

μ⎪tρj!−μijρzβ−x,t+μ!∑∑bβμh(zβ−x,t)−(it)⎬∑j!!ρj=0β⎪ρ=j

j

k

⎭

=e

ix,t

j−μij

μ!∑∑bβμhzβ−x,t j=k+1j!β进一步得到，对任何的

󰀀

t≤1/h，有

c1jhjt∑j!j=k+1

∞

j

g(t)e−ix,t≤μ!h−μ

bβμ∑β

≤const⋅h

−μ⋅hk+1t

k+1

󰀀

k+1

⋅expc1ht

−μ()

。

≤const⋅hk+1−μt

由（2-16）当

t→∞时，有g(t)≤ch+t

μ2.5 径向基函数插值法解偏微分方程

用径向基函数方法求解偏微分方程

[21]

Lu(x)=s(x), in⊂Rn， （2-19）

Bu(x)/∂Ω=f(x). （2-20）

首先介绍一下插值点

Θh={(xi)/i=1,N−M⊂Ω,(xi)/i=N−M+1,N⊂∂Ω}

uh(xi)=∑αjgj(x)+αN+1, x∈Ω=∪∂Ω, （2-21）

j=1N

N

把代入（2-19），（2-20）用插值点 得到下列插值系统

Luh(xi)=∑αjLgj(x)+αN+1L1=s(xi),i=1,L,N−M （2-22）

j=1N

Buh(xi)=∑αjBgj(x)+αN+1B1=

j=1

f(xi),i=N−M+1L,N.（2-23）

上面的L1和B1分别是L和B在点1处的

17

成都理工大学硕士学位论文

GL

GB

LgN(x1)⎤⎡Lg1(x1)L

⎥ （2-24） =⎢MMM⎢⎥

⎢⎣Lg1(xN−M)LLgN(xN−M)⎥⎦⎡Bg1(xN+M−1)LBgN(xN+M−1)⎤

⎥ （2-25） MMM=⎢⎢⎥

⎢LBgN(xN)⎥⎣Bg1(xN)⎦⎡GL=⎢⎢GB

T⎢p⎣

Lp⎤

(N+1)×(N+1)

Bp⎥∈R ⎥0⎥⎦

W

b=(s(x1),L,s(xN−M),f(xN−M+1),Lf(xN),0)T∈RN+1

可以把式子（2-22）和（2-23）写出

Wα=b 得到 α=W−1b （2-26）

把（2-26）带回（2-21）就可以得到

uh(xi)=g(xi)W−1b （2-27）

2.6 Multiquadric 的发展

Multi quadric是由R.LHardy于1968年提出来的一种径向基函数[14]。人们通常将Multi quadric简记为MQ据Hardy的综述文章所述：Multi

quadric-niharmonic(MQ-B)方法发现于1968年,Hardy总结了自MQ发现以20来年间所取得的成就。他列举了在大地测量学（geodesy）,地球物理学(geophysics)测绘学（surveying and mapping）,摄影测量学(photogrammetry),遥感与信号处理

(remote sensing and signal processing),地理学(geography),地质与采矿(geology and mining),数字地形模型(digital terrain models)及水文学(hydrology)诸方面的应用。他还将这20年细化为三个时期。即从1968年至1972年为研究基本的MQ方法的时期一这些方法包括配置模型(collocation mode)，最小二乘模型（least square ，密切模型(osculating mode ),区域分解及其它(domain decomposition and mode）

other matter)。1972年至1979年为过渡时期，这一时期主要表现为MQ在上述诸方面的广泛应用。第三个时期就是1980年至1988的MQ-B及其数学论证时期。

Hardy最后在论文的第五部分一Concluding remarks中提出的9条意见中的第4条中指出MQ可应用于并已应用于微分方程和积分方程。事实上,1968年至1988年这20年,MQ虽然在微分方程和积分方程中有应用，但人们主要将MQ用于插值。

R.Franke在其评论文章中指出就(accuracy)精度、（stability）稳定性、（efficiency）有效性、（memory requirement）内存需要和(ease to implementation)易于实现而言，

MQ在所有种散乱数据插值格式中首屈一指。W.R.Madych和S.A.Nelson证明了用MQ作为插值，总是得到（minimal semi-norm error）最小半范数误差。

18

第2章 径向基函数插值法及解偏微分相关知识

M.D.Buhmann,对Multoqiadric(MQ)：

φ(r)=r2+c2和逆MQ：

φ(r)=

1研究了奇数维Euclidean空间上的插值的收敛阶，得出：当n≥322r+c为奇数时，它的收敛阶分别为n-1阶和n+1阶。

M.D.Buhmann和C.A.Micchili研究了多变量整数网格上MQ的局部化性质 的改进问题。

由于Multiqudric(MQ)函数不是正定函数，因而，插值问题势必面临可解性的问题，为此，有的学者转而寻求解决拟插值的问题,该问题不存在可解性的问题，而且，只在拟插值构造得好，也可以获得理想的精度。

19

成都理工大学硕士学位论文

第3章 径向基函数插值法及解微分方程算例

通过上一章的介绍，我们已经对径向基函数插值法及其解偏微分方程有了一个很深的认识和理解，至于所提方法是否合理可言，这就需要我们通过经典例子对此方法验证。本章基于这个目的，利用MATLAB语言偏写的径向基函数插值程序对一维、二维函数拟合，并利用径向基函数无网格方法解经典的泊松方程及其解Hemholtz方程和对误差和节点数对误差影响分析。主要说明此方法的精确性和可行性，利用其函数图像和误差图像比较分析，更加形象的分析此方法的优点。

3.1径向基函数插值算例

3.1.1 一维插值算例

本算例用的径向基函数插值是全局径向基函数插值，通过在区域[0，10]内分别取20个、40个和100个均匀的节点对余弦函数y=cosx、余弦函数一阶导数利用全局径向基函数插值法对上述三y=−sinx和余弦函数二阶导数y=−cosx，

个函数进行拟合，分别在同一个直角坐标系中对插值解图与真实图比较，通过图像分别显示在区域[0,10]不同均匀节点误差，通过误差图分析取不同节点时余弦函数、余弦函数一阶导数和余弦函数二阶导数利用全局径向基函数插值分析各自的拟合效果。

y=cosx,dy=−sinx,d2y=−cosx 计算区域为[0，10]，余弦函数y=cosx及一阶导数y=−sinx和二阶导数y=−cosx的真实图及拟合图

算例1

像，计算过程中径向基插值函数用的是MQ函数（参数c=2）

20

第3章 径向基函数插值法及解微分方程算例

10.80.60.40.20-0.2-0.4-0.6-0.8-1051010.80.60.40.20-0.2-0.4-0.6-0.8-1051010.80.60.40.20-0.2-0.4-0.6-0.8-10510图3-1真实图像与插值拟合图像对比

注：蓝线是实图像 ，绿线用的拟合图。

图3-1从左到右分别给出余弦函数

y=cosx和余弦函数的一阶导函数

y=−sinx、余弦函数的二阶导数y=−cosx的真实图和插值图，利用径向基函

数插值在区域[0,10]取20个均匀离散点，从插值结果看误差主要集中在峰谷和两端点，在峰谷和两端点拟合效果不是很好，但在其它区域拟合效果很高。

0.80.710.60.50.50.40.3-0.50.2-10.10-1.5-4-5-300-1-211.532051005100510图3-2 相对误差分布 （纵坐标为%） 21

成都理工大学硕士学位论文

图3-2为区域[0,10]取20个均匀离散点余弦函数和余弦函数的一阶导函数、余弦函数的二阶导数的相对误差图像，从图3-2结果相对误差都非常低，但是余弦函数的一阶导数、二阶导数误差相对高些，尤其二阶导数的相对误差最高。

10.80.60.40.20-0.2-0.4-0.6-0.8-1051010.80.60.40.20-0.2-0.4-0.6-0.8-1051010.80.60.40.20-0.2-0.4-0.6-0.8-10510图3-3真实图像与插值拟合图像对比

注：蓝线是实图像 ，绿线用的拟合图

图3-3从左到右分别给出余弦函数

y=cosx和余弦函数的一阶导函数

y=−sinx、余弦函数的二阶导数y=−cosx的真实图和插值图，利用径向基函

数插值时在区域[0,10]取40个均匀离散点，很明显从插值结果看比(图3-1)区域[0,10]取20个均匀离散点拟合精度提高很多，误差主要集中在峰谷和两端点，其他区域拟合效果很高。

22

第3章 径向基函数插值法及解微分方程算例

0.80.70.60.50.40.30.20.1021.5164200.50-0.5-1-2-4-6-8-10-1.5-2-120510-140−3

0510510图3-4 相对误差分布，（纵坐标为10）

图3-4为在区域[0,10]取40个离散点的余弦函数和一阶导函数、二阶导数的相对误差图像，从图3-4结果相对误差都非常低，但是余弦函数的一阶导数、二阶导数误差相对高些，尤其二阶导数的相对误差最高。

10.80.60.40.20-0.2-0.4-0.6-0.8-1051010.80.60.40.2-0.50-1-0.2-0.4-0.6-0.8-10510-1.5-2-2.5-31.510.500510 图3-5 真实图像与插值拟合图像对比 (间距󰀀h=0.1)

注：蓝线是实图像 ，绿线用的拟合图

23

成都理工大学硕士学位论文

图3-5从左到右分别给出余弦函数

y=cosx和余弦函数的一阶导函数

y=−sinx、余弦函数的二阶导数y=−cosx的真实图和插值图，利用径向基函

数插值在区域[0,10]取100个均匀离散点，很明显从插值结果看比(图3-1)区域[0,10]取20和40个均匀离散点拟合精度提高很多，其插值图像和真实图像几乎重合，拟合效果十分高。

0.80.740.620.5-200.40.3-20.2-40.10-6-100-120-60-800-4020066040051005100−5

510图3-6 相对误差分布，(纵坐标为10)

图3-6为在区域[0,10]取100个离散点的余弦函数和一阶导函数、二阶导数的相对误差图像，从图3-6结果相对误差都非常低，但是余弦函数的一阶导数、二阶导数误差相对高些，尤其是二阶导数的相对误差最高。

间距 h=0.5 h=0.25 h=0.1

误差和y 5.83% 1.01% 0.28%

表3-1 误差分析

误差和dy 10.16% 2.78% 0.325%

误差和ddy 29.06% 3.472% 0.4426%

从上面的表和图可以得到径向基函数插值取节点越多函数图像导数图像拟合效果越好，当在区域[0,10]取100个离散点时，真实图像和拟合图像几乎都重合了，但是其导数的误差总和增加了，不过只是仅仅在几个某节点上误差大，整体的拟合效果是越来越精确，用全局径向基函数插值具有很高的精度，但是它不能很好地拟合函数的二阶导数，还有其结果严重依赖于插值函数(MQ)参数的选取。

24

第3章 径向基函数插值法及解微分方程算例

3.1.2 二维插值算例

f(x,y)=x⋅e

(−x2−y2)

,x,y∈[−2,2]节点80*80；径向基函数采用的是MQ

函数，（参数C=2）；插值图像和误差如图3-7和3-8；

图 3-7 径向基函数插值图像

图3-8 径向基函数插值误差值图像

25

成都理工大学硕士学位论文

图3-9 MATALAB 函数插值图像

图3-10 MATALAB 函数插值误差

根据上面插值函数结果及误差结果很明显看出，径向基函数插值精度十分高，比MATLAB函数插值[11]精度高很多，但是拟合精度也与径向基函数中的参数选取也有关，参数选取一般凭经验，有待进不步研究。下面用径向基函数法应用于解经典偏微分方程[13]。

26

第3章 径向基函数插值法及解微分方程算例

3.2径向基函数插值法解微分方程算例

3.2.1 解泊松方程

泊松方程的第一类边值问题

2⎧

⎛⎞⎪∇2u(x,y)+⎜∂u⎟=2y+x4,(x,y)∈[0,1]×[0,1]

⎪⎜∂y⎟

⎝⎠⎪⎪

u(0,y)=0,u(1,y)=y (3-1) ⎨

⎪2

(,0)0,(,1)ux=ux=x⎪

⎪⎪⎩

上述泊松方程的解析解为

u(x,y)=x2y (3-2) 在计算中对求解区域[0,1]×[0,1]用 N =20×20 个规则节点均匀划分，支撑域半径为scale=7，MQ 中参数c = 2.0。

图 3.11泊松方程在 y=0.5 上数值结果和精确值比较

图 3-11 比较了不同径向基函数（MQ、CSRBF）得到的 y = 0.5上各点解的近似解和解析解。计算结果表明，几种基函数近似解与解析解吻合得比较好。

27

成都理工大学硕士学位论文

图 3-12 精确解与 MQ 近似解比较

从图3-11和3-12计算中发现，只要每个径向基函数的参数取的合理计算精度是很高的，另外，即使将各紧支径向基函数的支撑域取得足够大，甚至各支撑域中包含所有的节点，紧支径向基函数插值的精度仍然低于全局径向基函数的插值精度。显然，强迫紧支径向基函数满足一致性条件可以显著提高近似的精度。

对于同一基函数，节点的数目、支撑域半径的大小以及参数c对计算结果都是有影响的。为了考察它们对数值结果的影响，针对问题二用 CSRBF 基函数计算了关于节点数目、支撑域半径大小(scale)的相对误差，本文主要用的径向基函数MQ 基函数，文献[18]给出了 MQ 基函数中参数c对结果的影响。

对数值结果影响很大，本文所用两种基函数中 CSRBF 较1) 不同的径向基函数，差，而 MQ的结果相对较好；

2) 对于同一基函数，节点取得越密，计算精度越高，但是计算量也相对的增加很多；

3) 支撑域半径的大小[8]，对计算结果有明显影响，半径越大支撑域内包含的节点数目就越多，它的精度也越高；但当支撑域半径达到一定程度时，增加节点数目对精度的影响变小。

3.2.2 解Hemholts方程

考虑二维区域Hemholtz方程第一边界值问题[13]

2⎧⎪Δu(x,y)+ku(x,y)=f(x,y),(x,y)∈Ω (3-3) ⎨

u(x,y)=u,(x,y)∈Γ⎪0⎩

本例中取MQ函数(c=2),k=2区域分为[0,8]×[0,8]，通过取不同的均匀节点数对径向基函数求解Hemholtz方程[24]和解析解误差分析比较，节点数对误差的形象，从中得到一定结论，所得结果仍然满足工程计算精度要求，精度非常高。

28

第3章 径向基函数插值法及解微分方程算例

Absolute error for a Helmholtz case using 10*10 nodes-6x 1042absolute error0-2-48642y002x648图3-11 Hemholtz方程径向基函数插值法误差（10*10节点） 图3-11是在区域取10×10均匀节点 ，相对误差单位是10−6，从上图相对误差结果是非常小。

Absolute error for a Helmholtz case using 50*50 nodesx 101-70.5absolute error0-0.5-18642y002x648图3-12 Hemholtz方程径向基函数插值法误差（50*50节点）

图3-12是在区域上取50×50均匀节点，相对误差单位是10−7，从上图相对误差结果是比区域取10×10均匀节点小很多。

29

成都理工大学硕士学位论文

Absolute error for a Helmholtz case using 100*100 nodes-6x 1010.5absolute error0-0.5-18642y002x648图3-13 Hemholtz方程径向基函数插值法误差（100*100节点）

图3-13是在区域上取100×100均匀节点，相对误差单位是10−6，从上图结果分析相对误差结果比50×50节点相对误差稍微增高了，因此取离散点是应采用适当的节点数。

如果采用适当的径向基函数的参数及分布点，无网格法算出来的结果精度相对很高，编程而且简单易懂。图3-11、3-12和3-13图像径向基函数插值方法的形函数构造基于区域内任意分布的节点，取节点很少时误差相对比较大，随着节点数越来越多其误差就越小，但是节点很多时误差又会增大，所以利用无网格方法计算时要取适当的节点数，选择适当的径向基函数的参数大小误差才会是相对很小。

径向基函数插值法解偏微分方程克服了有限元法有时需要重构网格的缺点它易于处理任意不规则的边界，克服了差分法对不规则区域形状而引起的降低精度的缺点。由数值算例可以看到该法给出的结果要比传统的有限差分法好,而计算量要小得多。特别对于非直线边界它具有很好的优势。上面通过径向基函数插值法解泊松方程和Hemholtz方程，上述数值试验得到此方法解偏微分方程计算精度十分高，下一章将此方法解偏微分方程应用到一个实际二维恒定渗流问题。

30

第4章 渗流问题数学模型及径向基函数插值法界二维渗流问题

第4章 渗流问题的数学模型及径向基函数插值法解

二维渗流

4.1描述渗流问题的基本方程

4.1.1 达西定律

渗流运动的基本规律，早在1855年由达西通过试验研究而总结出来，一般称为达西定律。其表达式为：

U=kJ (4-1)

或：v=kJ （4-2）

式中u为渗流流速；k为比例系数，也就是渗透系数，它取决于土的透水性能；J是渗透水流的水力坡降；v为断面平均渗流流速。达西定律表明对均匀沙土中的恒定均匀渗流，其水力坡度(单位流程的水头损失)与渗流流速的一次方成正比。因此也称为渗流线性定律。

达西定律经过后来的大量实践和研究，认为可将其推广运用到其它孔隙介质的恒定渗流和非恒定渗流中去。此时的达西定律表示式只能用式（4-1），而且对非恒定渗流的水力坡度和渗流流速应改写成如下形式：

u=−k

∂H

（4-3） ∂s∂HJ=− （4-4）

∂s这是液体做层流运动所遵 达西定律表明渗流的水头损失与流速呈线性关系，

循的规律。研究证明，达西定律只能适应于层流渗流(线性渗流)。

巴甫洛夫斯基给出达西定律应用范围的临界流速uk的经验公式为： uk

(0.75n+0.23)Nv

（4-5）

d 式中n为孔隙率；d为粒径，可用d10代替；v为水的运动粘性系数；N为常

=

数，约为7-9

附带指出，在水利工程中，大多数渗流运动是服从达西定律的。但是在堆石坝、堆石排水等大孔隙介质中，渗流为紊流，这时应采用非线性的渗流定律，巴甫洛夫斯基建议的公式为：

v=ktJ(cm/s) （4-6）

式中kt是渗流为紊流时的渗透系数。根据伊兹巴什的试验，对大粒径砾石

(d>5 cm) ，kt式可由下式确定：

31

成都理工大学硕士学位论文

kt

d⎠⎝

其中n为孔隙率；a为一常数，对完整块石，a为14；对破碎块石(n=0.4)，a为5。

a⎞=n⎛20−⎜⎟d(cm/s) （4-7）

4.1.2 渗流的连续性方程

渗流简化模型假设渗流区内的全部空间都被连续的水流所充满，因此，渗流的连续性和其它水流的连续性是一样的，两者有完全相同的连续性方程。对不可压缩的流体，有

∂ux∂uy∂uz++=0 （4-8） ∂x∂y∂z在各向同性土体的情况下，根据达西定律可得：

∂H

(4-9) ∂x∂H

（4-10） uy=−k∂y∂H

（4-11） uz=−k∂zp

式中，k单元土体的渗透系数;H渗透水头，即H=+y； p单元土体中心处

ρg的水压力；ρ流体的密度；g重力加速度；y单元土体中心处的位置水头。

ux=−k

将公式(4-9).、(4-10)、（4-11)代入(4-8)，则渗流的连续方程将变为下列形式：

∂2H∂2H∂2H

+2+2=0 （4-12） 2∂x∂y∂z

公式(4-12)表示三维无旋流的流态。

在二维平面渗流问题[15]的情况下，公式（4-12)变为：

∂2H∂2H

+2=0 （4-13） 2∂x∂y

对于均匀的各向异性土体的情况下，根据达西定律，渗流的流速可表示为：

∂H

（4-14） ∂x∂H

（4-15） uy=−kx

∂y∂H

（4-16） uz=−kx∂z式中kx、ky、kz分别为土体沿x、y、z方向的渗透系数。

将公式(4-14) 󰀀 (4-16)代入公式(4-8)，则得均匀各向异性土体情况下，

ux=−kx

三维渗流的连续方程为：

32

第4章 渗流问题数学模型及径向基函数插值法界二维渗流问题

∂2H∂2H∂2H

kx2+ky2+kz2=0 (4-17)

∂x∂y∂z

在二维平面渗流的情况下，均匀各向异性土体的渗流连续方程为：

∂2H∂2H

kx2+ky2=0 (4-18)

∂x∂y

4.1.3 渗流的运动方程

渗流时液体在介质的孔隙中流动受到孔隙周界的阻力，这种阻力是很大的，必须在渗流简化模型[10]中如实反映。由于土粒直径与渗流区尺度相比是很微小的，可以认为在分析水流作用力时，所取的脱离体中仍包含着足够多的土粒。这样，土粒对流动的阻力就可以认为是均匀分布在脱离体内。因此，可以把渗流阻力看作是作用在简化模型的液体内的一个体积力。

设

fx、fy和fz为渗流简化模型中单位质量液体的阻力f在各坐标轴方向的

分力，单位质量液体所受的其它质量力在各坐标轴方向的分力为X、Y、Z参照 描述水流运动的欧拉运动方程，可写出渗流的运动方程为：

∂u∂u∂u∂u⎫1∂ρ+fx=x+uxx+uyx+uzx⎪ρ∂x∂z⎪∂t∂x∂y∂uy∂uy∂uy∂uy⎪1∂ρ⎪

Y− +fy=+uy+uy+uz⎬ （4-19）

ρ∂y∂t∂x∂y∂z⎪1∂ρ∂uz∂uz∂uz∂uz⎪

⎪Z−+f=+ux+uy+uz

∂x∂y∂z⎪ρ∂zz∂t⎭

X−

阻力f的表达要利用达西定律，设f的方向与流向相反，液体沿流向移动距离

ds时，单位质量液体阻力所做的功为−fds，对单位重量液体而言，这个功为 −fds/g。若以H代表某一点的总水头(在渗流中是测压管水头)，则−dH是

单位重量液体所提供的能量，在数量上等于阻力所做的功。因此有：

fds

=dH （4-20） g对非恒定渗流，可写为

f=g∂H (4-21)

∂s由达西定律，有

∂Hu

=− （4-22） ∂sk

于是得f的表达式为

33

成都理工大学硕士学位论文

g

f=−u （4-23）

k写成分量形式，分别为：

fx=−gux,fy=−guy,fz=−guz （4-24）

hhh考虑到渗流流速极小，加速度中的位变加速度为二阶小量，故可以忽略，于是式(4-19)可简化为:

∂ux⎫1∂ρg

−u=X−

ρ∂xhx∂t⎪⎪

∂uy⎪1∂ρg⎪

−uy= Y−⎬ （4-25） ρ∂yh∂t⎪

∂uz⎪1∂ρg

⎪Z−−uz=

ρ∂zh∂t⎪

⎭

当质量力只有重力，用z轴垂直向上的坐标系，则Z=−g，X=Y=0由

ρ H=z+ γ有

∂p∂H=γ∂x∂x

，

∂p∂H

=γ∂y∂y

，

∂p∂H⎞

=γ⎛−1⎜⎟

∂z∂z⎝⎠

则式（4-25） 变为：

⎫1∂ux∂Hux++=0⎪g∂t∂xk⎪

⎪∂uy∂Huy1 ++=0⎪⎬ （4-26）

g∂t∂yk⎪

⎪1∂uz∂Huz++=0⎪g∂t∂zk⎪⎭

这就是非恒定渗流运动方程。

1888年，茹可夫斯基对于恒定渗流情况，将上式简化为： ∂H⎫

⎪∂x⎪∂H

uy=−k⎪⎬ （4-27）

∂y⎪∂H⎪

uz=−k⎪

∂z⎭

此即在恒定渗流中广泛应用的运动方程式。

ux=−k

4.1.4恒定渗流的流速势和拉普拉斯方程

对于恒定渗流而言，当土质为均质、各向同性时，渗透系数k是个常数。故在式(4-27)中，可以将k写在偏导数里面。

34

第4章 渗流问题数学模型及径向基函数插值法界二维渗流问题

ϕ=−kH (4-28)

∂ϕ⎫

⎪∂x⎪∂ϕ⎪

则得 uy=⎬ （4-29）

∂y⎪∂ϕ⎪uz=⎪

∂z⎭

ux=

上式表明必就是渗流的流速势。由此可得如下结论：在重力作用下，可以将均质各向同性土的恒定渗流看作是一种特殊的有势流动。

如把式(4-29)代入连续方程(4-8)，可得渗流的拉普拉斯方程: .

∂2ϕ∂2ϕ∂2ϕ+2+2=0 （4-30） 2∂x∂y∂z

∂2H∂2H∂2H

+2+2=0 （4-31） 2∂x∂y∂z

将式3-28)代入上式得

故渗流水头H也满足拉普拉斯方程。

由此可见，恒定渗流问题可归结为解拉普拉斯方程求渗流流速势(或渗流水头H)的问题。当H求得，亦即压强P己知，当少己知，由式(4-29)就可求渗 流流速。

由于渗流流速势ϕ和水头H有式(4-28)的关系，故等势面必然也是等水头

面。因为等势面上任一点的流速矢量与等势面垂直，故流速矢量也必与等水头面垂直。

4.2渗流数学模型基本方程的定解条件

4.2.1 初始条件和边界条件

以二维恒定渗流的连续方程:

∂ux∂uy+=0 (4-32) ∂x∂yux=−k

和运动方程

∂H⎫∂x⎪⎪

⎬ （4-33）

∂H⎪

uy=−k

∂y⎪⎭

为例，件所组成的方程组共有三个微分方程，其中包含

u,u

x

y

和H三个未知函

数，在一定的初始条件和边界条件下，即可求得渗流的流速场和水头场(或压强场)。

35

成都理工大学硕士学位论文

由于直接求解式(4-32)、(4-33)较为麻烦，大家通常选用的办法为求解由式(4-32), (4-33)变换后的二维恒定渗流流速势的拉普拉斯方程。由方程推导可知，上述方程均适用于渗流为有压与无压两种情况。对于渗流有压流与无压流，只是边界条件明显不同，二维恒定渗流基本方程定解的边界条件有如下几个方面：

1.不透水边界

对于无压渗流，不透水边界为下面不透水岩层，对于有压渗流，不透水层除下面不透水岩层外，还有上面不透水建筑物轮廓。不透水边界是一条流线，垂直于边界的流速分量必等于零，即

un或

=ucos(n,u)=0

∂H∂ϕ==0 ∂n∂n式中n为边界的法线方向。

2.透水边界

透水边界上各点水 透水边界指水流渗入的边界和低于下游水位的渗出边界，

头H相同，是一条等水头线(或等势线)，渗流流速必与透水边界垂直。

3．浸润面边界

对于无压渗流，浸润面即重力水区与毛细水区的分界面，该面上的压强等于大气压强，即相对压强p=0。该面上各点的水头H=y，故不是常数，它是随各点的垂直坐标位置y而定。在恒定渗流中，浸润面位置不变，浸润面是由流线组成的面。渗流流速在浸润面的法线方向上的分量为零，故也像不透水边界一样，有

∂H∂ϕ==0 ∂n∂n应该注意：浸润面本身的位置事先是不能给定的，它需待渗流问题解决的同时才能确定。对于有压渗流[10]，无浸润面边界。

4.渗出段边界

对于无压渗流，浸润面出口位置常高于下游水位，由此形成渗出段边。渗出段边界各点的压强等于大气压强，各点的水头H=y，即随各点的垂直坐标位置而变。对于有压渗流，由于下游水位一般高于不透水建筑物底面高程，无此渗出段边界。

4.2.2 几种渗流模型及边界条件

1、在平面渗流(二维渗流)的情况下，对于符合达西定律的各向异性连续体中的渗流，可以用下列微分方程来表示:

36

第4章 渗流问题数学模型及径向基函数插值法界二维渗流问题

∂⎛∂H⎞∂⎛∂H⎞

+⎜+Q=0 （4-34） kxky⎟⎜⎟⎜⎟⎜⎟∂x⎝∂x⎠∂y⎝∂y⎠

kx、ky，是x方向和y方向的渗流系数；Q是入渗或蒸发的流量；H是水头，

可以用下式表示:

H=

P

+y （4-35） ρg公式(4-35)的边界条件通常有三类，即:

(1)边界上水头H己知，即H=Ho；

(2)边界上流入或流出的水量己知，此时边界上应满足下列流量平衡方程式： kx

∂H∂Hlx+kyl+q=0 （4-36） ∂x∂yy

式中，lx，ly是边界表面的法相余弦。 （3）对于不透水边界，应满足：

∂H

=0，式中n为边界的外法线方向。 ∂n2、符合达西定律的各向异性连续介质中的三维空间稳定渗流，应满足下列渗流 的基本微分方程：

∂⎛∂H⎞∂⎛∂H⎞∂⎛∂H⎞kx+⎜ky+Q=0 (4-37) ⎟+⎜kz⎜⎟⎟∂x⎝∂x⎠∂y⎝∂y⎠∂z⎝∂z⎠

式中kx、ky、kz渗流介质在x、y、z方向的渗流系数；Q是单位边界表面的流出或流入量；H水头函数，即：

H=

P

+y ρg式中P计算点的压力；P水的密度；g重力加速度；y计算点的纵坐标高度。 公式((4-37)的边界条件为:

(1)边界面上的水头H为已知，即H=Ho

(2)已知边界面上渗流的流出或流入量，此时边界面上应满足下列流量平 衡方程式：

kx

∂H∂H∂Hlx+kyly+kzl+q=0 (4-38) ∂x∂y∂zz

式中lx、ly、lz是边界表面的法向余弦；Q是单位边界表面的流出或流入量。

(3) 对于不透水边界面，应满足下列条件：

∂H

=0 (4-39) ∂n式中n不透水边界的外法线方向

(4)在渗流的自由表面上应满足条件H=y. 3、在二维的非稳定渗流场中，渗流运动可用方程式(4-40)来描述，即:

37

成都理工大学硕士学位论文

μ为坝体土料的排水度，也称为给水度，即单位土体在饱和含水情况下水位降

落后在重力作用下，土的孔隙中自由排出的水量体积与土体积之比。当自由水面线的变化比较平缓时，上式中的H可用含水层的平均水头H来代替，则公式(4-40)可变为:

⎛∂H⎞∂⎛∂H⎞⎤∂HH⎡∂⎥ (4-41) =⎢⎜kx+⎜ky⎟⎟⎜⎟⎜⎟∂tμ⎢∂x∂x⎠∂y⎝∂y⎠⎥

⎣⎝⎦

−

⎛∂H⎞∂⎛∂H⎞⎤∂H1⎡∂⎥ (4-40) =H⎢⎜kx+⎜ky⎟⎟⎜⎟⎜⎟∂tμ⎢∂x∂x⎠∂y⎝∂y⎠⎥

⎣⎝⎦

上式的边界条件为:

(1)当边界面上的水头H己知时，H=Ho .

(2)己知边界面上渗流的流出或流入量，此时边界面上应满足下列流量平衡方 程式:

kx

∂H∂Hlx+kyl+q=0 ∂x∂yy

式中lx\\ly，边界表面的法向余弦q单位边界表面的流出或流入量。

(3)于不透水边界应满足下列条件:

∂H

=0 ∂n式中n是边界的外法线方向。

(4)渗流的自由表面上应满足条件：H=y.

4.3径向基函数计算渗流问题

上面介绍了渗流的模型是:设想渗流区内的全部土粒骨架不存在，整个渗流区的全部空间是被水所充满的连续流动，把土壤对水流的作用概化成作用与水流的力。由于渗流模型将渗流看作是连续空间内的连续介质的运动，因此可以按水力学的分类方法，将渗流分为恒定渗流和非恒定渗流，均匀渗流和非均匀渗流，渐变渗流和急变渗流，有压渗流和无压渗流以及空间(三维)渗流、平面(二维)渗流和线性(一维)渗流等。

渗流的渗流运动的基本规律，长期以来，一直受到许多水利科学工作者的广泛关注，并通过不断的研究与讨论，取得了卓有成效的进展，在工程生产实践中发挥了重要作用。特别是近几年来，随着计算机技术及数值计算方法的迅猛发展，越来越多的学者开始采用数学模型来解决渗流问题，并取得了显著突破。但因各种数值计算方法对渗流计算所得到的解为数值近似解，且各种方法所建立的理论有明显差别，因此对某一具体问题而言不同数值方法计算的结果在收敛性、稳定

38

第4章 渗流问题数学模型及径向基函数插值法界二维渗流问题

性等方面会有所区别。因此，开展二维渗流数值解法的对比研究在工程实践中具有重要意义[16]。

渗流问题数值解法常用的主要是有限差分法、有限元和有限体积法等，它们在理论上已比较成熟，是常用的解偏微分方程的方法，无网格法解偏微分方程是几十年新兴起的方法，在上面已经介绍了无网格法解偏微分方程方法之一径向基函数法，通过上面的数值试验已经验证了径向基函数法解偏微分方程可行性，只要径向基函数中参数取适当的值和在计算区域取适当的点的计算精度非常高，下面把它应用于解实际渗流问题。

本文我们渗流问题考虑是二维稳定渗流，在二维稳定渗流的计算中，本文直接径向基函数插值法解二维稳定渗流，其控制方程是拉普拉斯方程的形式。 计算实例

f(u)=k∇2u(x,y)+Q=0 x,y∈Ω （4-42）

u=u x,y∈Γ1 （4-43） n⋅k∇u(x)+a0u(x)=q x,y∈Γ2 （4-44）

∂2u(x,y)∂2u(x,y)2

+其中u(x)为水头，k为渗透系数，∇u(x,y)=,Γ1为定22∂x∂y

水头边界，Γ2为流量边界（在本文研究中，取隔水边界）。u和q分别是Γ1和Γ2上给定的水头及流量。n 为相应边界的外法向矢量。Ω为待求解区域。Q为子域内的流量，抽水取负.

本文主要是探讨径向基函数插值法渗流研究中的应用问题，因此选择稳 定流的渗流情况，渗流的控制方程和定解条件，方程变为LaplaceQ＝0，a0＝0，方程即为

−

−

∂2u(x,y)∂2u(x,y)

=0 in Ω （4-45） + ∇u(x,y)=22∂x∂y

u=f in ∂Ω （4-46） 边界极区域Ω满足下列条件

∂Ω:={x∈IR2:x=R(t),t∈[0,2π]}，R:[0,2π]→IR,2π−周期曲线 R(t)=ρ(t)(cos(t),sin(t)),t∈[0,2π]，0<ρ(t)≤R,t∈[0,2π]。

2

ρ(t)=

c2+sin(t/2)/c2+1,c=0.3216，在边界点

f(x,y)=y

径向基函数用MQ（参数C=2.0）

39

成都理工大学硕士学位论文

图4-2 节点分布图

图4-2 径向基函数无网格法数值解图

40

第4章 渗流问题数学模型及径向基函数插值法界二维渗流问题

图4-3 解析解

图 4-4 误差图

上面算例采用不规则节点，结果来看计算精度也十分高，采用规则和不规则节点对计算误差影响在文献[18]有详细的分析。

从上面简单的径向基函数插值法解简单的Laplace方程结果看误差是相对非常小，数值解和解析解十分吻合，所以此方法解这个简单渗流问题精度比较高。从图4-4结果看误差在边界尤其在点（1，1）附近相对高一些，但是本文只是通过二维稳定流的拉普拉斯方程来验证方法的可行性。实际工作中，大多数渗流问题是非稳定恒流，其模型也很复杂多为三维模型，所以引入各种渗流模型，对此方法进行改进研究。

41

成都理工大学硕士学位论文

本文研究工作及创新点

本文主要研究无网格法之一径向基插值法，下面是是本文的主要研究工作： （1）首先对径向基函数插值及解偏微分方程方法介绍，并且利用径向基函数插值法拟合余弦函数及其一导、二导的图像，并且在区域内取不同的点分析各自的误差，随着节点越多其各自的总体误差越来越小，但是其二导图像拟合时在某点误差相对大些；其误差主要集中在凸峰、凹峰及两端点上。

（2）利用此方法与MATLAB中插值函数分别对二维高斯函数进行插值拟合，结合显示此方法比MATLAB中插值函数法拟合精度高很多，其误差值主要集中在峰谷附近及边界附近。

（3）利用径向基函数插值法解泊松方程，取不同的径向基函数与解析解比较，结果显示数值解的精度都十分高。

（4）利用MQ径向基函数解二维区域Hemholtz方程，分别取不同均匀节点进行数值计算，从计算误差结果是节点相对多时误差非常小，但是节点太多误差反而会稍微有所增高，所以在利用此方法解偏微分方程时要取适当的点。 （5）利用径向基函数法解经典偏微分方程，从过分析不同的径向基函数对计算精度分析及取不同的离散点对误差影响，取经典MQ径向基函数和适当的离散点应用于解一个简单的渗流问题（Laplace方程），结果显示精度很高，将径向基函数插值法应用于渗流问题数值模拟中是非常有意义的尝试，因为这种方法的易用性和精度都较传统网格方法有了较大提高，但在一些复杂的实际工程问题有待进一步的研究和探讨。

通过上文研究及数值试验本文创新点主要体现在以下几点：

首先用径向基函数插值法这种无网格法通过取不同节点，应用于余弦函数和它的一阶导数和二阶导数的图像拟合，取得很好的拟合效果，表明它应用在余弦函数图像拟合时精度非常高； 创新点之二是将径向基函数插值法应用于解二维恒定渗流问题，取得很好的结果，表明径向基函数插值法无网格法应用到解渗流问题可行性；创新点之三通过数值试验提出径向基函数解偏微分方程时，表明了节点数要取适当，效果才最好，不是在计算区域取节点越多计算精度就越高。

42

结论和展望

结论和展望

本文只是用无网格法的其中之一径向基函数插值法，利用它求解了几个简单例子，由于时间比较紧只是研究了几个简单例子，利用全局径向基函数插值法拟合余弦函数和其一阶导数、二阶导数时候，在区域取100均匀节点函数精度非常高，其二阶导数精度相对低一些，但是从整体来看效果还是非常好；二维函数时径向基函数插值法比MATLAB网格插值精度比较，结果前者精度很高；利用径向基函数法解偏微分方程时要在计算区域取适当的节点，径向基函数的参数要去适当的值，计算节点的多少和参数的选取对计算精度都有一定的影响；最后，把径向基函数插值法解解偏微分方程通过二维稳定流的拉普拉斯方程来验证方法的可行性。大部分的径向基函数中都带有参数,而现在也没有形成一种理论来指导如何有效来选取参数,这也是径向基函数方法在实际应用中有待解决的问题.尽管人们都在尝试使用径向基函数解决实际问题,但是对于它的数学理论基础还很不完善,虽然许多专家学者在这方面做出了一些贡献,但仍有大量的数学问题还没有解决,而关于方面的数学理论的研究就更少了,这就需要我们关注并尝试解决这方面的问题.径向基函数插值法作为一种比较新的微分方程数值求解方法,目前还不完善,还必须考虑如何实现和传统数值方法相结合,使二者的优势互补.相信经过广大研究人员的不懈努力,这种方法会越来越完善,在实际应用中会发挥越来越大的作用。

本文只是通过二维稳定流的拉普拉斯方程来验证方法的可行性。实际工作中，可以引入各种渗流模型，尤其是不稳定二维和三维渗流问题，要进行并行计算，对此方法进行改进。

通过计算发现，影响计算精度的因素有很多。影响域的半径及域中节点的数目，节点的分布，基函数中参数的影响等。目前还没有确定这些因素的一般方法。本文采用 matlab 编程，计算效率较 c 与 fortran 语言编程低，改进程序仍有减少计算时间的可能。

无网格法是最近 10 年才新兴起来的一种数值模拟方法。虽然无网格法本 身仍在理论探讨中阶段，但因其在数值模拟过程不需要建立网格，无网格法在

应用研究中展现了巨大的潜力。

43

成都理工大学硕士学位论文

致 谢

时间过的飞快，转眼间研究生三年即将过去，学校、老师和同学给我留下了太多美好的回忆，在这毕业之际，借此向曾经给予我帮助和支持的人们表示最衷心的感谢。

感谢我的导师龚灏教授，在学习和科研方面给我莫大的帮助，给我提供了良好的学习条件和实践锻炼机会，他敏锐的洞察力、渊博的知识、严谨的治学态度、精益求精的工作作风精神给我留下了刻骨铭心的印象，这些都使我受益匪浅。恩师严谨的治学态度、对我的严格要求以及为人处世的坦荡胸襟将使我终身收益，是我一生学习的榜样。

在生活方面龚老师也给了我极大的帮助，使我没有任何的后顾之忧，安心学习。在我感到迷茫的时候，龚老师给了我许多宝贵的建议。对于龚老师对我各方面帮助和悉心教诲，我会永远铭记于心无限感激。在攻读硕士的这三年里，导师不仅为我创造了优越的科研和学习环境，使我得以在数学学科领域中自由翱翔，同时在思想上、人生态度和意志品质方面给予了谆谆教诲，这些教益必将激励着我在今后的人生道路上奋勇向前。

谨此向我的导师龚灏教授致以衷心的感谢和崇高的敬意！ 感谢周仲礼老师三年来在学习和生活上的帮助和鼓励。

感谢上海大学上海市应用数学与力学所长江学者钱跃竑教授和赵鹏师兄在我最困难的时候总能给予我安慰和鼓励，令我重拾信心成功考上上海大学博士。

感谢我的爸妈、大爷大娘和斌哥、斌嫂、弟弟和妹妹这些年来一如既往无尽的关心、爱护和支持，亲恩似海。

感谢我所有的师兄弟妹们，从他们身上我真切地感受到兄弟姐妹般的情谊以及相互学习、团结合作的乐趣，真诚的感谢你们。

感谢信息管理学院郭科教授及所有的老师们、在我研究生期间给予的学习和生活方面的帮助，你们传授给我的专业知识是我不断成长的源泉，使我学到很多东西，使我受益终生。再次的感谢你们。

还有许多帮助过我的人，他们对我的支持和帮助在我的记忆深处留下了深刻的印象。在此无法一一罗列，但对他们，我始终心怀感激。

最后，我要感谢所有曾经关心过我、鼓励过我、帮助过我的亲人、朋友、同学、师长，在此向他们致以深深的谢意。

44

参考文献

参考文献

[1] 李荣华,冯果忱.微分方程数值解法[M].北京:高等教育出版社,1980 [2] 陆金甫,关治.偏微分方程数值解法[M].北京:清华大学出版社,2004

[3] 张文生.科学计算中的偏微分方程有限差分法[M].北京:高等教育出版社,2006 [4] 吴宗敏.散乱数据拟合的模型、方法和理论[M] 科学出版社. [5]李庆扬，王能超，易大易。数值分析[M] 清华大学出版社。

[6]孙庆新,等. 数值分析[M] . 沈阳:东北工学院出版社,1990 :115 - 121. [7]李岳生,黄友谦. 数值逼近[M] . 北京:人民教育出版社,1978 :12 - 20. [8]张 雄，刘岩 无网格法[M] 北京：清华大学出版社，2004

[9]G.R.LIU Y.T.GU著 王建明 周学军 译 无网格法理论及程序设计 [M] 山东大学出

版社.

[10]刘鹤年 主编 流体力学 [M] 中国建筑工业出版社. [11]王沫然 MATLAB6.0与科学计算[M] 电子工业出版社 [12]吴望一 编著 流体力学 北京大学出版社

[13]Edward B. Magrab 等著 MATLAB原理与工程应用 电子工业出版社

[14]陈荣华 径向基函数拟插值理论及其在微分方程数值解中的应用[D] 复旦大学博士论

文

[15]向 波 非规则区域有限五点格式计算渗流问题[D] 武汉大学硕士

[16] 李晓春 小湾水电站坝肩岩体裂隙网络三维渗流网络与无网格法耦合模型渗流研究

[D]吉林大学

[17] 李明辉, 王桂杰。用径向基函数的方法求解常微分方程[J] 沈阳化工学院学报

第20 卷第1 期 2006.3 :68-71.

[18] 秦伶俐 黄文彬 周 喆.径向基函数在无单元方法中的应用[J] 中国农业大学学报

2004 ,9 (6) : 80～84.

[19] Chen, W. & Tanaka, M., New Insights into Boundary-only and Domain-type RBF Methods.

Int. J. Nonlinear Sci. & Numer. Simulation.2000,1(3), 145-151．

[20] Chen, W. & He, W., A note on radial basis function computing, Int. J. Nonlinear

Modelling Sci. & Engng. 2001, 59-65 ．

[21] Franke, C. and Schaback, R., Solving partial differential equation by collocation using radial

basis function, Appl. Math. Comput. 93(1998), 73-82.

[22]Niederreiter, H., Random Number Generation and Quasi-Monte Carlo Methods, SIAM,

CBMS 63, Philadelphia, PA, 1992.

[23] Chen, C.S., Golberg, M.A. and Hon, Y.C., The method of fundamental solutions and

quasi-monte-carlo method for diffusion equations, Int. J. Numer. Meth. Engng., 43(1998),

45

成都理工大学硕士学位论文

1421-1435.

[24]Chen, W. and Tanaka, M, New insights in boundary-only and domain-type RBF methods,

Inter. J. Nonlinear Sci Numer. Simulation, 3(2000), 145-152.

[25] Kansa, E.J. and Hon, Y.C., Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations, Comput. Math. Appls., 39(2000) 123-137.

[26]Duchon, J., Splines minimizing rotation invariant semi-norms in Sobolov spaces, in: “Constructive Theory of Functions of Several Variables”, Springer-Verlag, Berlin, 1976. [27] Powell, M.J.D., Recent advances in Cambridge on radial basis functions, The Second Inter.

Dortmunt Meeting on Approximation Theory, 1998.

[28] Fedoseyev, A.I., Friedman, M.J. and Kansa, E.J., Improved multiquadratic method for elliptic

partial differential equations via PDE collocation on the boundary, Comput. Math. Appl., (submitted), 2000.

[29]X.Zhang,K.Z.song,M.W.Lu,X.Liu Meshless methods based on collocation with rasia

function, Computational Mechanics 26(2000) 333-343

46