数值分析方法(讲义).doc

第十章 数值分析方法在生产实际中，常常要处理由实验或测量所得到的一批离散数据，数值分析中的插值与拟合方法就是要通过这些数据去确定某一类已经函数的参数，或寻求某个近似函数使之与已知数据有较高的拟合精度。插值与拟合的方法很多，这里主要介绍线性插值方法、多项式插值方法和样条插值方法，以及最小二乘拟合方法在实际问题中的应用。相应的理论和算法是数值分析的内容，这里不作详细介绍。1 数据插值方法及应用 在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样点，要求由此建立变量之间的函数关系或得到样点之外的数据。与此有关的一类问题是当原始数据精度较高，要求确定一个初等函数（一般用多项式或分段多项式函数）通过已知各数据点（节点），即，或要求得函数在另外一些点（插值点）处的数值，这便是插值问题。1、分段线性插值 这是最通俗的一种方法，直观上就是将各数据点用折线连接起来。如果那么分段线性插值公式为可以证明，当分点足够细时，分段线性插值是收敛的。其缺点是不能形成一条光滑曲线。例1、已知欧洲一个国家的地图，为了算出它的国土面积，对地图作了如下测量：以由西向东方向为x轴，由南向北方向为y轴，选择方便的原点，并将从最西边界点到最东边界点在x轴上的区间适当的分为若干段，在每个分点的y方向测出南边界点和北边界点的y坐标y1和y2，这样就得到下表的数据（单位：mm）。x7.010.513.017.534.040.544.548.056.0y1444547505038303034y24459707293100110110110x61.068.576.580.591.096.0101.0104.0106.5y1363441454643373328y2117118116118118121124121121x111.5118.0123.5136.5142.0146.0150.0157.0158.0y1326555545250666668y2121122116838182868568根据地图的比例，18 mm相当于40 km。根据测量数据，利用MATLAB软件对上下边界进行线性多项式插值，分别求出上边界函数，下边界函数，利用求平面图形面积的数值积分方法将该面积近似分成若干个小长方形，分别求出这些长方形的面积后相加即为该面积的近似解。式中，。这里线性插值和面积计算源程序如下：clear allx=7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0;y1=44 45 47 50 50 38 30 30 34 36 34 41 45 46 43 37 33 28 32 65 55 54 52 50 66 66 68;y2=44 59 70 72 93 100 110 110 110 117 118 116 118 118 121 124 121 121 121 122 116 83 81 82 86 85 68;newx=7:0.1:158;newy1=interp1(x,y1,newx,linear);newy2=interp1(x,y2,newx,linear);Area=sum(newy2- newy1)*0.1/182*1600最后计算的面积约为42414平方公里。 2、多项式插值 设有次多项式通过所有个点，那么就有可以证明当且时，这样的多项式存在且唯一。若要求得到函数表达式，可直接解上面方程组。若只要求得函数在插值点处数值，可用下列Lagrange插值公式多项式插值光滑但不具有收敛性，一般不宜采用高次多项式（如）插值。例2、在万能拉拨机中有一个园柱形凸轮，其底园半径R=300mm，凸轮的上端面不在同一平面上，而要根据动杆位移变化的需要进行设计制造。按设计要求，将底园周18等分，旋转一周。第个分点对应柱高，数据见下表。为了数控加工，需要计算出园周上任一点的柱高。凸轮高度的数据（单位：mm）分点0和1812345柱高502.8525.0514.3451.0326.5188.6分点67891011柱高92.259.662.2102.7147.1191.6分点121314151617柱高236.0280.5324.9369.4413.8458.3我们将园周展开，借助MATLAB软件画出对应的柱高曲线散点图(左下图)。clear;close;x=linspace(0,2*pi*300,19);y=502.8 ,525.0,514.3,451.0,326.5,188.6,92.2,59.6,62.2,102.7,147.1,191.6,236.0,280.5,324.9,369.4,413.8,458.3,502.8;plot(x,y,o);axis(0,2000,0,550); 可见，可以用三次多项式插值，下面给出借助MATLAB软件画出的柱高插值曲线图(右上图)。xi=0:2*pi*300;yi=interp1(x,y,xi,cubic);plot(xi,yi);3、样条插值这是最常用的插值方法。数学上所说的样条，实质上是指分段多项式的光滑连接。设有称分段函数为k次样条函数,若它满足(1) 在每个小区间上是次数不超过次的多项式；(2) 在上具有直到阶的连续导数。用样条函数作出的插值称为样条插值。工程上广泛采用三次样条插值。例3、某居民区的自来水是由一个园柱形的水塔提供。水塔高12.2米，直径17.4米。水塔由水泵根据塔中水位高低自动加水，一般每天水泵工作两次。按照设计，当水塔内的水位降至约8.2米时，水泵自动启动加水；当水位升至约10.8米时，水泵停止工作。现在需要了解该居民区用水规律，这可以通过用水率（单位时间的用水量）来反映。通过间隔一段时间测量水塔中的水位来估算用水率。下表是某一天的测量记录数据，测量了28个时刻（单位：小时）的水位（单位：米），但由于其中有3个时刻正遇到水泵在向水塔供水，而无水位记录（表中用符号/表示）。时刻00.9211.8432.9493.8714.9785.900水位9.6779.4799.3089.1258.9828.8148.686时刻7.0067.9828.9679.98110.92510.95412.032水位8.5258.3888.220/10.82010.500时刻12.95413.87514.98215.90316.82617.93119.037水位10.2109.9369.6539.4099.1808.9218.662时刻19.95920.83922.01522.95823.88024.98625.908水位8.4338.220/10.82010.59110.35410.180先通过体积公式，利用上表中的水位高，得到不同时刻水塔中水的体积。为提高精度，采用二阶差商来估算时刻的水流速度，即。 具体地，因为所有数据被水泵两次工作分割成三组数据，对每组数据的中间数据采用中心差商，前后两个数据不能够采用中心差商，改用向前或向后差商。中心差商公式 向前差商公式 向后差商公式 估算出水塔中水的流速（单位：立方米/小时）见下表。时刻00.9211.8432.9493.8714.9785.900流速54.51642.32038.08541.67933.29737.81430.748时刻7.0067.9828.9679.98110.92510.95412.032流速38.45532.12241.718/73.68676.434时刻12.95413.87514.98215.90316.82617.93119.037流速71.68660.19068.33359.21752.01156.62663.023时刻19.95920.83922.01522.95823.88024.98625.908流速54.85955.439/57.60257.76651.89136.464先用MATLAB画出水流速散点图。t=0 0.921 1.843 2.949 3.871 4.978 5.9 7.006 7.982 8.967 10.954 12.032 12.954 13.875 14.982 15.903 16.826 17.931 19.037 19.959 20.839 22.958 23.88 24.986 25.908; r=54.516 42.320 38.085 41.679 33.297 37.814 30.748 38.455 32.122 41.718 73.686 76.434 71.686 60.19 68.333 59.217 52.011 56.626 63.023 54.859 55.439 57.602 57.766 51.891 36.464;plot(t,r,b+); % (t,r)表示时间和流速title(流速散点图)；xlabel（时间(小时)）; ylabel(流速(立方米/小时)使用MATLAB软件中的三次样条插值命令得到用水率函数如下图所示。x0=t;y0=r;l,n=size (x0); dl=x0(n)-x0(1);x=x0(1):1/3600:x0(n); %被插值点ys=interp1 (x0,y0,x,spline); %样条插值输出plot (x,ys);title(样条插值下的流速图)；xlabel（时间(小时)）; ylabel(流速(立方米/小时) 2 数据拟合方法及应用 在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样点，要求由此建立变量之间的函数关系或得到样点之外的数据。与此有关的另一类问题是拟合问题。当原始数据有误差时，我们确定的初等函数并不要求经过数据点，而是要求在某种距离意义下的误差达到最小（通常考虑使各数据点误差平方和最小）。 假设已知函数（这里可以是多个未知参数）的一批有误差的数据要求据此确定参数，这样的问题称为数据拟合。最小二乘法就是求使得残差平方和达到最小。这里的建模原理实质上就是数理统计中的回归分析。 1、线性函数（1）建立回归方程若离散样点集中在一条直线的附近，这时可建立线性回归方程按最小二乘法得到的具体算法是记 则有（2）线性回归的显著性检验采用如下检验统计量： ，其中的若，则认为所建立的线性回归方程正确。（3）利用回归方程作区间预测当时，的置信度为的预测区间是：例1、为研究某一化学反应过程中温度x对产品得率y的影响，测得数据如下。x ( )100 110 120 130 140 150 160 170 180 190y ( % ) 45 51 54 61 66 70 74 78 85 89要求拟合出它们的函数关系。x=100 110 120 130 140 150 160 170 180 190;y=45 51 54 61 66 70 74 78 85 89;close; plot(x,y,o);由于这些点落在一条直线附近，可以用线性函数来拟合，按上面算法通过命令p=polyfit(x,y,1)求出。2、可线性化的函数 根据专业知识或离散样点的形状，有时可选择适当的非线性函数来拟合。为确定其中的未知参数，可通过变量转换，把非线性函数转换成线性函数，然后借助线性函数的方法来实现数据拟合。诸如下面所列的曲线函数都能做到线性化：（1） （2） （3）（4） （5） （6）以为例，我们只要作变换，就可化为线性函数。例2、电容器充电后，电压达到100伏，然后开始放电，测得时刻（秒）时电压（伏）如下表。t 0 1 2 3 4 5 6 7 8 9 10u100 75 55 40 30 20 15 10 10 5 5要求拟合电压u与放电时间t的函数关系。t=0 1 2 3 4 5 6 7 8 9 10;u=100 75 55 40 30 20 15 10 10 5 5;close; plot(t,u,o) 由于这些点落在曲线附近，可通过变量转换，化成线性函数。按上面算法通过如下命令来实现：T=t; U=log(u);p=polyfit(T, U,1)b=p(1)a=exp(p(2)求出，这样所拟合的函数3、多项式函数 若离散样点的形状表明既不能用线性函数来拟合，又不能用可以线性化的函数来拟合的话，从理论上讲，用一个多项式函数来拟合总是可行的。在实际应用中，最常用的是二次和三次多项式函数。下面通过一个例子来说明。例3、某种产品在生产过程中的废品率y与它所含的某种物质量x有关，现将试验所得16组数据记录列于下表。x34 36 37 38 39 39 39 40y 1.30 1.00 0.73 0.90 0.81 0.70 0.60 0.50x 40 41 42 43 43 45 47 48y 0.44 0.56 0.30 0.42 0.35 0.40 0.41 0.60要求拟合y与x的函数关系。x=34 36 37 38 39 39 39 40 40 41 42 43 43 45 47 48;y=1.30 1.00 0.73 0.90 0.81 0.70 0.60 0.50 0.44 0.56