廓模型与曲线拟合 一、.pdf

《廓模型与曲线拟合 一、.pdf》由会员分享，可在线阅读，更多相关《廓模型与曲线拟合 一、.pdf（13页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第二讲第二讲 量纲分析、轮廓模型与曲线拟合量纲分析、轮廓模型与曲线拟合 一、量纲分析一、量纲分析 1 单位与量纲单位与量纲 物理量可以分为两类：基本物理量（它们相互独立并可以通过自然规律的各种定律构成其他的物理量）长度，质量，时间；衍生物理量（由基本物理量和定律导出的物理量）速度，加速度等。国际通用的单位制（SI 制）由七个基本单位构成：长度L，质量M，时间T，电流强度I，温度，光强和物质的量：JN物理量 单位 符号 长度 米 m 质量 千克 kg 时间 秒 s 电流强度 安培 A 温度 开尔文 K 光强 坎德拉 cd 物质的量 摩尔 mol 一个物理量Q一般可以表示为基本物理量幂次之积，则称

2、这个乘幂之积的表达式为 该 物 理 量 对 选 定 的 一 组 基 本 量 的 量 纲 表 达 式，NJITMLQ=,称为量纲指数。例：速度=1LT，能量=22TML，压强=21TML。若某个物理量的量纲指数全为 0，则称该物理量为无量纲量，如圆周率，角度。注意：常数未必都无量纲（如万有引力常数），无量纲的量未必都是常数（如角度）。22 量纲分析法（Buckingham 定理，量纲齐次法则）量纲分析法（Buckingham 定理，量纲齐次法则）量纲齐次法则量纲齐次法则 设个物理量之间存在函数关系nnxxx,21L0),(21=nxxxL。有个 基 本 量 纲。的 量 纲 可 以 表 示 为。若

3、矩阵m,21mXXXLixniXxijmjji,2,1,1L=mnijB=)(的秩为r，则存在rn个无量纲量rn,21L使 得0),(1=rnL与0),(21xxL=nx相 互 等 价，其 中，而是方程组的基本解（称=miissix1si0=TBTB为量纲矩阵）。按照上述定理，量纲分析方法的一般步骤如下：（1）将与问题有关的有量纲的物理量记作，按照物理意义确定问题的基本量纲；nxxx,21L,21mXXXL（2）按照物理知识用基本量纲表示的量纲为。ixniXxijmjji,2,1,1L=（3）设为 无 量 纲 量，其 中=niiix1i待 定。的 量 纲 表 达 式 为，由 于=niiijii

4、jmjjnimjjXX1111)(无 量 纲，可 以 得 到。mjniiij,2,1,01L=（4）解线性方程组。若方程组的秩为mjniiij,2,1,01L=r，则有rn个基本解，记做。从而为无量纲量。rnsTsnss=,2,1,),(1LLrnsxniissi=,2,1,1L 例例 1 单摆的运动周期单摆的运动周期 考虑质量为的小球系在长度为l的线的一端，偏离平衡位置后小球在重力的作用下做往复摆动，忽略阻力，球摆动周期t的表达式。mmg在 该 问 题 中 出 现 了 四 个 物 理 量，量 纲 分 别 为。设 四 个 量 之 间 存 在 关 系，为无量纲量。，则有 gmlt,2,=LTgM

5、mLlTt0),(=gmltfvuyxgmlt=vyuvxvuyxLMTLTMLT+=22)(=+=0002vyuvx 易知该线性方程组的基础解系中有唯一的解向量，因此可以确定唯一的无量纲量，由上述定理T)1,0,1,2(glt12=0),(=gmltf与等价，按照隐函数存在定理，可以确定出函数，即0)()(12=glt0)(12=gltCglt=12glCt=，其中C待定。例例 2 船在水中航行受到的阻力船在水中航行受到的阻力 假设：（1）船长为l，吃水深度为h，船的航行速度为；v（2）不考虑风对船航行的影响（3）水是均匀的，密度一致，粘性系数一致 模型建立：（1）此问题中涉及的物理量有：阻

6、力，船长，吃水深度，速度，水的密度flhv，水的粘性系数，重力加速度。g（2）各物理量的量纲 112312,=MTLLTgMLLTvLhLlMLTf（注：对的量纲的说明，由关系xvp=）量纲矩阵为=210100201100011131111A。（3）假设gvhlf,之间存在关系0),(=gvhlf，按照量纲齐次原则，构造为无量纲量，为方程组的解。7654321yyyyyyygvhlf=621,yyyL0=Ay（4）求解方程组，注意到0=Ay3)(=Arank，则方程组有 4 个基本解。TTTTyyyy)1,0,0,2,0,1,0()0,1,1,1,0,1,0()0,0,1,2,0,2,1()0

7、,0,0,0,1,1,0(4321=（5）四个基本解对应的无量纲量分别为，。因此11=lhglv22=13=lv1224=vfl0),(=gvhlf等价于，由隐函数存在定理，存在0),(122121=vfllvglvlh使得),(12122=lvglvlhvlf。（注：在流体力学中称lgvFr=为 Froude 数，lv=Re为 Reynold 数。）3 无量纲化无量纲化 单位和量纲在建模过程中是一个需要注意的问题，在建立模型时，为了满足量纲齐次原则需要引入新的参量，这使得模型十分复杂；在建立和分析模型时，模型所描述的实际问题的内涵性质一般应该独立于度量单位的选择。因此在机理模型建立过程中如何

8、使得模型摆脱度量单位的影响，以得到更一般反映普遍规律的模型，下面例子将说明无量纲化方法在建模中的作用。例例 3 火箭发射模型火箭发射模型 模型假设：（1）火箭在地球的表面发射，忽略空气阻力；（2）火箭的初速度为，地球的半径为vr，地球表面的引力加速度为，火箭和地球的质量分别为；g21,mm（3）火箭在t时刻的高度为，)(ty0)0(=y。模型建立：由万有引力定律221221)(rymmkdtydm+=，vyy=)0(,0)0(。由假设 2，。在方程始终令，则有gy=)0(0=t22rmkg=，则模型可以简化为 vyyrygry=+=)0(,0)0(,)(22。在模型中有三个参量vgr,，两个变

9、量。这些量都是有量纲的，下面将利用无量纲化的方法化简模型，将变量化为无量纲的。所谓无量纲化就是将模型中的变量化为无量纲的变量。在该模型中，因此可以选择两个分别具有时间和长度量纲的常数，通过变换yt,LyTt=,ccyt,ccyyyttt/,/=生成两个无量纲变量yt,，在该变换下 22222)(,dtydytt dyddtdyytt dydcccc=，从而模型用新的无量纲变量表示为 ccccccccccyvtyyyryygtryygrytdtydytt dyd=+=+=)0(,0)0(,)1(1)()(222222222222（1）取，即vrtrycc/,=ryyrvtt/,/=，通过简单的求

10、导数运算，可以得到微分方程2)1(1+=yy，其中rgv2=为无量纲量，初始条件为1)0(,0)0(=yy。（2）令vrtrycc/,=，则有=+=)0(,0)0(,)1(12yyyy；（3）令，则有gvtgvycc/,/2=1)0(,0)0(,)1(12=+=yyyy。利 用 这 三 种 变 换 都 可 以 将 参 数 降 为 1 个，进 一 步 可 以 讨 论 对的 处 理（vsmrg/8000，从而1）。问题：在三个模型中能否作为一个非常小的数而视为 0？为什么？（从将0=简化后的模型进行合理性分析。）第一个变换中视为 0，得到0)1(12=+y，无意义；第二个变换中视为 0，得到0)0

11、(,0)0(,)1(12=+=yyyy，由01)0(=y和Taylor 展开式知道，在 0 点的附近，05.0)0(2100)(22=+ttyty，无意义；第三个变换中视为 0，得到1)0(,0)0(,1=yyy，方程有解tty+=22，代会原来的变量可以得到vtgty+=221，可以看成是原问题模型的近似解。二、轮廓模型二、轮廓模型 1量的比例关系 1量的比例关系 因为模型表达了不同量纲的量之间的转换规律，不同量纲的量的乘幂之间一定存在比例关系。所以在同一模型中，若量和的量纲分别为 和，则一定有。1X2XXX=1XX=2/21kXX=例 4（几何上的比例关系）例 4（几何上的比例关系）对于正

12、立方体：设棱长为，底面周长为aL=1aL42=，底面对角线长aL23=，立体对角线长aL34=；表面积，底面面积，对角面面积 216aS=22aS=232aS=；体积。31aV=结论：在简单的几何体中，相应部位的面积与相应部位长度的平方呈正比：，相应部位的体积与相应部位长度的立方呈正比：，相应部位的体积与相应部位面积的 3/2 次方呈正比：。2jiLS 3jiLV 2/3jiSV 对于一般的立体上述比例关系同样成立。例 5 生活中的长度、面积和体积例 5 生活中的长度、面积和体积（1）黑鲈的体重W和体长L W(ozs)17 16 17 23 26 27 41 49 L(in)12.50 12.

13、63 12.63 14.13 14.50 14.50 17.25 17.75 L3 1953 2015 2015 2821 3049 3049 5133 5592 W/L3 .0087 .0079 .0084 .008 .0085 .0089 .008 .0088 黑鲈鱼的体重与体长关系（）3LW（2）人的体重 W 和身高 L W(kg)12 17 22 35 48 54 66 75 L(cm)86 108 116 135 155 167 178 185 L3(103cm3)636 1260 1560 2460 3724 4657 5640 6332 W/L3 .0189 .0135 .014

14、1 .0142 .0129 .0116 .0117 .0118 人的体重 W 和身高 L 关系（）3LW（3）老虎的身长（不含头尾）与体重 问题：问题：如何根据老虎的体长（不包括头尾）估计其体重？模型假设：模型假设：1 将四足动物的躯干（不含头尾）视为质量为 m 的圆柱体，长度为l，截面面积，直径为 d；如图 s 2 把圆柱体的躯干看作一根支撑在四肢上的弹性梁，动物在体重作用下的最大下垂为f，即梁的最大弯曲，根据弹性力学弯曲度理论，有：23sdfl；3 以生物进化学的角度，可认为动物的相对下垂度l已达到一个最合适的数值，也即l为常数。模型建立：模型建立：ldvf42=，42ds=，又231sd

15、flk=为比例常数，从而1kflk2514=，422152144lkklkf=，。4klf=2.轮廓模型 2.轮廓模型 直接利用不同量纲的量之间的比例关系所得到的模型称之为轮廓模型。例 6 商品的包装与成本 例 6 商品的包装与成本 商 品 价格 含量 单价 价格 含量 单价 高露洁牙膏 15.7 元 190g 8.3 元/100g 5.8 元 60g 9.7 元/100g 诗芬洗发液 35.9 元 400ml 9 元/100ml 23.1 元 200ml 11.5 元/100ml 富丽饼干 8.8 元 450g 1.9 元/100g 3.0 元 150g 2 元/100g 奇宝饼 5.9 元

16、 250g 2.3 元/100g 4.3 元 150g 2.87 元/100g 建模分析为什么小包装的商品比大包装的要贵一些?假设:假设:（1）包装只计装包工时和包装材料。（2）不同规格的商品包装外观相似，包装材料相似，至少在价格上没有太大的差异。（3）不同规格的商品装包时工作效率相同。（4）不考虑利润及其他因素对商品价格的影响。符号说明：符号说明：W每件商品净重（产品的含量）)(WC每件商品的总成本 1C每件商品中产品的成本 2C每件商品装包工时投入 3C每件商品包装材料成本 S包装材料用量 WWCWA/)()(=商品单位重量的平均成本.分析：分析：321)(CCCWC+=；，WkC11=W

17、kC22=，3/2433WkSkC=模型：模型：3/1453/245)(,)(+=+=WkkWAWkWkWC 应用于价格预测：对于康尔乃奶粉的两个包装价格 32.4 元/400g；67.1 元/900g，联立方程可以得到；利用该方程可以预测其他规格的同产品的价格，如 预测:，实际价格为 115.9。3/23192.43791.5)(WWWC+=49.126)(,1800=WCW模型分析 模型分析（1）成本的降低率3/4431)()(=WkdWWdAWr是商品量的减函数；（2）支出的节省率3/1431)()(=WkWWrWS也是商品量的减函数。三、曲线拟合与插值 三、曲线拟合与插值 问题：问题：

18、给定一批数据点（输入变量与输出变量的数据），需确定满足特定要求的曲线或曲面。如果输入变量和输出变量都只有一个，则属于一元函数的拟合和插值；而若输入变量有多个，则为多元函数的拟合和插值（有点回归分析的意思）解决方案：解决方案：（1）若要求所求曲线（面）通过所给所有数据点，就是插值问题插值问题；（2）若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合数据拟合，又称曲线拟合或曲面拟合。注意：插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。注意：

19、插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。实例：实例：下面数据是某次实验所得，希望得到 X 和 f 之间的关系？x 1 2 4 7 9 12 13 15 17 f 1.5 3.9 6.6 11.7 15.6 18.8 19.6 20.6 21.1 曲线拟合问题最常用的解法最小二乘法的基本思路 第一步：确定拟合的函数类型),;(21maaaxfyL=，其中为待定系数。（函数类型的确定可以根据内在的规律确定，如果无现成的规则，则可以通过散点图，联系曲线的形状进行分析）

20、maaa,21L第二步：确定的最小二乘准则：要求个已知点与曲线 的距离的平方和最小。maaa,21Ln),(iiyx)(xfy=id=niiixfy12)(用 MATLAB 作拟合 1多项式拟合。作多项式拟合，可利用 mmmaxaxay+=L110a=polyfit(x,y,m)其中 x,y 为给出的数据，m 为多项式的次数。多项式在 x 处的值 y 可用以下命令计算：y=polyval（a，x）2用 MATLAB 作非线性最小二乘拟合 Matlab 的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit 和 lsqnonlin。两个命令都要先建立 M-文件 fun.m，在其中定义函数

21、 f(x)。（1）x=lsqcurvefit(fun,x0,xdata,ydata);（2）x=lsqcurvefit(fun,x0,xdata,ydata,options);（3）x=lsqcurvefit(fun,x0,xdata,ydata,options,grad);（4）x,options=lsqcurvefit(fun,x0,xdata,ydata,);（5）x,options,funval=lsqcurvefit(fun,x0,xdata,ydata,);（6）x,options,funval,Jacob=lsqcurvefit(fun,x0,xdata,ydata,);（7）x

22、=lsqnonlin（fun，x0）；（8）x=lsqnonlin（fun，x0，options）；（9）x=lsqnonlin（fun，x0，options，grad）；（10）x，options=lsqnonlin（fun，x0，）；（11）x，options，funval=lsqnonlin（fun，x0，）；插值插值 一维插值：一维插值：已知个节点1+nniyxii,1,0),(L=，求任意点处的函数值。常用的插值方法有拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite 插值和三次样条插值。*x*y?分段线性插值：将各数据点用折线连接起来?多项式插值：求一个多项式通过所有数据点

23、，可以假设出多项式的系数，最后通过求解方程得到每个系数（拉格朗日插值，用次多项式描述n1+n个点）?样条插值：分段多项式的光滑连接（三次样条插值）?牛顿插值：利用节点之间的各阶差商和差分构造多项式?Hermite 插值：对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值 MATLAB 命令：y=interp1(x0,y0,x,method)method 指定插值的方法，默认为线性插值。其值可为：nearest 最近项插值 linear 线性插值 spline 立方样条插值 cubic 立方插值。所有的插值方法要求 x0 是单调的。当 x0 为等距时可

24、以用快速插值法，使用快速插值法的格式为*nearest、*linear、*spline、*cubic。三次样条插值在 Matlab 中的实现 在 Matlab 中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。这个条件强迫第 1 个和第 2 个三次多项式的三阶导数相等。对最后一个和倒数第 2 个三次多项式也做同样地处理。Matlab 中三次样条插值也有现成的函数：y=interp1(x0,y0,x,spline)；y=spline(x0,y0,x)；pp=csape(x0,y0,conds),pp=csape(x0,y0,conds,v

25、alconds)，y=ppval(pp,x)。其中 x0,y0 是已知数据点，x 是插值点，y 是插值点的函数值。对于三次样条插值，我们提倡使用函数 csape，csape 的返回值是 pp 形式，要求插值点的函数值，必须调用函数 ppval。pp=csape(x0,y0)：使用默认的边界条件，即 Lagrange 边界条件。pp=csape(x0,y0,conds,valconds)中的 conds 指定插值的边界条件，其值可为：complete 边界为一阶导数，一阶导数的值在 valconds 参数中给出，若忽略 valconds参数，则按缺省情况处理。not-a-knot 非扭结条件 p

26、eriodic 周期条件 second 边界为二阶导数，二阶导数的值在 valconds 参数中给出，若忽略 valconds参数，二阶导数的缺省值为0,0。variational 设置边界的二阶导数值为0,0。对于一些特殊的边界条件，可以通过 conds 的一个21矩阵来表示，conds 元素的取值为 0，1，2。conds(i)=j 的含义是给定端点 的ij阶导数，即 conds 的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=2,1表示左边界是二阶导数，右边界是一阶导数，对应的值由 valconds 给出。例例 7 机床加工机床加工 待加工零件的外形根据工艺要求由一

27、组数据给出（在平面情况下），用程控铣床加工时每一刀只能沿),(yxx方向和方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的坐标。y),(yx表中给出的yx,数据位于机翼断面的下轮廓线上，假设需要得到x坐标每改变 0.1 时的坐标。试完成加工所需数据，画出曲线，并求出y0=x处的曲线斜率和范围内的最小值。1513 xyx 0 3 5 7 9 11 12 13 14 15 y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 要求用分段线性和三次样条计算。解 编写以下程序：x0=0 3 5 7 9 11 12 13 14 15;y0=0 1.2 1.7 2

28、.0 2.1 2.0 1.8 1.2 1.0 1.6;x=0:0.1:15;y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,spline);pp1=csape(x0,y0);y4=ppval(pp1,x);pp2=csape(x0,y0,second);y5=ppval(pp2,x);x,y1,y2,y3,y4,y5 subplot(2,2,1)plot(x0,y0,+,x,y2)title(Piecewise linear)subplot(2,2,2)plot(x0,y0,+,x,y3)title(Spline1)subplot(2,2,3)plot(x0,y0

29、,+,x,y4)title(Spline2)dx=diff(x);dy=diff(y3);dy_dx=dy./dx;dy_dx0=dy_dx(1)ytemp=y3(131:151);ymin=min(ytemp);index=find(y3=ymin);xmin=x(index);xmin,ymin 二维插值：二维插值：前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。若节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的高程（节点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这些点的高程（插值）。（1）网格节点：已知个节点 nm

30、njmizyxijji,2,1;,2,1),(LL=，构造一个二元函数通过全部已知节点，即),(yxfz=),(jiijyxfz=，再利用插值。),(yxfz=（2）散乱节点：已知 n 个节点 nizyxiii,2,1),(L=，构造一个二元函数通过全部已知节点，即),(yxfz=),(iiiyxfz=，再利用插值。),(yxfz=二维插值方法：（1）最邻近插值：二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。（2）分片线性插值：将四个插值点（矩形的四个顶点）处的函数值依次简记为：f(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,

31、yj+1)=f4 分两片的函数表达式如下：第一片(下三角形区域)：)()(),(23121jiyyffxxfffyxf+=第二片(上三角形区域)：)()(),(43141ijxxffyyfffyxf+=（3）双线性插值：双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：。其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。)(),(dcybaxyxf+=用 MATLAB 作网格节点数据的插值：z=interp2(x0,y0,z0,x,y,method)method 可以选择 nearest：最邻近插值；linear：双线性插值

32、；cubic：双三次插值；缺省时,双线性插值。如果是三次样条插值，可以使用命令 pp=csape(x0,y0,z0,conds,valconds)，z=fnval(pp,x,y)其中 x0,y0 分别为维和n维向量，z0 为mnm维矩阵，为矩阵，它的行数为zx的维数，列数为的维数，表示得到的插值，具体使用方法同一维插值。y例8 例8 在一丘陵地带测量高程，x和方向每隔100米测一个点，得高程如下表，试插值一曲面，确定合适的模型，并由此找出最高点和该点的高程。y x y 100 200 300 400 500 100 636 697 624 478 450 200 698 712 630 478

33、 420 300 680 674 598 412 400 400 662 626 552 334 310 解 编写程序如下：x=100:100:500;y=100:100:400;z=636 697 624 478 450;698 712 630 478 420;680 674 598 412 400;662 626 552 334 310;pp=csape(x,y,z)xi=100:10:500;yi=100:10:400 cz1=fnval(pp,xi,yi)cz2=interp2(x,y,z,xi,yi,spline)i,j=find(cz1=max(max(cz1)x=xi(i),y=

34、yi(j),zmax=cz1(i,j)用 MATLAB 作散点数据的插值计算 cz=griddata（x,y,z,cx,cy,method）method 可以选择nearest：最邻近插值；linear ：双线性插值；cubic：双三次插值；v4：Matlab 提供的插值方法；缺省时,双线性插值 例例 9 在某海域测得一些点（x,y）处的水深 z 由下表给出，在矩形区域（75,200）(-50,150)内画出海底曲面的图形。x 129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5 y 7.5 141.5 23 147 22

35、.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 解 编写程序如下：x=129,140,103.5,88,185.5,195,105,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;xi=75:1:200;yi=-50:1:150;zi=griddata(x,y,z,xi,yi,cubic)su

36、bplot(1,2,1)plot(x,y,*)subplot(1,2,2)mesh(xi,yi,zi)作业：作业：1 速度为的风吹在迎风面积为的风车上，空气密度为vS，用量纲分析方法确定风车获得的功率与P,Sv的关系。2 原子弹爆炸是巨大的能量从爆炸点以冲击波的形式向四周传播。据分析在时刻冲击波达到的半径tr与释放能量E，大气的密度，大气压强p有关（设时）。用量纲分析方法证明0=t0=r)()(32655/12EtpEtr=，为未定函数。3 质量的小球以速度v竖直上抛，阻力与速度成正比，比例系数。初始位置为，mk0=xx轴竖直向上，运动方程为vxxmgxkxm=+)0(,0)0(,0，方程的解可以表示为),;(kmgvtxx=。选择两种特殊尺度将问题无量纲化，并讨论当很小时求近似解的可能性。k4 雨滴匀速下降，空气阻力与雨滴的表面积和速度平方的乘积成正比，请确定雨速与雨滴质量的关系。