《地球物理反演》
——读书报告
专业所在院(系、部)核工程技术学院
研究生姓名郭猛猛
学号 2010070807
专业名称固体地球物理学
日期 2011年6月30日
《地球物理反演》读书报告
1、拉格朗日插值法
在数值分析中,拉格朗日插值法是以法国十八世纪数学家约瑟夫·路易斯·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律,而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测,在若干个不同的地方得到相应的观测值,拉格朗日插值法可以找到一个多项式,其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日(插值)多项式。数学上来说,拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。拉格朗日插值法最早被英国数学家爱德华·华林于1779年发现,不久后(1783年)由莱昂哈德·欧拉再次发现。1795年,拉格朗日在其著作《师范学校数学基础教程》中发表了这个插值方法,从此他的名字就和这个方法联系在一起。
定义:
对于给定的若n+1个点,对应于它们的次数不超过n的拉格朗日多项式只有一个。如果计入次数更高的多项式,则有无穷个,因为所有与相差的多项式都满足条件。
对某个多项式函数,已知有给定的k + 1个取值点:
其中x j对应着自变量的位置,而y j对应着函数在这个位置的取值。
假设任意两个不同的x j都互不相同,那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为:
其中每个为拉格朗日基本多项式(或称插值基函数),其表达式为:
拉格朗日基本多项式的特点是在x j上取值为1,在其它的点
上取值为0。
证明:
存在性:
对于给定的k+1个点:,拉格朗日插值法的思路是找到一个在一点x j取值为1,而在其他点取值都是0的多项式。这样,多
项式在点x j取值为y j,而在其他点取值都是0。而多项式
就可以满足
在其它点取值为0的多项式容易找到,例如:
它在点x j取值为:。由于已经假定x i两两互不相同,因此上面的取值不等于0。于是,将多项式除以这个取值,就得到一个满足“在x j取值为1,而在其他点取值都是0的多项式”:
这就是拉格朗日基本多项式。
唯一性:
次数不超过k的拉格朗日多项式至多只有一个,因为对任意两个次数不超过k的拉格朗日多项式:P1和P2,它们的差P1?P2在所有k+1个点上取值都是0,因此必然是多项式的倍数。因此,如果这个差P1?P2不等于0,次数就一定不小于k+1。但是P1?P2是两个次数不超过
k的多项式之差,它的次数也不超过k。所以P
?P2 = 0,也就是说P1 = P2。
1
这样就证明了唯一性。
优点与缺点:
拉格朗日插值法的公式结构整齐紧凑,在理论分析中十分方便,然而在计算中,当插值点增加或减少一个时,所对应的基本多项式就需要全部重新计算,于是整个公式都会变化,非常繁琐。这时可以用重心拉格朗日插值法或牛顿插值法来代替。此外,当插值点比较多的时候,拉格朗日插值多项式的次数可能会很高,因此具有数值不稳定的特点,也就是说尽管在已知的几个点取到给定的数值,但在附近却会和“实际上”的值之间有很大的偏差。这类现象也被称为龙格现象,解决的办法是分段用较低次数的插值多项式。
2、牛顿迭代法
牛顿迭代法(Newton's method)又称为牛顿-拉夫逊方法
(Newton-Raphson method),它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。多数方程不存在求根公式,因此求精确根非常困难,甚至不可能,从而寻找方程的近似根就显得特别重要。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。牛顿迭代法是求方程根的重要方法之一,其最大优点是在方程f(x) = 0的单根附近具有平方收敛,而且该法还可以用来求方程的重根、复根,此时线性收敛,但是可通过一些方法变成超线性收敛。另外该方法广泛用于计算机编程中。
设r是f(x) = 0的根,选取x0作为r初始近似值,过点(x0,f(x0))做曲线y = f(x)的切线L,L的方程为y = f(x0)+f'(x0)(x-x0),求出L与x轴交点的横坐标x1 = x0-f(x0)/f'(x0),称x1为r的一次近似值。过点(x1,f(x1))做曲线y = f(x)的切线,并求该切线与x轴交点的横坐标x2 = x1-f(x1)/f'(x1),称x2为r的二次近似值。重复以上过程,得r的近似值序列,其中x(n+1)=x(n)-f(x(n))/f'(x(n)),称为r的n+1次近似值,上式称为牛顿迭代公式。
解非线性方程f(x)=0的牛顿法是把非线性方程线性化的一种近似方法。把f(x)在x0点附近展开成泰勒级数f(x) = f(x0)+(x-x0)f'(x0)+(x-
x0)^2*f''(x0)/2! +… 取其线性部分,作为非线性方程f(x) = 0的近似方程,即泰勒展开的前两项,则有f(x0)+f'(x0)(x-x0)=0 设f'(x0)≠0则其解为x1=x0-f(x0)/f'(x0) 这样,得到牛顿法的一个迭代序列:x(n+1)=x(n)-
f(x(n))/f'(x(n))。
军人在进攻时常采用交替掩护进攻的方式,若在数轴上的点表示A,B 两人的位置,规定在前面的数大于后面的数,则是A>B,B>A交替出现。但现在假设军中有一个胆小鬼,同时大家又都很照顾他,每次冲锋都是让他跟在后面,每当前面的人占据一个新的位置,就把位置交给他,然后其他人再往前占领新的位置。也就是A始终在B的前面,A向前迈进,B跟上,A把自己的位置交给B(即执行B = A操作),然后A 再前进占领新的
位置,B再跟上……直到占领所有的阵地,前进结束。像这种两个数一前一后逐步向某个位置逼近的方法称之为迭代法。
迭代法也称辗转法,是一种不断用变量的旧值递推新值的过程,跟迭代法相对应的是直接法(或者称为一次解法),即一次性解决问题。迭代算法是用计算机解决问题的一种基本方法。它利用计算机运算速度快、适合做重复性操作的特点,让计算机对一组指令(或一定步骤)进行重复执行,在每次执行这组指令(或这些步骤)时,都从变量的原值推出它的一个新值。
利用迭代算法解决问题,需要做好以下三个方面的工作:
一、确定迭代变量。在可以用迭代算法解决的问题中,至少存在一个直接或间接地不断由旧值递推出新值的变量,这个变量就是迭代变量。
二、建立迭代关系式。所谓迭代关系式,指如何从变量的前一个值推出其下一个值的公式(或关系)。迭代关系式的建立是解决迭代问题的关键,通常可以使用递推或倒推的方法来完成。
三、对迭代过程进行控制。在什么时候结束迭代过程?这是编写迭代程序必须考虑的问题。不能让迭代过程无休止地重复执行下去。迭代过程的控制通常可分为两种情况:一种是所需的迭代次数是个确定的值,可以计算出来;另一种是所需的迭代次数无法确定。对于前一种情况,可以构建一个固定次数的循环来实现对迭代过程的控制;对于后一种情况,需要进一步分析出用来结束迭代过程的条件。
3、最小二乘法
1801年,意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后,由于谷神星运行至太阳背后,使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星,但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。
高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中。
法国科学家勒让德于1806年独立发现“最小二乘法”,但因不为世人所知而默默无闻。
勒让德曾与高斯为谁最早创立最小二乘法原理发生争执。
1829年,高斯提供了最小二乘法的优化效果强于其他方法的证明,因此被称为高斯-莫卡夫定理。(来自于wikipedia)
在我们研究两个变量(x, y)之间的相互关系时,通常可以得到一系列成对
的数据(x1, y1.x2, y2... xm , ym);将这些数据描绘在x -y直角坐标系中,若发现这些点在一条直线附近,可以令这条直线方程如(式1-1)。
Y计= a0 + a1 X (式1-1)
其中:a0、a1 是任意实数
为建立这直线方程就要确定a0和a1,应用《最小二乘法原理》,将实测值Yi与利用(式1-1)计算值(Y计=a0+a1X)的离差(Yi-Y计)的平方和
〔∑(Yi - Y计)2〕最小为“优化判据”。
令:φ = ∑(Yi - Y计)2 (式1-2)
把(式1-1)代入(式1-2)中得:
φ = ∑(Yi - a0 - a1 Xi)2 (式1-3)
当∑(Yi-Y计)平方最小时,可用函数φ对a0、a1求偏导数,令这两个偏导数等于零。
(式1-4)
(式1-5)
亦即:
m a0 + (∑Xi ) a1 = ∑Yi (式1-6)
(∑Xi ) a0 + (∑Xi2 ) a1 = ∑(Xi, Yi) (式1-7)
得到的两个关于a0、 a1为未知数的两个方程组,解这两个方程组得出:
a0 = (∑Yi) / m - a1(∑Xi) / m (式1-8)
a1 = [m∑Xi Yi - (∑Xi ∑Yi)] / [m∑Xi2 - (∑Xi)2 )] (式1-9) 这时把a0、a1代入(式1-1)中,此时的(式1-1)就是我们回归的元线性方程即:数学模型。
在回归过程中,回归的关联式是不可能全部通过每个回归数据点(x1, y1. x2, y2...xm,ym),为了判断关联式的好坏,可借助相关系数“R”,统计量“F”,剩余标准偏差“S”进行判断;“R”越趋近于 1 越好;“F”的绝对值越大越好;“S”越趋近于 0 越好。
R = [∑XiYi - m (∑Xi / m)(∑Yi / m)]/ SQR{[∑Xi2 - m (∑Xi / m)2][∑Yi2 - m (∑Yi / m)2]} (式1-10) *
在(式1-1)中,m为样本容量,即实验次数;Xi、Yi分别任意一组实验X、Y的数值。
此即最小二乘的一次完成算法,现代的递推算法,更适用于计算机的在线辨识。最小二乘是一种最基本的辨识方法,但它具有两方面的缺陷[1]:一是当模
型噪声是有色噪声时,最小二乘估计不是无偏、一致估计;二是随着数据的增长,将出现所谓的“数据饱和”现象。针对这两个问题,出现了相应的辨识算法,如遗
忘因子法、限定记忆法、偏差补偿法、增广最小二乘、广义最小二乘、辅助变量法、二步法及多级最小二乘法等。
每一种算法都有其优缺点,选择一种好的算法是对我们对方法本身理解透彻的与否。附:牛顿迭代法Matlab源程序:
1.定义函数
function y=f(x)
y=f(x);%函数f(x)的表达式
function y=z(x)
y=z(x);%函数z(x)的表达式
2.主程序
x=X;%迭代初值
i=0;%迭代次数计算
while i<= I;%迭代次数
y=x-y(x)/z(x);%牛顿迭代格式
if abs(y-x)>ε;%收敛判断
x=y;
else break
end
i=i+1;
end
fprintf('\n%s%.4f\t%s%d','x=',x,'i=',i) %输出结果