数学建模算法章节.ppt

1、武汉科技大学,数学建模,数学建模概述,数值计算方法建模,基本方法建模,数学规划方法建模,统计分析方法建模,第1章,第2章,第3章,第4章,第5章,目 录,图论方法建模,第6章,第2章 数值计算方法建模,2.1 非线性方程求解 非线性方程简介 次代数方程 超越方程,第2章 数值计算方法建模,非线性方程求解的MATLAB实现 x,fval,exitflag,output=fzero(f ,x0,options) r=roots(c) x,fval,exitflag,output=fsolve(f,x0,options),第2章 数值计算方法建模,例1 在无阻尼强迫振荡的研究中会经常遇到函数 . 试

2、求一点 ，满足 。 解：建立M函数文件h.m： function y=h(x) y=x*sin(x)-1; 在MATLAB指令窗中输入下面指令 x0=0,2;x,fval,exitflag=fzero(h,x0) 运行得结果为：x =1.1142，fval =2.2204e-016，exitflag =1，即所求非线性方程的解为1.1142.,第2章 数值计算方法建模,例2 求函数 的零点. 解：（1）为确定其零点的大体位置，先做出它的图形 ； （2）将图形放大得5个零点，利用ginput指令取其坐标； （3） 利用指令fzero 计算其相应精确坐标,图2.1 图形放大法,第2章 数值计算方法

3、建模,例3 求解非线性方程组 解：建立M文件f.m如下： function y=f(x) y(1)=x(1)2+x(2)2-4;y(2)=x(1)2-x(2)2-1; 在MATLAB指令窗中输入指令 x0=2,2;x,fval,exitflag=fsolve(f,x0,),第2章 数值计算方法建模,建模示例：贷款问题 问题： 某人从银行贷款购房，若他今年初贷款10万元，月利率0.5%，每月还1000元，试计算他每年末欠银行多少钱，多少时间才能还清？如果要10年还清，每月需还多少？,第2章 数值计算方法建模,1. 问题的分析 每月的还款金额应包括本金和利息两部分 。 2. 模型的建立 记第 个月

4、初此人欠银行 元，月利率为 ，每月还款 元，则 建立数学模型如下：,将上式依次递推可得,第2章 数值计算方法建模,令 得,将所给数据代入上式可得,=1110.3（元）,第2章 数值计算方法建模,课堂练习： 在一条20m宽的道路两侧，分别安装了一只2kW和一只3kW的路灯，它们离地面的高度分别为5m和6m（如下图）.在漆黑的夜晚，当两只路灯开启时，试求两只路灯连线的路面上最暗和最亮的点.（提 示：光源在 点的照度为 ，其中 为光源的功率， 为 光源到 的光线与水平面的夹角，为光源到 点的距离.）,第2章 数值计算方法建模,2.2 线性方程组的数值解法 线性方程组简介 阶线性方程组 矩阵形式,第2

5、章 数值计算方法建模,线性方程组求解的MATLAB实现 x=Ab A=sparse(r,c,v,m,n) full(A),第2章 数值计算方法建模,例1 求解下列线性方程组 解：在MATLAB指令窗中输入下面指令 ： A=1 1 0 0;2 -1 5 0;0 3 -4 2;0 0 2 -6;b=5;-9;19;2; x=Ab 运行得结果为： x=-2.0000,7.0000,0.4000,-0.2000,第2章 数值计算方法建模,例2 验证Hilbert矩阵是一个典型的病态矩阵 （a）用4阶Hilbert矩阵求解 的精确解（用分数表示所有的元素并进行精确计算）： （b）使用精度为4位有效数字的

6、算术计算求解,第2章 数值计算方法建模,解：（a） A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5 1/6;1/4 1/5 1/6 1/7; b=1;0;0;0; x1=Ab 运行得（a）的解为： x=16.0000, -120.0000, 240.0000,-140.0000 （b） A=1.0000 0.5000 0.3333 0.2500;0.5000 0.3333 0.2500 0.2000; 0.3333 0.2500 0.2000 0.1667;0.2500 0.2000 0.1667 0.1429; b=1;0;0;0;x2= Ab, co

7、nd(A) 运行得（b）的解为： x=18.7308,-149.6053,310.0628, -185.0881 cond(A)=19808,第2章 数值计算方法建模,例3 已知带状稀疏方程组 用稀疏矩阵和满矩阵分别求解，并对运行时间进行比较,第2章 数值计算方法建模,解：在MATLAB指令窗中输入下面指令 A1=sparse(1:500,1:500,1,500,500); A2=sparse(2:500,1:499,1,500,500); A3=sparse(3:500,1:498,1,500,500); A=12*A1-2*A2-2*A2+A3+A3; b=5*ones(500,1); t

8、ic;x=Ab;t1=toc AA=full(A); tic;xx=AAb;t2=toc y=sum(x),yy=sum(xx) 运行得结果为: t1=0.0320，t2=0.2007，y=250.0147，yy=250.0147,第2章 数值计算方法建模,建模示例：种群繁殖问题 问题： 种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群而言，为了保证稳定地收获，各个年龄的种群数量应维持不变.由于种群繁殖主要取决于雌性个体，所以下面种群数量均指其中的雌性. 已知某种群年龄为 ，第 年 年龄种群的数量为 ，繁殖率为 ，自然存活率为 ，对给定收获量 ，建立数学模型，使得各年龄的种群数量维

9、持不变，并就 ， ， 时，求各年龄种群的数量.,第2章 数值计算方法建模,1. 问题的分析 为了保证稳定地收获，需要维持各年龄的种群数量保持不变，根据这个条件可以建立不同种群数量的等式关系 2. 模型的建立 假设此种群的最高年龄为 ，则对第 +1年 年龄种群的数量 ，有 由于各年龄种群繁殖的后代均为年龄1的种群，所以有 要使得各年龄的种群数量维持不变，需满足 ，于是建立下面数学模型：,第2章 数值计算方法建模,令 ， ， 得线性方程组,第2章 数值计算方法建模,求解得x=8481.0, 2892.4, 1335.4, 6012.6，1405.1, 即各年各年龄种群的数量应为81,2892,13

10、35,6012,1405，才能保证种群数量的平衡 课堂练习： 某地区有三个重要产业：一个煤矿，一个发电厂和一条地方铁路开采一元钱的煤，煤矿需要支付0.25元的电费和0.25元的运输费；生产一元钱的电力，发电厂需要支付0.65元的煤费，0.05元的电费和0.05元的运输费；创收一元钱的运输费，铁路需要支付0.55元的煤费和0.10元的电费在某一周期内煤矿接到外地金额为50000元的订货，发电厂接到外地金额为25000元的订货，外地对地方铁路没有需求，问三个企业间一周内总产值多少时才能满足自身及外界需求？三个企业间相互支付多少金额？三个企业各创造多少新价值？,第2章 数值计算方法建模,2.3 插值

11、 插值简介 什么是插值 插值在数学建模中的应用,第2章 数值计算方法建模,插值的MATLAB实现 interp1(x,y,cx,method) interp2（x,y,z,cx,cy, method） griddata（x,y,z,cx,cy, method）,第2章 数值计算方法建模,例1 试用分段线性插值和三次样条插值分别对函数 （ ）进行插值，并与实际曲线图比较 解：在MATLAB指令窗中输入下面指令 ： x0=-5:5;y0=1./(1+x0.2); %产生节点（x0,y0) x=-5:0.1:5;y=1./(1+x.2); %产生插值点x y1=interp1(x0,y0,x); %

12、使用分段线性插值并作图 figure(1),plot(x,y,b,x,y1,k:),grid y2=interp1(x0,y0,x,spline); %使用三次样条插值并作图 figure(2),plot(x,y,b,x,y2,k:),grid 运行得图形见下图,第2章 数值计算方法建模,例2 在某山区（平面区域 内，单位：m）测得一些点的高度（单位：m）见下表，试作出该山区的地貌图,第2章 数值计算方法建模,解：编写M文件shanqu.m如下： x=1200:400:4800;y=1200:400:4800; z=1130 1250 1280 1230 1040 900 500 700 78

13、0 750; 1320 1450 1420 1400 1300 700 900 850 840 380;1390 1500 1500 1400 900 1100 1060 950 870 900; 1500 1200 1100 1350 1450 1200 1150 1010 880 1000;1500 1200 1100 1550 1600 1550 1380 1070 900 1050;1500 1550 1600 1550 1600 1600 1600 1550 1500 1500; 1480 1500 1550 1510 1430 1300 1200 980 850 750; 1450

14、 1470 1320 1280 1200 1080 940 780 620 460;1430 1440 1140 1110 1050 950 820 690 540 380;1400 1410 960 940 880 800 690 570 430 290; figure(1),mesh(x,y,z), x1,y1=meshgrid(1200:50:4800,1200:50:4800); z1=interp2(x,y,z,x1,y1,spline); figure(2),mesh(x1,y1,z1) 运行得得图形见下图。,第2章 数值计算方法建模,第2章 数值计算方法建模,例3 在某水道（平面

15、区域 单位：m）测得一些点的深度，数据见表2-3，已知某船只的吃水线为5米，试画出该水道的海底地貌图及船的禁入区 解：在MATLAB中编写M文件haiyum如下： 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;,第2章 数值计算方法建模,cx,cy=meshgrid(75:5:200,-90:5:150); cz=g

16、riddata(x,y,z,cx,cy,cubic); figure(1),mesh(cx,cy,cz);view(-60,30); figure(2), contour(cx,cy,cz,-5,-5,k) %绘制等高线运行 运行程序得图形见下图,第2章 数值计算方法建模,建模示例：零件加工问题 问题： 待加工零件的外形根据手艺要求由一组数据 给出（在平面情况下），用数控机床加工时刀具必须沿这些数据点前进，由于刀具每次只能沿 方向或 方向走非常小的一步，所以需要将已知数据加密，得到加工所要求的步长很小的 坐标 已知某机翼断面的下轮廓线上的部分数据，见表2-4现在需要得到 坐标 每改变0.1时

17、的坐标试完成加工所需数据，画出曲线，并求 范围内 的最小值,第2章 数值计算方法建模,1. 问题的分析 解决零件加工问题的关键在于计算加密点的纵坐标，此坐标显然可以通过插值得到，而最小值，则可以通过编程逐点搜索获得 2. 模型的建立及其求解 以所给数据为插值节点，以0. 1为步长，在 范围内作三次样条插值，相应的MATLAB指令为： 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),

18、grid 运行可得加密点的纵坐标，从而得下轮廓插值曲线见下图,第2章 数值计算方法建模,在 范围内搜索可得 的最小值为0.9828，相应的 值为13.8.,第2章 数值计算方法建模,课堂练习： 在一丘陵地带测量高程， 和 方向每隔100m测一个点，得高度见表2-7，试作出该丘陵地带的地形图，并求出该地带的最高点和该点的高度,表2-7 某丘陵测量数据表,第2章 数值计算方法建模,2.4 数据拟合 数据拟合简介 最小二乘思想 数据拟合在数学建模中的应用,第2章 数值计算方法建模,数据拟合的MATLAB实现 a=polyfit(x,y,m) x,resnorm,residual=lsqcurvefi

19、t(fun,x0,xdata,ydata) x,resnorm,residual=lsqnonlin(fun,x0,xdata,ydata),第2章 数值计算方法建模,例1 已知一组温度 和电阻 的数据见表2-8, 拟合电阻与温度之间的关系 ，并预测60时的电阻有多大 解：在MATLAB指令窗中输入下面指令 ： t=20.5 32.7 51.0 73.0 95.7;r=765 826 873 942 1032; aa=polyfit(t,r,1);a=aa(1),b=aa(2) y=polyval(aa,t); plot(t,r,k+,t,y,r),xlabel(t),ylabel(R) 运行

20、得a=3.3940,b=702.4918,R(60)= 906.0212，生成图形见下图。,表2-8温度和电阻数据表,第2章 数值计算方法建模,电阻和温度的关系式为,第2章 数值计算方法建模,例2 影响经济增长的主要因素有增加投资、增加劳动力以及技术革新等当科学技术发展不是很快时，经过简化假设与分析，可推导出经济学中著名的Cobb-Douglas生产函数： 其中， 分别表示产值、资金和劳动力， 可由经济统计数据确定现有美国马萨诸塞州1900-1925年上述三个经济指数的统计数据，见表2-9，试用数据拟合的最小二乘法求出上式中的参数。,第2章 数值计算方法建模,表2-9 美国马萨诸塞州1900-

21、1925年三个经济指数的统计数据,第2章 数值计算方法建模,解：（1）建立M函数文件jingji.m如下： function Q=jingji(x,y) Q=x(1)*(y(1,:).x(2).*(y(2,:).(1-x(2); 其中x为待求参数向量，y为已知自变量数据向量。 （2）建立M脚本文件nihe1.m Q=1.05,1.18,1.29,1.30,1.30,1.42,1.50,1.52,1.46,1.60,1.69, 1.81, 1.93 ,1.95, 2.01,2.00,2.09,1.96,2.20,2.12,2.16,2.08,2.24,2.56, 2.34,2.45; y=1.0

22、4,1.06,1.16,1.22,1.27,1.37,1.44,1.53,1.57,2.05,2.51,2.63, 2.74, 2.82, 3.24,3.24,3.61,4.10,4.36,4.77,4.75,4.54,4.58,4.58, 4.58, 4.58;1.05,1.08,1.18,1.22,1.17,1.30,1.39,1.47,1.31,1.43,1.58, 1.59, 1.66,1.68,1.65,1.62,1.86,1.93,1.96,1.95,1.90,1.58,1.67, 1.82, 1.60,1.61; x0=0.1,0.1; x,resnorm=lsqcurvefit

23、(jingji,x0,y,Q) 运行得结果为：x=1.0270,0.2187,0.7813 .于是得经济增长模型：,第2章 数值计算方法建模,建模示例：录像机计数问题 问题： 在老式的录像机上会有一个四位数字的计数器，用于计时一盘标明180分钟的录像带从头转到尾，用时184分钟，计数器读数从“0000”运转到“6061”此外，我们还在不同时间测试了这盘录像带的计数器读数，测试数据见表2-10现在计数器上显示为“4450”，问剩下的录像带还能否再记录1小时长的节目？,第2章 数值计算方法建模,表2-10 一盘录像带的测试数据,图2.11 录像机计数器工作原理示意图,录像带,轮盘,0000,磁头,

24、计数器,第2章 数值计算方法建模,1. 问题的分析 由录像机计数器的工作原理可知，计数器读数不是随着录像带的转动均匀增加，而是先快后慢增加开始时，轮盘是空的，读数为0000随着录像带的不断运动，轮盘的半径不断增加，由于录像带的运动速度（线速度）为常数，所以轮盘的转动越来越慢，从而导致计数器读数的增长也越来越慢在录像的过程中，缠绕在轮盘上的录像带的长度可通过录像的时间获得，另外此长度还与轮盘（包括缠绕的录像带部分）的半径密不可分，由此可以找出计数器读数与录像带转过时间的关系 2. 模型的假设 （1）录像带的线速度是常数 ； （2）计数器读数 与轮盘转的圈数 成正比，即 为比例系数；,第2章 数值

25、计算方法建模,（3）录像带的厚度是均匀的，设为 ，轮盘半径记为 ； （4）录像带各圈松紧均匀； （5）初始时刻 时， 。 3.模型的建立 由假设（1）可得时间 内录像带转过的总长度为： 再由假设（2）（4）可知，,第2章 数值计算方法建模,于是建立录像机计数问题的数学模型如下 简记为 4.模型的求解 t=0 20 40 60 80 100 120 140 160 184; n=0 1141 2019 2760 3413 4004 4545 5051 5525 6061; x0=0.1,0.1;x,resnorm=lsqcurvefit(jishu,x0,n,t) 其中M函数文件jishu.m为

26、： function t=jishu(x,n) t=x(1)*n+x(2)*(n.2);,第2章 数值计算方法建模,运行得录像机计数问题的数学模型： 5.模型的检验 6.模型的应用 =116.4（分钟），已录制时间不足120分钟，所以录完,第2章 数值计算方法建模,课堂练习： 某种医用薄膜允许一种物质的分子穿透它，即允许它从高浓度的溶液向低浓度的溶液扩散.在试制时需确定薄膜被这种分子穿透的能力.测定方法如下：用面积为 的薄膜将容器分成体积分别为 的两部分，在这两部分中分别注满该物质的两种不同浓度的溶液.此时该物质的分子就会从高浓度部分穿过薄膜向低浓度部分扩散.通过单位面积薄膜分子扩散的速度与薄

27、膜两侧溶液的浓度差成正比，比例系数 表征了薄膜被该物质分子穿透的能力，称为渗透率.定时测量容器中薄膜某一侧的溶液浓度值，便能确定其值。 参考数据： ，对容器 部分溶液浓度的测试结果见表2-13.,表2-13,第2章 数值计算方法建模,2.5 数值差分与数值微分 数值差分与数值微分简介 差商 三点公式,第2章 数值计算方法建模,数值微分的MATLAB实现 dx=diff(x) pp=spline(x,y),dy=ppd(pp);dyy=ppval(dy,xx); function dy=ppd(pp) breaks,coefs,m=unmkpp(pp); for i=1:m coefsm(i,:

28、)=polyder(coefs(i,:); end dy=mkpp(breaks,coefsm);,第2章 数值计算方法建模,例1 用三点公式计算 在 1.0,1.2,1.4处的导数值， 的值由表2-15给出 解：在MATLAB指令窗中输入指令： x=1.0,1.1,1.2,1.3,1.4; y=0.2500,0.2268,0.2066,0.1890,0.1736; diff3(x,y) 运行得各点的导数值为： -0.2470，-0.2170，-0.1890，-0.1650，-0.0014,表 2-15,第2章 数值计算方法建模,例2 基于正弦函数 的数据点，利用三点公式和三次样条插值分别求导

29、，并与解析所求得的导数进行比较 解: 编写M脚本文件bijiao.m如下： h=0.1*pi;x=0:h:2*pi;y=sin(x); dy1=diff3(x,y); pp=spline(x,y);dy=ppd(pp);dy2=ppval(dy,x); z=cos(x); error1=norm(dy1-z),error2=norm(dy2-z) plot(x,dy1,k:,x,dy2,r-,x,z,b) 运行得结果为：error1 =0.0666，error2 =0.0025，生成图形见图2.13,第2章 数值计算方法建模,图2.13 三点公式、三次样条插值与解析求导比较图,第2章 数值计算

30、方法建模,建模示例：湖水温度变化问题 问题： 湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡对某个湖的水温进行观测得数据见表2-16 试找出湖水温度变化最大的深度,表2-16 某湖的水温观测数据,第2章 数值计算方法建模,1. 问题的分析 湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度 2. 模型的建立与求解 记湖水的深度为 （m），相应的温度为 （），且有 ，并假定函数 可导 h=0 2.3 4.9 9.

31、1 13.7 18.3 22.9 27.2; T=22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1; hh=0:0.1:27.2; pp=spline(h,T);dT=ppd(pp);dTT=ppval(dT,hh); dTTmax,i=max(abs(dTT),hh(i) plot(hh,dTT, b ,hh(i),dTT(i), r. ),grid on,第2章 数值计算方法建模,运行得导函数绝对值的最大值点为： =11.4，最大值为1.6139,即湖水在深度为11.4m时温度变化最大，如图2.14所示（黑点为温度变化最大的点）,图2.14 湖水温度变化曲线图

32、,第2章 数值计算方法建模,课堂练习： 某居民小区有一个直径10m的圆柱形水塔，每天午夜24时向水塔供水，此后每隔2h记录水位，所得数据见表2-20计算小区在这些时刻每小时的用水量,表2-20 水位记录表,第2章 数值计算方法建模,2.5 数值积分 数值积分简介 数值积分公式 代数精度,第2章 数值计算方法建模,数值微分的MATLAB实现 trapz(x) quad(fun,a,b,tol) quadl(fun,a,b,tol),第2章 数值计算方法建模,例1 分别利用梯形公式、Simpson公式和Gauss-Lobatto法计算 ，并与其精确值比较 解：先对积分作符号运算，然后将其计算结果转

33、换为数值型，再将其与这三种方法求得的数值解比较，其MATLAB指令为： syms xx z0=simple(int(sqrt(1+xx2),0,1) z=double(z0);z=vpa(z,8) x=0:0.01:1;y=sqrt(1+x.2); z1=trapz(y)*0.01;z1=vpa(z1,8),err1=z-z1;err1=vpa(err1,8) z2=quad(sqrt(1+x.2),0,1);z2=vpa(z2,8),err2=z-z2;err2=vpa(err2,8) z3=quadl(sqrt(1+x.2),0,1);z3=vpa(z3,8),err3=z-z3;err3

34、=vpa(err3,8) 运行得精确值为 1.1477936，三种公式计算得数值积 分值分别为1.1477995，1.1477935和1.1477936，其相应误差分别为-.59e- 5，.1e-6和0.,第2章 数值计算方法建模,例2 人造地球卫星轨道可视为平面上的椭圆我国第一颗人造地球卫星近地点距地球表面439km，远地点距地球表面2384km，地球半径为6371km，求该卫星的轨道长度 解:卫星轨道椭圆的参数方程为 分别是长、短半轴，则椭圆的长度为 计算其数值解，编写M函数文件如下： function y=y(t) a=8755;b=6810; y=4*sqrt(a2*sin(t).2+b2*cos(t).2); 在MATLAB指令窗中输入以下