常州大学数值分析上机练习二参考解答

上机练习二参考解答

江苏工业学院数理学院 徐明华

2008-5-21

一、实验目的

1、 熟悉插值法的理论基础和方法基础

2、 了解软件Matlab 中有关函数插值的工具

二、实验内容

1、 对自己给定的函数()y f x =，选取适当的节点，编写MA TLAB 程序，分别给出利用Lagrange 插

值方法、Newton 插值方法确定的逼近()y f x =的插值逼近多项式，并将函数()y f x =、逼近多项式及插值余项的图形画在同一坐标系中。

要求：在数值实验中观察不同的被插函数()y f x =及不同的插值节点选取对插值余项的影响，并结合插值余项定理做适当的解释。

1） 算法介绍

已知列表函数

求一个次数不超过n 的多项式()n P x 逼近上述列表函数，满足条件：

(),0,1,2,,n i i P x y i n == （1）

其中

,0,1,2,,i x i n = 互不相同，称为插值节点；条件（1）称为插值条件；满足条件（1）的多项式()

n P x 称为插值逼近多项式。插值多项式的计算方法如下：

方法一：古典方法 （1） 方法介绍

方法一是一种很自然的想法。简要过程如下： 令

2012()n n n P x a a x a x a x =++++

由插值条件（1）可得待定系数01,,,n a a a 的一个线性方程组：

2010200020112111

2012n

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 ?++++=?++++=??

?

?++++=?

（2） 方程组（2）的系数矩阵的行列式为

2000

21112111n

n n n

n n

x x x x x x x x x

这是一个范德蒙（Vandermonde ）行列式，当,

0,1,2,,i x i n = 互不相同时不为零，因此，方程组

（2）有唯一解。即，我们可以通过解方程组（2）得到满足插值条件的插值多项式()n P x 。 （2） 方法说明

满足插值条件的多项式是存在并唯一的。即，不论用什么方法求出的满足插值条件的多项式本质

上都是一样的，只可能在表现形式上可能有所差异。 插值误差（余项）：

(1)01()

()()()()()()(1)!

n n n n f R x f x P x x x x x x x n ξ+=-=---+

由此插值误差可知，若函数()f x 的高阶导数有界，则增加插值节点个数，提高插值多项式的次

数，可能会降低插值误差。

Runge 现象：插值误差并不总是随着插值多项式的次数的增大而减小。

工作量：通过解方程组（2）确定插值多项式的工作量大约为3

(/3)O n 次乘除法。

根据上述算法分析，下面介绍确定插值多项式()n P x 的另外两种方法，达到克服存在问题，改进算法的目的。

方法二：Lagrange 插值方法 （1）算法简介

Lagrange 插值方法是一种构造性方法，它构造一个满足插值条件的多项式。具体如下： 令

0011()()()()n n n L x l x y l x y l x y =+++ （3） 其中

011011()()()()

()()()()()

i i n i i i i i i i n x x x x x x x x l x x x x x x x x x -+-+----=

----

满足

1,()0,i j j i

l x j i

=?=?≠?

因此，容易验证()n L x 满足插值条件：

(),0,1,,n i i L x y i n == 。

称由（3）式所确定的插值多项式为Lagrange 插值多项式。

说明：根据满足插值条件的插值多项式的唯一性，Lagrange 方法得到的多项式与方法一所得的插值多项式完全一致，可能差异仅仅在多项式的表达形式上。

（2）Lagrange 插值方法分析

优点：Lagrange 插值方法简单，且计算量比方法一少。

确定()n L x 的计算量大概是2

(2)O n 次乘除法；

缺点：（1）设根据列表函数

已经得到Lagrange 插值多项式()n L x ，若在列表函数中增加一个新的节点得到如下列表函数

则Lagrange 插值方法与方法一一样不能方便地利用已有的计算结果得到新的插值多项式1()n L x +

（2）与方法一一样可能会有Runge 现象。

为了克服Lagrange 插值方法存在的缺点(1), 下面我们介绍Newton 插值法。

方法三： Newton 插值方法

（1） 方法思想 构造如下形式的多项式

01020101()()()()()()

n n n N x a a x x a x x x x a x x x x -=+-+--+

+--

若()n N x 满足插值条件

(),0,1,2,,n i i N x y i n == ，

上述形式的插值多项式称为Newton 插值多项式。我们将会看到增加节点（11,n n x y ++）后，新的满足插值条件的插值多项式1()n N x +与()n N x 之间有如下关系

1101()()()()()n n n n N x N x a x x x x x x ++=+--- 。

由此可见在已有()n N x 的基础上确定1()n N x +时，仅仅需要计算1n a +，因此下面我们关心的问题主要是Newton 插值多项式中系数,0,1,2,k a k = ，的计算方法。为此引入差商概念

（2） 差商的概念

对于给定的函数()y f x =，定义 一阶差商：()()[,]i j i j i j

f x f x f x x x x -=

-

二阶差商：[,][,]

[,,]i j j k i j k i k

f x x f x x f x x x x x -=-

……

k 阶差商：

01112010[,,,][,,,]

[,,,]k k k k

f x x x f x x x f x x x x x --=

-

（3） 差商与Newton 插值多项式系数之间的关系

性质1：

010

011()

[,,,]()()()()

k

i k i i i i i k f x f x x x x x x x x x x x =-+=----∑

说明：由性质1可知差商与节点有关，与节点的顺序无关。

性质2：设n 次多项式

01020101()()()()()()

n n n N x a a x x a x x x x a x x x x -=+-+--+

+--

满足插值条件(),0,1,2,,n i i N x y i n == ，则

01[,,,],0,1,2,k k a f x x x k == ，这里000[]()a f x f x ==。

说明：01[,,,],0,1,2,k k a f x x x k == 仅与节点01,,,k x x x 有关，与新增节点无关。这就保证了Newton 插值多项式在增加新的节点时，已有的工作量能够被充分被利用。

（4） 差商表

下面给出差商的计算表，在实际计算过程中，可以逐行计算表中的数据，若新增加节点，则相关差商

可以在已得到的差商表的基础上，比较方便地得到。

0x 0()f x

1

x

1

()f x 0

1

[,]f x x 2x 2()f x 02[,]f x x 012[,,]f x x x

3x 3()f x 03[,]f x x

013[,,]f x x x 0123[,,,]f x x x x

4x

4()f x

04[,]f x x

014[,,]f x x x

0124[,,,]f x x x x

01234[,,,,]f x x x x x

2） 程序

这里仅给出Lagrange 插值方法和Newton 插值方法的Matlab 程序及相应的主调程序。

(1) Lagrange 插值方法的程序及其主程序

将三角函数sin y x =在[,ππ-]上用插值多项式逼近，插值节点分别取成

-3.1416 -1.8850 -0.6283 0.6283 1.8850 3.1416

和

-3.1416 -1.5708 0 1.5708 3.1416

按照实验要求的程序如下：

（a ）主程序

function test_lagrange() x0=-pi:0.4*pi:pi; ％x0=-pi:0.5*pi:pi; y0=sin(x0);

x=-pi:0.05*pi:pi;

[C,L,Y]=lagrange(x0,y0,x) y=sin(x); r=y-Y;

plot(x,Y ,'r--',x,y,'b-',x,r,'k-.','LineWidth',2);

legend('Lagrange polynomial','The original f(x) ','Error',0);

（b ）子程序

function [C,L,Y]=lagrange(x0,y0,x)

%Input -x0 is a row vector that contains a list of abscissas % -y0 is a row vector that contains a list of ordinates % -x is a vector that f(x) will be approximately % computed by using Lagrange polynomial

%Output Y is a vector that contains the approximate value of f(x). % -C is a matrix that contains the coefficients of the % Lagrange interpolatory polynomial.

% -L is a matrix that contains the Lagrange % coefficient of basic interpolating polynomials

if nargin < 2 | nargin > 3

error('Incorrect Number of Inputs'); end

if length(x0)~=length(y0)

error('The length of x0 must be equal to it of y0'); end

w=length(x0); n=w-1;

L=zeros(w,w);

% Form the Lagrange coefficient of basic interpolating polynomials for k=1:n+1 V=1;

for j=1:n+1 if k~=j

if abs(x0(k)-x0(j))

error('Divided by Zero,there are two nodes are the same'); end

V=conv(V,poly(x0(j)))/(x0(k)-x0(j)); end end

L(k,:)=V; end

% Determine the coefficients of the Lagrange % interpolating polynomial C=y0*L; if nargin==3

Y=polyval(C,x); end

(2) Newton 插值方法的程序

下述程序可以根据列表函数，给出Newton 插值多项式的系数、差商表，并能够根据给定的插值点x 计算被插函数在x 处的插值。

下面将函数2

1/(1)y x =+在[5,5-]上用Newton 插值多项式逼近，插值节点分别取成

-5 -4 -3 -2 -1 0 1 2 3 4 5

按照实验要求的程序如下： (a) 主程序 clc;

% test Newton interpolation polynomial %test_Newpoly.m x0=-5:1:5;

y0=1./(1+x0.*x0); x=-5:0.1:5;

[C,D,Y]=newpoly(x0,y0,x); y=1./(1+x.*x); r=y-Y;

plot(x,Y ,'r--',x,y,'b-',x,r,'k-.','LineWidth',2);

legend('Newton polynomial','The original f(x) ','Error',0);

（b ）子程序

function [C,D,Y]=newpoly(x0,y0,x)

%Input -x0 is a row vector that contains a list of abscissas % -y0 is a row vector that contains a list of ordinates % -x is a row vector that f(x) will be approximately % computed by using Newton interpolation polynomial %Output Y is a vector that contains the approximate value of f(x). % -C is a vector that contains the coefficients of the % Newton interpolatory polynomial. % -D is the divided-difference table. if nargin < 2 | nargin > 3

error('Incorrect Number of Inputs'); end

if length(x0)~=length(y0)

error('The length of x0 must be equal to it of y0'); end

n=length(x0); D=zeros(n,n); D(:,1)=y0';

% Use formula to form the divided-difference table for j=2:n for k=j:n

if abs(x0(k)-x0(k-j+1))

error('Divided by Zero,there are two nodes are the same'); end

D(k,j)=(D(k,j-1)-D(k-1,j-1))/(x0(k)-x0(k-j+1)); end end

% Determine the coefficients of the Newton interpolating polynomial C=D(n,n);

for k=(n-1):-1:1

C=conv(C,poly(x0(k))); m=length(C);

C(m)=C(m)+D(k,k); end

if nargin==3

Y=polyval(C,x); end

3） 实验结果

（1）用Lagrange 插值方法插值逼近sin y x =在[,ππ-]上的表达式

（a ）插值节点取成

-3.1416 -1.8850 -0.6283 0.6283 1.8850 3.1416

则结果如下：

图 1

（b ）插值节点取成

-3.1416 -1.5708 0 1.5708 3.1416

比上一实验减少一个节点，则结果如下如下：

图 2

(3) Newton 插值方法的程序

（a ）用Newton 插值多项式逼近函数2

1/(1)y x =+在[5,5-]上，插值节点分别取成

-5 -4 -3 -2 -1 0 1 2 3 4 5

的数值结果如下：

图 3

（b ）用Newton 插值多项式逼近函数2

1/(1)y x =+在[5,5-]上，插值节点分别取成

-5.0000 -4.5000 -4.0000 -3.5000 -3.0000 -2.5000 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 4.5000 5.0000

的数值结果如下：

图4

4） 结果分析

（1）用Lagrange 插值方法插值逼近sin y x =在[,ππ-]上的表达式，由于被插函数的高阶导数有界，根据插值余项定理可知，增加节点个数可以提高逼近精度，从上述两个实验结果图1和图2可以看出试验结果符合理论推测。

（2）用插值多项式逼近函数函数2

1/(1)y x =+，x ∈ [5,5]-，从图3及图4（注意图4的纵坐标单位与图3的纵坐标单位的差异）可知，此时，增加插值多项式的次数并不能够达到降低插值误差的目的，出现Runge 现象。

2、 熟悉Matlab\polyfun 下面的一些函数（可以在MA TLAB 命令窗口下执行help polyfun 得到

Matlab\polyfun 下面有哪些函数），并根据自己学习的需要给出其中一些函数的用法。 解：

1）在Matlab 命令窗口下输入：help polyfun 回车可得：

help polyfun

Interpolation and polynomials.

Data interpolation.

pchip - Piecewise cubic Hermite interpolating polynomial. interp1 - 1-D interpolation (table lookup). interp1q - Quick 1-D linear interpolation.

interpft - 1-D interpolation using FFT method. interp2 - 2-D interpolation (table lookup). interp3 - 3-D interpolation (table lookup). interpn - N-D interpolation (table lookup).

griddata - Data gridding and surface fitting.

griddata3 - Data gridding and hyper-surface fitting for 3-dimensional data. griddatan - Data gridding and hyper-surface fitting (dimension >= 2).

Spline interpolation.

spline - Cubic spline interpolation.

ppval - Evaluate piecewise polynomial.

Geometric analysis.

delaunay - Delaunay triangulation. delaunay3 - 3-D Delaunay tessellation. delaunayn - N-D Delaunay tessellation.

dsearch - Search Delaunay triangulation for nearest point. dsearchn - Search N-D Delaunay tessellation for nearest point. tsearch - Closest triangle search. tsearchn - N-D closest triangle search. convhull - Convex hull. convhulln - N-D convex hull. voronoi - Voronoi diagram.

voronoin - N-D V oronoi diagram.

inpolygon - True for points inside polygonal region. rectint - Rectangle intersection area. polyarea - Area of polygon.

Polynomials.

roots - Find polynomial roots.

poly - Convert roots to polynomial. polyval - Evaluate polynomial.

polyvalm - Evaluate polynomial with matrix argument. residue - Partial-fraction expansion (residues). polyfit - Fit polynomial to data. polyder - Differentiate polynomial.

polyint - Integrate polynomial analytically. conv - Multiply polynomials. deconv - Divide polynomials.

2）注释与说明

这部分内容包括一维、二维、高维插值（Data interpolation ）方面的程序，包括样条插值（spline interpolation ）方面的一些函数，还有一些与多项式（Polynomials ）有关的函数等。

3）举例

下面我们利用函数spline （样条插值函数）来处理2

1/(1)y x =+在[5,5-]上的逼近问题。

（a) 插值节点取成

-5 -4 -3 -2 -1 0 1 2 3 4 5

时的主程序及数值结果如下：

% test function spline x0=-5:1:5;

y0=1./(1+x0.*x0); x=-5:0.1:5;

Y=spline(x0,y0,x) y=1./(1+x.*x); r=y-Y;

plot(x,Y ,'r--',x,y,'b-',x,r,'k-.','LineWidth',2);

legend('Cubic Spline','The original f(x) ','Error',0);

图5

（b) 插值节点取成

-5.0000 -4.5000 -4.0000 -3.5000 -3.0000 -2.5000 -2.0000 -1.5000

-1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 2.5000

3.0000 3.5000

4.0000 4.5000

5.0000

时的主程序及数值结果如下：

% test function spline

x0=-5:0.5:5;

y0=1./(1+x0.*x0);

x=-5:0.1:5;

Y=spline(x0,y0,x)

y=1./(1+x.*x);

r=y-Y;

plot(x,Y,'r--',x,y,'b-',x,r,'k-.','LineWidth',2);

legend('Cubic Spline','The original f(x) ','Error',0);

4）实验结果说明

比较图3，4与图5，6可知，样条插值有效地避免了Runge现象，比较图5 与图6 可知，对上述问

题，增加插值节点可以减少样条插值的误差。

图6

3、熟悉网络课程中与本章内容有关的程序，并通过数值实验给出其中一些程序的用法。

课堂上已讲过，从略。

实验结果让我们进一步体会到各种插值方法存在的必要性，并提醒大家在使用各种方法时一定要注意

方法的使用条件。