《数学建模中的数据处理方法》由会员分享,可在线阅读,更多相关《数学建模中的数据处理方法(80页珍藏版)》请在金锄头文库上搜索。
1、数学建模中的数据处理方法,范筑军,主要内容,曲线插值与拟合 数值微分与积分 微分方程数值解 优化问题 回归分析 判别分析,曲线插值与拟合,一维插值 二维插值 曲线拟合,一维插值,对表格给出的函数,求出没有给出的函数值。 在实际工作中,经常会遇到插值问题。 下表是待加工零件下轮廓线的一组数据,现需要得到x坐标每改变0.1时所对应的y的坐标.,一维插值,下面是关于插值的两条命令(专门用来解决这类问题): y=interp1(x0,y0,x,method) 分段线性插值 y=spline(x0,y0,x) 三次样条插值 x0,y0是已知的节点坐标,是同维向量。 y对应于x处的插值。y与x是同维向量。
2、 method可选nearest(最近邻插值),linear(线性插值),spline(三次样条插值),cubic(三次多项式插值),一维插值,解决上述问题,我们可分两步:用原始数据绘图作为选用插值方法的参考.确定插值方法进行插值计算,一维插值(px_lc11.m),对于上述问题,可键入以下的命令: 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 plot(x0,y0) %完成第一步工作 x=0:0.1:15; y=interp1(x0,y0,x); %用分段线性插值完成第二步工作 plot(x,y)
3、y=spline(x0,y0,x); plot(x,y) %用三次样条插值完成第二步工作,练习,对y=1/(1+x2),-5x5,用n(=11)个节点(等分)作上述两种插值,用m(=21)个插值点(等分)作图,比较结果。(see:px_ex_lc1.m) 在某处测得海洋不同深度处水温如下表:求深度为500、1000、1500米处的水温。 (see:px_ex_lc2.m),二维插值,MATLAB中二维插值的命令是: z=interp2(x0,y0,z0,x,y,meth),二维插值,在一个长为5个单位,宽为3个单位的金属薄片上测得15个点的温度值,试求出此薄片的温度分布,并绘出等温线图。(数据
4、如下表),二维插值(px_lc21.m),temps=82,81,80,82,84;79,63,61,65,87;84,84,82,85,86; mesh(temps) %根据原始数据绘出温度分布图,可看到此图的粗造度。,二维插值,%下面开始进行二维函数的三阶插值。width=1:5; depth=1:3; di=1:0.2:3; wi=1:0.2:5;WI,DI=meshgrid(wi,di);%增加了节点数目ZI=interp2(width,depth,temps,WI,DI,cubic); % 对数据(width,depth,temps)进% 行三阶插值拟合。surfc(WI,DI,ZI
5、) contour(WI,DI,ZI),二维插值,曲线拟合,假设一函数g(x)是以表格形式给出的,现要求一函数f(x),使f(x)在某一准则下与表格函数(数据)最为接近。 由于与插值的提法不同,所以在数学上理论根据不同,解决问题的方法也不同。 此处,我们总假设f(x)是多项式。,曲线拟合,问题:弹簧在力F的作用下伸长x厘米。F和x在一定的范围内服从虎克定律。试根据下列数据确定弹性系数k,并给出不服从虎克定律时的近似公式。,曲线拟合,解题思路:可以用一阶多项式拟合求出k,以及近似公式。 在MATLAB中,用以下命令拟合多项式。 polyfit(x0,y0,n) 一般,也需先观察原始数据的图像,然
6、后再确定拟和成什么曲线。,曲线拟合(px_lc31.m),对于上述问题,可键入以下的命令: 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;plot(x,F,.) 从图像上我们发现:前5个数据应与直线拟合,后5个数据应与二次曲线拟合。于是键入 : a=polyfit(x(1:5),F(1:5),1); a=polyfit(x(5:9),F(5:9),2),曲线拟合,注意:有时,面对一个实际问题,究竟是用插值还是用拟合不好确定,还需大家在实际中仔细区分。同时,大家(包括学过计算方法的同学)注意去掌握相应的理
7、论知识。,数值微分与积分,数值积分 数值微分,数值积分,先看一个例子: 现要根据瑞士地图计算其国土面积。于是对地图作如下的测量:以西东方向为横轴,以南北方向为纵轴。(选适当的点为原点)将国土最西到最东边界在x轴上的区间划取足够多的分点xi,在每个分点处可测出南北边界点的对应坐标y1 ,y2。用这样的方法得到下表 根据地图比例知18mm相当于40km,试由上表计算瑞士国土的近似面积。(精确值为41288km2)。,数值积分,数值积分,解题思路:数据实际上表示了两条曲线,实际上我们要求由两曲线所围成的图形的面积。 解此问题的方法是数值积分的方法。具体解时我们遇到两个问题: 1。数据如何输入; 2。
8、没有现成的命令可用。,数值积分(px_wj11.m),对于第一个问题,我们可把数据拷贝成M文件(或纯文本文件)。 然后,利用数据绘制平面图形。键入 load mianji.txt A=mianji; plot(A(:,1),A(:,2),r,A(:,1),A(:,3),g),数值积分,数值积分,接下来可以计算面积。键入: a1=trapz(A(:,1)*40/18,A(:,2)*40/18); a2=trapz(A(:,1)*40/18,A(:,3)*40/18); d=a2-a1 d = 4.2414e+004,数值积分,至此,问题可以说得到了解决。 之所以说还有问题,是我们觉得误差较大。但
9、计算方法的理论给了我们更精确计算方法。只是MATLAB没有相应的命令。 想得到更理想的结果,我们可以自己设计解决问题的方法。(可以编写辛普森数值计算公式的程序,或用拟合的方法求出被积函数,再利用MATLAB的命令quad,quad8),数值微分,已知20世纪美国人口统计数据如下,根据数据计算人口增长率。(其实还可以对于后十年人口进行预测),数值微分,解题思路:设人口是时间的函数x(t).于是人口的增长率就是x(t)对t的导数.如果计算出人口的相关变化率 。那么人口增长满足 ,它在初始条件x(0)=x0下的解为 .(用以检查计算结果的正确性),数值微分,解:此问题的特点是以离散变量给出函数x(t
10、),所以就要用差分来表示函数x(t)的导数.,常用后一个公式。(因为,它实际上是用二次插值函数来代替曲线x(t))即常用三点公式来代替函数在各分点的导数值:,数值微分,MATLAB用命令diff按两点公式计算差分;此题自编程序用三点公式计算相关变化率.编程如下(diff3.m): for i=1:length(x)if i=1r(1)=(-3*x(1)+4*x(1+1)-x(1+2)/(20*x(1);elseif i=length(x)r(i)=(x(i+1)-x(i-1)/(20*x(i);elser(length(x)=(x(length(x)-2)-4*x(length(x)-1)+3
11、*x(length(x)/(20*x(length(x);end end r=r;,数值微分,保存为diff3.m文件听候调用.再在命令窗内键入 X=1900,1910,1920,1930,1940,1950,1960,1970,1980,1990; x=76.0, 92.0, 106.5, 123.2, 131.7, 150.7, 179.3, 204.0, 226.5, 251.4; diff3; 由于r以离散数据给出,所以要用数值积分计算.键入 x(1,1)*exp(trapz(X(1,1:9),r(1:9) 数值积分命令:trapz(x),trapz(x,y),quad(fun,a,b
12、)等.,微分方程数值解(单摆问题),单摆问题的数学模型是在初始角度不大时,问题可以得到很好地解决,但如果初始角较大,此方程无法求出解析解.现问题是当初始角为100和300时,求出其解,画出解的图形进行比较。,微分方程数值解(单摆问题),解:若0较小,则原方程可用 来近似.其解析解为(t)= 0cost, . 若不用线性方程来近似,那么有两个模型:,微分方程数值解(单摆问题),取g=9.8,l=25, 100=0.1745, 300=0.5236.用MATLAB求这两个模型的数值解,先要作如下的处理:令x1=,x2=,则模型变为,微分方程数值解(单摆问题),再编函数文件(danbai.m) fu
13、nction xdot=danbai(t,x) xdot=zeros(2,1); xdot(1)=x(2);xdot(2)=-9.8/25*sin(x(1);,微分方程数值解(单摆问题),在命令窗口键入() t,x=ode45(danbai,0:0.1:20,0.1745,0); t,y=ode45(danbai,0:0.1:20,0.5236,0); plot(t,x(:,1),r,t,y(:,1),k);,优化问题,线性规划有约束极小问题 非线性规划有约束极小问题 非线性无约束极小问题 非线性最小二乘问题 二次规划,线性规划有约束极小问题,模型 用命令 x, fval= linprog(f
14、,A,b,A1,b1,lb,ub),线性规划有约束极小问题,Find x that minimizes f(x)=-5x1-4x2-6x3 subject to x1-x2+x320 3x1+2x2+4x342 3x1+2x230 0x1, 0x2,0x3,线性规划有约束极小问题,线性规划有约束极小问题,解问题把问题极小化并将约束标准化,线性规划有约束极小问题,键入c=-2,-3,5;a=-2,5,-1; b=-10;a1=1,1,1;b1=7;LB=0,0,0; x,y=linprog(c,a,b,a1,b1,LB) 得当X=(6.4286,0.5714,0.0000)时, z=-14.57
15、14最大.,线性规划有约束极小问题,解问题,线性规划有约束极小问题,解:键入 c=-2,-1,1;a=1,4,-1;2,-2,1; b=4;12;a1=1,1,2;b1=6; lb=0;0;-inf;ub=inf;inf;5; x,z=linprog(c,a,b,a1,b1,lb,ub) 得当X=(4.6667,0.0000,0.6667)时, z=-8.6667最小.,非线性规划有约束极小问题,模型:MATLAB求解此问题的命令是: x,fval,exitflag,output,lambda,grad,hessian=fmincon(fun,x0,A,b,A1,b1,LB,UB,nonlcon,options,p1,p2,) fun是目标函数的m_文件名.nonlcon是约束函数C(x)和C1(x)的m_文件名.文件输出为C,C1.,非线性规划有约束极小问题,求解最优化问题,非线性规划有约束极小问题,第1步: 建立目标函数和非线性约束的m_文件. function y=e1511(x)% 目标函数的m_文件 y=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);function c1,c2=e1511b(x)% 非线性约束的m_文件 c1=1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10; c2=0;,