1、微分方程组应用微分方程组应用 微分方程组在工程技术中的应用是非常广泛的，涉及的领域有机械、电工技术、通讯、医学等许多领域。下面给出几个例子，用以阐述关于微分方程组的实际应用及其建模思想。1 1、两自由度的振动问题两自由度的振动问题 常系数线性方程组在工程技术与科学研究中有很多应用，不少问题都归结为它的求解问题。两个自由度的振动问题。为了消除不需要的振动，两个自由度的振动问题。为了消除不需要的振动，常常在振动系统中设置减振器。常常在振动系统中设置减振器。其中1m2m1k2k1k2kcF1x2x1m2m1m2m的受力情况。是原机械部件的质量；是减振器的质量；和是两个弹簧，它们的弹性系数（或称为刚度
2、）也分别用和是减速器表示；（假定阻力与速度成正比）的阻尼系数；是强迫力；和分别表示和距它们的平衡位置的位移。体和下面来建立这个系统的运动方程。先分别考虑物2m1k2k1m1x2m2x2k21xx2k221kxx2222212d xmkxxdt 律，有方程(1)物体的受力情况假定弹簧和都满Hook定律。当物体位移时，物体同时位移同时位移这时，弹簧变形（拉长或压缩）的长度为。因此，这时弹簧的弹性力是（力的方向与位移方向相反），由牛顿第二定(2)物体 的受力情况1)沿位移方向的外力：2)阻尼力：（方向与速度方向相反）；3)这时物体 受到两个弹簧的作用：1m0sinFFwt1dxcdt1m弹簧 的弹性
3、力 ；弹簧 的弹性力 。由牛顿第二定律，有因此，上述运动系统满足微分方程组1k11k x2k221kxx21110112212sind xdxmFwtck xkxxdtdt211112122022222212sin0d xdxmckkxk xFwtdtdtd xmkxxdt上述方程组在变换 及 之下，就变成了一个常系数线性非齐次方程组 (4.6.1)11dxydt22dxydt11112212101112222212sindxydtdykkkcxxyFwtdtmmmdxydtdykxxdtm 为了简单，只考虑无阻尼自由振动的情形，即 。于是，(4.6.1)变成 (4.6.2)00,0cF111
4、12212112222212dxydtdykkkxxdtmmdxydtdykxxdtm 它的特征方程为 设 ，则(4.6.3)可写成解之，得122124212120(4.6.3)kkkk kmmm m2z12212212120kkkk kzzmmm m2122122121,2121212142kkkkkkk kzmmmmm m容易看出，因此，可令 ，这时因此，方程(4.6.3)最后可写成故(4.6.3)的特征根全是纯虚根。(4.6.1)解的第一个函数可写成形如120,0zz 2212,zz 22122122122121212142kkkkkkk kmmmmm m2222011234sincos
5、sincosxCtCtCtCt11122sinsinxAtAt2m2k2m2k部件的运动频率。这个事实可以用来防止机器或这是两个简谐运动的叠合。每一个简谐运动的角频（即与）均与减振器的参数与有关。因此，调整与可以改变原设备与外力发生共振现象，以及减轻外力的干扰等。3 3、人造卫星的轨迹方程、人造卫星的轨迹方程我们知道，人造卫星在最后一段时间运载火箭熄灭之后，即进入它的轨道，轨道的形状因发射角度和发射速度的不同，而分别出现椭圆、抛物线或双曲线。下面来讨论这些问题。地球与人造卫星是相互吸引的，但因二者的质量相差很大，因此，可假设地球是不动的，又因人造卫星的体积与地球相比是很小的，故可把它看作质点。
6、为简单起见，我们不考虑太阳、月亮和其它星球的作用，并略去空气阻力。在上面的假设下，可把问题归结为：从地球表面上一点 ，以倾角 ，初速度 射出一质量为 的物体，如图4.5所示，求此物体运动的轨道方程。图4.5A0vm,p x yAOyyxoyyx0t txy22mMFfxy 根据万有引力定律，地球对卫星的引力大小为过发射点和地心的直线作轴，轴与发射方向所成的平面为平面，平面通过地心，取垂直于轴且过地心的直线为轴，取开始发射时间为，经过时间后，卫星位于点，下面建立和所满足的方程。其方向指向地心，其中 是引力系数，是地球质量，是地球与卫星间的距离。如图4.6所示，这个引力在 轴方向上的分力分别为f2
7、0326.685 10/fkmkg sM245.98 10Mkg22xy,x y卫星在 轴上所获得的分加速度分别为 和 。由牛顿第二定律，得到卫星的运动方程为：32223222cossinxyfmMxFFxyfmMyFFxy ,x y22d xdt22d ydt232222232222(4.6.5)d xfmMxmdtxyd yfmMymdtxy 0t 0v0t 00,0 xyR6370Rkm地球半径,x y轴上的分量分别为当时，卫星在地表面以倾角，初速度射出，所以，在时，。卫星的初速度在因此，初始条件为 （4.6.6）下面利用首次积分法来求方程组(4.6.5)满足初始条件(4.6.6)的解。
8、将方程(4.6.5)两端消去 后，以 乘第一个方程，以 乘第二个方程。0000cos,sinttdxdyvvdtdt 000000,0,cos,sinttxyRdxvdtdyvdtmyx然后相减，得因为故有对上式两边积分，得首次积分为22220d xd yyxdtdt2222ddydxd yd xxyxydtdtdtdtdt1dydxxycdtdt0ddydxxydtdtdt再以 乘方程组(4.6.5)中第一个方程，以 乘第二个方程，然后两式相加，得由于dxdtdydt22322222dx d xdy dyfMdxdyxydt dtdt dtdtdtxy 2222222ddxdydx d xd
9、y dydtdtdtdt dtdt dt及从而得积分得上式另一首次积分1322222222dfMfMdxdyxydtdtdtxyxy 22222ddxdydfMdtdtdtdtxy22212222dxdyfMcdtdtxy于是，原方程组(4.6.5)降为较低阶的方程组 12221222(4.6.7)2dydxxycdtdtdxdyfMcdtdtxycos,sinxryrcossin4.6.8sincosdxdrdrdtdtdtdydrdrdtdtdt将它们代入(4.6.7)得作极坐标变换，并求导 两式联立消去 ，得 212222(4.6.9)2drcdtdrdfMrcdtdtrddt21222
10、(4.6.10)drfMccdtrr rr t rr t t rr tt 这就是卫星运动轨道的极坐标参数方程。若将参这里得到一个仅含有一个未知数得一阶微分方程，若由此解出，代入(4.6.9)中第一式，便可确定，由此得到数消去，便得出卫星运动轨道的极坐标方程。2212212drrfMccdcrr21rdtdc2112222111cosfMrcccfMcc利用分离变量求该方程的解得为此，由(4.6.9)的第一式求得,并代入(4.6.10)得1222221111cosfMcfMcrccc2212 12,1cc cpefMfM11cosecrpp，则上式化为整理得令1cosprec,p e c12,c
11、 c c0t 0,02rR0000sin,costtdrdvvdtdtR 知或(4.6.11)这就是所求的卫星运动轨道的极坐标方程，其中有三个任意常数(或)，它们可由初始条件(4.6.6)确定。注意到当时，因此，。且由(4.6.6)及(4.6.8)把它们代入(4.6.9)及(4.6.11)得 (4.6.12)我们知道，(4.6.11)式是圆锥曲线的极坐标方程10220cos21sincRvfMcvRpRce 当 时，轨道是圆；当 时，轨道是椭圆；当 时，轨道是抛物线；当 时，轨道是双曲线 下面进一步讨论卫星发射的初速度与卫星轨道形状的关系因为0e 01e1e 1e 22 121c cefM22
12、 1221c cefM 12,c c22224220022coscossin1RvR vefMfM222210coscR vpfMfM故上故将(4.6.12)式中的代入上式，并整理得(4.6.13)注意到(4.6.12)式及式可化为因此，当 时，得由此得222221tanppeRR0e 1,tan0pR2222200cosR vR vpRfMfM20fMvRRMf22062.76/vkm s07.9/vkm s07.9/vkm s卫星的轨道是一个圆.即把地球半径，质量及引力常数的具体数值代入上式，并计算得即称为第一宇宙速度第一宇宙速度，此时1e 2020fMvR02 fMvR2011.2/vk
13、m s当时，由(4.6.13)得即所求速度是第一宇宙速度的倍，即，称为第二宇宙速度第二宇宙速度，它的轨道是抛物线。01e1e 200220,fMfMvvRR1e 02 fMvR是双曲线（一支）。当时，因，由(4.6.13)式可知这表明，当初速度小于第二宇宙速度时，卫星轨道是一个椭圆。当时，由(4.6.13)式可得因此，当初速度大于第二宇宙速度时，它的轨道由电路和机械装置组装成的永磁体扩音器模型如图4.7所示。4 4、扩音器振动模型、扩音器振动模型图4.7 一个时变电源电压 驱动着一个音圈能换器，能换器的转化系数为 ，通过能换器使扬声器的振动膜发生振动。由音圈组成的能换器，本质上是一个在永磁场内
14、的自由运动。当变化的电流通过音圈时，音圈在永磁体的磁力和电流产生的磁力相互作用下进行运动。用 代表能换器与扬声器相互作用力，是能换器电阻，是能换器的感应系数.E tTfRLmCkT,dxeTfTidt ex x压定律，可得的微分方程组质量为的扬声器的振动是一个具有阻尼的弹簧系统振动。阻尼系数为，弹簧的弹性系数为，能换器转化系数与有关变量间的相互关系为这里是音圈两端的电压降；是音圈位移；代表音圈的速度，由牛顿第二定律及回路电 22d xdxmckxTidtdtdxdiTLRiE tdtdtdxydt(4.6.14)令，把(4.6.14)化为一阶微分方程组 dxydtdykcTxyidtmmmdi
15、TRyiE tdtLL xAxg t(4.6.15)方程组(4.6.15)的向量形式为，这里设 则 0100,0,0 xkcTAg txymmmE tiTRLL 1,1,1,1,2,1,cosmkcTRLE twt 0100111,0012cosAg twt 矩阵 的特征方程为特征根为 ，其相应特征向量分别为A3234201231,1,1ii 1231111,1,11iiii xA x cossinsincossincossincosttttttttteetetteettetteetet g t 因此，可求得齐次线形方程组的基解矩阵为为求(4.6.15)的通解，需求出(4.6.15)的一个特写
16、成解，我们用待定系数法来求此特解。把因此，方程组(4.6.15)有特解将 代入方程 得sincoscossincoswwt awwt bwt Aawt Abwt a 0cos,01g twt 这里cossinpxwt awt b px xAxg t sinwtcoswt(4.6.16)waAbwbAaaa22Aw E bwaab比较上述方程两边和的系数得把(4.6.16)的消去得由此可求得和分别为xx2264264222264264242426426424324444324,444422224444w wwwwwwwwwwwwbawwwwwwwww wwwwwwww xAxg t cossin
17、xt Cwt awt b 的通解为方程这里 是任意常向量。进一步，可求得原方程中音圈的位移 和驱动能换器的电流 分别为Cix123226421234242642cossin32 cos4 sin44sincos22 cos22 sin44ttttttxc ec etc etwwtw wwtwwwic ec etc etwwwtw wwwtwww5、人体含铅量问题、人体含铅量问题铅是如何进入体内的铅是如何进入体内的?消化道消化道 通常是铅吸收的主要途径，儿童由这一途径吸收的铅所占比例约为85。铅在肠道内通过主动转运和被动扩散两种方式由小肠吸收人血液。主动转运占吸收总量的80以上。主动转运依赖铅与
18、肠粘膜上一种转运蛋白质结合，由转运蛋白作载体将铅转运人血液。被动扩散则是铅由肠腔向血液自然扩散。肠腔中铅浓度越高、血液中铅浓度越低，被动扩散的量就越大。呼吸道呼吸道 空气中含铅的灰尘经鼻孔吸人呼吸道后，一部分被鼻毛、气管纤毛和支气管纤毛、细支气管纤毛阻挡，最后以痰液的形式返回口腔，儿童将其咽人消化道，再由消化道吸收入血。另一部分特别微小的铅尘到达肺泡，沉积在肺泡里，然后由吞噬细胞等吸收，进入血液。皮肤皮肤 当我们接触有机铅的时候，铅会被皮肤直接吸收。铅中毒会引起神经、消化、循环系统紊乱，表现为贫血、腹痛、高血压等。血铅高组的儿童的总智商、操作智商、语言智商分别比低血铅组落后14分、14分和13
19、分，而每升血液中的铅浓度上升100微克，儿童的身高将降低13厘米。低剂量的铅接触可以对人体的红细胞、肾脏、免疫系统、骨髓和中枢神经的功能产生不良影响，而所有这些影响发生前都可能没有明显的临床症状。铅中毒会影响婴幼儿最初站立、行走和说话的年龄，也可能引起孩子注意力涣散、记忆力减退、理解力降低及学习困难等。儿童体内的铅水平可分从血、骨、齿、尿、发检测到，儿童血铅水平分为5级。级：血铅值低于10微克分升，身体处于相对安全状态；级：血铅值为10微克19微克分升，属于轻度铅中毒。影响造血、神经传导和认知能力，儿童易出现头昏、烦躁、注意力涣散、多动、厌食、腹胀、轻度贫血；级：血铅值达到20微克44微克分升
20、，为中度铅中毒。引起缺钙、缺锌、缺铁、免疫力低下，运动不协调，视力和听力受损，学习困难、智力下降，生长发育迟缓，贫血、腹绞痛、食癖、反应迟钝等；级：血铅值达到45微克69微克分升，就是重度铅中毒。可出现性格改变、易激怒、攻击性行为、运动失调、贫血、腹绞痛、高血压、心律失常和痴呆等；级：当血铅值大于70微克分升时，为极重度铅中毒，可导致脏器损害、铅性脑病、瘫痪、昏迷等。问题和模型问题和模型 铅是一种人体所必须的微量元素，但体内铅含量过多时就会引起铅中毒.铅是一种重金属元素，通过食物、饮料、空气进入人体，经过呼吸和消化系统后进入血液，再经过血液循环慢慢进入人体和骨头中.铅可以经过人体的排泄系统、通
21、出出汗、剪头发、剪指甲排出体外.我们根据铅在人体内的变化情况将人体分为血液、组织、骨头3个仓室，铅在这三个仓室中的转化关系如图所示x1(t)表示t时刻血液中的含铅量，x2(t)表示t时刻组织中的含铅量，x3(t)表示t时刻骨头中的含铅量.假设在单位时间内从环境经过消化、吸收系统进入血液的铅为L，从血液进入组织、骨头的铅分别为a31*x1(t)和a21*x1(t)，从组织、骨头再进入血液的铅分别为a13*x2(t)、a12*x2(t)，血液和骨头向外界排出的铅分别为a01*x1(t)和a02*x3(t).示意图示意图 (1)根据前面的假设建立人体血液、组织和骨头中含铅量所满足的微分方程组,并求其
22、通解.(2)取时间单位为天，对一个志愿者在一种环境中的生活情况进行测定得a21=0.011，a12=0.012,a31=0.0039，a13=0.000035,a01=0.021,a02=0.016,L=49.3.假设该志愿者开始时体内的含铅量为0，求他体内血液、组织和骨头中含铅量随时间变化的关系，画出这些解曲线的图，并求当x+时这些函数的极限.(4)假设该志愿者在此环境中生活了365天后搬到了一个无铅的环境中去(不再有外界的铅进入体内，即L=0)，再讨论他体内血液、组织和骨头中含铅量随时间变化的关系，并画出0t1460时这些含铅量曲线的图形.简化图简化图模型,0)0(,0)0(,0)0()(
23、)()()()()()()()()()()(321313131321202121231321213121011xxxtxatxadttdxtxaatxadttdxbtxatxatxaaadttdx模型模型 x=Ax+b,数据数据Lead Transfer Coefficients Lead Transfer Coefficients(Rabinowitz,et al.)Units:days-1a a2121=0.011=0.011a a1212=0.012=0.012from blood to tissue and backa a3131=0.0039=0.0039a a1313=0.0000
24、35=0.000035from blood to bone and backa a0101=0.021=0.021a a0202=0.016=0.016excretion from blood and tissue时间分为两段:0,365和365,1460#定义常数定义常数a01:=0.021;a02:=0.016;a12:=0.012;a13:=0.000035;a21:=0.011;a31:=0.0039;L:=49.3;#定义矩阵和向量定义矩阵和向量A:=matrix(3,3,-(a01+a21+a31),a12,a13,a21,-(a02+a12),0,a31,0,-a13);b:=m
25、atrix(3,1,L,0,0);#定义方程定义方程equn1:=diff(x1(t),t)=-(a01+a21+a31)*x1(t)+a12*x2(t)+a13*x3(t)+L;equn2:=diff(x2(t),t)=a21*x1(t)-(a02+a12)*x2(t);equn3:=diff(x3(t),t)=a31*x1(t)-a13*x3(t);#解初始值问题解初始值问题dsolve(equn1,equn2,equn3,x1(0)=0,x2(0)=0,x3(0)=0,x1(t),x2(t),x3(t);结果太复杂结果太复杂,利用迭加原理、特征值和特利用迭加原理、特征值和特征向量法征向量
26、法#非齐次方程的特解非齐次方程的特解with(linalg):with(plots):xe:=evalm(-(inverse(A)&*b);#特征值和特征向量特征值和特征向量eigenvals(A);eigenvects(A);lambda:=-.3061796847e-4,-.1980315266e-1,-.4410122938e-1;v1:=.1124436946e-1,.442226656e-2,10.00746814;v2:=-.5975248926,-.8018660787,.1178839056;v3:=-.8256380196,.5640574390,.7307156328e-1
27、;#特征向量组成矩阵特征向量组成矩阵 P:=augment(v1,v2,v3);#得到解的表达式得到解的表达式x1:=c1*P1,1*exp(lambda1*t)+c2*P1,2*exp(lambda2*t)+c3*P1,3*exp(lambda3*t)+xe1,1;x2:=c1*P2,1*exp(lambda1*t)+c2*P2,2*exp(lambda2*t)+c3*P2,3*exp(lambda3*t)+xe2,1;x3:=c1*P3,1*exp(lambda1*t)+c2*P3,2*exp(lambda2*t)+c3*P3,3*exp(lambda3*t)+xe3,1;#在解的表达式用
28、在解的表达式用t=0t=0代入代入x10:=simplify(subs(t=0,x1);x20:=simplify(subs(t=0,x2);x30:=simplify(subs(t=0,x3);#求解常数求解常数solve(x10=0,x20=0,x30=0,c1,c2,c3);assign(%);#作图作图plot(x1,x2,x3,t=0.365,color=red,blue,green,thickness=2);plot1:=%:#重新显示解重新显示解x1;x2;x3;#求极限求极限limit(x1,t=infinity);limit(x2,t=infinity);limit(x3,t
29、=infinity);#得到得到365365天后解的表达式天后解的表达式xx1:=cc1*P1,1*exp(lambda1*t)+cc2*P1,2*exp(lambda2*t)+cc3*P1,3*exp(lambda3*t);xx2:=cc1*P2,1*exp(lambda1*t)+cc2*P2,2*exp(lambda2*t)+cc3*P2,3*exp(lambda3*t);xx3:=cc1*P3,1*exp(lambda1*t)+cc2*P3,2*exp(lambda2*t)+cc3*P3,3*exp(lambda3*t);#计算解在计算解在t=365t=365时的值时的值 x1365:=
30、simplify(subs(t=365,x1);x2365:=simplify(subs(t=365,x2);x3365:=simplify(subs(t=365,x3);xx1365:=simplify(subs(t=365,xx1);xx2365:=simplify(subs(t=365,xx2);xx3365:=simplify(subs(t=365,xx3);#求解常数求解常数solve(xx1365=x1365,xx2365=x2365,xx3365=x3365,cc1,cc2,cc3);assign(%);#重新显示解重新显示解xx1;xx2;xx3;#求极限求极限limit(xx1,t=infinity);limit(xx2,t=infinity);limit(xx3,t=infinity);#作图作图plot(xx1,xx2,xx3,t=365.1460,color=red,blue,green,thickness=2);plot2:=%:#将两个图画在一个图上将两个图画在一个图上display(plot1,plot2);
