数值计算方法教案-插值方法.doc

复习：
1.数值计算方法的含义
2．误差及误差限
3.误差与有效数字
4.数值计算中应注意的问题
第二章 插值方法
一．插值的含义
问题提出：
已知函数在n+1个点上的函数值，求任意一点的函数值。
说明：函数可能是未知的；也可能是已知的，但它比较复杂，很难计算其函数值。
解决方法：
构造一个简单函数来替代未知（或复杂）函数，则用作为函数值的近似值。
二、泰勒（Taylor）插值
1.问题提出：
已知复杂函数在点的函数值，求附近另一点的函数值。
2.解决方法：
构造一个代数多项式函数，使得与在点充分逼近。
泰勒多项式为：

显然，与在点，具有相同的i阶导数值（i=0,1,…,n）。
3.几何意义为：
与都过点；
与在点处的切线重合；
与在点处具有相同的凹凸性；
其几何意义可以由下图描述，显然函数能相对较好地在点逼近。

4.误差分析（泰勒余项定理）：
，其中在与之间。
5.举例：
已知函数，求。
分析：本题理解为，已知“复杂”函数在=100点的函数值为，求的附近一点+15的函数值。
解：
（1）构造1次泰勒多项式函数：。
其中，，，则有：

故有
误差分析：

函数在[100,115]区间绝对值的极大值为，
则有：

于是近似值10.75有三位有效数字。
几何意义：显然，也过点（100,10），且就是函数在点（100,10）处的切线，如下图所示。
（2）构造2次泰勒多项式函数：
。
把，及代入，有
。
分析误差

函数在[100,115]区间绝对值的极大值为，则有

于是近似值10.721875有四位有效数字。
运行文件taylor.m：
%已知函数f(x)=x^(1/2)，求f(115)
%一次泰勒插值
subplot(1,2,1);
f=inline('x^(1/2)');
p1=inline('5+0.05*x');
fplot(f,[-50,300]);
hold on
fplot(p1,[-50,300]);
plot(115,10.75,'*')
line([115,115],[0,10.75])

%二次泰勒插值
subplot(1,2,2);
p2=inline('10+1/20*(x-100)-1/4000/2*(x-100)^2');
fplot(f,[-30,300]);
hold on
fplot(p2,[-30,300]);
plot(115,10.72,'*')
line([115,115],[0,10.72])
可以得到以下图形：

6.泰勒插值存在的问题：
1.函数必须存在n+1阶导函数，即使存在n+1阶导数，计算的工作量也比较大；
2.要求h为个小量，若h较大，则计算的误差就很大。

三．拉格朗日（Lagrange）插值
1.问题提出：
已知函数在n+1个点上的函数值，求任意一点的函数值。
说明：函数可能是未知的；也可能是已知的，但它比较复杂，很难计算其函数值。
2.解决方法：
构造一个n次代数多项式函数来替代未知（或复杂）函数，则用作为函数值的近似值。
设，构造即是确定n+1个多项式的系数。
3.构造的依据：
当多项式函数也同时过已知的n+1个点时，我们可以认为多项式函数逼近于原来的函数。根据这个条件，可以写出非齐次线性方程组：

其系数矩阵的行列式D为范德萌行列式：

故当n+1个点的横坐标各不相同时，方程组系数矩阵的行列式D不等于零，故方程组有唯一解。即有以下结论。
结论：当已知的n+1个点的横坐标各不相同时，则总能够构造唯一的n次多项式函数，使也过这n+1个点。
4.几何意义

5．举例：
已知函数，求。
分析：本题理解为，已知“复杂”函数，当x=81,100,121,144时，其对应的函数值为：y=9,10,11,12，当x=115时，求函数值。
解：
（1）线性插值：过已知的（100,10）和（121,11）两个点，构造1次多项式函数，于是有

则
。
（2）抛物插值：构造2次多项式函数，使得它过已知的（100,10）、（121,11）和（144,12）三个点。于是有2次拉格朗日插值多项式：
则有
10.72275550536420
6.拉格朗日n次插值多项式公式：

其中称为基函数（k=0,1,….,n），每一个基函数都是关于x的n次多项式，其表达式为：

拉格朗日公式特点：
1．把每一点的纵坐标单独组成一项；
2.每一项中的分子是关于x的n次多项式，分母是一个常数；
3.每一项的分子和分母的形式非常相似，不同的是：
分子是，而分母是
7.误差分析（拉格朗日余项定理）
，
其中在所界定的范围内。
针对以上例题的线性插值，有

函数在[100,115]区间绝对值的极大值为，
则有：

于是近似值有三位有效数字。
针对以上例题的抛物线插值，有

函数在[100,115]区间绝对值的极大值为，则有

于是近似值10.72275550536420有四位有效数字。
8.拉格朗日插值公式的优点
公式有较强的规律性，容易编写程序利用计算机进行数值计算。
9. 拉格朗日插值通用程序
程序流程图如下：

文件lagrange.m如下：
%拉格朗日插值
close all
n=input('已知的坐标点数n=?');
x=input('x1,x2,...,xn=?');
y=input('y1,y2,...,yn=?');
xx=input('插值点=？');
syms t %定义t为符号量
p=0;
for k=1:n
l=1;
for j=1:k-1
l=l*(t-x(j))/(x(k)-x(j));
end
for j=k+1:n
l=l*(t-x(j))/(x(k)-x(j));
end
p=p+l*y(k);
end
p=inline(p); %把符号算式p变为函数形式
fplot(p,[min(min(x),xx)-1,max(max(x),xx)+1]); %画多项式函数
hold on
p(xx) %显示插值点
plot(x,y,'o',xx,p(xx),'*'); %画已知点和插值点
在MATLAB命令窗口输入：
lagrange
然后有以下对话过程和结果，
已知的坐标点数n=?6
x1,x2,...,xn=?[1,3,5,7,9,11]
y1,y2,...,yn=?[-1,20,0,-1,12,3]
插值点=？8
ans =
5.67187500000000
有以下图形：

10.作业
1，已知函数sin(x)过以下数据点：
x
0.79
1.0
1.6
sin(x)
0.710353
0.841471
0.999574
请用线性插值和抛物插值，计算sin(0.63)的值，并分析误差。
四．牛顿（Newton）插值
复习：
（1）问题提出：已知函数在n+1个点的值(x0,y0),(x1,y1),….(xn,yn),求当x=x’时，y’的值。
（2）解决方法：构造n次多项式函数，使它也过已知的n+1个点。
（3）拉格朗日公式：，
（4）拉格朗日公式的优点：结构规律性强，便于编写程序。
（5）拉格朗日插值的缺点：无承袭性（继承性）
若算出3点的抛物插值精度不够，再进行4点的3次多项式插值时，必须从头算起，前面算出的3点抛物插值的计算结果不能利用。
而泰勒插值却是具有承袭性的，如线性插值的结果不精确，那么再加上一项，就变成了泰勒抛物插值，如：
泰勒1次插值：
泰勒2次插值：。
而牛顿插值就是具有承袭性的插值公式
1.差商的概念
设n+1个点互不相等，则定义：
和两点的一阶差商为：
,三点的二阶差商为：
,四点的三阶差商为：
……
n+1个点的n阶差商为：

差商具有对称性：；
2.牛顿插值解决的问题与拉格朗日插值解决的问题相同
只是表述 n次多项式的公式不同。
3.牛顿插公式的推导
根据差商的概念，有：
…………………是两点的一阶差商；
……是三点的二阶差商；
……

把以上各式从后向前逐次代入，可以得到：

其中

以上的表达式称为牛顿插值公式，可以证明，n次牛顿插值多项式与n次拉格朗日插值多项式完全相同，只是表达形式不同。
故，拉格朗日余项定理与牛顿余项定理相同：
，
其中在所界定的范围内。
则有公式：
4.牛顿插值差商表
xi
yi
一阶差商
二阶差商

n阶差商
*
x0
y0

1
x1
y1
f[x0,x1]

(x-x0)
x2
y2
f[x1,x2]
f[x0,x1,x2]

(x-x0)(x-x1)
x3
y3
f[x2,x3]
f[x1,x2,x3]

(x-x0)…(x-x2)
…
…

xn-1
yn-1

xn
yn
f[xn-1,xn]
f[xn-2,xn-1,xn]
…
f[x0,…,xn]
(x-x0)…(x-xn-1)

5.举例
例1：已知函数f(x)当x=-2,-1,0,1时，其对应函数值为f(x)=13,-8,-1,4。求f(0.5)的值。
解：根据已知点，填写以下差商表：

xi
yi
一阶差商
二阶差商
三阶差商
*
-2
13

1
-1
-8
-21

(x+2)
0
-1
7
14

(x+2)(x+1)
1
4
5
-1
-5
(x+2)(x+1)x
则四点三次牛顿插值多项式为：

故，=3.625
可以在MATLAB下运行程序newton01.m：
p3=inline('13-21*(x+2)+14*(x+2)*(x+1)-5*(x+2)*(x+1)*x');
fplot(p3,[-2.5,2.5]);
hold on
xi=[-2,-1,0,1];
yi=[13,-8,-1,4];
plot(xi,yi, '*');
plot(0.5,p3(0.5),'o');
可以得到以下图形：

例2：已知函数f(x)当x=-2,-1,0,1,2时，其对应函数值为f(x)=13,-8,-1,4,1。求f(0.5)的值。
解：该题目与例1相比，就是多了一个点，所以和例1的差商表相比，只需多一列，多一行：
xi
yi
一阶差商
二阶差商
三阶差商
四阶差商
*
-2
13

1
-1
-8
-21

(x+2)
0
-1
7
14

(x+2)(x+1)
1
4
5
-1
-5

(x+2)(x+1)x
2
1
-3
-4
-1
1
(x+2)(x+1)x(x-1)
而5个点的4次牛顿插值多项式是在的基础上多增加1项：

则可以在MATLAB下运行程序newton02.m：
p4=inline('13-21*(x+2)+14*(x+2)*(x+1)-5*(x+2)*(x+1)*x+(x+2)*(x+1)*x*(x-1)');
fplot(p4,[-2.5,2.5],'r');
hold on
xi=[-2,-1,0,1,2];
yi=[13,-8,-1,4,1];
plot(xi,yi, '*');
plot(0.5,p4(0.5),'o');
可以得到以下图形：

6.牛顿插值的优点
（1）具有承袭性质
（2）利用差商表，计算多点插值，比拉格朗日公式计算方便。
7.牛顿插值算法的通用程序
以下是程序流程图：

MATLAB的通用程序newton.m为：
%牛顿插值
close all
n=input('已知的坐标点数n=?');
x=input('x1,x2,...,xn=?');
y=input('y1,y2,...,yn=?');
xx=input('插值点=？');
% 计算差商：f[x1,x2],f[x1,x2,x3],...,f[x1,x2,...,xn]
f=y;
for i=1:n-1 % 计算第i阶差商
for k=n:-1:i+1
f(k)=(f(k)-f(k-1))/(x(k)-x(k-i));
end
end
syms t %定义t为符号量
p=f(1);
for k=2:n
l=1;
for j=1:k-1
l=l*(t-x(j));
end
p=p+l*f(k);
end
p=inline(p); %把符号算式p变为函数形式
fplot(p,[min(min(x),xx)-1,max(max(x),xx)+1]); %画多项式函数
hold on
p(xx) %显示插值点
plot(x,y,'o',xx,p(xx),'*'); %画已知点和插值点

在MATLAB命令窗口输入：
newton
然后有以下对话过程和结果，
已知的坐标点数n=?6
x1,x2,...,xn=?[1,3,5,7,9,11]
y1,y2,...,yn=?[-1,20,0,-1,12,3]
插值点=？8
ans =
5.67187500000000
有以下图形：

8.作业
（1）过(0,6)，(1,7)，(2,20)，(3,81)，(4,250)五个点做多项式函数p(x)，并求p(-2)的值。
（2）给出下列函数表，已知函数f(x)是一个多项式函数，试求其次数及x的最高幂的系数。
x
0
1
2
3
4
5
f(x)
-7
-4
5
26
65
128
（3）请写出下面数列中？的值
① 2,5,9,15,23,？
② 2,8,15,29,50,？,125
五 埃尔米特（Hermite）插值
1．问题提出
已知函数在n+1个点上的函数值及一阶导函数值，求任意一点的函数值。
2.解决方法：
构造一个2n+1次代数多项式函数，使得

即，多项式函数也过这n+1个点，且函数f(x)和在这n+1个点上具有相同的切线。
3. 埃尔米特插值公式：
当节点横坐标各不相同时，存在唯一的n+1次代数多项式函数：

其中，
，
4. 举例
例1．求满足下列条件的埃尔米特插值多项式。

1
2

2
3

1
-1
解：根据埃尔米特插值公式有：

把表中值代入，得：
例2．已知函数满足下列数据表：

1
2

1/3
0.2

0
-0.14
构造3次埃尔米特插值多项式。
解：根据埃尔米特插值公式可以构造为：

在Matlab命令窗口输入：
f=inline('x/(x^3+2)');
p3=inline('19/150*x^3-16/25*x^2+9/10*x-4/75');
fplot(f,[0,3]);
hold on
fplot(p3,[0,3],'r');
plot([1,2],[1/3,0.2],'*');
绘出如下图形

例3.求二次多项式满足，，。其中为已知常数。
解：设，根据已知条件有
，，，
于是基函数一定含有因子，基函数一定含有因子，基函数一定含有因子。
设，则有

解得：
，
则有：
，，

六 分段插值
1．龙格(Runge)现象（高次多项式插值的缺陷）
针对函数选取6个节点:
xi：[-5,-3,-1,1,3,5]；
yi：[1/26,0.1,0.5,0.5,0.1,1/26]
可以构造5次多项式函数
若选项11个节点
xi：[-5,-4,-3,-2,-1,0,1,2,3,4,5]；
yi：[1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26]
可以构造10次多项式函数
利用用拉格朗日插值的通用程序（或牛顿插值的通用程序）可以画出f(x),P5(x)和P10(x)的图形。程序runge.m如下：
%用拉格朗日插值公式分析龙格现象
close all
n=6;
x=[-5,-3,-1,1,3,5];
y=[1/26,0.1,0.5,0.5,0.1,1/26];
x6=x;
y6=y;
syms t %定义t为符号量
p=0;
for k=1:n
l=1;
for j=1:k-1
l=l*(t-x(j))/(x(k)-x(j));
end
for j=k+1:n
l=l*(t-x(j))/(x(k)-x(j));
end
p=p+l*y(k);
end
p5=inline(p); %把符号算式p变为函数形式
f=inline('1/(x^2+1)');
fplot(f,[-5,5]); %画原来的函数
hold on
fplot(p5,[-5,5],'g'); %画5次多项式函数

n=11;
x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];
y=[1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26];
syms t %定义t为符号量
p=0;
for k=1:n
l=1;
for j=1:k-1
l=l*(t-x(j))/(x(k)-x(j));
end
for j=k+1:n
l=l*(t-x(j))/(x(k)-x(j));
end
p=p+l*y(k);
end
p10=inline(p); %把符号算式p变为函数形式
fplot(p10,[-5,5],'r'); %画10次多项式函数
legend('f (x)','P_5(x)','P_1_0(x)')
plot(x6,y6,'*'); %画6个已知节点
plot(x,y,'o'); %画10个已知节点
plot([-5,5],[0,0],'k'); %画坐标轴
plot([0,0],[-0.5,2],'k');
运行该程序，可以绘制出如下图形：

从图中可以看出，随着节点的增加，采用高次多项式插值，可以在某些区域较好的逼近原来的函数（如在[-2,2]区间）；但在高次多项式的两端出现了激烈震荡的现象，这就是所谓的龙格现象。
从该图可以看出，在附近时，与f(x)偏离很远。例如，而。这就说明用高次插值多项式来近似f(x)的效果并不好，因而通常不用高次插值。
2．分段线性插值
当节点较多时，可以采用分段线性插值，公式如下：
，
以上公式即为两点的线性拉格朗日插值公式。
例如针对龙格现在的函数选取11个节点:
xi：[-5,-4,-3,-2,-1,0,1,2,3,4,5]；
yi：[1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26]
可以构造10个1次多项式函数，即分段线性函数。
在MATLAB命令窗口输入：
f=inline('1/(x^2+1)');
fplot(f,[-5,5]); %画原来的函数
hold on
x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];
y=[1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26];
plot(x,y,'o')
plot(x,y,'r')
可以得到以下图形：

显然和龙格现象相比，分段线性插值函数比和都能更好的逼近原函数f(x)。
3．三次样条插值
（1）样条函数的概念
分段线性插值在节点处没有连续的一阶导函数，其光滑性较差。对于飞机的机翼的型线及船舶型往往要求有二阶光滑度（即在节点处要求二阶导函数连续）。
样条函数的概念来源于工程设计的实践。所谓“样条”（spline）是早期工程设计中的一种绘图工具，它是富有弹性的细长条。绘图时，用压铁迫使样条通过指定的型值点，并保证样条的光滑外形。在绕度不大的情况下，样条的曲线即为三次样条函数。
（2）几何意义

（3）构造三次样条函数的理论分析
如上图所示，通过已知的六个点，构造5个三次多项式函数分别是：红色、蓝色、黑色、紫色和绿色5根曲线。为确定一根曲线，就需要确定4个待定系数，所以总共需要4*5=20个待定系数。
另外，分析需要的约束条件。
每一根函数都要过已知的左右两个点，则有5*2=10个约束条件。
此外，每两个相邻曲线在相邻点处要求充分光滑，即在连接点处左右两个函数在该点具有1次和2次的导函数连续，图中有4个“中间点”，故又有4*2=8个约束条件。
若在整个图形的两端在加2个约束条件，整个3次样条函数就确定了。如：
①左右两端点上的1阶导函数已知；
②左右两端点上的2阶导函数已知，如（称为自然边界条件）；
③若原来的函数f(x)是以xn-x0为周期的周期函数，则y0=yn，且。
（4）用MATLAB函数interp1进行三次样条函数的插值
例1.对龙格现象中的函数进行11个点的三次样条插值：
x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];
y=[1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26];
xi=-5:0.01:5;
yi=interp1(x,y,xi,'spline');
plot(xi,yi,x,y,'o')
hold on
f=inline('1/(x^2+1)');
fplot(f,[-5,5],'r')
可以绘出下图：

红色为原来的函数f(x)，蓝色为通过曲线f(x)上的11个点而构成的三次样条函数。从图中可以发现三次样条函数很好的地描述了函数f(x)。
例2.随机构造15个点，用牛顿法（或拉格朗日法）构造14次代数多项式函数，然后再根据这15个点构造3次样条函数。
现在在MATLAB命令窗口输入：
x=1:15;
y=round(10*randn(1,15));
然后运行程序newton，有以下对话过程：
已知的坐标点数n=?15
x1,x2,...,xn=?x
y1,y2,...,yn=?y
插值点=？2
ans =
0
xi=1:0.01:15;
yi=interp1(x,y,xi,'spline');
plot(xi,yi,'r')
可以绘出如下图形

显然图中，可以看到龙格现象，如果，在另一个图中重新画3次样条函数：
close
plot(xi,yi,x,y,'o')
可以得到：

七、小结
本章学习了
1．泰勒插值
2．拉格朗日插值
3.牛顿插值
4.埃尔米特插值
5.龙格现象
6.分段线性插值
7.分段三次插值(3次样条函数)

作业：
（1）已知单调连续函数y=f(x)的下列数据
xi
-1.1
0.0
1.2
2.1
yi
-2.2
-1.l
1.0
2.1
用插值法计算，x为多少时，f(x)=0。
提示：把xi和yi“颠倒”理解。
（2）用插值法计算矩阵的特征多项式。
提示：为3次多项式函数，故让分别取0,1,2,3时，求出的函数值，再构造3次多项式函数。