第 章 实例及数学模型一般地说,为了定量地解决一个实际问题,从中抽象,归纳出来的数学表述就是数学模型.数学模型可以描述为,对现实世界的一个研究对象,为了一个特定目的,做出必须的简化假设,根据对象的内在规律,运用适当的数学工具,得到的一个数学表述.而数学建模包括模型的建立,求解,分析和检验的全过程.从实际问题到数学模型,又从数学模型的求解结果回到现实对象.数学建模面临的实际问题多种多样,建模的目的不同,分析的方法不同,采用的数学工具不同,所得模型的类型也不同.1 初等数学模型实例1 商品市场占有率问题有R 和S 两家公司经营同类产品,这两家公司相互竞争.每年R 公司保持1/4的顾客,而3/4转移向S 公司;每年S 公司保持有2/3的顾客,而1/3转移向R 公司.当产品开始制造时R 公司占有3/5的市场份额,而S 公司占有2/5的市场份额.试问两年以后,两家公司所占有的市场份额变化怎样?5年以后会怎样?10年以后呢?是否有一组初始市场份额分配数据使以后每年的市场分配成为稳定不变?一 问题及分析模型建立根据两家公司每年顾客转移的数据资料,形成下面的转移矩阵⎪⎪⎪⎪⎭⎫⎝⎛=32433141A 又根据产品制造之初市场的初始分配数据可得如下向量⎪⎪⎪⎪⎭⎫ ⎝⎛=52530X所以一年后,市场分配为⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎭⎫⎝⎛==52533243314101AX X 两年以后,市场分配为022X A X =设n 年后市场分配的份额为n X ,则有01X A AX X nn n ==-设数据b a ,为R 公司和S 公司的初始市场份额,则有1=+b a为了使以后每年的市场份额分配不变,根据顾客转移的规律,有⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎭⎫⎝⎛b a b a 32433141 即⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎭⎫⎝⎛--0031433143b a 这是一个齐次方程组问题,如果方程组有解,则应该在非零解的集合中选取正数解作为市场份额稳定的初始份额.由上面的分析得该问题的数学模型为求两个线性方程组,即0X A X n n =⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭⎫⎝⎛⎪⎪⎪⎪⎭⎫⎝⎛--0031433143b a 二 模型的求解可以用[x1,x2]=solve(s1,s2,v1,v2)来求方程组0=AX 的解. 也可以用命令rref(A ),化A 为上三角阵,再求解. 计算程序为A=[1/4 1/3;3/4 2/3]; X0=[3/5;2/5]; X2=A^2*X0 X5=A^5*X0 X10=A^10*X0 运算结果为X2 =0.3097 0.6903X5 =0.3077 0.6923 X10 =0.3077 0.6923为了求a 和b 作为R 公司和S 公司稳定的初始市场份额,用命令rref 来求解齐次方程组 计算程序为format rat; rref(A-eye(2)) 运算结果为 ans =1 -4/90 0 由此得化简后的方程为094=-b a结合约束条件1=+b a可以得到134=a 139=b 这就是使市场稳定的两家公司的初始份额.2 微积分方法模型实例问题分析及模型建立模型的求解3 微分方程模型本节以实例分析并建立微分方程模型,对模型作了解析解以及MA TLAB 数值解.而当函数以离散数据形式表示时,函数的数值微分就得借助差分来计算,差分是微分的近似.故本节还分析了简单的差分方程模型.实例1 温度冷却由物理学知道,物体冷却的速度与当时的物体温度和周围环境温度之差成正比.今100℃的沸水注入杯里,放在室温为20℃的环境冷却,5min 后测得水温为60℃.求水温u 与时间t 的函数关系.一 问题分析及模型的建立设比例系数为)0(>k k ,根据题意可得微分方程)20(--=u k dtdu60,10050====t t u u二 模型的求解此为简单的一阶可分离变量微分方程,可得解析解5)5.0(8020t u +=另外,还可用MA TLAB 程序求其解析解和数值解. 解析解的程序为dsolve('Du+k*(u-20)=0','u(0)=100','t') %dsolve 为求常微分方程的符号解函数运算结果为u =20+80*exp(-k*t)再由605==t u ,可得52ln =k ,即5)5.0(8020tu +=数值解的程序为f=inline('-0.2*log(2)*(u-20)','t','u');[t,u]=ode45(f,[0, 100],100); %ode45为龙格库塔法求微分方程的数值解plot(t,u) %绘制0到100分钟的温度随时间变化的图形图 温度随时间变化从图可看出温度随时间逐渐趋于20℃.实例2 动物种群的相互竞争与相互依存的模型在生物的种群关系中,一种生物以另一种生物为食的现象,称为捕食.一般说来,由于捕食关系,当捕食动物数量增长时,被捕食动物数量就逐渐下降,捕食动物由于食物来源短缺,数量也随之下降,而被捕食动物数量却随之上升.这样周而复始,捕食动物与被捕食动物的数量随时间变化形成周期性的震荡.田鼠及其天敌的田间种群消长动态规律也是如此.实验调查数据表明:无论是田鼠还是其天敌的数量都呈周期性的变化,天鼠与天敌的作用系统随时间序列推移,田鼠密度逐渐增加,其天敌随之增加,但时间上落后一步.由于天敌密度增加,则田鼠密度降低,而田鼠密度的降低,则其天敌密度亦减少,如此往复循环,从而形成一定的周期.试用数学模型来概括这一现象,并总结出其数量变化的近似公式.一 问题分析及模型的建立设)(t x 和)(t y 分别表示t 时刻田鼠与其天敌的数量,如果单独生活,田鼠的增长速度正比于当时的数量,即x dt dxλ= 而田鼠的天敌由于没有被捕食对象,其数量减少的速率正比于当时的数量,即y dtdyμ-= 现在田鼠与其天敌生活一起,田鼠一部分遭到其天敌的消灭,于是以一定的速率α减少,减少的数量正比于天敌的数量,因此有x y dtdx)(αλ-= 类似地,田鼠的天敌有了食物,数量减少的速率μ减少β,减少的量正比于田鼠的数量,因此有y x dtdy)(βμ--= 上述公式,最后两个方程联合起来称为V olterra-Lot 方程,这里μλβα,,,均为正数,初始条件为00)0(,)0(y y x x ==现在通过实验调查所得到的数据如表,此数据为每隔两个月田间调查一次,得到的田鼠及其天敌种群数量的记录,数量的单位经过处理.试建立合理的数学模型. 表 田鼠种群数量记录29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5表 田鼠天敌种群数量记录1.6 1.3 1.1 1.2 1.1 1.3 1.82.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9 1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3二 模型的求解V olterra-Lotok 方程的解析解即y x ,的显示解难求出,因此公式的参数方程不宜直接用Matlab 函数来拟合解,可用如下的方法来求其近似解.V olterra-Lotok 可转化为⎩⎨⎧+-=-=dtx y d dty x d )(ln )(ln βμαλ 在区间],[1i i t t -上积分,得i i i i i S t t x x 111)(ln ln αλ--=--- i i i i i S t t y y 211)(ln ln βμ+--=---这里,⎰-=ii t t i ydt S 11,⎰-=ii t t i xdt S 22, m i ,,1 =于是得到方程组⎩⎨⎧==222111B P A B P A这里⎪⎪⎪⎪⎪⎭⎫⎝⎛------=-im m m S t t S t t S t t A 1121211011 ⎪⎪⎪⎪⎪⎭⎫⎝⎛------=-m m m S t t S t t S t t A 212212012 ⎪⎪⎭⎫ ⎝⎛=αλ1P ⎪⎪⎭⎫⎝⎛-=βμ2PT m m x x x x B )ln ,,(ln1011-= T m m y y y yB )ln ,,(ln 101-= 因此方程组参数的最小二乘解为 111111)(B A A A P T T-= 22122)(B A A A P TT -=由于)(t x 和)(t y 均为未知,因此21,S S i 用数值积分方法的梯形公式解 )(21111--+-≈=⎰-i i i i t t i y y t t ydt S ii )(21121--+-==⎰-i i i i t t x x t t xdt S ii 这样就可求得参数的近似值.模型参数求解的程序为 clear all,clcX=[29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 ...34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 ... 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5];Y=[1.6 1.3 1.1 1.2 1.1 1.3 1.8 2.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9 ... 1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3];N=[X;Y]; T=[0:2:60]; for i=1:30A(i,1)=T(i+1)-T(i);A(i,[2 3])=((T(i+1)-T(i))/2)*[-(N(1,i+1)+N(1,i)),-(N(2,i+1)+N(2,i))]; B(i,[1 2])=[log(N(1,i+1)/N(1,i)),log(N(2,i+1)/N(2,i))]; end;A1=A(:,[1 3]);P1=inv((A1'*A1))*A1'*B(:,1) A2=A(:,[1 2]);P2=inv((A2'*A2))*A2'*B(:,2)上述结果代入V olterra-Lotok 方程,用MA TLAB 函数ode45求方程在时间[0,60]的数值解.作图可看到田鼠及其天敌数量的周期震荡.求方程Volterra-Lotok 的数值解的程序为定义函数vlok 为 [vlok.m]function dydt=vlok(T,Y)dydt=[(0.8765-0.5468*Y(2))*Y(1);(-0.1037+0.0010*Y(1))*Y(2)]; clear all, clcX=[29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 ...34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 ... 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5];Y=[1.6 1.3 1.1 1.2 1.1 1.3 1.8 2.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9 ... 1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3]; N=[X,Y]; T=[0:2:60];[t,Y]=ode45(@vlok,[0:0.5:60],[29.7 1.6]); plot(t,Y(:,1)/100,'k'); hold on;plot(t,Y(:,2),'-.k');title('田鼠及其天敌的V olterra-Lotok 模型拟合曲线'); xlabel('时间');ylabel('数量(只/每百)'); gtext('田鼠'); gtext('天敌');legend('田鼠','天敌');legend('田鼠','天敌');图 田鼠及其天敌的模拟曲线实线和虚线分别为田鼠和天敌的实际值,田鼠的数量为y 坐标乘以100.上机实验研究种群竞争模型设有甲乙两个种群,当它们独自生存时数量演变服从Logistic 规律,即 )1()(11n xx r t x -=⋅)1()(22n y y r t y -=⋅这里)(),(t y t x 分别为甲乙种群的数量,21,r r 为它们的固有增长率,21,n n 为它们的最大容量.当两个种群在同环境中生存时,它们之间的关系是为了争夺同资源而进行竞争.考查由于乙消耗有限的资源对甲的增长产生影响,可以合理地修改甲的方程为 )1()(2111n y s n x x r t x --=⋅这里1s 的含义为:对于供养甲的资源而言,单位数量乙(相对2n )的消耗为单位数量甲(相对1n )消耗的1s 倍.类似地,甲的存在也影响了乙的增长,乙的方程应改写为 )1()(2122n yn x s y r t y --=⋅对2s 可作相应解释.当给定种群的初始值0)0(x x = 0)0(y y =及参数212121,,,,,n n s s r r 后,公式确定了两个种群数量的变化规律.方程无解析解,一般用数值解法研究该问题,试用数值解法研究下面的问题:设0,,5.0,00,100212121========y x s s n n r r ,计算)(),(t y t x 并绘出它们的图形,求出时间t 充分大以后)(),(t y t x 的变化趋势.实例 线性差分方程模型Florida 沙丘鹤属于濒危物种,据报道,生态学家估计它在较好的环境下,每年平均增长率仅为1.94%,而在中及较差的环境下,每年平均增长率则分别为-3.24%,-3.82%,即它将逐渐减少.如果在某地的保护区内开始有100只鹤,建立描述其数量变化规律的模型,并作数值计算.而人工孵化为挽救这个濒危物种的措施之一,如果每年人工孵化5只鹤放入该保护区,那么在中环境下沙丘鹤的数量的变化规律? 问题分析及模型的建立记第k 年沙丘鹤的数量为k x ,正常环境下的平均增长率为r ,记r a +=1,则第1+k 年鹤的数量为k k ax x =+1 r a +=1 ,1,0=k在人工孵化条件下,设每年孵化的数量为b ,则 b ax x k k +=+1 模型的求解在较好,中,较差的环境下,以0382.0,0324.0,0194.0--=r 以及1000=x 代入,用MATALAB 计算并作图,程序为function y=exf(x0,n,r) %exf 的函数M 文件 a=1+r;x=x0; %赋初值 for k=1:nx(k+1)=a*x(k); end xk=(0:20)';y1=exf(100,20,0.0194); y2=exf(100,20,-0.0324); y3=exf(100,20,-0.0382);round([k,y1,y2,y3]); %对结果四舍五入取整 plot(k,y1,k,y2,':',k,y3,'--'),gtext('r=0.0194'),gtext('r=-0.0324'),gtext('r=-0.0382')结果分析讨论时间充分长的变化趋势∞→=0x a x kk 11>+=r a 0>r 00→=x a x kk 11<+=r a 0<r在人工孵化的情况,当11<+=r a 即 0<r ,得到 abx a a b x a a a b x a x k kk kk -=→--+=++++=-111]1[01用5%,24.3=-=b r 代入上式即可.一阶线性常系数差分方程的解,平衡点及其稳定性方程形式为b ax x k k +=+1这里b a ,为已知常数.令x x x k k ==+1得到代数方程b ax x +=的根abx -=1称为差分方程的平衡点. 差分方程的解可表为abca x kk -+=1 ,1,0=k这里c 由初始值0x 确定.如果∞→k 时x x k →,称平衡点x 为稳定的,否则平衡点x 为不稳定的. 差分方程的平衡点稳定1<⇔a如果第1+k 时段变量1+k x 不仅取决于第k 时段变量k x ,而且与以前时段变量有关,这得用高阶差分方程刻划.实例 高阶线性常系数差分方程模型一年生植物春季发芽,夏天开花,秋季产种,没有腐烂,风干,被人为获去的那些种子可以活过冬天,其中的一部分能在第二年春季发芽,开花,产种,其中的另一部分虽未能发芽,但如又能活过一个冬天,则其中一部分可在第三年春季发芽,开花,产种,如此继续.一年生植物只能活1年,而近似地认为,种子最多可以活过两个冬天,试建立数学模型研究这种植物数量的变化规律,及它能够一直繁殖下去的条件. 问题分析及模型的建立记一棵植物秋季产种的平均数为c ,种子能够活过一个冬天的比例为b ,一岁的种子能在春季发芽的比例为1a ,未能发芽但又能活过一个冬天的比例为b ,两岁的种子能在春季发芽的比例为2a .设21,,a a c 固定,而b 可在一定范围内变化.记第k 年的植物数量为k x ,按照种子最多可以活过两个冬天的假定,k x 与1-k x 和2-k x 有关,由1-k x 决定的部分为11-k bcx a ,而由2-k x 决定的部分则为212)1(--k bcx a b a .如果今年(0=k )种下(并成活)的植物数量为0x ,可以得到011bcx a x = 21211)1(---+=k k k bcx a b a bcx a x 2=k 记bc a p 1-=,bc a b a q )1(12--=,则001=+px x 021=++--k k k qx px x 2=k 模型的求解设20.018.0,25.0,5.0,1021-====b a a c 以及1000=x ,用MA TLAB 计算的程序为 function y=exf(x0,n,b) c=10;a1=0.5;a2=0.25; p=-a1*b*c;q=-a2*(1-a1)*c*b^2; x(1)=x0;x(2)=-p*x(1);for k=3:nx(k)=-p*x(k-1)-q*x(k-2);endxk=(0:19)';y1=exf(100,20,0.18);y2=exf(100,20,0.19);y3=exf(100,20,0.20);round([k,y1,y2,y3]),plot(k,y1,k,y2,':',k,y3,'--'),gtext('b=0.18'),gtext('b=0.19'),gtext('b=0.20')运算结果为结果分析:可以看到,对于不同的b ,k x 的变化规律有较大差别.设二阶差分方程有形如k k x λ=的解,即02=++q p λλ此方程称为差分方程的特征方程,根为2422,1q p p -±-=λ 称为差分方程的特征根,方程的解可表为k k k c c x 2211λλ+=这里常数21,c c 由初始条件10,x x 确定.本例用 bc a b a q bc a p )1(,121--=-= 5.0,101==a c 25.02=a 得到b 23052,1±=λ 当20.0,19.0,18.0=b 时,),(21λλ分别为)0477.0,0477.1(),0453.0,9953.0(),0430.0,9430.0(---18.0=b ,用90,10010==x x 代入可得64.951=c 36.42=c ,于是k k k x )043.0(36.4)943.0(64.95-+= ,1,0=k可以看出,当12,1<λ时0→k x ,当12,1>λ时∞→k x .得到植物能够繁殖的条件为191.0>b高阶线性常系数差分方程的解,平衡点及其稳定性方程为b x a x a x a x a k n k n n k n k =+++++--++11110特征方程为01110=++++--n n n n a a a a λλλ差分方程的解为x c c x k n n k k +++=λλ 11令 x x x k n k ===+ 得到的平衡点,n c c ,,1 由初始值n x x ,,1 确定,当所有的特征根的模小1时,平衡点为稳定的.实例 状态转移方程组模型随着计算机通信网络系统特别是Internet 网络的应用日益广泛,计算机网络可靠性分析以及提高系统的可靠性意义重大.研究和分析具有实用性的高可靠性计算机通信网络系统,是国际上非常活跃的一个研究方向,计算机随时可能发生三种状态,无故障,间歇故障和永久故障.因此,计算机一般处于三种工作状态,无故障工作,带故障工作和不工作,这三种状态之间的转移过程为试建立该系统的状态转移模型,并进行可靠性分析.问题分析及模型的建立该问题属于状态转移问题,用马尔科夫状态转移原理,用)(),(21t P t P 和)(3t P分别表示系统处于无故障工作,带故障工作和不工作三种状态的概率,则有状态转移方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧+++=++-=+++-=)()()()2()()()]([)()2()()()(]2)2[()(213212211t P t P dt t dP t P t P dtt dP t P t P dt t dP t p t p t p t t t p λλλλλλγλγλλλ 初始条件为0)0(,0)0(,1)0(321===P P P ,参数取值p λ 510-至410-, t λ 410-至310- γ 0.01至0.1模型的求解此为带参数的微分方程组模型.求解此方程组模型有两种方法:用特征根法求解析解;用数值解法求数值解.假定模型的参数取值为01.0,10,1045===--γλλt p ,则解析法的程序为lp=10^(-5);lt=10^(-3);gm=0.01;A=[-(lp+lt),gm,0;lt/2,-(gm+lp+lt),0;lp+lt/2,lp+lt,0];[l,v]=eig(A)e=inv(l)*[1;0;0]运算结果为数值解法的程序为%微分方程组的M 函数文件function xdot=eqs0(t,p,flag,lp,lt,gm)A=[-(lp+lt),gm,0;lt/2,-(gm+lp+lt),0;lp+lt/2,lp+lt,0];P=[p(1);p(2);p(3)];xdot=a*p;命令函数文件ts=[0 1000];p0=[1;0;0];lp=10^(-5);lt=10^(-3);gm=0.01[t,p]=ode23(‘eqs0’,ts,p0,[],lp,lt,gm)plot(t,1-p(:,3));xlabel(‘时间t(小时)’);ylabel(‘可靠度R(t)’);title(‘参数取值lp=0.0001;lt=0.001;gm=0.01’);grid on;运算结果为计算结果表明:在时间h t 1000=的情况,]0047.0,9369.0[)](),(),([321 t p t p t p显示系统工作概率为,4 插值与拟合模型实例问题分析及模型建立模型的求解5 最优化模型实例问题分析及模型建立模型的求解水箱的水流问题在许多供水单位由于没有测量流入或流出水箱的设备,而只能测量水箱的水位.试通过测得的某时刻水箱水位的数据,估计在任意时刻t 流出水箱的流量)(t f .某社区有一供水水箱,在居民用水过程,当水箱的水位下降到最低水位l 时,水泵就自动向水箱输水直到最高水位H ,此期间不能测量水泵的供水量,因此,当水泵正在输水时不容易建立水箱水位和用水量之间的关系.水泵每天输水一次或两次,每次约2小时.已知该水箱是一个高为40ft(英尺),直径为57ft 的圆柱体,表为该居民区一天水箱水位的数据,当水位降至27.00ft 时水泵开始工作,水位增到35.50ft 时,停止输水(1ft=0.3048m).表 社区某天水箱水位时间/s 水位/0.01ft 时间/s 水位/0.01ft 0 9.68 12.95 10.020.92 9.45 13.88 9.941.84 9.31 14.98 9.652.95 9.13 15.90 9.413.87 8.98 16.83 9.184.98 8.81 17.93 8.925.90 8.69 19.05 8.667.00 8.52 19.96 8.457.93 8.39 20.84 8.228.97 8.22 22.02 水泵供水9.98 水泵供水 22.96 水泵供水10.93 水泵供水 23.88 10.5910.95 10.82 24.99 10.3512.03 10.50 25.91 10.18问题分析及模型的建立由于水箱是正圆柱体,横截面积为常数,所以在水泵不工作时段,流量容易根据水位和时间的变化计算出来,但怎样估计水泵供水时段的流量比较困难.水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到.因为水泵不工作时段的流量作为用于插值或拟合的原始数据,因此水泵不工作时段的流量越精确越好.这些流量可以用两种方法来计算:(1)对表的数据用数值微分计算出各时段的流量,从而拟合其他时刻或连续时间的流量;(2)先用表的数合数据拟合水位 时间函数,再求导数就可以得到连续时间的流量.有了每个时刻的流量,就可以计算水箱的总流量.水泵不工作时段的用水量可以由测量记录直接得到,由表的数据可以直接计算出;水泵工作时水箱的流量通过拟合出来的流量函数计算出.这样就可以计算出水箱的总流量.模型假设流量看做是时间的连续函数,为了便于计算,不妨将流量定义为单位时间流出的水的高度,即水位对时间变化率的绝对值,水箱的截面积为)()24.3.057(422m S ⨯⨯=π在计算总流量时将上面得到的结果乘以S 即可.流量只取决于水位差,与水位本身无关,即流出的水的流速正比于水面高度的平方根.题目给出水箱的最低和最高水位分别是8.1648m 和10.7352m(设出口的水位为0),计算得sqrt(10.7352/8.1648),大约为1,故可以忽略水位对流速的影响.流量估计方法用一种比较简单的方法计算水箱流量与时间的关,将表的数据分为5段,按时间t 排列:第一段:0—8.97;第二段:9.38—10.93;第三段:10.95—20.85;第四段:22.02—22.96;第五段:23.88—25.91.再对每一段的数据做如下的处理:设某段数据为{),(,),,(),,(1100n n y x y x y x },邻近数据中点的平均流速用公式流速=(左端点的水位-右端点的水位)/区间长度即 ii i i i i x x y y x x v --=++++111)2( 计算;每一段数据首尾点的流速用下面公式计算)/()43()(022100x x y y y x v -+-=)/()43()(221-----+-=n n n n n n x x y y y x v根据上面的公式,可以计算出时间与流速之间的数据如表表 时间与流速之间的数据表时间/h 流速/(cm/h) 时间/h 流速/(cm/h)0 29.89 11.50 29.850.46 22.05 12.49 31.521.38 18.47 13.52 29.032.395 16.22 14.52 26.503.52 16.29 15.50 26.094.52 15.30 16.47 24.795.45 13.05 17.38 23.676.45 15.45 18.49 23.507.465 13.98 19.52 25.208.45 16.45 20.50 23.858.97 19.29 20.85 22.259.98 水泵供水22.02 水泵供水10.93 水泵供水22.96 水泵供水10.95 30.50 25.91 13.15用两种计算方法建立模型(1)插值方法由表,对水泵不工作时段采取插值方法,可以得到任意时刻的流速,从而知道任意时刻的流量,这里分别采用拉格朗日(Lagrange)插值法,分段线性插值法和样条插值法作插值.对水泵工作时段2应用前后期的流速作插值,由于第5段水泵不工作时的数据太少,将其与水泵工作时段4合并一起作插值.这样就总共得对4段数据作插值(第1,3未供水时段,第2供水时段,第4,5时段的混合时段).(2)曲线拟合法拟合水位时间函数.根据表的测量记录知,一天有两次供水时段和三次未供水时段,分别对1,3未供水时段的测量数据直接作多项式拟合,可以得到水位函数,再由水位时间函数确定流量时间函数,这样也可以求出一天总用水的估计.模型求解的MA TLAB程序插值法以第一段未供水时数据为例分别用拉格朗日,线性多项式,样条插值方法计算出流量函数和用水量.由于MA TLAB没有直接提供拉格朗日插值法的命令函数,这里先给出用MATLAB语言实现的拉格朗日插值法的函数lglrcz.mfunction Y=lglrcz(X0,Y0,X)n=length(X0)m=length(X);for I=1:mz=X(i)s=0;for k=1:np=1.0for j=1:nif j~=kp=p*(z-X0(j))/(X0(k)-X0(j));endends=p*Y0(k)+sendY(i)=s;End%由表可得到t=[0 0.46 1.38 2.395 3.45 4.525 5.50 6.45 7.465 8.45 8.97];v=[29.89 20.75 18.45 16.22 16.25 15.32 13.05 15.45 13.98 16.25 19.25];t0=0:0.1:8.097lr=lglrcz(t,v,t0); %拉格朗日插值法lrjf=0.1*trapz(lr)fdcz=interp1(t,v,t0); %分段线性插值法fdczjf=0.1*trapz(fdcz)scz=interp1(t,v,t0,’spline’); %样条插值法sczjf=0.1*trapz(scz)plot(t,v,’*’,t0,lr,’r’,t0,fdcz,’g’,t0,scz,’b’)gtext(‘lglr’);gtext(‘fdxx’);gtext(‘syt’);其运算结果为一般情况,样条插值方法具有比较好的性质,大多数情况下都采用该方法.另外,其他时段的处理方法与第一段未供水时段的处理方法类似,只给出结果.表各时段和一天的总用水量(用水高度)第一未供水时段第二供水段第三供水时段第四混合段全天拉格朗日插值法145.622 258.866 54.2689 92.1335 550.6922分段线性插值法147.145 258.9697 49.6055 76.4688 532.1866样条插值法145.687 258.6557 53.35 81.7699 539.4652拟合法拟合水位时间函数t,分别为已输入的时刻和水位测量记录(由表得到,水泵供水的4个时刻不输入),第一设h未供水时段各时刻的水位可以由下面程序实现,如图所示t=[0 0.92 1.84 2.95 3.87 4.98 5.90 7.00 7.93 8.97 10.95 12.03 12.9513.88 14.98 15.90 16.85 17.93 19.04 19.96 20.85 23.88 24.99 25.66];h=[9.68 9.48 9.32 9.13 8.98 8.81 8.69 8.52 8.39 8.22 10.82 10.50 10.22 9.949.65 9.41 9.18 8.92 8.66 8.45 8.22 10.59 10.35 10.18];c1=polyfit(t(1:10),h(1:10),3);tp1=0:0.1:8.9x1=polyval(c1,tp1);pot(tp1,x1);变量X1存放了以0.1为步长计算出的各个时刻的水位高度.第二未供水时段时间水位图可以由下面程序实现,如图所示c2=polyfit(t(10:20),h(10:20),3);tp2=10.9:0.1:20.9X2=-polyval(c2,tp2);Plot(tp2,X2)确定流量时间函数c1=polyfit(t(1:10),h(1:10),3);c2=polyfit(t(10:20),h(10:20),3);a1=polyder(c1);a2=polyder(c2);tp1=0:0.01:8.97tp2=10.95:0.01:20.85X13=-polyval(a1,tp1);X013=-polyval(a1,[0:0.01:8.97]);Wgsysll=100*trapz(tp1,X013);X4=-polyval(a1,[7.93,8.97]);X23=-polyval(a2,tp2);X0=-polyval(a2,[10.95:0.01:20.85]);Wgsys=100*trapz(tp2,X0);X00=-polyval(a2,[10.95,12.03]);X=-polyval(a2,[19.96,20.85]);Plot(tp1,X13*100);Plot(tp2,X*100).结果如图第二供水段的流量则用前后时期的流量做拟合得到.为使流量函数在11,9==t t 连续,只取四个点,用三次多项式拟合得到第二供水时段的时间 流量如图.实现的程序为dygsdsj=[7.93 8.97 10.95 12.03];dygsdls=[X0,X];nhjg=polyfit(dygsdsj,dygsdls,3);nhsj=7.93:0.1:12.03nhlsjg=polyval(nhjg,nhsj);gssjl=8.97:0.01:10.95gsl=polyval(nhjg,[8.97:0.01:10.95]);gsysll=100*trapz(gssjl,gsl);plot(nhsj,100*nhlsjg)在第四供水时段之前取85.20,96.19=t 两点的流量,用第五未供水时段的三个记录做差分得到两个流量数据22.52,18.52,再用这四个数据做三次多项式拟合得到第四供水时段与第五未供水时段的时间 流量函数,如图,程序为t3=[19.96 20.85 t(22),t(23)];ls3=[X*100,22.52,18.52];nd=polyfit(t3,ls3,3);tp3=19.96:0.01:25.91;X=polyval(nd,tp3);Gsj=20.85:0.01:25;Gs2=polyval(nd,[20.85:0.01:25]);Gsys=trapz(gssj2,gs2);Plot(tp3,X);一天总用水量的估计分别对供水的两个时段和不供水的两个时段积分(流量对时间)并求和得到一天的总用水量ft(总用水高度,单位为cm).各时段用水量如表约为526.89352表各时段用水量及一天总用水量(单位:cm)时段落第一未供水时段第二供水时段第三未供水时段第四混合时段全天用水ft2用水高度145.65 260.66 46.60 73.9625 526.8925 微分方程模型最优化方法模型实例截断切割问题某公司经常得从一个长方体中加工出一个尺寸,位置预定的长方体(这两个长方体的对应表面是平行的),通常采用截断切割的加工方式,这里”截断切割”是指物体沿某个切割平面分成两部分.因此在一般情况下,得经过6次截断切割,分别截去原长方体的前,后,左,右,上,下的6个方向多余的部分.设水平切割单位面积的费用是垂直切割单位面积的r倍,且当先后两次垂直切割的平面不平行时,因调整刀具需额外费用e.如果截去各方向多余小块的先后顺序不同,则加工费用不同.试设计确定最优加工次序的方法,此处的最优是指加工费用最少(由工艺需求,与水平工作台接触的长方体底面是事先指定的).用下面实例验证所设计的方法:需加工长方体与成品长方体的长,宽,高分别为10,14.5,19和3,2 4,二者左面,前面,底面之间的距离分别为,6,7,9(单位:cm),垂直切割费用为1元/cm2,r 和e 的数据有四组.0,1==e r 0,5.1==e r 0,7==e r 152,5.1≤≤=e r 问题分析:这是一个优化问题,求切割顺序,使加工费用最低.决策变量为切割顺序,用),,(61x x X =表示切割顺序,i x 表示第i 次切割,可以分别表示前,后,左,右,上,下的切割,61,,x x 互不相同,可以取6,,1 的任意全排.目标函数:加工费用由切割费用和刀具调整费用构成.问题的已知条件有需加工长方体与成品长方体对应表面平行切割费用与切割面的面积成正比,具体地说就是垂直切割费用为1元/cm2,水平切割费用为r 元/cm2,且仅当先后两次垂直切割的切割面不平行时,才需调整刀具,调整刀具的费用为e .水平工作台接触的长方体底面是事先指定的.不考虑第一次切割前的刀具调整费用.数学模型设需加工的长方体的长为a ,宽为b ,高为c 为常数,需加工长方体与成品长方体两者的前,后,左,右,上,下面之间的距离为212121,,,,,c c b b a a 也为常数.可变参数有:水平切割费用r 元/平方厘米,调整刀具的费用e .在切割方式X ,对应的加工费用可表示为),,(r e X f .可得组合优化模型求切割方式min X X =,使加工费用),,(r e X f 达到最小,即),,(min r e X f SX ∈ 这里},,6,,1),,({61j i x x x x x X S j i i ≠≠===由于集合S 为有限集,只有720!6=种切割方式,当r e ,取定,切割顺序给定,很容易算出加工费用.可以依次求出各切割方式下的切割费用,比较最小者,就可得到最小费用的加工顺序. 解程序为情形一 0=e先用穷举法求出720种切割方式的费用,存放在数组c ,再用函数min(c)和find(c==minc)求最小费用及其对应的切割方式.%jieduan e=0% a0 三维向量,各分量为需加工长方体的长,宽,高% a1 三维向量,各分量为成品长方体的长,宽,高% d1 三维向量,各分量为需加工与成品长方体两者的前面,左面,底面之间的距离.% r 水平切割单位面积的切割费用% minc 最小费用%min X 列数为6的矩阵,各行为最小费用对应的切割顺序a0=[10, 14.5,19];a1=[3,2 4];d1=[6,7,9];r=1;d2=a0-a1-d1; d=[d1 d2];d=d([1,4,2 5,3,6]);p=0%可行的加工顺序表.For I=1:6For j=1:6,if(j-i)~=0,For k=1:6,if(k-I)*(k-j)~=0,For l=1:6,if(l-i)*(l-k)~=0,For l=1:6, if(m-i)*(m-j)*(m-k)*(m-l)~=0;For n=1:6If(n-i)*(n-j)*(n-k)*(n-l)*(n-m)~=0,P=p+1;X(p,:)=[I,j,k,l,m,n];End,end,end,end,end,end,end,end,end.%加工顺序表X 对应的切割费用表f=[1,1,2 2 3,3 ];for p=1:720o=X(p,:);const=0;a=a0;foe I=1:6j=o(i);a3=a;a3(f(j))=[];if f(j)==3const=cost+r*a3(1)*a3(2);elseconst=const+a3(1)*a3(2);enda(f(j))=a(f(j)-d(j));endc(p)=cost;end.%求最小费用及其对应的加工顺序minc=min(c),find(c==minc);minx=x(ans,:);运算结果为.因此,当0,1==e r 时,最优加工顺序为 下 前 左 上 后 右 或 下 前 上 左 后 右切割费用为374元情形二 0≠e加工费用是由切割费用和调整刀具的费用两者组成,即e z r Xf r e X f ⨯+=),0,(),,(这里,z 为加工顺序是X 时的调整刀具次数,全体切割顺序按调整刀具次数划分为三类,同类的刀具调整费用是相同的.可以先分别求出在0=e 时,每一类的最小费用及相应的加工顺序,它们就是各类的最优加工顺序.再用每一类的最小切割费用加上相应的刀具调整费用,得到加工总费用.各类的最优加工顺序进行比较,就可得整体的最优加工顺序.0>e 时,%jieduan e>0function[min,minx1,minx2 minx3]=cutordel(a0,a1,d,r)minc=[inf,inf,inf];minx1=[];minx2=[];minx3=[];k1=0;k2=0;k3=0;v1=[1 3 6];%三类可行的加工顺序表x1,x2 x3及相应的切割费用表.%c1 c2 c3for i1=1:6,ol=v1(i1);v2=v;v2(i)=[];for i2=:5,o2=v2(i2);v3=v2;v3(i2)=[];for i3=1:4,o3=v3(i3);v4=v3;v4(i3)=[];for i4=1:3,o4=v4(i4);v5=v4;v5(i4)=[];for i5=,o5=v5(i5);o6=(3-i5);x=[o1,o2 o5 o6];c=cost(x,a0,a1,d1,r);z=adjustnum(x);switch zcase 1k1=k1+1;x1(k1,:)=x;c1(k1)=c;case 2k2=k+1;x2(k2:)=x;c2(k2)=c;case 3k3=k3+1;x3(k3,:)=x;c3(k3)=c;end,end,end,end,end.Minc=[min(c1),min(c2),min(c3)];Find(c1==minc(1));Minx1=x1(ans,:);find(c2==minc(2));Minx2=x2(ans,:);find(c3=minc(3));Minx3=x3(ans,:).求切割顺序是x 时,切割费用的子函数const 为function c=const(x,a0,a1,d1,r)c=0;d2=a0-a1-d1;a=a0;for p=1:6switch x(p)case 1c=c+a(2)*a(3);a(1)=a(1)-d1(1);case 2c=c+a(2)*a(3);a(1)=a(1)-d2;case 3c=c+a(1)*a(3); a(2)=a(2)-d(2);case 4c=c+a(1)*a(3);a(2)=a(2)-d2(2);case 5c=c+r*a(1)*a(2);a(3)=a(3)-d1(3);case 6c=c+r*a(1)*a(2);a(3)=a(3)-d2(3);endend.%求加工顺序x的调整刀具次数的子函数adjustnum(x)为function z=adjustnum(x)z=-1;v0=0;for p=:6if x(p)<5if x(p)<3v=1;elsev=2;endif(v0-v)~=0z=z+1;v0=v;endendend在MA TLAB输入命令a0=[10 14.5 19];a1=[3 2 4];d1=[6 7 9];r=1.5;[minc minx1 minx2 minx3]=cutrode(a0,a1,d1,r)运算结果为因此,每一类的最小费用分别为:C1(e)=+e 此时调整一次刀具。
