Maple-与微分方程模型

Maple与微分方程模型,主要内容学习、了解数学软件Maple；用Maple解决常微分方程问题；对一些实际问题建立微分方程模型，用Maple软件求解；一些实际问题的动画演示,参考资料周义仓、靳祯、秦军林，常微分方程及其应用，科学出版社，2003年6月。李强，Maple8基础教程，中国水利水电出版社，2004DavidBarrow等，SolvingDifferentialEquationswithMapleV,Brooks/ColePublishingCompany,1998.0/resource/maple.htm0/maple/tour/index.htm,Maple9教程“四部曲”Maple9GettingStartedGuide,834KB/sheaweb/files/m9GettingStartedGuide.pdfMaple9LearningGuide,4.49MB/sheaweb/files/m9LearningGuide.pdfMaple9IntroductoryProgrammingGuide,7.08MB/sheaweb/files/m9IntroductoryProgrammingGuide.pdf4.Maple9AdvancedProgrammingGuide,3.63MB/sheaweb/files/m9AdvancedProgrammingGuide.pdfMaple应用,1.Maple软件包介绍1.1引言Maple是Waterloo大学开发的集数值计算、符号运算和图形显示于一体的软件包。它可以把人们从烦琐的代数运算中解脱出来，使得我们将主要精力集中与分析和解决问题的思路与方法中去。数学运算分为数值运算、符号运算、逻辑运算等。数值运算：如函数求值、方程求根、矩阵的特征值、特征向量等，数值计算是一项简单烦琐的工作，计算机的出现解决了数值计算的困难，使得大规模、高速度的计算可以实现，使得数学在天气预报、油藏模拟、航天等领域发挥了巨大的作用。,符号运算一种更加智能化的计算，所处理的符号可以代表各种数、代数式、数学结构（集合），代数运算是在一些代数规则下所进行的数学处理，主要是寻求一个简单完美的公式解答。用计算机进行符号运算是数学与计算机领域中一个新的发展方向，是人们用计算机代替大脑进行数值运算后用计算机代替大脑进行代数运算的成功实践，是数学走向机械化的重要步骤，这也使得计算机更加智能化。目前数学软件包很多，如Mathematica;Maple;Matlab;MathCAD;SAS;Lindo,Maple的历史1980年Wterloo大学符号计算研究小组成立；1985年：Maple3.3商业版发行；1992年：Windows操作系统下的MapleVRelease2面世；1994年，MapleVRelease3发行；1996年，MapleVRelease4发行；1997年，MapleVRelease5发行；1998年，MapleVRelease5.2,5.2发行；2000年，Maple6发行；2001年，Maple7发行；现在Maple10发行,Maple提供了一整套问题解决的数学环境它支持广泛、大量的数学操作如：数值分析、符号代数、图形图像。Maple是一个广泛应用的符号计算软件，它拥有强大的功能，主要表现在集符号运算、数值计算、可视化和程序设计于一体，这些功能是通过Maple提供的线性代数程序包、微积分方程程序包、统计程序包、偏微分方程程序包以及画图程序包等实现的。在和别的数学软件相比，Maple在处理微分方程中问题的功能特别好。,在线帮助Maple允许通过如下方式获得信息：在由建菜单选择HelponXXXX,或ctrl+F1得到鼠标指针位置文字的帮助。在help菜单选择Introduction使用帮助导航器,获得层级结构导航，单击感兴趣的标题直到出现相应帮助。在help菜单选择TopicSearch进行主题搜索。在工作表中可使用?命令获得明确主题的帮助页。例如?plot打开plot命令的帮助页。在help菜单选择FullTextSearch可以进行全文关键字检索。在help菜单选择History可以方便的回访已访问的帮助页。,1.2数值计算Maple可视为功能强大的计算器。请输入下面的指令：32*1213;200!;ifactor(%);expand(%);(230)/(320)*sqrt(3);evalf(%);evalf(Pi,500);Sum(1+i)/(1+i4),i=1.100);sum(1+i)/(1+i4),i=1.100);evalf(%);Sum(1/k2,k=1.infinity);value(%);Product(i2+3*i-11)/(i+3),i=0.10);value(%);evalf(%,50);,1.3代数运算展开、分解、化简表达式expr:=(x+y)15;expand(expr);factor(%);simplify(cos(x)5+sin(x)4+2*cos(x)2-2*sin(x)2-cos(2*x);normal(x3-y3)/(x2+x-y-y2);定义变量expr1:=(41*x2+x+1)2*(2*x-1);expr2:=expand(expr1);eval(expr2,x=1);top:=expr2;bottom:=expand(3*x+5)*(2*x-1);answer:=normal(top/bottom);,表达式变形my_expr:=(a*x2+b)/(x*(-3*x2-x+4);my_expr:=(a*x2+b)/(x*(-3*x2-x+4);定义函数f:=x-x2+1/2;f(2);f(a+b);使用unapply命令将表达式转化为函数g:=unapply(x2+1/2,x);解方程（组）eqn:=x3-1/2*a*x2+13/3*x2=13/6*a*x+10/3*x-5/3*a;solve(eqn,x);为验根我们计算方程在特殊点x的值eval(eqn,x=1/2*a);,定义如下方程组;eqn1:=a+2*b+3*c+4*d+5*e=41;eqn2:=5*a+5*b+4*c+3*d+2*e=20;eqn3:=3*b+4*c-8*d+2*e=125;eqn4:=a+b+c+d+e=9;用变量e来表示其他未知数a，b，c，d得到一组解solve(eqn1,eqn2,eqn3,eqn4,a,b,c,d);解不等式ineq:=x+y+4/(x+y)0,x-y0,x-yx*sin(a*x)+b*x2;#定义函数diff(f(x),x);#求导int(f(x),x);#积分int(f(x),x=1.2);#积分expr:=(2*x+3)/(7*x+5);#定义表达式Limit(expr,x=infinity);#显示求表达式的极限limit(expr,x=infinity);#给出表达式的极限expr:=sin(4*x)*cos(x):#定义表达式approx1:=series(expr,x=0);#展开为级数poly1:=convert(approx1,polynom);#转化为多项式plot(expr,poly1,x=-1.1,y=-2.2,title=cat(convert(expr,string),vs.SeriesApproximation);#比较展开式与原函数的图像,expr:=sin(4*x)*cos(x):#定义表达式approx1:=series(expr,x=0);#展开为级数poly1:=convert(approx1,polynom);#转化为多项式plot(expr,poly1,x=-1.1,y=-2.2,title=cat(convert(expr,string),vs.SeriesApproximation);#比较展开式与原函数的图像,1.6线性代数使用with(linalg)命令载入线性代数工具包restart;with(linalg);A:=matrix(3,3,1/2,-1/3,2,-5,14/3,9,0,11,-5/6);#定义矩阵det(A);#计算行列式值inverse(A);#计算逆矩阵B:=matrix(3,3,1/2,0,-2,sin(theta),1,phi2,0,phi-1,3/4);C:=multiply(A,B);#求矩阵A、B的积并存于Cdet(C);#计算行列式值eigenvals(A);#特征值eigenvects(A);#特征向量vandermonde(s,t,u,v,w);#Vandermonde矩阵A:=matrix(3,3);print(A);,1.7微分方程解微分方程命令格式:dsolve(equn,y(x);dsolve(equn,conds,y(x);equn为方程conds为条件例ode1:=diff(y(x),x)-y(x)-cos(x);ans1:=dsolve(ode1,y(x);ans2:=dsolve(ode1,y(0)=1,y(x);ode3:=diff(y(x),x)-y(x)2+y(x)*sin(x)-cos(x);ans3:=dsolve(ode1);,解微分方程组命令格式:dsolve(sysODE,ICs,funcs)sysODE为方程组ICs为条件组例restart;sys:=diff(x(t),t)=2*x(t)+y(t),diff(y(t),t)=3*x(t)+4*y(t);ans1:=dsolve(sys,x(t),y(t);ans2:=dsolve(sys,x(0)=0,y(0)=1,x(t),y(t);assign(%);x(t);plot(x(t),y(t),t=0.1);,问题：用逐次迭代法求下列初始值问题的解:解.该问题等价的积分方程为:利用Maple去进行这些重复性的迭代:y0:=1;y1:=1+int(x2+y02,x=0.x);y2:=1+int(x2+y12,x=0.x);y3:=1+int(x2+y22,x=0.x);y4:=1+int(x2+y32,x=0.x);,画向量场及积分曲线DEtoolsphaseportrait#画向量场及积分曲线；(diff(y(x),x)=-y(x),y(x),#定义微分方程y=-y;x=-2.2,#指出x的范围y(-2)=2,y(-2)=1,y(-2)=-2,#给出3个初始值dirgrid=17,17,#定义网格点密度arrows=LINE,#定义线段类型axes=NORMAL);#定义坐标系类型DEtoolsdfieldplot(diff(y(x),x)=x2-y(x),y(x),x=-2.2,y=-2.2,dirgrid=9,9,arrows=LINE,axes=NORMAL);,画向量场及积分曲线DEtoolsphaseportrait#画向量场及积分曲线；(diff(y(x),x)=x2-y(x),y(x),#定义微分方程x=-2.2,#指出x的范围y(-2)=1.3,y(-2)=-2,#给出3个初始值dirgrid=33,33,#定义网格点密度arrows=LINE,#定义线段类型axes=NORMAL);#定义坐标系类型,常微分方程的数值解原理：函数的近似展开在此基础上可以得到不同精度的算法,Euler折线及改进的Euler折线法printlev1:=0;h:=0.1;x0:=0;y0:=0.5;z0:=0.5;f1:=(x,y)-1+(y-x)2;f2:=(x,y)-2*(x-y)+2*(y-x)*(1+(y-x)2);fornfrom0to9dox|(n+1):=h*(n+1);y|(n+1):=y|n+h*f1(x|n,y|n);z|(n+1):=z|n+h*f1(x|n,z|n)+h2*f2(x|n,y|n)/2;print(x|(n+1),y|(n+1),z|(n+1);od;,整理数据data1:=x0,y0,x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7,x8,y8,x9,y9,x10,y10;data2:=x0,z0,x1,z1,x2,z2,x3,z3,x4,z4,x5,z5,x6,z6,x7,z7,x8,z8,x9,z9,x10,z10;作图比较plot(data1);plot(data2);plot(data1,data2,color=blue,red);,求精确解dsolve(diff(y(x),x)=1+(y(x)-x)3,y(0)=0.5,y(x);求精确解在个点的值w0:=0.5;h:=0.1;fornfrom0to9dox|(n+1):=h*(n+1);w|(n+1)=subs(x=x|(n+1),-1/2*(-8*x+4*x2)/(4-2*x)+1/(sqrt(4-2*x);od;整理数据data3:=x0,w0,x1,w1,x2,w2,x3,w3,x4,w4,x5,w5,x6,w6,x7,w7,x8,w8,x9,w9,x10,w10;作图比较plot(data3);plot(data1,data2,data3,color=blue,red,green);,1.8Maple编程Maple提供了简单有效的编程语言.1.8.1条件语句ifif条件then语句组fiif条件then语句组else语句组fiif条件then语句组elif条件then语句组fiif条件then语句组elif条件then语句组else语句组fix:=10:ifx2thenx2fi;ifx0theny:=1elsey:=0fi;,1.8.2循环语句for,while,dofor变量名from初值by步长to终值do语句组odwhile条件do语句组odforifrom6by10to50doprint(i)od;s:=0;forifrom11by2whilei0dot:=t2+1;ift5thenbreakfi;od;1.8.4过程Maple的绝大多数函数是以过程的形式存在，一个Maple语言的过程定义如下:lc:=proc(s,u,t,v)descriptionformalinearcombinationofthearguments;s*u+t*vendproc;lc(3,x,4,y);,L:=proc(n:nonnegint)ifn=1then1elifn=2then1elseL(n-1)+L(n-2)fiendproc:L(1);L(2);L(3);L(4);L(5);L(6);,线性系统的奇点：稳定结点,with(DEtools):DE931:=diff(x(t),t)=-x(t),diff(y(t),t)=-2*y(t);DEplot(DE931,x(t),y(t),t=-10.10,x(0)=0,y(0)=2,x(0)=0,y(0)=-2,x(0)=-2,y(0)=0,x(0)=2,y(0)=0,x(0)=1,y(0)=0.01,x(0)=1,y(0)=0.04,x(0)=1,y(0)=0.1,x(0)=1,y(0)=0.25,x(0)=1,y(0)=1,x(0)=1,y(0)=-0.01,x(0)=1,y(0)=-0.04,x(0)=1,y(0)=-0.1,x(0)=1,y(0)=-0.25,x(0)=1,y(0)=-1,x(0)=-1,y(0)=0.01,x(0)=-1,y(0)=0.04,x(0)=-1,y(0)=0.1,x(0)=-1,y(0)=0.25,x(0)=-1,y(0)=1,x(0)=-1,y(0)=-0.01,x(0)=-1,y(0)=-0.04,x(0)=-1,y(0)=-0.1,x(0)=-1,y(0)=-0.25,x(0)=-1,y(0)=-1,x=-8.8,y=-8.8,stepsize=0.05,dirgrid=21,21,color=red,linecolor=blue,axes=BOXED,title=LinearSystem:StableNode,arrows=SLIM);,线性系统的奇点：鞍点,with(DEtools):DE932:=diff(x(t),t)=-x(t),diff(y(t),t)=2*y(t);DEplot(DE932,x(t),y(t),t=-10.10,x(0)=0,y(0)=2,x(0)=0,y(0)=-2,x(0)=-2,y(0)=0,x(0)=2,y(0)=0,x(0)=1,y(0)=1,x(0)=2,y(0)=2,x(0)=3,y(0)=3,x(0)=4,y(0)=4.2,x(0)=5,y(0)=6,x(0)=1,y(0)=-1,x(0)=2,y(0)=-2,x(0)=3,y(0)=-3,x(0)=4,y(0)=-4.2,x(0)=5,y(0)=-6,x(0)=-1,y(0)=1,x(0)=-2,y(0)=2,x(0)=-3,y(0)=3,x(0)=-4,y(0)=4.2,x(0)=-5,y(0)=6,x(0)=-1,y(0)=-1,x(0)=-2,y(0)=-2,x(0)=-3,y(0)=-3,x(0)=-4,y(0)=-4.2,x(0)=-5,y(0)=-6,x=-8.8,y=-8.8,stepsize=0.05,dirgrid=21,21,color=red,linecolor=blue,axes=BOXED,title=LinearSystem:Saddle(UnstableEquilibrium),arrows=SLIM);,线性系统的奇点：稳定焦点,with(DEtools):DE933:=diff(x(t),t)=-x(t)/10-y(t),diff(y(t),t)=x(t)-y(t)/10;DEplot(DE933,x(t),y(t),t=-10.40,x(0)=8,y(0)=0,x(0)=6,y(0)=0,x=-8.8,y=-8.8,stepsize=0.05,dirgrid=21,21,color=red,linecolor=blue,axes=BOXED,title=LinearSystem:StableFocus(Spiral),arrows=SLIM);,线性系统的奇点：中心,with(DEtools):DE934:=diff(x(t),t)=-y(t),diff(y(t),t)=x(t);DEplot(DE934,x(t),y(t),t=-10.10,x(0)=1,y(0)=0,x(0)=2,y(0)=0,x(0)=3,y(0)=0,x(0)=4,y(0)=0,x(0)=5,y(0)=0,x(0)=6,y(0)=0,x(0)=7,y(0)=0,x(0)=8,y(0)=0,x=-8.8,y=-8.8,stepsize=0.05,dirgrid=21,21,color=red,linecolor=blue,axes=BOXED,title=LinearSystem:Center,arrows=SLIM);,线性系统的奇点：退化结点,with(DEtools):DE935:=diff(x(t),t)=-x(t),diff(y(t),t)=5*x(t)-y(t);DEplot(DE935,x(t),y(t),t=-10.10,x(0)=0,y(0)=2,x(0)=0,y(0)=-2,x(0)=2,y(0)=-8,x(0)=4,y(0)=-8,x(0)=6,y(0)=-8,x(0)=8,y(0)=-8,x(0)=8,y(0)=5,x(0)=-3,y(0)=-8,x(0)=-4,y(0)=13,x(0)=-8,y(0)=-8,x(0)=-8,y(0)=5,x=-8.8,y=-18.18,stepsize=0.05,dirgrid=21,21,color=red,linecolor=blue,axes=BOXED,title=LinearSystem:ImproperNode,arrows=SLIM);,稳定极限环,DEtoolsphaseportrait(diff(x(t),t)=y(t)-0.05*x(t)*(x(t)2+y(t)2-1)*(x(t)2+y(t)2-4),diff(y(t),t)=-x(t)-0.05*y(t)*(x(t)2+y(t)2-1)*(x(t)2+y(t)2-4),x(t),y(t),t=-10.10,x(0)=1,y(0)=0,x(0)=2,y(0)=0,x(0)=0.5,y(0)=0,x(0)=1.5,y(0)=0,x(0)=2.5,y(0)=0,x=-,y=-,dirgrid=17,17,stepsize=0.05,arrows=SLIM,linecolor=black);,p22:=proc(A)p:=-(A1,1+A2,2):q:=A1,1*A2,2-A1,2*A2,1:t:=p2-4*q:s1:=is(q,positive):s2:=is(q,negative):s3:=is(t,positive):s4:=is(t,0):s5:=is(t,negative):s6:=is(p,positive):s7:=is(p,0):s8:=is(p,negative):if(s1)and(s3)and(s6)thenprintf(稳定结点);elif(s1)and(s3)and(s8)thenprintf(不稳定结点);elif(s1)and(s4)and(s6)thenprintf(稳定的临界或退化结点);elif(s1)and(s4)and(s8)thenprintf(不稳定的临界或退化结点);elif(s1)and(s5)and(s6)thenprintf(稳定焦点);elif(s1)and(s5)and(s8)thenprintf(不稳定焦点);elif(s1)and(s7)thenprintf(中心);elif(s2)thenprintf(鞍点);elseprintf(无法断定);endif;end;,运行这个过程,A:=matrix(2,2,1,2,-4,3);P22(A);不稳定焦点A:=matrix(2,2,-9,2,4,-3);P22(A);不定结点A:=matrix(2,2,-1,2,4,-3);p22(A);鞍点,2.Maple应用举例2.1最速降线问题在数学史上最著名的问题之一就是最速降线问题：在给定点P(高)、Q(低)之间，找一条曲线使一质点沿该曲线无摩擦下滑的时间最短。建立坐标系求曲线y(x)y(x)光滑y(0)=0,y(a)=b;下滑时间最小先计算在小段时间再积分,总时间为测试:restart;a:=1;b:=1;y:=x;g:=9.8;dy:=diff(y,x);Times:=int(sqrt(1+dy2)/sqrt(2*g*y),x=0.a);time=0.63887restart;a:=1;b:=1;y:=sqrt(x);g:=9.8;dy:=diff(y,x);times:=int(sqrt(1+dy2)/sqrt(2*g*y),x=0.a);evalf(%);time=0.58440,restart;a:=1;b:=1;y:=sqrt(1-x2);g:=9.8;dy:=diff(y,x);times=int(sqrt(1+dy2)/sqrt(2*g*y),x=0.a);with(student);middlesum(sqrt(1+dy2)/sqrt(2*g*y),x=0.a,2000);Evalf(%);Time=0.57353,理论分析f(u)在u=0取得最小值,h(0)=0,h(a)=0,h(x)光滑.f(0)=0.经过推导后得到方程利用dsolve(y(x)*(diff(y(x),x)2+1)=c,y(x);求解后给出的结果太复杂,我们令,2.2人体含铅量问题铅是一种人体所必须的微量元素，但体内铅含量过多时就会引起铅中毒.铅是一种重金属元素，通过食物、饮料、空气进入人体，经过呼吸和消化系统后进入血液，再经过血液循环慢慢进入人体和骨头中.铅可以经过人体的排泄系统、通出出汗、剪头发、剪指甲排出体外.我们根据铅在人体内的变化情况将人体分为血液、组织、骨头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)根据前面的假设建立人体血液、组织和骨头中含铅量所满足的微分方程组,并求其通解.(2)取时间单位为天，对一个志愿者在一种环境中的生活情况进行测定得a21=0.011，a12=0.012,a31=0.0039，a13=0.000035,a01=0.021,a02=0.016,L=49.3.假设该志愿者开始时体内的含铅量为0，求他体内血液、组织和骨头中含铅量随时间变化的关系，画出这些解曲线的图，并求当x+时这些函数的极限.(3)假设该志愿者在此环境中生活了365天后搬到了一个无铅的环境中去(不再有外界的铅进入体内，即L=0)，再讨论他体内血液、组织和骨头中含铅量随时间变化的关系，并画出0t1460时这些含铅量曲线的图形.,铅是一种重金属元素，在尘土、空中飘浮物、汽车尾气、油漆、一些彩釉陶瓷制品以及塑料所含的增塑剂中含量较高，主要通过口腔-消化道进入生物体，对其神经系统和血液系统有较强的毒性，会破坏其正常功能，其中以神经毒性为主。众所周知，人的思想和行为受神经系统控制，一旦神经系统遭到破坏，其后果是可想而知的。儿童正处于神经系统发育完善期，是铅污染的最大受害者。儿童好奇心强，喜欢用手去触摸各种物品，而且一般有“手-口接触（吮吸）”频繁的行为特点，这就给铅进入体内创造了条件。由于儿童正处于生长发育阶段，许多器官的功能都不完善，大脑的发育尚未完成，“血-脑屏障（一种保护性组织，有防止血液中物质随意进入脑细胞的功能）”不能阻止铅离子（2+）进入大脑。虽然，环境因素造成的铅危害一般不会达到工业性中毒的程度，但低量、长期的铅接触会造成铅在体内积累。虽然血铅含量超标并不等于铅中毒，但研究表明，过量的铅接触会造成神经元能量代谢降低，使儿童，尤其是学龄前儿童的神经发育指数达不到应有的标准，造成记忆力和思维能力下降，影响智力的发展这种损害是不可逆的。,儿童体内的铅水平可分从血、骨、齿、尿、发检测到，儿童血铅水平分为5级。级：血铅值低于10微克分升，身体处于相对安全状态；级：血铅值为10微克19微克分升，属于轻度铅中毒。影响造血、神经传导和认知能力，儿童易出现头昏、烦躁、注意力涣散、多动、厌食、腹胀、轻度贫血；级：血铅值达到20微克44微克分升，为中度铅中毒。引起缺钙、缺锌、缺铁、免疫力低下，运动不协调，视力和听力受损，学习困难、智力下降，生长发育迟缓，贫血、腹绞痛、食癖、反应迟钝等；级：血铅值达到45微克69微克分升，就是重度铅中毒。可出现性格改变、易激怒、攻击性行为、运动失调、贫血、腹绞痛、高血压、心律失常和痴呆等；级：当血铅值大于70微克分升时，为极重度铅中毒，可导致脏器损害、铅性脑病、瘫痪、昏迷等。,铅中毒会引起神经、消化、循环系统紊乱，表现为贫血、腹痛、高血压等。血铅高组的儿童的总智商、操作智商、语言智商分别比低血铅组落后分、分和分，而每升血液中的铅浓度上升微克，儿童的身高将降低厘米。低剂量的铅接触可以对人体的红细胞、肾脏、免疫系统、骨髓和中枢神经的功能产生不良影响，而所有这些影响发生前都可能没有明显的临床症状。铅中毒会影响婴幼儿最初站立、行走和说话的年龄，也可能引起孩子注意力涣散、记忆力减退、理解力降低及学习困难等。铅通过各种方式进入人体后，可使人的智力下降，学习、工作成绩低落；蓄积到一定程度时会使人出现精神障碍、噩梦、失眠、头痛等慢性中毒症状；严重者还可有乏力、食欲不振、恶心、腹胀、腹痛或腹泻等。,图,简化图,模型,模型x=Ax+b,数据,时间分为两段: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:=matrix(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);,结果太复杂,利用迭加原理、特征值和特征向量法#非齐次方程的特解with(linalg):with(plots):xe:=evalm(-(inverse(A),#得到解的表达式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;#在解的表达式用t=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=infinity);#得到365天后解的表达式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=365时的值x1365:=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;