回归分析合金强度y与其中含碳量x有密切关系,如下表根据此表建立y(x)。并对结果作可信度进行检验、判断x对y影响是否显著、检查数据 中有无异常点、由x的取值对y作出预测。回归分析解:在x-y平面上画散点图,直观地知道y与x大致为线性关系。用命令polyf it(x,y,1)可得y=140.6194x+27.0269。作回归分析用命令[b,bint,r,rint,ststs]=re gress(y,x,alpha)可用help查阅此命令的具体用法残差及置信区间可以用rcoplot(r,rint)画图回归 分析设回归模型为y=β0+β1x,在MATLAB命令窗口中键入下列命令进行回归分析(px_reg11.m)x=0.1:0. 01:0.18;x=[x,0.2,0.21,0.23]'';y=[42,41.5,45,45.5,45,47.5,49,55,50 ,55,55.5,60.5]'';X=[ones(12,1),x];[b,bint,r,rint,stats]=regress( y,X,0.05);b,bint,stats,rcoplot(r,rint)数值积分解题思路:数据实际上表示了两条曲线,实 际上我们要求由两曲线所围成的图形的面积。解此问题的方法是数值积分的方法。具体解时我们遇到两个问题:1。数据如何输入;2。没有 现成的命令可用。数值积分(px_wj11.m)对于第一个问题,我们可把数据拷贝成M文件(或纯文本文件)。然后,利用数据绘制平 面图形。键入?loadmianji.txtA=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-a1d=4.2414e+ 004数值积分至此,问题可以说得到了解决。?之所以说还有问题,是我们觉得误差较大。但计算方法的理论给了我们更精确计算方法。只 是MATLAB没有相应的命令。?想得到更理想的结果,我们可以自己设计解决问题的方法。(可以编写辛普森数值计算公式的程序,或用拟合 的方法求出被积函数,再利用MATLAB的命令quad,quad8)数值微分已知20世纪美国人口统计数据如下,根据数据计算人口增 长率。(其实还可以对于后十年人口进行预测)数值微分解题思路:设人口是时间的函数x(t).于是人口的增长率就是x(t)对t的导数 .如果计算出人口的相关变化率。那么人口增长满足,它在初始条件x(0)=x0下的解 为.(用以检查计算结果的正确性)数值微分解:此问题的特点是以离散变量给出函数x(t),所以就 要用差分来表示函数x(t)的导数.常用后一个公式。(因为,它实际上是用二次插值函数来代替曲线x(t))即常用三点公式来代替函数在 各分点的导数值:数值微分MATLAB用命令diff按两点公式计算差分;此题自编程序用三点公式计算相关变化率.编程如下(diff 3.m):fori=1:length(x)ifi==1r(1)=(-3x(1)+4x(1+1 )-x(1+2))/(20x(1));elseifi~=length(x)r(i)=(x(i+1) -x(i-1))/(20x(i));elser(length(x))=(x(length(x)-2) -4x(length(x)-1)+3x(length(x)))/(20x(length(x)));endend r=r;数值微分保存为diff3.m文件听候调用.再在命令窗内键入X=[1900,1910,1920,1930,1940,1 950,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),qua d(‘fun’,a,b)等.微分方程数值解(单摆问题)单摆问题的数学模型是在初始角度不大时,问题可以得到很好地解决,但如 果初始角较大,此方程无法求出解析解.现问题是当初始角为100和300时,求出其解,画出解的图形进行比较。微分方程数值解(单摆问题 )解:若θ0较小,则原方程可用来近似.其解析解为θ(t)=θ0cosωt, .?若不用线性方程来近似,那么有两个模型:微分方程数值解(单摆问题)取g=9.8,l=25,100 =0.1745,300=0.5236.用MATLAB求这两个模型的数值解,先要作如下的处理:令x1=θ,x2=θ’,则模型变为 微分方程数值解(单摆问题)再编函数文件(danbai.m)functionxdot=danbai(t,x)xdot=zer os(2,1);xdot(1)=x(2);xdot(2)=-9.8/25sin(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,A,b,A1,b1,lb,ub)线性规 划有约束极小问题Findxthatminimizesf(x)=-5x1-4x2-6x3subjecttox1-x 2+x3≦203x1+2x2+4x3≦423x1+2x2≦300≦x1,0≦x2,0≦x3线性规划有约束极小问题线性 规划有约束极小问题解问题把问题极小化并将约束标准化线性规划有约束极小问题键入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.5714最大.线性规划有约束极小问题解问题线 性规划有约束极小问题解:键入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,l b,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_文件.functiony=e1511 (x)%目标函数的m_文件y=exp(x(1))(4x(1)^2+2x(2)^2+4x(1)x(2)+2x(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;非线性规划有约束极小问题第2步:运行程序.键入x0=[-1,1 ];a1=[1,1];b1=0;[x,f,exitflab,output]=fmincon(‘e1511’,x0,[],[],a 1,b1,[],[],’e1511b’)得结果.输出结果的意义:经过4次迭代(iterations:4)收敛到了(exitfa g=1)最优解x(1)=-1.2247,x(2)=1.2247,目标函数最优值为1.8951.非线性无约束极小问题用命 令x=fmin(''f'',x0)。或用命令x=fminu(''f'',x0),或用命令x=fmins(''f'',x0)。非线性最小 二乘问题用命令x=leastsq(''f'',x0),或用命令x=curvefit(''f'',x0)。二次规划用命 令x=qp(H,c,A,b)。关于这些命令的详细使用规则和例子,用借助help进行查阅。回归分析前面我们曾学过拟合。但从统计 的观点看,对拟合问题还需作回归分析。例如:有描述问题甲和问题乙的两组数据(x,y)和(x,z)。设x=[1,2,3,4];y=[1 .0,1.3,1.5,2.02.3];z=[0.6,1.95,0.9,2.85,1.8]。如果在平面上画出散点图,那么问题甲的四个 点基本在一条直线上而问题乙的四个点却很散乱。如果都用命令polyfit(x,y,1),polyfit(x,z,1)来拟合,将得到同 一条直线。回归分析自然对问题甲的信任程度会高于对问题乙的信任程度。所以有必要对所得结果作科学的评价分析。回归分析就是解决这种问 题的科学方法。下面结合三个具体的例子介绍MATLAB实现回归分析的命令。数学建模中的数据处理方法范筑军主要内容曲线插值与 拟合数值微分与积分微分方程数值解优化问题?回归分析判别分析曲线插值与拟合一维插值二维插值曲线拟合一维插值对表 格给出的函数,求出没有给出的函数值。在实际工作中,经常会遇到插值问题。下表是待加工零件下轮廓线的一组数据,现需要得到x坐标每改 变0.1时所对应的y的坐标.一维插值下面是关于插值的两条命令(专门用来解决这类问题):y=interp1(x0,y0,x, ’method’)分段线性插值y=spline(x0,y0,x)三次样条插值 x0,y0是已知的节点坐标,是同维向量。y对应于x处的插值。y与x是同维向量。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)y=spline(x0,y0,x '');plot(x,y)%用三次样条插值完成第二步工作练习对y=1/(1+x2),-5≤ x≤5,用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个点的温度值,试求出此薄片的温度分布,并绘出等温线图。(数据如下表)二维插值(px_lc21 .m)temps=[82,81,80,82,84;79,63,61,65,87;84,84,82,85,86];mesh(te mps)%根据原始数据绘出温度分布图,可看到此图的粗造度。二维插值%下面开始进行二维函数的三阶插值。wid th=1:5;depth=1:3;di=1:0.2:3;wi=1:0.2:5;[WI,DI]=meshgrid(wi,d i);%增加了节点数目ZI=interp2(width,depth,temps,WI,DI,''cubic'');%对数据( width,depth,temps)进%行三阶插值拟合。surfc(WI,DI,ZI)contour(WI,DI,Z I)二维插值曲线拟合假设一函数g(x)是以表格形式给出的,现要求一函数f(x),使f(x)在某一准则下与表格函数(数据)最为 接近。由于与插值的提法不同,所以在数学上理论根据不同,解决问题的方法也不同。此处,我们总假设f(x)是多项式。曲线拟合问题 :弹簧在力F的作用下伸长x厘米。F和x在一定的范围内服从虎克定律。试根据下列数据确定弹性系数k,并给出不服从虎克定律时的近似公式。 曲线拟合解题思路:可以用一阶多项式拟合求出k,以及近似公式。在MATLAB中,用以下命令拟合多项式。polyfit(x0, y0,n)一般,也需先观察原始数据的图像,然后再确定拟和成什么曲线。曲线拟合(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)曲线拟合注意:有时,面对一个实际问题,究竟是用插值还是用拟合不好确定,还需大家在实际中仔细区分。同时,大家(包括学过计算方法的同学)注意去掌握相应的理论知识。数值微分与积分数值积分?数值微分数值积分先看一个例子:现要根据瑞士地图计算其国土面积。于是对地图作如下的测量:以西东方向为横轴,以南北方向为纵轴。(选适当的点为原点)将国土最西到最东边界在x轴上的区间划取足够多的分点xi,在每个分点处可测出南北边界点的对应坐标y1,y2。用这样的方法得到下表根据地图比例知18mm相当于40km,试由上表计算瑞士国土的近似面积。(精确值为41288km2)。数值积分 |
|