Maple 简介

1. Maple 简介 周义仓 西安交通大学理学院应用数学系Tel. 82668741（O）中2-2210 Email: zhouyc@mail.xjtu.edu.cn

2. 主要内容和参考文献 • 1.1 主要内容 • 参考资料 • Maple 软件包介绍 • 引言，数值计算 ，代数运算 ， • 图形 ，微积分，线 性 代 数 ， • 微分方程， Maple 编程 • Maple 应用举例 • 最速降线问题，人体含铅量

3. 1.2 参考资料 • Maple与微分方程（复印）； • 周义仓、靳祯、秦军林，常微分方程及其应用，科学出版社，2003年6月。 • 精英科技，孙非，Maple6实例教程，中国电力出版社，2001。 • 李强，Maple8基础教程，中国水利水电出版社，2004 • David Barrow等， Solving Differential Equations with Maple V, Brooks/Cole Publishing Company, 1998. • http://202.112.84.70/resource/maple.htm • http://202.204.240.30/maple/tour/index.htm

4. Maple9教程“四部曲” • Maple 9 Getting Started Guide , 834KB http://bolide.lamost.org/sheaweb/files/m9GettingStartedGuide.pdf • Maple 9 Learning Guide, 4.49MB http://bolide.lamost.org/sheaweb/files/m9LearningGuide.pdf • Maple 9 Introductory Programming Guide, 7.08MB http://bolide.lamost.org/sheaweb/files/m9IntroductoryProgramming • Guide.pdf 4. Maple 9 Advanced Programming Guide , • 3.63MB http://bolide.lamost.org/sheaweb/files/m9AdvancedProgrammingGuide.pdf • Maple应用 • http://www.maplesoft.com/applications/index.aspx

5. Maple 软件包介绍 • 2.1 引言 • Maple是Waterloo大学开发的集数值计算、符号运算和图形显示于一体的软件包。它可以把人们从烦琐的代数运算中解脱出来，使得我们将主要精力集中与分析和解决问题的思路与方法中去。 • 数学运算分为数值运算、符号运算、逻辑运算等。 • 数值运算： 如函数求值、方程求根、矩阵的特征值、特征向量等，数值计算是一项简单烦琐的工作，计算机的出现解决了数值计算的困难，使得大规模、高速度的计算可以实现，使得数学在天气预报、油藏模拟、航天等领域发挥了巨大的作用。

6. 符号运算（代数运算）：这是一种更加智能化的计算，所处理的符号可以代表各种数、代数式、数学结构（集合），代数运算是在一些代数规则下所进行的数学处理，主要是寻求一个简单完美的公式解答。符号运算（代数运算）：这是一种更加智能化的计算，所处理的符号可以代表各种数、代数式、数学结构（集合），代数运算是在一些代数规则下所进行的数学处理，主要是寻求一个简单完美的公式解答。 用计算机进行符号运算是数学与计算机领域中一个新的发展方向，是人们用计算机代替大脑进行数值运算后用计算机代替大脑进行代数运算的成功实践，是数学走向机械化的重要步骤，这也使得计算机更加智能化。 目前数学软件包很多，如 Mathematica Maple Matlab MathCAD SAS Lindo ……

7. Maple 的历史 • 1980年Wterloo大学符号计算研究小组成立； • 1985年：Maple3.3商业版发行； • 1992年：Windows操作系统下的Maple V Release2面世； • 1994年，Maple V Release3发行； • 1996年，Maple V Release4发行； • 1997年，Maple V Release5发行； • 1998年，Maple V Release5.2,5.2发行； • 2000年，Maple 6发行； • 2001年，Maple 7发行； • 200?年，Maple 8发行； • 200?年，Maple 9发行； • 200?年，Maple 9.5发行；

8. Maple提供了一整套问题解决的数学环境 它支持广泛、大量的数学操作如：数值分析、符号代数、图 形图像。 Maple是一个广泛应用的符号计算软件，它拥有强大的功能，主要表现在集符号运算、数值计算、可视化和程序设计于一体，这些功能是通过Maple提供的线性代数程序包、微积分方程程序包、统计程序包、偏微分方程程序包以及画图程序包等实现的。 在和别的数学软件相比，Maple在处理微分方程中问题的功能特别好。

9. 在线帮助 • Maple 允许通过如下方式获得信息 ： • 在由建菜单选择 Help on XXXX,或 ctrl+F1得到鼠标指针位置文 字的帮助 。 • 在 help菜单选择 Introduction使用帮助导航器, 获得层级结构导航，单击感兴趣的标题直到出现相应帮助。 • 在help菜单选择Topic Search进行主题搜索。 • 在工作表中可使用?命令获得明确主题的帮助页。 • 例如 ?plot打开plot命令 的帮助页。 • 在 help菜单选择Full Text Search可以进行全文关键字检索。 • 在 help菜单选择History可以方便的回访已访问 的帮助页。

10. 2.2 数值计算 Maple可视为功能强大的计算器。请输入下面的指令： 32*12^13; 200!; ifactor(%); expand(%); (2^30)/(3^20)*sqrt(3); evalf(%); evalf(Pi,500); Sum((1+i)/(1+i^4),i=1..100); sum((1+i)/(1+i^4),i=1..100); evalf(%); Sum( 1/k^2, k=1..infinity ); value(%); Product( ((i^2+3*i-11)/(i+3)), i=0..10 ); value(%); evalf( % , 50 );

11. 2.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( (x^3-y^3)/(x^2+x-y-y^2) ); • 定义变量 • expr1 := (41*x^2+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 );

12. 表达式变形 • my_expr := (a*x^2+b)/(x*(-3*x^2-x+4)); • my_expr := (a*x^2+b)/(x*(-3*x^2-x+4)); • 定义函数 • f := x -> x^2+1/2 ; • f(2); • f(a+b); • 使用 unapply命令将表达式转化为函数 • g := unapply(x^2 + 1/2, x); • 解方程（组） • eqn:=x^3-1/2*a*x^2+13/3*x^2=13/6*a*x+10/3*x-5/3*a; • solve( eqn, {x} ); • 为验根我们计算方程在特殊点x的值 • eval( eqn , x=1/2*a );

13. 定义如下方程组; • 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) < 10: • solve( ineq, {x} );

14. 2.4.图 形 • 可以对显式、隐式 、参数型函数及数据集作图。 • 首先 用 with命令载入两个作图工具包, 列出可用的函数。 • with(plots); with(plottools); • 2D作图举例 • plot( tan(x), x=-2*Pi..2*Pi, y=-4..4); • 隐函数作图 • implicitplot( { x^2+y^2=1, y=exp(x) }, x=-Pi..Pi, y=-Pi..Pi,scaling=CONSTRAINED ); • 线性不等数组的图解 • inequal( { x+y > 0, x-y <= 1, y = 2 }, x=-3..3, y=-3..3,optionsfeasible=(color=red), optionsopen=(color=blue,thickness=2), optionsclosed=(color=green, thickness=3),optionsexcluded=(color=yellow) );

15. 极坐标 • plot(cos(16*t),t=-Pi..Pi,coords=polar,scaling=constrained); • plot(sin(t^(-2)),t=-5*Pi..5*Pi,coords=polar); • plot(cos(sin(t^(-5))),t=-10*Pi..10*Pi,coords=polar); • 3D图象 • plot3d( x*exp(-x^2-y^2), x=-2..2, y=-2..2, axes=BOXED,title=`A Surface Plot`); • 球坐标 • plot3d(1,t=0..2*Pi,p=0..Pi,coords=spherical); • plot3d(1,t=0..2*Pi,p=0..Pi,coords=spherical, scaling=constrained); • plot3d(x*sin(y),x=-1..3*Pi,y=0..2*Pi,coords=spherical); • 动画 • plots[animate](x*t,x=-1..1,t=1..30,numpoints=5,frames=100); • animate3d(cos(t*x)*sin(t*y),x=-Pi..Pi,y=-Pi..Pi,t=1..2000); • animate([(1+cos(n*t/180*Pi))*cos(n*t/180*Pi), (1+cos(n*t/180*Pi))*sin(n*t/180*Pi),t=0..1],n=1..360);

16. 2.5 微积分 f := x -> x*sin(a*x) + b*x^2; #定义函数 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. Series Approximation" ) ); #比较展开式与原函数的图像

17. 2.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,phi^2, 0,phi-1, 3/4] ); C := multiply( A, B ); #求矩阵 A、B的积并存于C det(C); #计算行列式值 eigenvals(A); #特征值 eigenvects(A); #特征向量 vandermonde( [s, t, u, v, w] ); #Vandermonde矩阵 A := matrix( 3, 3 ); print(A);

18. 2.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);

19. 解微分方程组 • 命令格式: • 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);

20. 第21页 用逐次迭代法求下列初始值问题的解: 解. 该问题等价的积分方程为: 利用Maple去进行这些重复性的迭代: y0:=1; y1:=1+int(x^2+y0^2,x=0..x); y2:=1+int(x^2+y1^2,x=0..x); y3:=1+int(x^2+y2^2,x=0..x); y4:=1+int(x^2+y3^2,x=0..x);

21. 第25页 DEtools[phaseportrait] #画向量场及积分曲线； ([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); #定义坐标系类型 DEtools[dfieldplot] ([diff(y(x),x)=x^2-y(x)],y(x), x=-2..2,y=-2..2, dirgrid=[9,9], arrows=LINE, axes=NORMAL);

22. 第26页 DEtools[phaseportrait] #画向量场及积分曲线； ([diff(y(x),x)= x^2-y(x)],y(x), #定义微分方程 x=-2..2, #指出x的范围 [[y(-2)=1.3],[y(-2)=-2]], #给出3个初始值 dirgrid=[33,33], #定义网格点密度 arrows=LINE, #定义线段类型 axes=NORMAL); #定义坐标系类型

23. 第88-91: 常微分方程的数值解 原理：函数的近似展开 在此基础上可以得到不同精度的算法

24. 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); for n from 0 to 9 do x||(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)+h^2*f2(x||n,y||n)/2; print (x||(n+1),y||(n+1),z||(n+1)); od;

25. ＃整理数据 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]);

26. ＃求精确解 dsolve({diff(y(x),x)=1+(y(x)-x)^3,y(0)=0.5},y(x)); ＃求精确解在１０个点的值 w0:=0.5;h:=0.1; for n from 0 to 9 do x||(n+1):=h*(n+1); w||(n+1)=subs(x=x||(n+1), -1/2*(-8*x+4*x^2)/(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]);

27. 2.8 Maple 编程 Maple提供了简单有效的编程语言. 2.8.1 条件语句 if if 条件 then 语句组 fi if 条件 then 语句组 else 语句组 fi if 条件 then 语句组 elif 条件then 语句组 fi if 条件 then 语句组 elif 条件 then 语句组else 语句组 fi x:=10:if x>2 then x^2 fi; if x<0 then y:=-1 elif x>0 then y:=1 else y:=0 fi;

28. 2.8.2 循环语句for,while,do for 变量名 from 初值 by 步长 to 终值 do 语句组 od while 条件 do 语句组 od for i from 6 by 10 to 50 do print(i) od; s := 0; for i from 11 by 2 while i < 18 do s := s+ i; od;

29. 2.8.3 控制流程语句 break t:=1: while t>0 do t:=t^2+1; if t>5 then break fi; od; 2.8.4过程 Maple的绝大多数函数是以过程的形式存在，一个Maple语言的过程定义如下: lc := proc( s, u, t, v ) description "form a linear combination of the arguments"; s * u + t * v end proc; lc(3,x,4,y);

30. L:=proc(n::nonnegint) if n=1 then 1 elif n=2 then 1 else L(n-1)+L(n-2) fi end proc: L(1); L(2); L(3); L(4); L(5); L(6);

31. p22:=proc(A) p:=-(A[1,1]+A[2,2]): q:=A[1,1]*A[2,2]-A[1,2]*A[2,1]: t:=p^2-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))then printf("稳定结点"); elif((s1)and(s3)and(s8))then printf("不稳定结点"); elif((s1)and(s4)and(s6))then printf("稳定的临界或退化结点"); elif((s1)and(s4)and(s8))then printf("不稳定的临界或退化结点"); elif((s1)and(s5)and(s6))then printf("稳定焦点"); elif((s1)and(s5)and(s8))then printf("不稳定焦点"); elif((s1)and(s7))then printf("中心"); elif(s2)then printf("鞍点"); else printf("无法断定"); end if; end;

32. 3.Maple 应用举例 3.1最速降线问题 在数学史上最著名的问题之一就是最速降线问题：在给定 点 P(高 )、Q(低 )之间，找一条曲线使一质点沿该曲线无摩擦下滑 的时间最短 。 建立坐标系 求曲线y(x) y(x)光滑 y(0)=0, y(a)=b; 下滑时间最小 先计算在 小段时间 再积分

33. 总时间为 测试: restart;a:=1;b:=1;y:=x;g:=9.8; dy:=diff(y,x); times=int(sqrt(1+dy^2)/sqrt(2*g*y),x=0..a); Time=0.63887 restart;a:=1;b:=1;y:=sqrt(x);g:=9.8; dy:=diff(y,x); times=int(sqrt(1+dy^2)/sqrt(2*g*y),x=0..a); evalf(%); Time=0.58440

34. restart;a:=1;b:=1;y:=sqrt(1-x^2);g:=9.8; dy:=diff(y,x); times=int(sqrt(1+dy^2)/sqrt(2*g*y),x=0..a); with(student); middlesum(sqrt(1+dy^2)/sqrt(2*g*y),x=0..a,2000); Evalf(%); Time=0.57353

35. 理论分析 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)); 求解后给出的结果太复杂,我们令

36. 3.2 人体含铅量（第260页 19题） 铅是一种人体所必须的微量元素，但体内铅含量过多时就会引起铅中毒.铅是一种重金属元素，通过食物、饮料、空气进入人体，经过呼吸和消化系统后进入血液，再经过血液循环慢慢进入人体和骨头中.铅可以经过人体的排泄系统、通出出汗、剪头发、剪指甲排出体外. 我们根据铅在人体内的变化情况将人体分为血液、组织、骨头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)，再讨论他体内血液、组织和骨头中含铅量随时间变化的关系，并画出0≤t≤1460时这些含铅量曲线的图形.

37. 铅是一种重金属元素，在尘土、空中飘浮物、汽车尾气、油漆、一些彩釉陶瓷制品以及塑料所含的增塑剂中含量较高，主要通过口腔---消化道进入生物体，对其神经系统和血液系统有较强的毒性，会破坏其正常功能，其中以神经毒性为主。铅是一种重金属元素，在尘土、空中飘浮物、汽车尾气、油漆、一些彩釉陶瓷制品以及塑料所含的增塑剂中含量较高，主要通过口腔---消化道进入生物体，对其神经系统和血液系统有较强的毒性，会破坏其正常功能，其中以神经毒性为主。 众所周知，人的思想和行为受神经系统控制，一旦神经系统遭到破坏，其后果是可想而知的。儿童正处于神经系统发育完善期，是铅污染的最大受害者。儿童好奇心强，喜欢用手去触摸各种物品，而且一般有“手--口接触（吮吸）”频繁的行为特点，这就给铅进入体内创造了条件。由于儿童正处于生长发育阶段，许多器官的功能都不完善，大脑的发育尚未完成，“血---脑屏障（一种保护性组织，有防止血液中物质随意进入脑细胞的功能）”不能阻止铅离子（Ｐｂ2+）进入大脑。虽然，环境因素造成的铅危害一般不会达到工业性中毒的程度，但低量、长期的铅接触会造成铅在体内积累。虽然血铅含量超标并不等于铅中毒，但研究表明，过量的铅接触会造成神经元能量代谢降低，使儿童，尤其是学龄前儿童的神经发育指数达不到应有的标准，造成记忆力和思维能力下降，影响智力的发展——这种损害是不可逆的。

38. 儿童体内的铅水平可分从血、骨、齿、尿、发检测到，儿童血铅水平分为5级。Ⅰ级：血铅值低于10微克／分升，身体处于相对安全状态；Ⅱ级：血铅值为10微克－19微克／分升，属于轻度铅中毒。影响造血、神经传导和认知能力，儿童易出现头昏、烦躁、注意力涣散、多动、厌食、腹胀、轻度贫血；Ⅲ级：血铅值达到20微克－44微克／分升，为中度铅中毒。引起缺钙、缺锌、缺铁、免疫力低下，运动不协调，视力和听力受损，学习困难、智力下降，生长发育迟缓，贫血、腹绞痛、食癖、反应迟钝等；Ⅳ级：血铅值达到45微克－69微克／分升，就是重度铅中毒。可出现性格改变、易激怒、攻击性行为、运动失调、贫血、腹绞痛、高血压、心律失常和痴呆等；Ⅴ级：当血铅值大于70微克／分升时，为极重度铅中毒，可导致脏器损害、铅性脑病、瘫痪、昏迷等。

39. 铅中毒会引起神经、消化、循环系统紊乱，表现为贫血、腹痛、高血压等。血铅高组的儿童的总智商、操作智商、语言智商分别比低血铅组落后１４分、１４分和１３分，而每升血液中的铅浓度上升１００微克，儿童的身高将降低１．３厘米。铅中毒会引起神经、消化、循环系统紊乱，表现为贫血、腹痛、高血压等。血铅高组的儿童的总智商、操作智商、语言智商分别比低血铅组落后１４分、１４分和１３分，而每升血液中的铅浓度上升１００微克，儿童的身高将降低１．３厘米。 低剂量的铅接触可以对人体的红细胞、肾脏、免疫系统、骨髓和中枢神经的功能产生不良影响，而所有这些影响发生前都可能没有明显的临床症状。铅中毒会影响婴幼儿最初站立、行走和说话的年龄，也可能引起孩子注意力涣散、记忆力减退、理解力降低及学习困难等。 铅通过各种方式进入人体后，可使人的智力下降，学习、工作成绩低落；蓄积到一定程度时会使人出现精神障碍、噩梦、失眠、头痛等慢性中毒症状；严重者还可有乏力、食欲不振、恶心、腹胀、腹痛或腹泻等。

40. 图

41. 简化图

42. 模型

43. 模型 x' = Ax + b,

44. 数据

45. 时间分为两段: [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)});

46. 结果太复杂,利用迭加原理、特征值和特征向量法结果太复杂,利用迭加原理、特征值和特征向量法 #非齐次方程的特解 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]; #特征向量组成矩阵 P:=augment(v1,v2,v3);

47. #得到解的表达式 x1:=c[1]*P[1,1]*exp(lambda[1]*t)+c[2]*P[1,2]*exp(lambda[2]*t)+c[3]*P[1,3]*exp(lambda[3]*t)+xe[1,1]; x2 := c[1]*P[2,1]*exp(lambda[1]*t)+c[2]*P[2,2]*exp(lambda[2]*t)+c[3]*P[2,3]*exp(lambda[3]*t)+xe[2,1]; x3 := c[1]*P[3,1]*exp(lambda[1]*t)+c[2]*P[3,2]*exp(lambda[2]*t)+c[3]*P[3,3]*exp(lambda[3]*t)+xe[3,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},{c[1],c[2],c[3]}); assign(%);

48. #作图 plot([x1, x2, x3], t=0..365, color=[red,blue,green], thickness=2); plot1:=%:

49. #重新显示解 x1;x2;x3; #求极限 limit(x1,t=infinity); limit(x2,t=infinity); limit(x3,t=infinity); #得到365天后解的表达式 xx1:=cc[1]*P[1,1]*exp(lambda[1]*t)+cc[2]*P[1,2]*exp(lambda[2]*t)+cc[3]*P[1,3]*exp(lambda[3]*t); xx2 := cc[1]*P[2,1]*exp(lambda[1]*t)+cc[2]*P[2,2]*exp(lambda[2]*t)+cc[3]*P[2,3]*exp(lambda[3]*t); xx3 := cc[1]*P[3,1]*exp(lambda[1]*t)+cc[2]*P[3,2]*exp(lambda[2]*t)+cc[3]*P[3,3]*exp(lambda[3]*t);

50. # 计算解在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}, {cc[1],cc[2],cc[3]}); assign(%); #重新显示解 xx1;xx2;xx3; #求极限 limit(xx1,t=infinity); limit(xx2,t=infinity); limit(xx3,t=infinity);