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

复习：

1.数值计算方法的含义

2.误差及误差限

3.误差与有效数字

4.数值计算中应注意的问题

第二章插值方法

一.插值的含义

问题提出：

函数y=/（式）在n+1个点七,x，…,X”上的函数值,求任意一点£的函数值

”0。

说明：函数），=/（"可能是未知的；也可能是的，但它比拟复杂，很难计算其函数值

”0。

解决方法：

构造一个简单函数P（x）来替代未知（或复杂）函数y=/（x）,那么用P（x'）作为函数

值/（X'）的近似值。

二、泰勒（Taylor］插值

1.问题提出：

复杂函数y=/（x）在点的函数值/（x（J，求与附近另一点马+〃的函数值/（%）+〃）。

2.解决方法：

构造一个代数多项式函数己（可，使得勺（力与/（力在1=与点充分逼近。

泰勒多项式为：

月（/）=/（/）+，（%）（工7。）+2^（%7。）2+-+2^^（工7。）"

2!n\

显然，8（x）与/（x）在x=/点，具有相同的i阶导数值（i=0』,…,n）。

3.几何意义为：

匕（x）与/（冗）都过点（％，/&））；

乙W与/（汇）在点（，％，/（、））处的切线重合；

《（X）与/（X）在点■）,/（戈0））处具有相同的凹凸性；

其几何意义可以由以下图描述，显然函数/3（力能相对较好地在/点逼近f（x）0

4.误差分析［泰勒余项定理〕：

.（力-/（力|=）‘其中g在与与x之间。

5.举例：

函数/（6=«,求/（115）。

分析：此题理解为，“复杂”函数/"）=«在/=100点的函数值为/（$）=10,求公

的附近一点与+15的函数值/（/+15）。

解：

<

⑴构造1次泰勒多项式函数［（式）：^（X）=/（Xo）+/（A^）（X-X0）o

其中/（%）=〃100）=10,r（x）="n）=r（ioo）q,那么有：

4（x）=5+0.05%

故有了（115）。6（115）=10.75

误差分析：

出（115）—5（115）1=华⑴5—100）2

函数/〃⑴在［100/15］区间绝对值的极大值为了〃（100）=2.5、10'4,

那么有：

|^（115）-/（115）|<0.028125<0.05

于是近似值10.75有三位有效数字。

几何意义:显然，:3也过点(100,10),且.(/)就是函数/(力=石在点(100,10]

处的切线，如以下图所示。

(2)构造2次泰勒多项式函数E(x)：

.(x)=/(Xo)+./(Xo)(To)+，

把/(100)=10,/(100)=-!-及/〃(100)=2.5x107代入，有

20

/(115)(115)=10.721875o

分析误差

函数二(“在[100/15]区间绝对值的极大值为尸(100)=3.75、10飞，那么有

|^(115)-/(115)|<0.002109375<0.005

于是近似值10.721875有四位有效数字。

运行文件taylor.m：

%函数f(x)=xNl/2),求f(115)

%一次泰勒插值

subplot(1,2,1);

f=inline(,xA(l/2),);

p1=inlinc('5+0.05*x');

fplot(f,[-5(),30()]);

holdon

fplot(p1,(-50,300]);

plol(115,10.75,'*')

line([l15,115],[0,10,75])

%二次泰勒插值

subplot(l,2,2);

p2=inline('10+1/20*(x-l00)-1/4000/2*(x-l00)A2,);

fplut(f,[-3O,3OO]);

holdon

fplot(p2,[-30,300]);

plot(115,10.72；*')

line([l15,115],[0,10,72])

可以得到以以下图形：

6.泰勒插值存在的问题：

1.函数必须存在升1阶导函数，即使存在n+1阶导数，计算的工作量也比拟大;

2.要求h为个小量，假设h较大，那么计算的误差就很大。

三.拉格朗日〔Lagrange〕插值

1.问题提出：

函数N=/（R）在n+1个点玉上的函数值为,y,,y„,求任意一点V的函数值

说明：函数y=/（x）可能是未知的；也可能是的，但它比拟复杂，很难计算其函数值

八0。

2.解决方法：

构造一个n次代数多项式函数4（外来替代未知（或复杂）函数y=/（耳，那么用月,（一）

作为函数值/（V）的近似值。

设C（x）=q）+4x+%/+…，构造2（%）即是确定n+1个多项式的系数

%，q，/，，a八o

3.构造的依据:

当多项式函数?(x)也同时过的n+1个点时，我们可以认为多项式函数匕(“逼近于原

来的函数/(可。根据这个条件，可以写出非齐次线性方程组:

%。A)++…+。,吊”=%

n

《)+q内+/x：+—+aflx]=>i

旬十+%工；+,••+=X,

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

故当n+l个点的横坐标八Z各不相同时，方程组系数矩阵的行列式D不等于零，故

方程组有唯一解。即有以下结论。

结论：当的n+1个点的横坐标飞,公各不相同时，那么总能够构造唯一的n次多项式

函数匕(x),使号(力也过这n+1个点。

4.几何意义

5.举例：

函数/(x)=«,求

分析：此题理解为，“复杂”函数/(力=«,当x=81/00」21,144时，其对应的函数

值为：y=9,10,11,12,当x=115时，求函数值/(115)。

解：

(1)线性插值：过的(100,10)和(121,11)两个点，构造1次多项式函数出力，于

是有

1(x)=x—121-x—100.

x1()Hx11

100-121-----121-100

那么

/(115)«(115)=10.71428571428572。

(2)抛物插值：构造2次多项式函数6(x)，使得它过的(100,10)、(⑵,11)和］144,12)

三个点。于是有2次拉格朗日插值多项式：

月(加上吗工/。+占/四、xU+7上。。哈

2、7(100-121)(100-144)(121-100)(121-144)(144-100)(144-121)

那么有

/(115)«^(115)=10.72275550536420

6.拉格朗日n次插值多项式公式:

P(v}=(x-xja-w)(17“)

“&T|)(%)T2)…(/—怎)°

(x-x0)(x-x2)(l〃)

+

(%-%0)(%一丹)("X")>1

(—)(-—%)(--J"

匕(X)=4)(X)乳+/I(x)y++/“(X)”=WX(X)果

Jt=o

其中田力称为基函数(k=0,l,…n),每一个基函数都是关于x的n次多项式，其表达

式为：

♦(〃)=门T―

j*k

拉格朗日公式特点：

1.把每一点的纵坐标,“单独组成一项；

2.每一项中的分子是关于x的n次多项式，分母是一个常数；

3.每一项的分子和分母的形式非常相似，不同的是：

分子是(了—,)，而分母是(迎-…)

7.误差分析〔拉格朗日余项定理〕

其中4在所界定的范围内。

针对以上例题的线性插值，有

|/］(115)-/(115)|=1^^1(115-100)(115-121)

函数/(x)在［100/15］区间绝对值的极大值为/〃(100)=2.5x107,

那么有：

|^(115)-/(115)|<0.01125<0.05

于是近似值/(115卜爪115)=10.71428571428572有三位有效数字。

针对以上例题的抛物线插值，有

|/^(115)-/(115)|=R^^(115-100)(115-121)(115-144)

函数/历⑴在［1()0［15］区间绝对值的极大值为―(18))=3.75、1()2,那么有

|^(115)-/(115)|<0.00163125<0.(X)5

于是近似值“115卜4(115)=10.72275550536420有四位有效数字。

8.拉格朗日插值公式的优点

公式有较强的规律性，容易编写程序利用计算机进行数值计算。

9.拉格朗日插值通用程序

程序流程图如下：

y

文件lagrange.m如下:

%拉格朗n插值

closeall

n=inputC的坐标点数n=?');

x=input('xl,x2,...,xn=？')；

y=input(,yl,y2,…,yn=？');

xx二input('插值点=？');

symst%定义t为符号量

P=0;

fork=l:n

1=1；

forj=l:k-l

l=l*(t-x(j))/(x(k)-x(j));

end

forj=k+l:n

l=l*(t-x(j))/(x(k)-x(j));

encl

p=p+l*y(k);

end

p=inline(p);%把符号算式p变为函数形式

fplot(p,[min(min(x),xx)T,max(max(x),xx)+l]);%画多项式函数

holdon

p(xx)%显示插值点

plot(x,y,'o',xx,p(xx),';%画点和插值点

在MATLAB命令窗口输入：

lagrange

然后有以下对话过程和结果，

的坐标点数n二?6

xl,x2,...,xn=?[l,3,5,7,9,11]

yl,y2,...,yn=?[-l,20,0,-1,12,3]

插值点二?8

ans

有以以下图形:

10.作业

1,函数sin(x)过以下数据点:

X0.791.01.6

sin(x)0.7103530.8414710.999574

请用线性插值和抛物插值，计算sin(0.63)的值，并分析误差。

四.牛顿〔Newton〕插值

复习：

⑴问题提出：函数在n+1个点的值(xO,yO),(xl,yl),.…(xn,yn),求当x=X，时，4的值。

(2)解决方法：构造n次多项式函数勺(耳，使它也过的n+1个点。

⑶拉格朗日公式:2(4)=方乙(“瓦，4(x)=n7

A=0y=oZ—

出

(4)拉格朗日公式的优点：结构规律性强，便于编写程序。

〔5〕拉格朗日插值的缺点：无承袭性〔继承性〕

假设算出3点的抛物插值精度不够，再进行4点的3次多项式插值时，必须从头算起,

前面算出的3点抛物插值的计算结果不能利用。

而泰勒插值却是具有承袭性的，如线性插值的结果不精确，那么再加上一项，就变成

了泰勒抛物插值，如：

泰勒1次插值：《(％)=/(/)+/(%)(%-%)

泰勒2次插值：鸟=+（犬7。）+务1（犬7。）2。

4•

而牛顿插值就是具有承袭性的插值公式

1.差商的概念

设n+1个点如和…，怎互不相等,那么定义:

/(七)-/(外)

天和X.（Z丰j）两点的一阶差商为：/[%勺]=

X「Xj

七,Xj,五三点的二阶差商为：/X）,七]=

茗一/

/[药，"£卜/[乙—]

天,弓,演,七四点的三阶差商为：/[西月,%,为]

XT

n+1个点/小,…,小的n阶差商为:

/[xXy•x]=/[为‘"卜/卜占'，5]

差商具有对称性：/[xpxy]=/[x;,x.]；/[七,孙占卜/[x/c]

2.牛顿插值解决的问题与拉格朗日插值解决的问题相同

只是表述n次多项式尺（力的公式不同。

3.牛顿插公式的推导

根据差商的概念，有：

/(x)=/C%)+/[x/o](X7。)/[X,内）]是X,,%两点的一阶差商；

/卜毛]=/卜2[]+/[乂与闻（X-%）...是心飞,西三点的二阶差商;

/[戈，/0,一•，*|]=/&，与，…，％]+/[%，飞，不…，']（工一人）

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

/(x)=/(x0)+/[x0,x,](x-x0)+/[x0,x,,x2](.r-x0)(x-x,)+

…+/[不。，内，・・，七』(工一工。)(1一%)・・,(工一七1)

+/[x,x0,x1,.-,x,J(x-x0)(x-xl)...(x-x„)

J(X)=%(%)+居(x)

其中匕(%)=/(%)+/[x"l(xf)+/闻玉闯(工一豌)(工一9)+

..•+/[^0^1,...,A,J(x-x0)(x-x1).-.(x-xZI_I)

尼（x）=/[x为，5，…，x”]（x-Xo）（x-X），（工一乙）

以上勺（6的表达式称为牛顿插值公式，可以证明，n次牛顿插值多项式与n次拉格朗

日插值多项式完全相同，只是表达形式不同。

故，拉格朗日余项定理与牛顿余项定理相同:

R；（力=%（力-〃刈=/jo小，…,七』口（工-占）=,

Jt=o[〃十1卜k=o

其中4在X（）,X|,…,所界定的范围内。

那么有公式：小,"，…闯=皆令

4.牛顿插值差商表

xiyi一阶差商二阶差商n阶差商

xOyo1

xlyif[xO,xl](x-xO)

x2y2f[xl,x2]f[xO,xl,x2](x-xO)(x-xl)

x3y3f[x2,x3]f[xl,x2,x3](x-xO)…(x-x2)

••••♦•

xn-1yn-1

xnynf[xn-1,xn]f[xn-2,xn-1,xn]♦♦・f[xO,,,,,xn](x-xO)•••(x-xn-1)

5.举例

例1：函数f（x）当X=-2「1,O,1时，其对应函数值为f（x）=13,-8,-l,4。求f（0.5）的

值。

解：根据点，填写以下差商表：

xiyi一阶差商二阶差商三阶差商*

-2131

-1-8-21(x+2)

0-1714(x+2)(x+1)

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

那么四点三次牛顿插值多项式4(“为：

/J(x)=13—21(x+2)+14(x+2)(x+l)-5(x+2)(x+l)x

故，〃0.5卜4(0.5)=13—21(0.5十2)十14(0.5十2)(0.5十1)-5(0.5十2)(0.5十1)0.5=3.625

可以在MATLAB下运行程序newtonOl.m：

p3=inlineC13-21*(x+2)+14*(x+2)*(x+1)-5*(x+2)*(x+1)*x*);

f?lot(p3,[-2.5,2.5]);

holdon

xi=[-2,-1,0,1];

yi=[13,-8,-1,4];

plot(xi,yi,';

plot(0.5,p3(0.5),'o');

可以得到以以下图形：

601-------1--------!--------1--------1--------1--------1--------1--------1--------1--------1

20卜\T

、

•20-]

-40--

例2：函数f(x)当x=-2,7,0,1,2时，其对应函数值为f(x)=13,-8,T,4,1。求f(0.5)的

值。

解：该题目与例1相比，就是多了一个点，所以和例1的差商表相比，只需多一列，多一

行：

xiyi一阶差商二阶差商二阶差商四阶差商*

-2131

-1-8-21(x+2)

0-1714(x+2)(x+1)

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

21-3-4-11(x+2)(x+1)x(x-1)

而5个点的4次牛顿插值多项式乙(x)是在《(x)的根底上多增加1项：

=13-21(x+2)+14(x+2)(x+1)-5(JT+2)(X+l)x+(x+2)(^+l)x(x-l)

那么

/(O.5)«/>(0.5)=13-21(0.54-2)+14(0.54-2)(0.5+1)-5(0.5+2)(0.5+1)0.54-(0.5+2)(0.54-1)0.5(0.5-1)

=2.6875

可以在MATLAB卜.运行程序r.ewton02.m：

p4=inline('13-21本(x+2)+14*(x+2)木(x+l)-5*(x+2)木(x+1)*x+(x+2)*(x+1)*x*(x_l)');

fplot(p4,[-2.5,2.5],'r');

holdon

xi=[-2,-l,0,1,2];

yi=[13,-8,-1,4,1]；

plot(xi,yi,';

plot(0.5,p4(0.5),'o');

可以得到以以下图形：

6,牛顿插值的优点

(1)具有承袭性质

(2)利用差商表，计算多点插值，比拉格朗口公式计算方便。

7.牛顿插值算法的通用程序

以下是程序流程图：

开始

I

MATLAB的通用程序newton.m为:

%牛顿插值

closeall

n=input('的坐标点数n=？')；

x二input('xl,x2,...,xn=?’)；

y=input('yl,y2,…,yn=?');

xx=input('插值点=?');

%计算差商：f[xl,x2],f[xl,x2,x3],...,f[xl,x2,...,xn]

f=y；

fol-i=l：n-l%计算笫i阶差商

fork=n:-l:i+1

f(k)=(f(k)-f(k-1))/(x(k)-x(k-i));

end

end

symst%定义t为符号量

P=f(1)；

fork=2:n

1=1；

forj=l:k-l

1=1*(t-x(j));

end

p=p+l*f(k);

end

p=inline(p);%把符号算式p变为函数形式

fplot(p,[min(min(x),xx)-l,max(max(x),xx)+1]);%画多项式函数

holdon

p(xx)/显示插值点

plot(x,y,'o',xx,p(xx),';用画点和插值点

在MATLAB命令窗口输入：

newton

然后有以下对话过程和结果，

的坐标点数n二?6

xl,x2,...,xn=?[l,3,5,7,9,11]

yl,y2,...,yn=?[-l,20,0,-1,12,3]

插值点二？8

ans=

有以以下图形:

40

8.作业

(1)过(0,6),(1,7),(2,20),(3,81),(4,250)五个点做多项式函数p(x),并求p(-2)

的值。

(2)给出以下函数表，函数f(x)是一个多项式函数，试求其次数及x的最高募的系数。

X012345

f(x)-7-452665128

(3)请写出下面数列中？的值

①2,5,9,15,23,?

②2,8,15,29,50,?,125

五埃尔米特〔Hermite〕插值

1.问题提出

函数y=/(x)在n+l个点/,内,…，毛上的函数值此加〃及一阶导函数值

为‘，，；，…，)；求任意一点x的函数值/(x)o

2.解决方法：

构造一个2n+l次代数多项式函数旦川(犬)，使得

&+G)=y,A1

2=0,1,…，〃

即，多项式函数8e(x)也过这n+1个点，且函数f(x)和鸟用(X)在这n+1个点上具有

相同的切线。

3.埃尔米特插值公式:

当节点横坐标各不相同时，存在唯一的n+1次代数多项式函数总川（x）：

乌13=汽「1-2（X7£）〃（尤）”+之3娱

k=0L」k=Q

其中,

〃Y—X"1

；=0Xk-Xjj=oxk~Xj

j"j*k

4.举例

例1.求满足以下条件的埃尔米特插值多项式。

玉12

23

%

91-1

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

也可也….）力［2门

y

3一/卜「与

+（…MV'

把表中值代入，得：75（x）=-2x3+8x2-9x+5

例2.函数f（©=—二满足以下数据表:

X+2

12

1/30.2

X

90-0.14

构造3次埃尔米特插值多项式A（X）。

解：根据埃尔米特插值公式可以构造6"）为:

D/19316

PAx)、=——x------x

-150251075

在Matlab命令窗口输入:

f二inline('x/(x-3+2)');

p3=in1ine(*19/150*x3-16/25*x2+9/10*x-4/75,);

fplot(f,[0,3]);

holdon

fplot(p3,[0,3],'r');

plot([1,2],[1/3,0.2],'*');

绘出如以下图形

例3.求二次多项式满足£（%）=%,2（X）=X，6'（%）=%'。其中Xo，N，No，X，No'为

常数。

叱Y—X

解：设巴（力=4）（大）为+/"%）,+4（/）为'，根据条件有4（x）=n

7=0xk~xj

，0（工0）=1(%)=。4)(犬。)=。

」0（X1）=0».0=1,«4(A；)=0,

，0’（与）=。4(%)=i

于是基函数/。（“一定含有因子（x-xj,基函数“X）一定含有因子（工-%）2,基函数

4（同一定含有因子（不一与）（，一内）。

设，o（x）=（X-X］）（0¥+〃），那么有

，o(%)=(%-%)("+。)=1

(or()+人)+〃(％)-xj=0

解得:

1,b二XTq

&(毛一百)(王一/)

(x0-x1)(x,-x0)

那么有:

皿「*雅音，叱。")=

六分段插值

1.龙格(Runge)现象[高次多项式插值的缺陷〕

针对函数/(x)=选取6个节点：

xi：[-5,-3,-1,1,3,5]；

yi：[1/26,0.1,0.5,0.5,0.1,1/26]

可以构造5次多项式函数1(x)

假设选项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次多项式函数E)(x)

利用用拉格朗日插值的通用程序(或牛顿插值的通用程序)可以画出f(x),P5(x)和

P10(x)的图形。程序rungc.m如下:

%用拉格朗日插值公式分析龙格现象

closeall

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;

symst%定义t为符号量

P=0；

fork=l:u

1=1；

forj=l:k-l

1=1*(t-x(j))/(x(k)-x(j));

end

forj=k+l:n

l=l*(t-x(j))/(x(k)-x(j));

end

p=p+l*y(k);

end

p5=inline(p);为把符号算式P变为函数形式

f=inlineCl/(x"2+D，);

fplot(f,[-5,5]);%画原来的函数

holdon

fplot(p5,[-5,5],'g');%画5次多项式函数

n=ll;

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];

symst%定义t为符号量

P=0;

fork=l:n

1=1；

forj=l:k-l

l=l*(t-x(j))/(x(k)-x(j));

end

forj=k+l:n

l=l*(t-x(j))/(x(k)-x(j));

end

p=p+l*y(k);

encl

plO=iiiline(p);%把符号算式p变为函数形式

fplot(plO,[-5,5]/^);%画10次多项式函数

legend('f(x)','P_5(x)','P_l_0(x)')

plot(x6,y6,';%画6个节点

plot(x,y,'o');%画10个节点

plot([-5,5],[0,0],'k');%画坐标轴

plot([0,0],[-0.5,2],'。);

运行该程序，可以绘制出如以下图形：

从图中可以看出，随着节点的增加，采用高次多项式插值，可以在某些区域较好的逼

近原来的函数（如在［-2,2］区间）；但在高次多项式的两端出现了剧烈震荡的现象，这就是

所谓的龙格现象。

从该图可以看出，在x=±5附近时,.（力与f（x）偏离很远。例如为（4.8）=1.8044,

而/（4.8）=0.0416o这就说明用高次插值多项式月（到来近似f（x）的效果并不好，因而通常

不用高次插值。

2.分段线性插值

当节点较多时，可以采用分段线性插值，公式如下：

x-x;

0（力=——-

%一%%1f

以上公式即为两点Q，x）,（心,为）的线性拉格朗日插值公式。

例如针对龙格现在的函数/（"=廿口选取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('l/(x-2+l)');

fplot(f,[-5,5]);%画原来的函数

holdon

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)

可以得到以以下图形：

显然和龙格现象相比，分段线性插值函数夕（司比6（"和儿（X）都能更好的逼近原函数

f（x）O

3.三次样条插值

（1）样条函数的概念

分段线性插值在节点处没有连续的一阶导函数，其光滑性较差。对于飞机的机翼的型

线及船舶型往往要求有二阶光滑度（即在节点处要求二阶导函数连续）。

样条函数的概念来源于工程设计的实践。所谓“样条”（spline）是早期工程设计中的

一种绘图工具，它是富有弹性的细长条。绘图时，用压铁迫使样条通过指定的型值点，并

俣证样条的光滑外形。在绕或