数学建模算法章节

1、武汉科技大学武汉科技大学武汉科技大学武汉科技大学数学建模概述数学建模概述 数值计算方法建模数值计算方法建模 基本方法建模基本方法建模 数学规划方法建模数学规划方法建模 统计分析方法建模统计分析方法建模第第1 1章章第第2 2章章第第3 3章章第第4 4章章第第5 5章章目目 录录图论方法建模图论方法建模第第6 6章章武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 2.1 非线性方程求解非线性方程求解 非线性方程简介非线性方程简介 u 次代数方程次代数方程 u 超越方程超越方程 n00111axaxaxannnn0)2sin(xex武汉科技大学武汉科技大学第第2章章 数值

2、计算方法建模数值计算方法建模 非线性方程求解的非线性方程求解的MATLAB实现实现 x,fval,exitflag,output=fzero(f ,x0,options) r=roots(c) x,fval,exitflag,output=fsolve(f,x0,options)武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例1 在无阻尼强迫振荡的研究中会经常遇到函数 .试求一点 ，满足 。 解：建立M函数文件h.m： function y=h(x) y=x*sin(x)-1; 在MATLAB指令窗中输入下面指令 x0=0,2;x,fval,exitflag=fze

3、ro(h,x0) 运行得结果为：x =1.1142，fval =2.2204e-016，exitflag =1，即所求非线性方程的解为1.1142.xxxhsin)(2 , 0 x1)(xh武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例2 求函数 的零点. 解：（1）为确定其零点的大体位置，先做出它的图形 ； （2）将图形放大得5个零点，利用ginput指令取其坐标； （3） 利用指令fzero 计算其相应精确坐标20.1( )(sin)0.5tf tt et图图2.1 图形放大法图形放大法 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v

4、例例3 求解非线性方程组v 解：建立解：建立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,)1422212221xxxx武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 建模示例：贷款问题建模示例：贷款问题u 问题：问题： 某人从银行贷款购房，若他今年初贷款某人从银行贷款购房，若他今年初贷款10万元，月利率万元，月利率0.5%，每月还每月还1000元，试计算他每年

5、末欠银行多少钱，多少时间才能还元，试计算他每年末欠银行多少钱，多少时间才能还清？如果要清？如果要10年还清，每月需还多少？年还清，每月需还多少？ 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 1. 问题的分析问题的分析 每月的还款金额应包括本金和利息两部分每月的还款金额应包括本金和利息两部分 。 2. 模型的建立模型的建立 记第记第 个月初此人欠银行个月初此人欠银行 元，月利率为元，月利率为 ，每月还款，每月还款 元，元，则则 建立数学模型如下：建立数学模型如下： ttxpQ100000, 2 , 1),(0111xntpxQxxxtttn 将上式依次递推可得将上式依次

6、递推可得 QppxpQQpxpQxpxnnnnn1)1 ()1 ( )1 ()1 ()1 (1121武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 令令 得得 01nxppxpQnn1)1 ()1 (1将所给数据代入上式可得将所给数据代入上式可得Q=1110.3（元） 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模& 课堂练习：课堂练习： 在一条在一条20m宽的道路两侧，分别安装了一只宽的道路两侧，分别安装了一只2kW和一只和一只3kW的路的路灯，它们离地面的高度分别为灯，它们离地面的高度分别为5m和和6m（如下图）（如下图）.在漆黑的夜晚，

7、当在漆黑的夜晚，当两只路灯开启时，试求两只路灯连线的路面上最暗和最亮的点两只路灯开启时，试求两只路灯连线的路面上最暗和最亮的点.（提（提 示：光源在示：光源在 点的照度为点的照度为 ，其中，其中 为光源的功率，为光源的功率， 为为 光源到光源到 的光线与水平面的夹角，为光源到的光线与水平面的夹角，为光源到 点的距离点的距离.） Q2sinrPkIPQrQxsxOy1P2PQ1r2r1h2h12武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 2.2 线性方程组的数值解法线性方程组的数值解法 线性方程组简介线性方程组简介 u 阶线性方程组阶线性方程组 u 矩阵形式矩阵形式

8、n11 11221121 1222221 122 nnnnnnnnnna xa xa xba xa xa xba xa xa xbbAx 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 线性方程组求解的线性方程组求解的MATLAB实现实现 x=Ab A=sparse(r,c,v,m,n) full(A) 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例1 求解下列线性方程组求解下列线性方程组 解：在解：在MATLAB指令窗中输入下面指令指令窗中输入下面指令 ： A=1 1 0 0;2 -1 5 0;0 3 -4 2;0 0 2 -6;b=5;-

9、9;19;2; x=Ab 运行得结果为：运行得结果为： x=-2.0000,7.0000,0.4000,-0.2000262 19243 -9 52x5 4343232121xxxxxxxxx武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例2 验证验证Hilbert矩阵是一个典型的病态矩阵矩阵是一个典型的病态矩阵 （a）用）用4阶阶Hilbert矩阵求解矩阵求解 的精确解（用分数表示所有的元素的精确解（用分数表示所有的元素并进行精确计算）：并进行精确计算）： （b）使用精度为）使用精度为4位有效数字的算术计算求解位有效数字的算术计算求解 AXb71615141615

10、14131514131214131211A1000b AXb武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模1429. 01667. 02000. 02500. 01667. 02000. 02500. 03333. 02000. 02500. 03333. 05000. 02500. 03333. 05000. 00000. 1A1000b v 解：（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.

11、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, cond(A) 运行得（b）的解为： x=18.7308,-149.6053,310.0628, -185.0881 cond(A)=19808武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例3 已知带状稀疏方程组 用稀疏矩阵和满矩阵分别求解，并

12、对运行时间进行比较5 122 5 2122 52122. . .5 2122 5 2122 5 21225 21250049949850049949849750049949849749665432543214321321xxxxxxxxxxxxxxxxxxxxxxxxxxxxx武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 解：在解：在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,

13、500,500); A=12*A1-2*A2-2*A2+A3+A3; b=5*ones(500,1); tic;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章章 数值计算方法建模数值计算方法建模 建模示例：种群繁殖问题建模示例：种群繁殖问题u 问题问题： 种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群而言，为

14、了保证稳定地收获，各个年龄的种群数量应维持不的种群而言，为了保证稳定地收获，各个年龄的种群数量应维持不变变.由于种群繁殖主要取决于雌性个体，所以下面种群数量均指其中由于种群繁殖主要取决于雌性个体，所以下面种群数量均指其中的雌性的雌性. 已知某种群年龄为已知某种群年龄为 ，第，第 年年 年龄种群的数量为年龄种群的数量为 ，繁殖率为繁殖率为 ，自然存活率为，自然存活率为 ，对给定收获量，对给定收获量 ，建立数学模型，建立数学模型，使得各年龄的种群数量维持不变，并就使得各年龄的种群数量维持不变，并就 ， ， 时，求各年龄种群的数量时，求各年龄种群的数量. nk, 1 tktkxkbkskh5n,20

15、0,400,500, 6 . 0, 4 . 0, 3, 5, 0321324143521hhhssssbbbbb100,10054hh武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 1. 问题的分析问题的分析 为了保证稳定地收获，需要维持各年龄的种群数量保持不变，根据为了保证稳定地收获，需要维持各年龄的种群数量保持不变，根据这个条件可以建立不同种群数量的等式关系这个条件可以建立不同种群数量的等式关系 2. 模型的建立模型的建立 假设此种群的最高年龄为假设此种群的最高年龄为 ，则对第，则对第 +1年年 年龄种群的数量年龄种群的数量 ，有，有 由于各年龄种群繁殖的后代均为年龄

16、由于各年龄种群繁殖的后代均为年龄1的种群，所以有的种群，所以有 要使得各年龄的种群数量维持不变，需满足要使得各年龄的种群数量维持不变，需满足 ，于是建立，于是建立下面数学模型：下面数学模型： ntkktx, 11, 2 , 1, 2 , 1,1, 1nkthxsxktkkktnktkktxbx11 , 1kttkxx, 1武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 令令 ， ， 得线性方程组得线性方程组 111112211nnnnnkkkhxsxhxsxxbx1-00 s00001000111 -n211321nnbssbbbbAnnxxxxx1211210nnhhh

17、ccAx 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 求解得求解得x=8481.0, 2892.4, 1335.4, 6012.6，1405.1, 即各年各年即各年各年龄种群的数量应为龄种群的数量应为81,2892,1335,6012,1405，才能保证种群数量的平，才能保证种群数量的平衡衡& 课堂练习：课堂练习： 某地区有三个重要产业：一个煤矿，一个发电厂和一条地方铁某地区有三个重要产业：一个煤矿，一个发电厂和一条地方铁路开采一元钱的煤，煤矿需要支付路开采一元钱的煤，煤矿需要支付0.25元的电费和元的电费和0.25元的运输费；元的运输费；生产一元钱的电力，发

18、电厂需要支付生产一元钱的电力，发电厂需要支付0.65元的煤费，元的煤费，0.05元的电费和元的电费和0.05元的运输费；创收一元钱的运输费，铁路需要支付元的运输费；创收一元钱的运输费，铁路需要支付0.55元的煤费元的煤费和和0.10元的电费在某一周期内煤矿接到外地金额为元的电费在某一周期内煤矿接到外地金额为50000元的订货，元的订货，发电厂接到外地金额为发电厂接到外地金额为25000元的订货，外地对地方铁路没有需求，元的订货，外地对地方铁路没有需求，问三个企业间一周内总产值多少时才能满足自身及外界需求？三个企问三个企业间一周内总产值多少时才能满足自身及外界需求？三个企业间相互支付多少金额？三

19、个企业各创造多少新价值？业间相互支付多少金额？三个企业各创造多少新价值？ 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 2.3 插值插值 插值简介插值简介 u 什么是插值什么是插值u 插值在数学建模中的应用插值在数学建模中的应用 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 插值的插值的MATLAB实现实现 interp1(x,y,cx,method) interp2（x,y,z,cx,cy, method） griddata（x,y,z,cx,cy, method） 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例

20、1 试用分段线性插值和三次样条插值分别对函数试用分段线性插值和三次样条插值分别对函数 （ ）进行插值，并与实际曲线图比较进行插值，并与实际曲线图比较v 解：在解：在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); %使用分段线性插值并作图使用分段线性插值并作图 figure(1),plot(x,y,b,x,y1,k:),grid y2=interp1(x0,y0,x,spline); %

21、使用三次样条插值并作图使用三次样条插值并作图 figure(2),plot(x,y,b,x,y2,k:),grid 运行得图形见下图运行得图形见下图 211xy55x武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例2 在某山区（平面区域在某山区（平面区域 内，单位：内，单位：m）测得一些点的高度（单位：测得一些点的高度（单位：m）见下表，试作出该山区的地貌图）见下表，试作出该山区的地貌图 8004800, 04800 xyxy12001200160016002000200024002400280028003200320036003600400040004400440

22、048004800120012001130113012501250128012801230123010401040900900500500700700780780750750160016001320132014501450142014201400140013001300700700900900850850840840380380200020001390139015001500150015001400140090090011001100106010609509508708709009002400240015001500120012001100110013501350145014501200120

23、011501150101010108808801000100028002800150015001200120011001100155015501600160015501550138013801070107090090010501050320032001500150015501550160016001550155016001600160016001600160015501550150015001500150036003600148014801500150015501550151015101430143013001300120012009809808508507507504000400014501

24、4501470147013201320128012801200120010801080940940780780620620460460440044001430143014401440114011401110111010501050950950820820690690540540380380480048001400140014101410960960940940880880800800690690570570430430290290武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 解：编写解：编写M文件文件shanqu.m如下：如下： x=1200:400:4800;y=1

25、200:400:4800; z=1130 1250 1280 1230 1040 900 500 700 780 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 150

26、0; 1480 1500 1550 1510 1430 1300 1200 980 850 750; 1450 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

27、,y1,z1) 运行得得图形见下图。运行得得图形见下图。 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例3 在某水道（平面区域在某水道（平面区域 单位：单位：m）测得一些）测得一些点的深度，数据见表点的深度，数据见表2-3，已知某船只的吃水线为，已知某船只的吃水线为5米，试画出该水道的米，试画出该水道的海底地貌图及船的禁入区海底地貌图及船的禁入区 v 解解：在MATLAB中编写M文件haiyum如下： x=129 140 103.5 88 185.5 195 105 157.5 107.5 77

28、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;15090,20075yxxyz129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.57.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5 4 8 6 8 6 8 8 9 9 8 8 9 4 9武汉科技大学武汉科技大学第第2

29、章章 数值计算方法建模数值计算方法建模 cx,cy=meshgrid(75:5:200,-90:5:150); cz=griddata(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章章 数值计算方法建模数值计算方法建模 建模示例：建模示例：零件加工问题零件加工问题 u 问题：问题： 待加工零件的外形根据手艺要求由一组数据待加工零件的外形根据手艺要求

30、由一组数据 给出（在平面情给出（在平面情况下），用数控机床加工时刀具必须沿这些数据点前进，由于刀具每次况下），用数控机床加工时刀具必须沿这些数据点前进，由于刀具每次只能沿只能沿 方向或方向或 方向走非常小的一步，所以需要将已知数据加密，得方向走非常小的一步，所以需要将已知数据加密，得到加工所要求的步长很小的到加工所要求的步长很小的 坐标坐标 已知某机翼断面的下轮廓线上的部分数据，见表已知某机翼断面的下轮廓线上的部分数据，见表2-4现在需要得到现在需要得到 坐标坐标 每改变每改变0.1时时 的坐标试完成加工所需数据，画出曲线，并求的坐标试完成加工所需数据，画出曲线，并求 范围内范围内 的最小值的

31、最小值 ( , )x yxy( , )x yxy1513 xy 0 3 5 7 9 11 12 13 14 15 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 1. 问题的分析问题的分析 解决零件加工问题的关键在于计算加密点的纵坐标，此坐标显然可以解决零件加工问题的关键在于计算加密点的纵坐标，此坐标显然可以通过插值得到，而最小值，则可以通过编程逐点搜索获得通过插值得到，而最小值，则可以通过编程逐点搜索获得 2. 模型的建立及其求解模型的建立及其求解 以所给数据为插值节点，以以所给数据为插值节点，以

32、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),grid 运行可得加密点的纵坐标，从而得下轮廓插值曲线见下图运行可得加密点的纵坐标，从而得下轮廓插值曲线见下图 150 x武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模在在 范围内搜索可得范围内搜索

33、可得 的最小值为的最小值为0.9828，相应的，相应的 值为值为13.8. 1513 xyx武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模&课堂练习：课堂练习： 在一丘陵地带测量高程，在一丘陵地带测量高程， 和和 方向每隔方向每隔100m测一个点，得高度测一个点，得高度见表见表2-7，试作出该丘陵地带的地形图，并求出该地带的最高点和该，试作出该丘陵地带的地形图，并求出该地带的最高点和该点的高度点的高度 xy表表2-7 某丘陵测量数据表某丘陵测量数据表xy10010020020030030040040010010063663669769762462447847820

34、0200698698712712630630478478300300680680674674598598412412400400662662626626552552334334武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 2.4 数据拟合数据拟合 数据拟合简介数据拟合简介 u 最小二乘思想最小二乘思想u 数据拟合在数学建模中的应用数据拟合在数学建模中的应用 niiiniimyxfccc121210)(),(min武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 数据拟合的数据拟合的MATLAB实现实现 a=polyfit(x,y,m) x,resn

35、orm,residual=lsqcurvefit(fun,x0,xdata,ydata) x,resnorm,residual=lsqnonlin(fun,x0,xdata,ydata) 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例1 已知一组温度已知一组温度 和电阻和电阻 的数据见表的数据见表2-8, 拟合电阻与温度之间的关拟合电阻与温度之间的关系系 ，并预测，并预测60时的电阻有多大时的电阻有多大 v 解：在解：在MATLAB指令窗中输入下面指令指令窗中输入下面指令 ： t=20.5 32.7 51.0 73.0 95.7;r=765 826 873 942

36、 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) 运行得运行得a=3.3940,b=702.4918,R(60)= 906.0212，生成图形见下图。，生成图形见下图。 tRRatb表表2-8温度和电阻数据表温度和电阻数据表tR（ ） 20.5 32.7 51.0 73.0 95.7 （ ） 765 826 873 942 1032 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模电阻和温度的关系式为电阻和温度的关系式为3.3940

37、702.4918Rt武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例2 影响经济增长的主要因素有增加投资、增加劳动力以及技术革新影响经济增长的主要因素有增加投资、增加劳动力以及技术革新等当科学技术发展不是很快时，经过简化假设与分析，可推导出经济学等当科学技术发展不是很快时，经过简化假设与分析，可推导出经济学中著名的中著名的Cobb-Douglas生产函数：生产函数： 其中，其中， 分别表示产值、资金和劳动力，分别表示产值、资金和劳动力， 可由经济统计数据确可由经济统计数据确定现有美国马萨诸塞州定现有美国马萨诸塞州1900-1925年上述三个经济指数的统计数据，见年上

38、述三个经济指数的统计数据，见表表2-9，试用数据拟合的最小二乘法求出上式中的参数。，试用数据拟合的最小二乘法求出上式中的参数。 (, ),0,1Q K LaK L 1,Q K L, ,a 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模表表2-9 美国马萨诸塞州美国马萨诸塞州1900-1925年三个经济指数的统计数据年三个经济指数的统计数据 QQKKLL1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 19121.05 1.18 1.29 1.30 1.30 1.42 1.50 1.52 1.46 1.6

39、0 1.69 1.81 1.931.04 1.06 1.16 1.22 1.27 1.37 1.44 1.53 1.57 2.05 2.51 2.63 2.741.05 1.08 1.18 1.22 1.17 1.30 1.39 1.47 1.31 1.43 1.58 1.59 1.661913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 19251.95 2.01 2.00 2.09 1.96 2.20 2.12 2.16 2.08 2.24 2.56 2.34 2.452.82 3.24 3.24 3.61 4.10 4.3

40、6 4.77 4.75 4.54 4.58 4.58 4.58 4.581.68 1.65 1.62 1.86 1.93 1.96 1.95 1.90 1.58 1.67 1.82 1.60 1.61tt武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 解：（解：（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脚本文件脚本文件ni

41、he1.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.04,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

42、,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(jingji,x0,y,Q) 运行得结果为：运行得结果为：x=1.0270,0.2187,0.7813 .于是得经济增长模型：于是得经济增长模型：武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模7813. 02187. 00270. 1),(LKLKQ 建模示例：录像机计数问题建

43、模示例：录像机计数问题 问题：问题： 在老式的录像机上会有一个四位数字的计数器，用于计时一盘标在老式的录像机上会有一个四位数字的计数器，用于计时一盘标明明180分钟的录像带从头转到尾，用时分钟的录像带从头转到尾，用时184分钟，计数器读数从分钟，计数器读数从“0000”运转到运转到“6061”此外，我们还在不同时间测试了这盘录像带的计数器此外，我们还在不同时间测试了这盘录像带的计数器读数，测试数据见表读数，测试数据见表2-10现在计数器上显示为现在计数器上显示为“4450”，问剩下的录，问剩下的录像带还能否再记录像带还能否再记录1小时长的节目？小时长的节目？ 武汉科技大学武汉科技大学第第2章章

44、 数值计算方法建模数值计算方法建模表表2-10 一盘录像带的测试数据一盘录像带的测试数据tntn（分钟）（分钟）0 0101020203030404050506060707080809090（读数）（读数）00000000061706171141114116011601201920192403240327602760309630963413341337153715（分钟）（分钟）100100110110120120130130140140150150160160170170184184（读数）（读数）40044004428042804545454548034803505150515291529

45、1552555255752575260616061图图2.11 录像机计数器工作原理示意图录像机计数器工作原理示意图录像带录像带轮盘轮盘0000磁头磁头计数器计数器武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 1. 问题的分析问题的分析 由录像机计数器的工作原理可知，计数器读数不是随着录像带的转动由录像机计数器的工作原理可知，计数器读数不是随着录像带的转动均匀增加，而是先快后慢增加开始时，轮盘是空的，读数为均匀增加，而是先快后慢增加开始时，轮盘是空的，读数为0000随随着录像带的不断运动，轮盘的半径不断增加，由于录像带的运动速度着录像带的不断运动，轮盘的半径不断增加，由

46、于录像带的运动速度（线速度）为常数，所以轮盘的转动越来越慢，从而导致计数器读数的（线速度）为常数，所以轮盘的转动越来越慢，从而导致计数器读数的增长也越来越慢在录像的过程中，缠绕在轮盘上的录像带的长度可通增长也越来越慢在录像的过程中，缠绕在轮盘上的录像带的长度可通过录像的时间获得，另外此长度还与轮盘（包括缠绕的录像带部分）的过录像的时间获得，另外此长度还与轮盘（包括缠绕的录像带部分）的半径密不可分，由此可以找出计数器读数与录像带转过时间的关系半径密不可分，由此可以找出计数器读数与录像带转过时间的关系 2. 模型的假设模型的假设 （1）录像带的线速度是常数）录像带的线速度是常数 ； （2）计数器读

47、数）计数器读数 与轮盘转的圈数与轮盘转的圈数 成正比，即成正比，即 为比例为比例系数；系数；vmnkknm,武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 （3）录像带的厚度是均匀的，设为）录像带的厚度是均匀的，设为 ，轮盘半径记为，轮盘半径记为 ； （4）录像带各圈松紧均匀；）录像带各圈松紧均匀； （5）初始时刻）初始时刻 时，时， 。 3.模型的建立模型的建立 由假设（由假设（1）可得时间）可得时间 内录像带转过的总长度为：内录像带转过的总长度为： 再由假设（再由假设（2）（）（4）可知，）可知， dr0t0ntvtL miidrL1)(2武汉科技大学武汉科技大学第第

48、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为：为： function t=jishu(x,n) t=x(1)*n+x(2)*(n.2);222nvdknvdkrkt2bnant武

49、汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 运行得录像机计数问题的数学模型：运行得录像机计数问题的数学模型： 5.模型的检验模型的检验 6.模型的应用模型的应用 =116.4（分钟），已录制时间不足（分钟），已录制时间不足120分钟，所以录完分钟，所以录完262106113. 2104530. 1nntt武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模&课堂练习：课堂练习： 某种医用薄膜允许一种物质的分子穿透它，即允许它从高浓度的溶液某种医用薄膜允许一种物质的分子穿透它，即允许它从高浓度的溶液向低浓度的溶液扩散向低浓度的溶液扩散.在试制时需确

50、定薄膜被这种分子穿透的能力在试制时需确定薄膜被这种分子穿透的能力.测定方测定方法如下：用面积为法如下：用面积为 的薄膜将容器分成体积分别为的薄膜将容器分成体积分别为 的两部分，在这的两部分，在这两部分中分别注满该物质的两种不同浓度的溶液两部分中分别注满该物质的两种不同浓度的溶液.此时该物质的分子就会此时该物质的分子就会从高浓度部分穿过薄膜向低浓度部分扩散从高浓度部分穿过薄膜向低浓度部分扩散.通过单位面积薄膜分子扩散的通过单位面积薄膜分子扩散的速度与薄膜两侧溶液的浓度差成正比，比例系数速度与薄膜两侧溶液的浓度差成正比，比例系数 表征了薄膜被该物质表征了薄膜被该物质分子穿透的能力，称为渗透率分子穿

51、透的能力，称为渗透率.定时测量容器中薄膜某一侧的溶液浓度值，定时测量容器中薄膜某一侧的溶液浓度值，便能确定其值。便能确定其值。 参考数据：参考数据： ，对容器，对容器 部分溶液浓度的测试部分溶液浓度的测试结果见表结果见表2-13. S,ABV VK321000cm ,10cmABVVSB表表2-13 时时 间间10010020020030030040040050050060060070070080080090090010001000浓度浓度mg/cm3（10-5）4.544.544.994.995.355.355.655.655.905.906.106.106.266.266.396.396.

52、506.506.596.59武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 2.5 数值差分与数值微分数值差分与数值微分 数值差分与数值微分简介数值差分与数值微分简介 u 差商差商u 三点公式三点公式 hxfhxfxf)()()(*000()()()2f xhf xhfxhhyyyxf243)(2100hyyyxfnnnn234)(12武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 数值微分的数值微分的MATLAB实现实现 dx=diff(x) pp=spline(x,y),dy=ppd(pp);dyy=ppval(dy,xx); function

53、 dy=ppd(pp) breaks,coefs,m=unmkpp(pp); for i=1:m coefsm(i,:)=polyder(coefs(i,:); end dy=mkpp(breaks,coefsm);武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例1 用三点公式计算用三点公式计算 在在 1.0,1.2,1.4处的导数值，处的导数值， 的值由表的值由表2-15给出给出 v 解：在解：在MATLAB指令窗中输入指令：指令窗中输入指令： x=1.0,1.1,1.2,1.3,1.4; y=0.2500,0.2268,0.2066,0.1890,0.1736;

54、 diff3(x,y) 运行得各点的导数值为：运行得各点的导数值为： -0.2470，-0.2170，-0.1890，-0.1650，-0.0014 ( )yf xx( )f x表表 2-15 x)(xf1.0 1.1 1.2 1.3 1.4 1.0 1.1 1.2 1.3 1.4 0.2500 0.2268 0.2066 0.1890 0.17360.2500 0.2268 0.2066 0.1890 0.1736武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 例例2 基于正弦函数基于正弦函数 的数据点，利用三点公式和三次样条插的数据点，利用三点公式和三次样条插值分别

55、求导，并与解析所求得的导数进行比较值分别求导，并与解析所求得的导数进行比较 v 解解: 编写编写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

56、 sinyx武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模图图2.13 三点公式、三次样条插值与解析求导比较图三点公式、三次样条插值与解析求导比较图武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 建模示例：湖水温度变化问题建模示例：湖水温度变化问题 问题：问题： 湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低这种现象会影响水的对流和混合过程，使得下层越往下水的温度越低这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡对某个湖的水温进行观测得数据见

57、表水域缺氧，导致水生鱼类死亡对某个湖的水温进行观测得数据见表2-16 试找出湖水温度变化最大的深度试找出湖水温度变化最大的深度 表表2-16 某湖的水温观测数据某湖的水温观测数据深度（深度（m m）0 02.32.34.94.99.19.113.713.718.318.322.922.927.227.2温度（温度（）22.822.822.822.822.822.820.620.613.913.911.711.711.111.111.111.1武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 1. 问题的分析问题的分析 湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便

58、转湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度温度变化最大的深度 2. 模型的建立与求解模型的建立与求解 记湖水的深度为记湖水的深度为 （m），相应的温度为），相应的温度为 （），且有），且有 ，并，并假定函数假定函数 可导可导 h=0 2.3 4.9 9.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,

59、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 onhT)(hTT )(hT武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模 运行得导函数绝对值的最大值点为：运行得导函数绝对值的最大值点为： =11.4，最大值为，最大值为1.6139,即湖即湖水在深度为水在深度为11.4m时温度变化最大，如图时温度变化最大，如图2.14所示（黑点为温度变化最大所示（黑点为温度变化最大的点）的点） h图图2.14 湖水温度变化曲线图湖水温度变

60、化曲线图 武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模&课堂练习：课堂练习： 某居民小区有一个直径某居民小区有一个直径10m的圆柱形水塔，每天午夜的圆柱形水塔，每天午夜24时向水塔供水，此后时向水塔供水，此后每隔每隔2h记录水位，所得数据见表记录水位，所得数据见表2-20计算小区在这些时刻每小时的用水量计算小区在这些时刻每小时的用水量 表表2-20 水位记录表水位记录表时刻时刻/h24681012141618202224水位水位/cm305298290265246225207189165148130114武汉科技大学武汉科技大学第第2章章 数值计算方法建模数值计算方法建模v 2.5 数值积分数值积