数值积分第二章插值

计算方法江西理工大学理学院学而时习之不亦说乎

第二章插值与曲线拟合引言:已知熔盐在423～473K的密度和粘度如下表所示，估计450K时的密度和粘度。温度(K)密度(kg/m3)粘度(Pa·S)4231976177.584331967146.514431959122.794531951104.6463194390.26473193478.79

这样，对于函数在区间[a,b]上的各种计算，就用对插值函数的计算取而代之。构造插值函数需要关心下列问题：（1）、插值函数是否存在；（2）、插值函数是否惟一；（3）、如何表示插值函数；（4）、如何估计被插函数与插值函数的误差。

由于代数多项式具有简单和一些良好的特性，故常用代数多项式作为插值函数。

在实际实用中，人们不采用待定系数法，因为：

(1)、计算复杂；

(2)、不容易计算误差。

3、几何意义：（见flash5-1）

§2.1Lagrange插值一、先讨论最简单情形，只有两个节点，的插值多项式(即线性插值)解上述方程组得:

设插值多项式为，且满足插值条件:将

代入插值多项式得：令：即有：则：

是

和

的线性组合。显然，满足：

1i=j0i≠j(i=0,1j=0,1)

是过两点的直线，从几何上看就是过两点的直线来近似代替这种插值称为线性插值。称为线性插值基函数。二、抛物线插值（二次插值） 已知函数的三个点

设插值多项式为且满足插值条件解上述方程组,求出,,,并代入整理后可得令

显然，二次插值多项式可以写成，，的线性组合。即

并且:

满足条件:1i=j0i≠j(i=0,1,2,j=0,1,2)

则称之为二次插值的基函数。 从几何上看，二次插值就是用过三点，

的抛物线来近似代替曲线。因此三点二次插值又称为抛物线插值。三、n+1个结点的插值函数

(2)n+1个结点的Lagrange插值多项式四、插值余项与误差估计：关于余项的几点讨论:例题:化简整理后得:当n=1时，线性插值的余项为当n=2时，抛物插值的余项为五、Lagrange算法:step1:输入插值节点数n,插值点序列(xi,yi),i=0,1,…n,

待计算的函数点x;step2:fori=0ton{{2.1forj=0ton/*对于给定的x计算基函数li(x)*/if(j!=i)temp=temp*(x-xj)/(xi-xj);}2.2fx=fx+temp*yi;}step3:输出Ln(x)的计算结果.一．差商(均差)§2.2Newton插值公式1.定义3.差商的计算(造差商表)----实用!

各级差商的计算可按下表进行2.2Newton插值公式有了差商的概念，前面介绍的线性插值公式可表示为：

(2.2.1)(2.2.1)称为一次Newton插值多项式,记为N1(x),即xi0阶差商一阶差商二阶差商三阶差商四阶差商…x0f(x0)x1f(x1)f[x0,x1]x2f(x2)f[x1,x2]f[x0,x1,x2]x3f(x3)f[x2,x3]f[x1,x2,x3]f[x0,x1,x2,x3]x4f(x4)f[x3,x4]f[x2,x3,x4]f[x1,x2,x3,x4]f[x0,x1,x2,x3,x4]x5f(x5)f[x4,x5]f[x3,x4,x5]f[x2,x3,x4,x5]f[x1,x2,x3,x4,x5]

由上面推导可以看出,Nn(x)至多是n次多项式。另外，由(2.2.2)(2.2.3)

由Newton插值多项式的定义可以看出有如下的递推公式

(2.2.4)

由(2.2.2)式知：每增加一个插值节点，Newton插值多项式只增加一项。这是Newton插值多项式最大的优点,另外,Newton插值多项式的计算量也比Lagrange插值多项式小.

因为,满足插值条件的插值多项式是惟一的,所以,Newton插值余项与Lagrange插值余项应该相等,即

据此可以得到导数与差商的关系三、Newton插值算法:step1:输入插值节点数n,插值点序列(xi,yi),要计算的插值点u.step2:形成差商表g[k],k=0,1,…,n.step3:置初值t=1;newton=f(x0);step4:fori=0ton{t=(u-x(i-1))*t;/*形成(x-x0)…(x-xi-1)*/newton=newton+t*g[i];}step5:输出f(x)=Nn(u)=newton;

clc;formatlong;%显示15位数据

x0=[0.40.50.70.8];%输入初始数

y0=[-0.916291-0.693147-0.3566750.223144];x=0.6;%插值点

n=max(size(x0));y=y0(1);%迭代初始值

disp(y);s=1;dx=y0;

%每次循环用dx记录差商表的一列fori=1:n-1%构造差商表

%%n为节点个数

dx0=dx;

forj=1:n-i

dx(j)=(dx0(j+1)-dx0(j))/(x0(i+j)-x0(j));

end

df=dx(1);

s=s*(x-x0(i));

y=y+s*df;%计算

disp(y);end解：先造差商表

由Newton公式得四次插值多项式为：例2.3在公路建设中,无论在路基或路面结构层施工中,压实度是施工质量控制一个极其重要的指标。要准确测定压实度以及控制施工过程中的材料含水量,必须先通过室内试验确定材料的最大干密度和最佳含水量,而最大干密度和最佳含水量是路基施工密实度控制的关键。实际工作中,通常将材料击实试验的结果通过作图法求得材料的最佳含水量和最大干密度。设有一组实际的击实试验数据如下表所示含水率%干密度g/cm^3含水率%干密度g/cm^311.61.72617.71.81513.51.77619.61.77015.71.825根据表一得如下差商表xf(x)一阶差商二阶差商三阶差商四阶差商11.61.72613.51.7760.02631615.71.8250.022273-0.00098617.71.815-0.005-0.006494-0.00090319.61.770-0.023684-0.0047910.0002790.000148由差商表和牛顿插值公式得到四次插值多项式若要确定最佳含水量可以通过求

牛顿基本插值公式对结点是否等距没有限制。当结点为等距时，牛顿插值公式可进行简化。为此引入差分概念。§2.3等距结点插值公式一.差分1.定义设叫步长.记,为函数在点处以h为步长的一级前向差分,记为称即:一般地,m级差分可以定义为:

(2.3.1)(2.3.2)二、差分表

1.前向差分xkfk=f(xk)△fk△2fk△3fk△4fkx0f0x1f1△f0x2f2△f1△2f0x3f3△f2△2f1△3f0x4f4△f3△2f2△3f1△4f0由式(2.3.1)和式(2.3.2)定义的差分,通常称为向前差分.而称：(2.3.3)(2.3.4)2.后向差分表:xkfk=f(xk)▽fk▽2fk▽3fk▽4fkx0f0x1f1▽f1x2f2▽f2▽2f2x3f3▽f3▽2f3▽3f3x4f4▽f4▽2f4▽3f4▽4f4三、差分的性质性质1:各阶差分均可表示成函数值的线性组合mjmj例如：

性质2:向前差分和向后差分的关系为:(2.3.5)性质3:差分与差商之间有如下关系:四、等距节点的插值公式当插值节点是等距离的时候,插值公式可以用差分表示,设将式(2.3.5)代入式(2.2.2)得等距节点插值公式:(2.3.6)在上式中令x=x0+th,则式(2.3.6)变成则上式为Newton向前插值公式,其余项为:(2.3.7)也可以用向后差分公式表示Newton插值公式,令x=xn+th,x∈[x0,xn],则有上式称为Newton向后插值公式,其余项为(2.3.8‘)(2.3.8)例:已知y=f(x)=sin(x)的函数表如下X0.40.50.60.7sinx0.389420.479430.564640.64422分别用Newton向前、向后插值公式求sin0.57891的近似值。解:x0=0.4x1=0.5x2=0.6x3=0.7，写出差分表如下:xifi=sin(xi)一阶差分二阶差分三阶差分0.40.389420.50.479430.090010.60.564640.08521-0.004800.70.644220.07958-0.00563-0.00083Newton向前插值公式为上面差分表中最后一行数据就是向后差分值§2.4Hermite插值多项式一、Hermite插值多项式设x0,x1…xn为区间[a,b]上的n+1个各互异节点,且。要求一个2n+1次的插值多项式H2n+1,使得(2.4.1)构造两组插值基函数每个基函数为2n+1次多项式,且满足如下条件其中:1j=k0j≠k(2.4.2)可以证明其中j=0,1,2,…,n因此Hermite插值多项式为:特别：当n=1时，三次Hermite插值多项式在应用中非常重要，列出其详细计算公式，取插值节点xk，xk+1，三次Hermite插值多项式H3(x)满足:相应的插值基函数为:因为n=1,故j=k,k+1于是有:二、Hermite插值余项定理:设x0，x1，…，xn为[a，b]上n+1个各互异节点，H2n+1(x)为f(x)的过该组节点的2n+1次Hermite插值多项式，若f(x)∈C2n+1[a，b],f(2n+2)(x)存在,则对任意x∈[a，b]，插值余项为:特别对于三次Hermite插值余项为:例:设函数f(x)在x0=1.3,x1=1.6,x2=1.9的函数值及导数值如下:xif(xi)(xi)1.30.6200860-0.52202321.60.4554022-0.5698591.90.2818186-0.5811571应用Hermite插值求f(1.5)的近似值。解首先来计算Lagrange基本插值多项式及其导数其次,计算多项式最后得到例2:给定函数值表如下:属非标准的插值问题§2.5分段低次插值一、高次插值的病态性质:

在代数插值中,过分地提高插值的次数会带来一些新的问题,即随着插值节点的增加,插值多项式并不一定能很好的逼近,有时差异甚至很大,如函数:

当用去逼近时,就出现所谓的龙格(Runge)现象(见动画演示).

因此,高次插值多项式(7-8次以上)在实际中很少使用.二、分段线性插值设在区间[a,b]上给定n+1个节点,以及各节点的函数值。作一插值函数使满足:

我们称函数P(x)为[a,b]上关于数据(xi,yi)(i=0,1,…,n)的分段线性插值函数.

由Lagrange线性插值公式,可写出分段表达式xi0123451.000000.500000.200000.100000.058820.03846i=1,2,3,4>>x0=-5:0.1:5;>>y1=1./(1+x0.^2);>>x=-5:1:5;//坐标>>y=1./(1+x.^2);>>y2=interp1(x,y,x0);>>plot(x0,y1)>>holdon>>plot(x0,y2,'*m‘)二、Hermite分段插值多项式以上介绍的分段线性插值只能保持节点处函数连续因而光滑性较差,下面介绍的分段三次Hermite插值多项式是插值区间上的光滑函数.

设给定函数表:

在每个小区间[xi,xi+1](i=0,1,…,n－1)上作三次Hermite插值,可得到一个分段三次Hermite插值多项式H(x),使H(x)满足:§2.6三次样条插值一、三次样条插值函数的基本概念定义:如果函数S(x)在区间[a,b]上满足:

(2.6.1)和(2.6.2)共表示有4n-2个条件,因此还需要两个条件才能确定S(x),通常实际问题要求三次样条插值函数满足一定的边界条件,一般使用的边界条件有以下三类.

第一类边界条件:(2.6.3)第二类边界条件(自然边界条件):(2.6.4)第三类边界条件:

设f(x)是周期函数,不妨设以xn－x0为周期,这时边界条件满足:

这样三次样条插值函数S(x)满足(2.6.1),(2.6.2)和(2.6.3)－(2.6.5)中的某一类边界条件时,共计4n个条件,可以确定4n个待定参数.三、样条插值函数的建立讨论用节点处的二阶导数来建立插值函数,设:其中:Aj，Bj为积分常数,根据插值条件:S(xj)=yj，S(xj+1)=yj+1，得(2.6.8)

(2.6.8)式中的Mj(j=0,1,…,n)未知,为了确定Mj对S(x)求导数其中:即:(2.6.13)(2.6.14)写成矩阵形式例题:p44例7

根据已知数据作出差商表:ixifif[xi,xi+1]f[xi,xi+1,xi+2]027.74.11284.32/32294.1-0.2-0.66673303.0-1.9-0.4500解方程组得:M0=-23.531,M1=0.395M2=0.830,M3=-9.115最后由2.6.8式得三次样条插值函数的表示式(见教材)

可以证明，对于三次样条插值函数，当插值节点逐渐加密时，不但样条函数收敛于被插函数本身，其导数也同样收敛于被插函数的导数，显然其性质优于多项式插值。三次样条有明确的力学背景：样条曲线可以看成是弹性细梁受到集中载荷作用而产生的挠度曲线，在挠动不大的情况下，这种挠度曲线在数学上恰好表现为三次样条函数，集中载荷的作用点就是三次样条函数的节点。误差与收敛的讨论见教材P46(2.6.3)§2.7数据拟合的最小二乘法一、最小二乘法定义:设(xi,yi)(i=1,2,…,n)是给定的一组数据,在某个函数类H中寻求一个函数p(x),使p(x)满足：p(x)与yi的差的平方和

称这种求近似函数的方法为数据拟合的最小二乘法.用最小二乘法解决实际问题包含两个基本环节：2、按最小二乘法原则求最小二乘解。

对于函数类H，通常取次数较低的多项式或其它较简单的函数集合。1、根据所给数据点的变化趋势与问题的实际背景确定函数类H,即确定P(x)所具有的形式。二、多项式拟合

即对给定的数据(xi,yi)(i=1,2,…,n),求一个m次多项式(m<n)Pm(x)=a0+a1x+…+amxm使得即:

上式称为正则方程组.可以写成如下矩阵形式:

通过解正则方程组可求出参数a0,a1,…,am(可以证明方程组存在惟一解),进而得到已给数据组(xi,yi)(i=1,2,…,n)的最小二乘m次拟合多项式.三、例题例1已知数据如表1所示,求它的拟合函数编号123456789xi1345678910yi1054211234表1解：根据已知数据画出散点图如图1所示。图1

由图1可以看出它近似一条抛物线,故可设拟合曲线为建立正则方程组,计算下面各个值编号xiyixiyixi2xi2yixi3xi41110101101123515945278134416166464256452102550125625561636362161296671749493432401782166412851240968932781243729656191044010040010001000053321473811025301725317写出a0,a1,a2的正规方程组:例2由下表的数据求一个形如1234567815.320.527.436