MATLAB数值算法数学建模

数学建模-MATLAB的使用简介一、 曲线插值与拟合二、 数值微分与积分三、 微分方程数值解四、 优化问题五、 回归分析1. 一维插值对表格给出的函数，求出没有给出的函数值。在实际工作中，经常会遇到插值问题。例1：表1是待加工零件下轮廓线的一组数据，现需要得到x坐标每改变0.1时所对应的y的坐标.x035791112131415y01.21.72.02.12.01.81.21.01.6下面是关于插值的两条命令(专门用来解决这类问题)：y=interp1(x0,y0,x) 分段线性插值y=spline(x0,y0,x) 三次样条插值其中x0,y0是已知的节点坐标，是同维向量。y对应于x处的插值。y与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.6plot(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),-5x5，用n（=11）个节点（等分）作上述两种插值，用m（=21）个插值点（等分）作图，比较结果。解：键入并运行如下命令n=11;m=21;x=-5:10/(m-1):5;y=1./(1+x.2);xo=-5:10/(n-1):5;yo=1./(1+xo.2);y1=interp1(xo,yo,x);y2=spline(xo,yo,x);plot(x,y,r,x,y1,b,x,y2,k)练习：在某处测得海洋不同深度处水温如下：深度44671495014221634水温7.044.283.402.542.13求深度为500、1000、1500米处的水温。解：输入程序：D=446,714,950,1422,1634;T=7.04,4.28,3.40,2.54,2.13;Di=500,1000,1500;Ti=interp1(D,T,Di)MATLAB的命令interp1(X,Y,Xi,method)用于一元插值.其中Method可选nearest(最近邻插值),linear(线性插值),spline(三次样条插值),cubic(三次多项式插值)2. 二维插值MATLAB中二维插值的命令是：z=interp2（x0,y0,z0,x,y,meth） 例2：在一个长为5个单位，宽为3个单位的金属薄片上测得15个点的温度值，试求出此薄片的温度分布，并绘出等温线图。（数据如下表）y x 1 2 3 4 5 1 82 81 80 82 84 2 79 63 61 65 87 3 84 84 82 85 86程序：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)contour(WI,DI,ZI) 3曲线拟合假设一函数g(x)是以表格形式给出的，现要求一函数f(x)，使f(x)在某一准则下与表格函数（数据）最为接近。由于与插值的提法不同，所以在数学上理论根据不同，解决问题的方法也不同。此处，我们总假设f(x)是多项式。例3：弹簧在力F的作用下伸长x厘米。F和x在一定的范围内服从虎克定律。试根据下列数据确定弹性系数k,并给出不服从虎克定律时的近似公式。x1247912131517F1.53.96.611.715.618.819.620.621.1解题思路：可以用一阶多项式拟合求出k，以及近似公式。在MATLAB中，用以下命令拟合多项式。polyfit（x0,y0,n）一般，也需先观察原始数据的图像，然后再确定拟和成什么曲线。对于上述问题,可键入以下的命令: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)得注意：有时，面对一个实际问题，究竟是用插值还是用拟合不好确定，还需大家在实际中仔细区分。同时，大家（包括学过计算方法的同学）注意去掌握相应的理论知识。4数值积分先看一个例子：例4现要根据瑞士地图计算其国土面积。于是对地图作如下的测量：以西东方向为横轴，以南北方向为纵轴。（选适当的点为原点）将国土最西到最东边界在x轴上的区间划取足够多的分点xi，在每个分点处可测出南北边界点的对应坐标y1 ，y2。用这样的方法得到下表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根据地图比例知18mm相当于40km，试由上表计算瑞士国土的近似面积。（精确值为41288km2）。解题思路：数据实际上表示了两条曲线，实际上我们要求由两曲线所围成的图形的面积。解此问题的方法是数值积分的方法。具体解时我们遇到两个问题：1。数据如何输入；2。没有现成的命令可用。解：对于第一个问题，我们可把数据考备成M文件（或纯文本文件）。然后，利用数据绘制平面图形。键入load mianji.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）5数值微分实际的例：已知20世纪美国人口统计数据如下，根据数据计算人口增长率。（其实还可以对于后十年人口进行预测）年份1900191019201930194019501960197019801990人口10676.092.0106.5123.2131.7150.7179.3204.0226.5251.4解题思路：设人口是时间的函数x(t).于是人口的增长率就是x(t)对t的导数.如果计算出人口的相关变化率。那么人口增长满足，它在初始条件x(0)=x0下的解为.（用以检查计算结果的正确性）解：此问题的特点是以离散变量给出函数x(t),所以就要用差分来表示函数x(t)的导数.常用后一个公式。（因为，它实际上是用二次插值函数来代替曲线x(t)）即常用三点公式来代替函数在各分点的导数值：MATLAB用命令diff按两点公式计算差分;此题自编程序用三点公式计算相关变化率.编程如下:for i=1:length(x) if i=1 r(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); else r(length(x)=(x(length(x)-2)-4*x(length(x)-1)+3*x(length(x)/(20*x(length(x); endendr=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)等.6：微分方程数值解（单摆问题）单摆问题的数学模型是在初始角度不大时,问题可以得到很好地解决,但如果初始角较大,此方程无法求出解析解.现问题是当初始角为100和300时,求出其解,画出解的图形进行比较。解:若0较小,则原方程可用来近似.其解析解为(t)= 0cost,. 若不用线性方程来近似,那么有两个模型: 取g=9.8,l=25, 100=0.1745, 300=0.5236.用MATLAB求这两个模型的数值解,先要作如下的处理:令x1=,x2=,则模型变为再编函数文件function 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);参考书:数学实验,高等教育出版社7解优化问题线性规划有约束极小问题模型用命令x, fval= linprog(f,A,b,A1,b1,lb,ub)例1：Find x that minimizes f(x)=-5x1-4x2-6x3subject tox1-x2+x3203x1+2x2+4x3423x1+2x2300x1, 0x2,0x3First, enter the coefficients:f = -5; -4; -6A = 1 -1 1 3 2 4 3 2 0;b = 20; 42; 30;lb = zeros(3,1);Next, call a linear programming routine:x,fval,exitflag,output,lambda = linprog(f,A,b,lb)；Entering x, fval,lambda.ineqlin, and lambda.lower getsx = 0.0000 15.0000 3.0000fval = -78.0000和其它信息。例2：解问题把问题极小化并将约束标准化键入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最大.例3: 解问题解:键入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最小.非线性规划有约束极小问题模型1用命令x=constr(f ,x0)。ExamplesFind values of x that minimize f(x)=-x1x2x3, starting at the point x = 10; 10; 10 and subject to the constraints 0x1+2x2+2x372.-x1-2x2-2x30,x1+2x2+2x372,第一步：编写M文件function f,g=myfun(x)f=-x(1)*x(2)*x(3);g(1)=-x(1)-2*x(2)-2*x(3);g(2)=x(1)+2*x(2)+2*x(3)-72;第二步：求解在MATLAB工作窗中键入x0=10,10,10;x=constr(myfun,x0)即可模型2：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;第2步: 运行程序.键入x0=-1,1;a1=1,1;b1=0;x,f,exitflab,output=fmincon(e1511,x0,a1,b1,e1511b)得结果.输出结果的意义:经过4次迭代(iterations:4)收敛到了(exitfag=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进行查阅。参考书:科学计算技术与MATLAB,科学出版社8回归分析前面我们曾学过拟合。但从统计的观点看，对拟合问题还需作回归分析。例如：有描述问题甲和问题乙的两组数据(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实现回归分析的命令。1 合金强度y与其中含碳量x有密切关系，如下表x0.100.110.120.130.140.150.160.170.180.200.210.23y42.041.545.045.545.047.549.055.050.055.055.560.5根据此表建立y(x)。并对结果作可信度进行检验、判断x对y影响是否显著、检查数据中有无异常点、由x的取值对y作出预测。2 将17至19岁的运动员每两岁一组分为7组，每组两人测量其旋转定向能力，以考察年龄(x)对这种运动能力(y)的影响。现得到一组数据如下表年龄17192123252729第一人2048251326153002612031935第二人243528112633142692257213试建立关系y(x)，并作必要的统计分析。3 某厂生产的某产品的销售量与竞争对手的价格x1和本厂的价格x2有关。下表是该产品在10个城市的销售记录。x1120140190130155175125145180150x210011090150210150250270300250y(个)10210012077469326696585试建立关系y(x1，x2)，对结果进行检验。若某城市本厂产品售价160（元），对手售价170（元），预测此产品在该城市的销售量。解1：在x-y平面上画散点图，直观地知道y与x大致为线性关系。用命令polyfit(x,y,1)可得y=140.6194x+27.0269。作回归分析用命令b,bint,r,rint,ststs=regress(y,x,alpha) 可用help查阅此命令的具体用法残差及置信区间可以用rcoplot(r,rint)画图设回归模型为 y=0+1x,在MATLAB命令窗口中键入下列命令进行回归分析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)得结果和图b = 27.0269 140.6194bint = 22.3226 31.7313111.7842 169.4546 stats =0.9219 118.0670 0.0000结果含义为0=27.0269 1=140.61940的置信区间是 22.3226,31.73131的置信区间是 111.7842,169.4546R2=0.9219 F=118.0670, p10-4.R是衡量y与x的相关程度的指标，称为相关系数。R越大，x与y关系越密切。通常R大于0.9才认为相关关系成立。F是一统计指标p是与F对应的概率，当 p0.05时，回归模型成立。此例中 p=0 10-40.05，所以，所得回归模型成立。观察所得残差分布图，看到第8个数据的残差置信区间不含零点，此点视为异常点，剔除后重新计算。此时键入：X(8,:)=;y(8)=;b,bint,r,rint,stats=regress(y,X);b,bint,stats,rcoplot(r,rint)得：b = 27.0992 137.8085bint = 23.8563 30.3421 117.8534 157.7636stats =0.9644 244.0571 0.0000可以看到：置信区间缩小；R2、F变大，所以应采用修改后的结果。解2：在x-y平面上画散点图，直观地知道y与x大致为二次函数关系。设模型为y=a1x2+a2x+a3此问题可以利用命令polyfit(x,y,2)来解，也可以象上题一样求解。下面介绍用命令polytool来解。首先在命令窗口键入x=17:2:29;x=x,x;y=20.48,25.13,26.15 30,26.1,20.3,19.35,24.35,28.11,26.3,31.4,26.92,25.7,21.3;polytool（x,y,2）得到一个交互式窗口窗口中绿线为拟合曲线、红线为y的置信区间、可通过移动鼠标的十字线或通过在窗口下方输入来设定x值，窗口左边则输出与x对应的y值及y的置信区间。通过左下方的Export下拉菜单可输出回归系数等。更详细的解释可通过help查阅。解3.这是一个多元回归问题。若设回归模型是线性的，即设y=0+1x1+2x2那么依然用regress(y,x,alpha)求回归系数。键入x1=120,140,190,130,155,175,125,145,180,150;x2=100,110,90,150,210,150,250,270,300,250;y=102,100,120,77,46,93,26,69,65,85;x=ones(10,1),x1,x2;b,bint,r,rint,stats=regress(y,x);b,bint,stats,得b = 66.5176 0.4139 -0.2698bint = -32.5060 165.5411 -0.2018 1.0296 -0.4611 -0.0785stats =0.6527 6.5786 0.0247p=0.0247,若显著水平取0,01，则模型不能用；R2=0.6527较小;0，1的置信区间包含零点。因此结果不理想。于是设模型为二次函数。此题设模型为纯二次函数：y=0+1x1+2x2+11x12+22x22MATLAB提供的多元二项式回归命令为rstool(x,y,model,alpha).其中alpha为显著水平、model在下列模型中选一个：linear（线性）purequadratic（纯二次）interaction（交叉）quadratic（完全二次）对此例，在命令窗中键入x(:,1)=;rstool(x,y,purequadratic)得到一个对话窗：其意义与前面的对话窗意义类似。若要回答“本厂售价160，对手售价170，预测该市销售量”的问题，只需在下方窗口中分别肩入160和170，就可在左方窗口中读到答案及其置信区间。下拉菜单Export向工作窗输出数据具体操作为：弹出菜单，选all，点击确定。此时可到工作窗中读取数据。可读数据包括：beta（回归系数） rmse（剩余标准差） residuals （残差）本题只要键入beta,rmse,residuals得beta = -312.5871 7.2701 -1.7337 -0.0228 0.0037rmse = 16.6436residuals = 6.6846 -12.6703 -0.2013 6.4855 -19.6533 7.9989 -11.4737 5.4303 -4.9932 22.3926大家不妨用此命令选其它模型作一个比较。注：关于各种模型，可查阅有关理论参考书:数学实验,高等教育出版社参考书:科学计算技术与MATLAB,科学出版社MATLAB中还包括神经网络工具箱，小波分析工具箱，在网上还可以下载遗传算法工具箱，有兴趣的同学可以借这次机会，结合学习MATLAB，好好学习一下相关理论知识。最后，祝大家学习，竞赛都取得成功。谢谢大家。16