工程实践中必不可少的数学方法(数据处理） 第五章 数值分析法 工程实践中必不可少的数学方法(数据处理） 用连续的观点处理离散问题 理论与经验的有机结合

5.1曲线拟合法 经验公式 含 义 根据一组实验观测数据确定自 变量x与因变量y的一个最“逼近” 的函数关系式

几何解释 找一条“最佳” 的线，使 之与 靠得最近 作 用 用连续函数分析方 法进行建模讨论。

一个实际问题 程序控制的铣床精密工件加工工件的表面是 光滑的，走刀方向是逐段线性的，实际上是利 用逐段线性函数近似连续光滑函数。

曲 线 拟 合 步 骤 确定经验公式形式 确定经验公式中的系数 检验经验公式有效性

利用已知的结论 描点作图法 如何确定经验公式形式 多项式近似 曲改直

利用已知结论 相关的定理、定律 前人比较成熟的成果，公认的结论 普遍采用的公式 根据经验的假设，假想（要验证）

多项式近似是工程中十分常见的方法， 它首先需要我们确定多项式的次数， 一般可以用差分法、差商法来估计。 （逼近论 类似的思想：泰勒公式）

差分与差商概念 一阶向前差分 二阶向前差分 ………… m阶向前差分 一阶差商 二阶差商 m阶差商

一般地，等距节点用差分，不等距节点用差商 作用： 作为微分与导数的近似估计，便于确 定多项式的阶数

差商与导数的联系 （微分中值定理） 若y=f(x)在[a，b]上m次可导，且 则 若结点为等距分割点时，有 ，h为结点距，且 因此对n阶多项式有 常数。据此,我们可以根据数据的差分来确定多项式的次数。

问题 氨蒸汽的压力和温度关心 用低温状态下的压力等有关数据进行外 推。表5.1给出了氨蒸汽的一组温度和压 力数据。 考虑到测量设备等的限制，我们希望利 用低温状态下的压力等有关数据进行外 推。表5.1给出了氨蒸汽的一组温度和压 力数据。

目的 希望利用低温状态下的压力等有关数据进行 外推。能否从所列的数据中计算75度 氨蒸汽的 压力？

表5.1 温度( ) 20 25 30 35 40 45 50 55 60 压力(kN/m2) 805 985 1170 1365 1570 1790 2030 2300 2610 根据5.1的数据,可以绘制图5.2。

根据表5.1的数据可以得到差分表5.2 3阶或4阶多项式比较合适 p 20 805 25 985 180 30 1170 185 5 20 805 25 985 180 30 1170 185 5 35 1365 195 10 5 40 1570 205 10 0 -5 45 1790 220 15 5 5 50 2030 240 20 5 0 55 2300 270 30 10 5 60 2610 310 40 10 0 3阶或4阶多项式比较合适

曲线改直 是工程中又一常用的判断曲线形式的方法,许多常见的函数都可以通过适当的变换转化为线性函数。 （1）幂函数 （2）指数函数 （3）抛物函数

电弧电流I与电压降V之间的经验公式的确定。 试确定V与I的经验公式的形式。

>> x=[0.5 1 2 4 8 12 ]; >> y=[160 120 94 75 62 56]; >> x1=log(x);y1=log(y-30); >> y2=log(y-50); >> subplot(2,2,1); plot(x,y,'o'); title('I-V图') subplot(2,2,2); plot(x,y1,'.'); title('I-log(V-30)图') subplot(2,2,3); plot(x,y2,'.'); y2=log(y-50); title('I-log(V-50)图') subplot(2,2,4); plot(x1,y1,'.'); title('log(I)-log(V-30)图') >>

下图就是先通过全对数图确定经验公式形式， 再借助其他方法获得的经验公式图形与实际 数据比较结果。 利用原始数据作图，或作全对数图，办对数 图，可以帮助我们进一步寻找经验公式的形式， 下图就是先通过全对数图确定经验公式形式， 再借助其他方法获得的经验公式图形与实际 数据比较结果。

经验函数形式

I 0.5 1 2 4 8 12 LnI -0.301 0 0.301 0.602 0.903 1.079 V 160 120 94 75 62 56 ln(V-30) 2.114 1.954 1.806 1.653 1.505 1.415 lnI与ln(V-30)呈线性关系

x=[0.5 1 2 4 8 12 ]; y=[160 120 94 75 62 56]; x1=log(x);y1=log(y-30); plot(x1,y1,'ro',x1,y1); hold on polyfit(x1,y1,1) text(0.5,4.5,'log(V-30)=-0.5.4log(I)+4.509') hold off

x=[0.5 1 2 4 8 12 ]; y=[160 120 94 75 62 56]; x1=log(x);y1=log(y-50); z1=-0.8847*x1+4.249; plot(x1,y1,'ro',x1,z1); hold on text(0.5,4.,'log(V-50)=-0.8847log(I)+4.249') hold off

图示法 均值法 确定系数常用方法 差分法 最小二乘法 插值法

图示法、均值法主要用于线性关系，且精度较低； 差分可以用于低阶多项式，精度较低； 最小二乘法与插值法在处理实际问题时较常用。

最小二乘法 基 本 思 想

几个最常见的形式

设因变量 与自变量 有如下关系: 是 n 次观测值

记 A 称为观测矩阵

最小二乘法原理 最小二乘法原理就是找一组数值 ，使 为最小。 数学原理 高等数学中介绍的多元函数求极值 S取极值的必要条件

最小二乘法应用 例如用最小二乘法求出的例1的经验公式为 或者

x=[20 25 30 35 40 45 50 55 60]; >> y=[805 985 1170 1365 1570 1790 2030 2300 2610]; >> y1=0.007138,-0.511111,48.887205,-27.121212; >> y2=0.000130536,0.01375,0.6844*,19.76,223.79; >> plot(x,y,'o',x,y1,'r-',x,y2,'b-.') >> hold on >> legend('原始数据','3阶多项式','4阶多项式 ') >> hold off

5.2 插值法建模 目的 无论是从理论和计算的角度,还是从应用的角度看,多项式都是最简单的函数,因此,多项式插值是最基本的插值方法,不妨设 已知n+1个节点(xi,yi),i=0,1,2,…,n,其中,xi互不相同, 求一个过这些已知点的函数曲线,据此推断出任一点x*处的值。 无论是从理论和计算的角度,还是从应用的角度看,多项式都是最简单的函数,因此,多项式插值是最基本的插值方法,不妨设 Ln(x)是n次多项式,记作

对于节点(xi,yi) 应有 将(6.4)式代入(6.3)得

易见 是范得蒙行列式, , 拉格朗日插值多项式 在求插值多项式时,通常的做法不是解方程组(6.5),而是先构造一组基函数: （6.6）

(6.7) 令 (6.8)

（5.9） 如果 则 图 5.4 与 的比较

利用拉格朗日插值法需要注意的是，随着节点数的增加,拉格朗日插值多项式的次数也增加,高次插值多项式有时会出现很大的振荡行为。 理论上,不能保证 处处收敛到g(x)。 请看下面的一个例子 。图5.4给出了 与 的图形比较。从图5.4中可以看出，高次插值多项式并没有获得理想的近似结果。高次插值多项式的这些缺陷，促使人们寻求简单的低次多项式插值。

分段线性插值 简单地说,分段线性插值就是将每两个相邻的节点用直线连接起来,如此形成的一条折线就是分段线性插值函数,记作 ,它满足 ,且 在每个小区间 上是线性函数。 如果已知y=g(x)在节点 的值 ，则 L(x)可以表示为 (5.10)

其中 (5.11) 易见Ln(x)满足插值条件。 记 利用（5.9）式可以得到 (5.12) (5.12)式给出了分段线性插值的误差估计。

分段二次插值 分段二次插值的基本思想是利用三个相邻节点的数据分别作二次插值,形成二次曲线段,然后将这些曲线段连接起来,形成的曲线,就是分段二次插值曲线。 如果给定g(x)在节点 的值 对应n的奇偶情况,分别讨论如下: （1）设n为偶数 对于数据 ,作二次插值函数 (5.13)

（2）当n为奇数时,在小区间 上作二次插 值函数 ，其它的区间上作法与n为偶数相同。 令 称 为 在这n个节点的分段二次插值函数。 可以证明, (5.14) 其中, 。

分段二次插值函数有很好的收敛性,计算也 比较简单,还可以根据函数g(x)的具体情况 在不同的小区间上采用不同的插值公式。 分段插值函数虽然在节点处连续,但是不一 定可导,因此光滑性较差。下面介绍的三次 样条插值将克服这一缺点。

样条插值函数 它的基本思想就是对节点分段插值，但是需 要将这些分段函数曲线光滑地连接起来，形 成一条光滑曲线。

最常用的是三次样条插值。 定义1 设在区间[a，b]上给定n+1个插值节点 及其函数g（x）相应的值 。 如果函数 满足： (1) ; (2) 在每个小区间 上是三次多项式; (3) 在[a,b]上有连续的二阶导数。 则称为[a，b]区间上的三次样条插值函数。

其中系数待定,且要满足下列条件: （i）插值条件 （ii）连接条件 插值条件与连接条件共给出了4n-2个等式,而要求出S(x)需要确定4n个参数,因此还需要附加两个条件,通常称为边界条件,常用的有如下两种:

（iii） 边界条件 (5.15) 或者 (5.16) 利用插值条件、连接条件、边界条件可以唯一确定S(x)。直接计算S(x)一般计算量很大，下面我们不加证明地给出一种简单的构造三次样条插值函数的方法。 设在区间[a，b]上给定n+1个插值节点 及其函数 相应的值 。 记

(5.17) (5.18) 式中 是未知的,可通过下面方法求出来。 对于由(5.15)式给出的边界条件 取 由下列矩阵方程求出

(5.19) 对于由(5.16)式给出的边界条件，取 其它的未知参数由下列方程求出 (5.20)

以上两种形式方程组的系数矩阵都为严格对 角占优的三对角方程组,存在唯一解，求 出 后，代入(5.18)式,就可以得到S(x)的具体表达式。

x=0:10; y=[2,4,5,7,6,7,5,8,9,12,17]; x0=0:0.1:10; y1=interp1(x,y,x0); y2=interp1(x,y,x0,'spline'); plot(x,y,'ro',x0,y1,'-.',x0,y2) legend('原始数据','分段线性','3次样条')

实例 美国某州的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天所用的总水量。许多社区没有测量流入或流出当地水塔的水量的装置，他们只能代之以每小时厕量水塔中的水位，其精度不超过0.5% 。更重要的是，当水塔中的水位下降到最低水位L时，水泵就启动向水塔输水直到最高水位H才停止，但也不能测量水泵的供水量。当水泵正在输水时不容易建立水塔水位和水泵工作时间的关系。水泵每天输水一次或两次

当水泵正在输水时不容易建立水塔水位和水泵工作时间的关系。水泵每天输水一次或两次每次大约两小时。 试估计任何时刻（包括水泵正在输水时间）从水塔流出的流量，并估计一天的总用水量。表5.5给出了某个小镇一天中真实的数据。

背景信息 经了解,该水塔是一个高为40英尺,直径为57英尺的正圆柱。通常当水塔水位降至约为27英尺时水泵开始启动，当水位上升至35.5英尺时水泵停止工作。 (注:1英尺=0.305米,1加仑=4.546升=4.546×10-3 立方米)

表5.5 某小镇某天水塔水位统计 时间(秒) 水位(0.01英尺) 时间(秒) 水位(0.01英尺) 0 3175 46636 3350 3361 3110 49953 3260 6635 3054 53936 3167 10619 2994 57254 3087 13937 2947 60574 3012 17921 2892 64554 2927 21240 2850 68535 2842 25223 2795 71854 2767 时间(秒) 水位(0.01英尺) 时间(秒) 水位(0.01英尺) 28543 2752 75021 2697 32284 2697 79254 水泵开动 35932 水泵开动 82649 水泵开动 39332 水泵开动 85968 3475 39435 3550 89953 3397

模 型 假 设 根据流体力学中的Torricell定律,从水塔中流出水的最大流速正比于水位高度的平方根。 假设水泵供水速度是常数 从水塔流出的水流速度比水泵供水速度小， 不存在断水情况 每天从水塔中流出的水流速度随时间的变化可以用光滑曲线近似

当水泵不供水时，根据水位随时间的变化可以确定从水塔中流出水流的速度 。 建立数学模型 当水泵不供水时，根据水位随时间的变化可以确定从水塔中流出水流的速度 。 (5.21)

时间中点 平均流速 时间中点 平均流速 时间中点 平均流速 0 无数据 9.47 无数据 18.49 12.2 利用差商确定 ， 表示函数 关于 的函数 表5.6 观测时间中点与该时间区间平均水流速度 时间中点 平均流速 时间中点 平均流速 时间中点 平均流速 (小时) (×103加仑/小时) (小时) (×103加仑/小时) (小时) (×103加仑/小时) 0 无数据 9.47 无数据 18.49 12.2 0.46 11.2 10.45 无数据 19.50 12.9 1.38 9.7 10.94 无数据 20.40 12.6 2.4 8.6 11.49 15.6 21.43 无数据 3.41 8.1 12.49 16.4 22.49 无数据 4.43 9.3 13.42 15.5 23.42 无数据 5.44 7.2 14.43 13.4 24.43 11.2 6.45 7.9 15.44 13.8 25.32 3.5 7.47 7.4 16.37 12.9 8.45 8.4 17.38 12.2

(2) 用三次样条曲线拟合流速曲线 利用流速的光滑性假设,采用光滑的三次样条插值近似表示流速曲线,得到如图5.5的插值曲线。 (3) 计算总用水量 通过三次样条插值,可以获得水塔的流速V与时间的关系式V(t),据此, 我们可以求出这个小 镇一天的总的用水量 具体算法此处不再详细说明, 留待读者自己考虑。 图5.5 三次样条插值曲线

机床加工问题 用程控铣床加工机翼断面的下轮廓线时 每一刀只能沿x方向和y方向走非常小的一步。 表3-1给出了下轮廓线上的部分数据 试完成加工所需的数据，画出曲线.

用MATLAB作一维插值计算 yi=interp1(x，y，xi，'method') xi处的插值结果 插值节点 被插值点 插值方法 ‘nearest’ ：最邻近插值‘linear’ ： 线性插值； ‘spline’ ： 三次样条插值； ‘cubic’ ： 立方插值。 缺省时： 线性插值。 插值方法要求x是单调的，且xi不能够超过x的范围。

求解机床加工问题 x0=[0 3 5 7 9 11 12 13 14 15 ]; y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ]; x=0:0.1:15; y=interp1(x0,y0,x,'spline'); plot(x0,y0,'k+',x,y,'r') grid on x0=[0 3 5 7 9 11 12 13 14 15 ]; y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ]; x=0:0.1:15; y=interp1(x0,y0,x,'spline'); plot(x0,y0,'k+',x,y,'r') grid on

用MATLAB作网格节点数据的插值(二维) z=interp2(x0,y0,z0,x,y,’method’) 被插值点的函数值 插值节点 被插值点 插值方法 ‘nearest’ 最邻近插值 ‘linear’ 双线性插值 ‘cubic’ 双三次插值 缺省时, 双线性插值 要求x0,y0单调；x，y可取为矩阵，或x取行向量，y取为列向量，x,y的值分别不能超出x0,y0的范围。

用MATLAB作散点数据的插值计算 插值节点 被插值点的函数值 被插值点 插值方法 要求cx取行向量，cy取为列向量。 cz =griddata（x，y，z，cx，cy，‘method’） 被插值点的函数值 被插值点 插值方法 插值节点 ‘nearest’ 最邻近插值 ‘linear’ 双线性插值 ‘cubic’ 双三次插值 'v4‘ Matlab提供的插值方法 缺省时, 双线性插值 要求cx取行向量，cy取为列向量。

航行区域的警示线 某海域上频繁地有各种吨位的船只经过。 为保证船只的航行安全，有关机构在低潮时对水深进行了测量,下表是测量数据： 表3 水道水深的测量数据 x 129.0 140.0 103.5 88.0 185.5 195.0 105.5 y 7.5 141.5 23.0 147.0 22.5 137.5 85.5 z 4 8 6 8 6 8 8 x 157.5 107.5 77.0 81.0 162.0 162.0 117.5 y -6.5 -81.0 3.0 56.5 -66.5 84.0 -33.5 z 9 9 8 8 9 4 9

航行区域的警示线 其中（x, y）为测量点，z为（x, y）处的水深 (英尺)，水深z是区域坐标（x, y）的函数z= z (x, y）， 船的吨位可以用其吃水深度来反映，分为 4英 尺、4.5英尺、5英尺和 5.5英尺 4 档。 航运部门要在矩形海域（75，200）×（－50，150）上为不同吨位的航船设置警示标记。 请根据测量的数据描述该海域的地貌，并绘制不同吨位的警示线，供航运部门使用。

Matlab求解 x=[129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5]; y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5]; z=[-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9]; cx=75:0.5:200; cy=-70:0.5:150; cz=griddata(x,y,z,cx,cy','cubic'); meshz(cx,cy,cz), xlabel('X'),ylabel('Y'),zlabel('Z') figure(2),contour(cx,cy,cz,[-5 -5]); grid on, hold on plot(x,y,'+') xlabel('X'),ylabel('Y') x=[129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5]; y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5]; z=[-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9]; cx=75:0.5:200; cy=-70:0.5:150; cz=griddata(x,y,z,cx,cy','cubic'); meshz(cx,cy,cz),rotate3d xlabel('X'),ylabel('Y'),zlabel('Z') %pause figure(2),contour(cx,cy,cz,[-5 -5]);grid on, hold on plot(x,y,'+') xlabel('X'),ylabel('Y')

自变量不单调 x=[0 30 50 70 80 90 120 148 170 180 202 212 230 248 268 271 280 290 300 312 320 340 360 372 382 390 416 430 478 440 420 380 360 340 320 314 280 240 200]; y=[80 64 47 42 48 66 80 120 121 138 160 182 200 208 212 210 200 196 188 186 200 184 188 200 202 240 246 280 296 308 334 328 334 346 356 360 392 390 400];