插值法综述《计算方法》学习报告.讲义

插值法综述

一、插值法及其国内外研究进展

1.插值法简介

插值法是一种古老的数学方法，它来自生产实践，早在一千多年前，我国科学家在研究历法上就应用了线性插值与二次插值，但它的基本理论却是在微积分产生之后才逐渐完善的，其应用也日益增多，特别是在计算机广泛使用之后，由于航空、机械加工、自动控制等实际问题的需要，使插值法在实践和理论上都显得更为重要，并得到了空前的发展。

2.国内外研究进展

● 插值法在预测地基沉降的应用

● 插值法在不排水不可压缩条件下两相介质的两重网格算法的应用 ● 拉格朗日插值法在地震动的模拟研究中的应用 ● 插值法在结构抗震可靠性分析中的应用

● 插值法在应力集中应变分布规律实验分析中的应用 3.代表性文献

● 不等时距GM(1%2c1)模型预测地基沉降研究 秦亚琼 武汉理工大学学报

(交通科学与工程版) 2008.2

● 不排水不可压缩条件下两相介质的两重网格算法 牛志伟 岩土力学 2008.3 ● 基于拉格朗日插值法的地震动的模拟 白 可 山西建筑 2010.10 ● 响应表面法用于结构抗震可靠性分析 张文元 世界地震工程 1997 ● 小议应力集中应变分布规律的实验方法 查珑珑 淮海工学院学报

（自然科学版）2004.6

二、插值法的原理

【原理】

设有n+1个互不相同的节点（i x ，i y ） （i=0，1，2，...n ）则存在唯一的多项式： 2012()...(1)n

n n L x a a x a x a x =++++

使得()(0,1,2,...)(2)n j j L x y j n ==

证明：构造方程组

20102000201121112012......(3)...n n n

n n n

n n n n a a x a x a x y a a x a x a x y a a x a x a x y

?++++=?++++=???

?++++=?

令：0

01

1111n

n n n

n x x x x A x x ???

??

?=????????

01

n a a X a ??????=??????01n y y Y y ??????=??

???? 方程组的矩阵形式如下：(4)AX Y

=

由于1

10

()0n

n i j i j A x x -===-≠∏∏所以方程组（4）有唯一解。

从而2012()...n n n L x a a x a x a x =++++唯一存在。 三、常用插值法

3.1 Lagrange 插值法

3.1.1 Lagrange 插值法的一般提法

给定))(,(i i x f x ),,1,0(n i =，多项式

∑∏

∑=≠=

=????

?

??--==n i n

i

j j j i j i n

i i i n x x x x y x l y x 000)()(? 称为)(x f 关于n x x x ,,,10 的n 次Lagrange 插值多项式。

3.1.2 Lagrange 插值多项式的构造 已知n+1个节点(,)(0,1,

,j j x y j n =其中j x 互不相同，不妨设

01),n a x x x b =<<<=要求形如：

1

110()n n n n n P x a x a x a x a --=++++ 的插值多项式。

若n 次多项式()(0,1,,)j l x j n =在n+1个节点01n x x x <<

<上满足条件：

1,

;()(,0,1,

,)0,

.

j k k j l x j k n k j =?==?

≠?

就称这n+1个n 次多项式01(),(),,()n l x l x l x 为节点01,,

,n x x x 的n 次插值基函数。

3.1.3 Lagrange 插值法的程序设计 f[x_]:=Exp[x]

A=Table[{x,f[x]},{x,0,0.8,0.2}]//N

g1=ListPlot[Table[A],Prolog->AbsolutePointSize[18]]; Interpolation[A,InterpolationOrder->3] g2=Plot[%[x],{x,0,0.8}] Show[g1,g2] N[%%%[0.12],20] N[%%%%[0.72],20] N[f[0.12],20] N[f[0.72],20]

3.1.4 Lagrange 插值法典型例题及其解法

5===，构造二次拉格朗日插值多项式。 （1

）计算

（2）估计误差并与实际误差相比较。 解

（1）以插值点(27,3), (64,4), (125,5)代入插值公式，得

222()()2000x x j x y l x y i i i

x x i i j i j j i

φ??

?-

?==∑∑∏

?-===

? ?≠?

?

=

(64)(125)(27)(125)(27)(64)

345(2764)(27125)(6427)(64125)(12527)(12564)

x x x x x x ------?+?+?------

(100)

2

(10064)(100125)(10027)(100125)(10027)(10064)

345(2764)(27125)(6427)(64125)(12527)(12564)

?≈------=?+?+?------ 4.68782=

(2) 由误差公式有

(3)()

()(27)(64)(125)3!

f R x x x x ξ=---

记810(3)(3)

3()(),()27

f x f x x f x -==在[27，125]上是单调递减函数。

(3)(3)

5()(27) 5.6450310f x f

-≤≈?

(3)()

(100)(10027)(10064)(100125)0.6181316f R ξ∴≤---≈

2(100)0.04623?=。

3.1.5 Lagrange 插值法误差估计

(1)0

()()()()(),(,)(1)!n n

n n j j f R x f x L x x x a b n ξξ+==-=-∈+∏ (1)

1()n n f

M ξ++≤10

()(1)!n

n n j j M R x x x n +=?≤-+∏ 3.2 Newton 插值法

3.1.1 Newton 插值法的一般提法

]

,,,[)())((],,[))((],[)()()(1011021010000n n n x x x f x x x x x x x x x f x x x x x x f x x x f x N

----++--+-+=

称为Newton 插值多项式。

3.1.2 Newton 插值多项式的构造

()(-)(-)(-)(-)(-)0102010-1

N x a a x x a x x x x L a x x L x x n n n =++++

由插值条件()(0,1,

,)N x f j n n j j == 当0x x =时,()000

N x a f n ==.

当1

x x =时,1()()1010N x a a x x f n =+-=,推得 1000110

f f

a f

a x x

-==-. 当2

x x =时，()()()()20120220212N x a a x x a x x x x f n =+-+--=

2010

2010221

f f f f

a x x

---=

- 则引入记号：[][]10[][,]0110f x f x f f x f x x k k x x -==- [,][,]0201[,,]01221f x x f x x f x x x x x -=-2010

201021f f f f x x x x x x

-----=

- 依次递推可得a k 的一般表达式：[,,...,][,,...]012011[,,...,]011

f x x x x f x x x k k k f x x x k x x

k k ---=--

3.1.3 Newton 插值法的程序设计

{x[0],x[1],x[2],x[3],x[4]}={10,11,12,13,14}; y[k_]:=Log[x[k]] Table[y[k],{k,0,4}]//N; MatrixForm[%]

f[i_,j_]:=(y[j]-y[i])/(x[j]-x[i]) Table[f[i,i+1],{i,0,3}]//N; MatrixForm[%]

f[i_,j_,k_]:=(f[j,k]-f[i,j])/(x[k]-x[i]) Table[f[i,i+1,i+2],{i,0,2}]//N; MatrixForm[%]

f[i_,j_,k_,l_]:=(f[j,k,l]-f[i,j,k])/(x[l]-x[i]) Table[f[i,i+1,i+2,i+3],{i,0,1}]//N; MatrixForm[%]

f[i_,j_,k_,l_,m_]:=(f[j,k,l,m]-f[i,j,k,l])/(x[m]-x[i]) Table[f[i,i+1,i+2,i+3,i+4],{i,0,0}]//N; MatrixForm[%]

A={{y[0],y[1],y[2],y[3],y[4]},{0,f[0,1],f[1,2],f[2,3],f[3,4]}, {0,0,f[0,1,2],f[1,2,3],f[2,3,4]},{0,0,0,f[0,1,2,3],f[1,2,3,4]}, {0,0,0,0,f[0,1,2,3,4]}};

Transpose[A]//N; MatrixForm[%] a[0]=y[0]; a[1]=f[0,1]; a[2]=f[0,1,2]; a[3]=f[0,1,2,3]; a[4]=f[0,1,2,3,4];

N[x]=Sum[a[k]*Product[(x-x[m]),{m,0,k-1}],{k,0,4}]//N Expand[%]

3.1.4 Newton 插值法典型例题及其解法 已知函数()f x 的函数表如下：

求四次牛顿插值多项式，并由此求()596.0f 的近似值。

分析 表中给出六对数据，故最高可构造五次多项式。但由于0.596接近于

40.00=x ，因此可取前五对数据来做差商表。 解 构造差商表如下：

故四次牛顿插值多项式为

()()()()+--+-+=55.04.028000.04.011600

.141075.04x x x x P ()()()()()?--+---55.04.003134.065.055.04.019733

.0x x x x x

()()80.065.0--x x

于是()596.0f ≈()596.04P ＝0.631 95。

3.1.5 Newton 插值法误差估计

(1)

1(1)

()()

()()()(1)!

n n t t t n h f

R x f x N x th n n n

n ξ++++=-+=+ 其中(,).0x x n

ξ∈

四、插值法的比较

Lagrange 插值是利用基函数方法构造的插值多项式，在理论上十分重要，但计算不太方便。基函数方法是将插值问题划归为特定条件下容易实现的插值问题，本质上是广义的坐标系方法。Newton 插值在计算插值多项式及求解函数近似值都比较方便且计算量相对较小，是求函数近似值常用的方法，尤其是等距节点的差分插值公式最为常用。 五、插值法在结构工程专业的应用案例

1.案例叙述

应用Matlab 求解水道测量数据问题

摘要：水道测量数据问题是一个给定数据散乱、随机分布的二维离散数据的插值问题。本文以水道测量数据问题为例，应用Matlab 软件提供的求解三维网格点数据的函数，对求解决给定数据散乱、随机分布的二维离散数据插值问题，给出了一个简便易行的方法。

问题重述

水道测量数据问题是1986年美国大学生数学建模竞赛的A 题，由加州海军研究生院数学系的Richard Franke 提供。问题如下：

在某海域测得一些点（x, y ）处的水深z （单位：英尺）由表1给出，水深数据是在低潮时测得的。船的吃水深度为5英尺，问在矩形区域 (75，200) ? (-50，150) 里的哪些地方船要避免进入。

表1 水道水深测量数据（单位：英尺）

2.案例解法

●假设与问题分析

由题目给出的信息是很少的，除了14个位置的水深之外一无所知。显然，题目要求我们找出水深不到5英尺的区域。为了讨论方便，下面三个假设是合理的：（1）所给数据是精确的；

（2）讨论区域的海底曲面是光滑的，更确切地说，可以认为曲面的一阶、二阶导数是连续的。因为我们可以认为讨论区域为浅水海域，由于长期的海水水流作用，形成的是以砾石或沙为主要组成部分的海底，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形。

（3）水深是一个按区域来划分的变量，在某个位置的水深与其周围区域的水深是相互依赖的，但这种依赖作用随距离的增大而减小。就我们讨论的问题来说，每一个给定数据点影响周围的每一个未知点，一个给定数据点离未知点越近，作用就越大。

根据假设，海底曲面是连续光滑的，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形，因而很自然的想法就是用某种光滑的拟合曲面去逼近已知的14个数据点或以14个已知的数据点为基础，利用二维插值补充一些点的水深，以求得水深不超过5米的区域。

●问题求解

题目中给定的14个已知数据点，是一组散乱、随机分布的二维离散数据集合，一般首先采用改进的Shepard方法，从给定的数据恢复出规则分布点上的数据，然后再应用双三次样条插值或其它的二维数据插值方法来处理。然而，利用Matlab中求解三维网格点数据的函数griddata，却可直接对散乱、随机分布的二维离散数据进行插值。

函数griddata的调用格式为：

[xi, yi, zi] = griddata(x, y, z, xi, yi, 'v4')，

其返回与向量x、y和z所描述的数据点集合相匹配的表面f(x, y) 上网格点的z坐标矩阵zi。函数griddata在点 (xi, yi) 处对表面函数f(x, y) 进行插值，从而得到zi的值。

在此，我们采用求解三维网格点数据的函数griddata，对题目中给出的二维离散数据集合进行插值，作出矩形区域 (75，200) ?(-50，150) 范围内的海底地形图、水深不超过5米的危险区域的海底地貌图、矩形区域 (75，200) ?(-50，150) 范围内的海底等高线图以及水深不超过5米的危险区域的平面图，并求出水深不超过5米的危险海域范围为：[113.75, 200]?[0, 100]。

问题求解的Matlab程序及运行结果附后。

3.求解案例的程序设计

●求解水道测量数据问题的Matlat程序

●clear;

x = [129, 140, 108.5, 88, 185.5, 195, 105.5, 157.5, 107.5, 77, 81, 162, 162, 117.5];

y = [7.5, 141.5, 28, 147, 22.5, 137.5, 85.5, -6.5, -81, 3, 56.5, -66.5, 84,

-38.5];

z = [-4, -8, -6, -8, -6, -8, -8, -9, -9, -8, -8, -9, -4, -9];

nx = 100;

px = linspace(75, 200, nx);

ny = 200;

py = linspace(-50, 150, ny);

[xi, yi] = meshgrid(px, py);

[xi, yi, zi] = griddata(x, y, z, xi, yi, 'v4');

figure(1), meshc(xi, yi, zi+5);

title('(75, 200), (-50, 150) 范围内的海底地形图');

rotate3d

figure(2), contour(xi, yi, zi);

title('(75, 200), (-50, 150) 范围内的海底等高线图');

grid

figure(3), contour(xi, yi, zi, [-5 -5]);

title('水深不超过5米的危险区域的平面图');

grid

[a, b] = find(zi>=-5);

amin = min(a);

amax = max(a);

bmin = min(b);

bmax = max(b);

xmin = 75+((200-75)/100)*bmin

xmax = 75+((200-75)/100)*bmax

ymin = -50+((150+50)/200)*amin

ymax = -50+((150+50)/200)*amax

[i1, j1] = find(zi<-5);

for k = 1:length(i1)

zi(i1(k), j1(k)) = -5;

end

figure(4), meshc(xi, yi, zi);

title('水深不超过5米的危险区域的海底地貌图');

rotate3d

运行结果如下：

xmin = 113.75，xmax = 200，ymin = 0，ymax = 100。

4.方法的推广与探究

本文虽然针对的只是水道测量数据问题，但给出的方法却具有一般性：即利用

Matlab中给出的求解三维网格点数据的函数griddata，可以解决给定数据散乱、随机分布的二维离散数据插值问题，而且方法简便易行。用于水道测量数据问题的求解，其精度与 [2]—[4] 给出的求解方法相当。

六、参考文献