Maple_与微分方程模型

1、Maple Maple 与微分方程模型与微分方程模型主要内容 学习、了解数学软件Maple； 用Maple解决常微分方程问题； 对一些实际问题建立微分方程模 型， 用Maple软件求解； 一些实际问题的动画演示 参考资料周义仓、靳祯、秦军林，常微分方程及其应用，科学出版社，2003年6月。李强，Maple8基础教程，中国水利水电出版社，2004David Barrow等， Solving Differential Equations with Maple V, Brooks/Cole Publishing Company, 1998.0/resource/m

2、aple.htm0/maple/tour/index.htmMaple9教程“四部曲” Maple 9 Getting Started Guide , 834KB /sheaweb/files/m9GettingStartedGuide.pdfMaple 9 Learning Guide, 4.49MB /sheaweb/files/m9LearningGuide.pdf Maple 9 Introductory Programming Guide, 7.08MB

3、/sheaweb/files/m9IntroductoryProgrammingGuide.pdf 4. Maple 9 Advanced Programming Guide , 3.63MB /sheaweb/files/m9AdvancedProgrammingGuide.pdfMaple应用http:/ Maple 软件包介绍软件包介绍1.1 引言引言 Maple是是Waterloo大学开发的集数值计算、符号运算和图大学开发的集数值计算、符号运算和图形显示于一体的软件包。它可以把人们从烦琐的代数运

4、算形显示于一体的软件包。它可以把人们从烦琐的代数运算中解脱出来，使得我们将主要精力集中与分析和解决问题中解脱出来，使得我们将主要精力集中与分析和解决问题的思路与方法中去。的思路与方法中去。数学运算分为数值运算、符号运算、逻辑运算等。数学运算分为数值运算、符号运算、逻辑运算等。 数值运算：数值运算： 如函数求值、方程求根、矩阵的特征值、特如函数求值、方程求根、矩阵的特征值、特征向量等，数值计算是一项简单烦琐的工作，计算机的出征向量等，数值计算是一项简单烦琐的工作，计算机的出现解决了数值计算的困难，使得大规模、高速度的计算可现解决了数值计算的困难，使得大规模、高速度的计算可以实现，使得数学在天气预

5、报、油藏模拟、航天等领域发以实现，使得数学在天气预报、油藏模拟、航天等领域发挥了巨大的作用。挥了巨大的作用。 符号运算符号运算 一种更加智能化的计算，所处理的符号可以代表各种数、代一种更加智能化的计算，所处理的符号可以代表各种数、代数式、数学结构（集合），代数运算是在一些代数规则下所进数式、数学结构（集合），代数运算是在一些代数规则下所进行的数学处理，主要是寻求一个简单完美的公式解答。行的数学处理，主要是寻求一个简单完美的公式解答。 用计算机进行符号运算是数学与计算机领域中一个新的发展用计算机进行符号运算是数学与计算机领域中一个新的发展方向，是人们用计算机代替大脑进行数值运算后用计算机代替方向

6、，是人们用计算机代替大脑进行数值运算后用计算机代替大脑进行代数运算的成功实践，是数学走向机械化的重要步骤，大脑进行代数运算的成功实践，是数学走向机械化的重要步骤，这也使得计算机更加智能化。这也使得计算机更加智能化。 目前数学软件包很多，如目前数学软件包很多，如 Mathematica; Maple; Matlab; MathCAD; SAS; LindoMaple 的历史的历史1980年年Wterloo大学符号计算研究小组成立；大学符号计算研究小组成立；1985年：年：Maple3.3商业版发行；商业版发行；1992年：年：Windows操作系统下的操作系统下的Maple V Release2

7、面世；面世；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发行；发行； 现在 Maple 10发行发行 Maple提供了一整套问题解决的数学环境 它支持广泛、大量的数学操作如：数值分析、符号代数、图 形图像。 Maple是一个广泛应用的符号计算软件，它拥有强大的功能，主要表现在集符号运算、数值计算、可视化和程序设

8、计于一体，这些功能是通过Maple提供的线性代数程序包、微积分方程程序包、统计程序包、偏微分方程程序包以及画图程序包等实现的。 在和别的数学软件相比，Maple在处理微分方程中问题的功能特别好。在线帮助在线帮助 Maple 允许通过如下方式获得信息 ： 在由建菜单选择 Help on XXXX,或 ctrl+F1得到鼠标指针位置文 字的帮助 。 在 help菜单选择 Introduction使用帮助导航器, 获得层级结构导航，单击感兴趣的标题直到出现相应帮助。 在help菜单选择Topic Search 进行主题搜索。在工作表中可使用?命令获得明确主题的帮助页。例如 ?plot 打开plot命

9、令 的帮助页。 在 help菜单选择Full Text Search可以进行全文关键字检索。 在 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(%); Produc

10、t( (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

11、, 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

12、/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

13、-y 0, x-y x*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, polyn

14、om ); #转化为多项式 plot( expr, poly1, x=-1.1, y=-2.2, title = cat( convert(expr, string), vs. Series Approximation ) ); #比较展开式与原函数的图像 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(

15、expr, string), vs. Series Approximation ) ); #比较展开式与原函数的图像 := approx14 x383x342130 x5()O x6 := poly14 x383x342130 x5 := approx14 x383x342130 x5100391260 x724660190720 x961256599979200 x11()O x12 := poly14 x383x342130 x5100391260 x724660190720 x961256599979200 x111.6 线线 性性 代代 数数 使用with(linalg)命令载入线 性

16、代数工具包 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的积并存于C det(C); #计算行列式值 eigenvals(A); #特征值 eigenvects(A); #特征向量 vandermonde( s

17、, 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);解微分方程组命令

18、格式: 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+y0

19、2,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);1)0(,13yyyxdxxyxy03)(1 1)( := y01 := y1113x3x := y21163x7215x516x423x3x2xy31159535x15412285x1312268x1218451975x114525x1029911340 x941630 x847315x7 := 2990 x6815x556x443x3x2xy41x243x3232315x71615x5910 x676x48699832432430

20、0 x134937794638512875x152622593742200 x12 := 20093155925x111264156700 x10401911340 x913332520 x85800099340540200 x145244178513505000 x2044916481249567318000 x192307204145972927000 x1845901107270163000 x2423505261304776501875x234282079722431710000 x2215983441804526222500 x21842329638512875x1713357767

21、140864824000 x1611890355320 x2811763214119435204875x2733115846728625x26653698348271573350000 x25821210236775x291109876902975x31x画向量场及积分曲线画向量场及积分曲线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=

22、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)=

23、1.3,y(-2)=-2, #给出给出3个初始值个初始值dirgrid=33,33, #定义网格点密度定义网格点密度arrows=LINE, #定义线段类型定义线段类型axes=NORMAL); #定义坐标系类型定义坐标系类型常微分方程的数值解原理：函数的近似展开在此基础上可以得到不同精度的算法)()( 2)( )()(202000hoxfhxhfxfhxfEuler折线及改进的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

24、 n from 0 to 9 dox|(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;作图比

25、较作图比较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;for n from 0 to 9 dox|(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,x

26、5,w5,x6,w6,x7,w7,x8,w8,x9,w9,x10,w10;作图比较作图比较 plot(data3); plot(data1,data2,data3,color=blue, red,green);1.8 Maple 编程Maple提供了简单有效的编程语言.1.8.1 条件语句 if if 条件 then 语句组 fiif 条件 then 语句组 else 语句组 fiif 条件 then 语句组 elif 条件then 语句组 fiif 条件 then 语句组 elif 条件 then 语句组else 语句组 fix:=10:if x2 then x2 fi; if x0 the

27、n y:=1 else y:=0 fi;1.8.2 循环语句for,while,do for 变量名 from 初值 by 步长 to 终值 do 语句组 odwhile 条件 do 语句组 od for i from 6 by 10 to 50 do print(i) od; s := 0; for i from 11 by 2 while i 0 dot:=t2+1;if t5 then break fi; od;1.8.4过程 Maple的绝大多数函数是以过程的形式存在，一个Maple语言的过程定义如下:lc := proc( s, u, t, v ) description form

28、a linear combination of the arguments; s * u + t * vend proc;lc(3,x,4,y);L:=proc(n:nonnegint) if n=1 then 1 elif n=2 then 1 else L(n-1)+L(n-2)fiend proc: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

29、(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

30、(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= Linear System: Stable Node ,arrows=SLIM);线性系统的奇点：鞍点with(DEtools): DE932 :=diff(x(t),t)=-x(t),diff(y(t),t)=

31、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,

32、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=Linear System : Saddle(UnstableEquilibrium), arrows=SLIM);线性系统的奇点：稳定焦点with(DEtools):

33、 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=Linear System : Stable Focus(Spiral) , arrows=SLIM);线性系统的奇点：中心with(DEtools):DE934 :=di

34、ff(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=Linear System : Center ,arro

35、ws=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

36、=-8.8,y=-18.18,stepsize=0.05, dirgrid=21,21,color=red, linecolor=blue,axes=BOXED,title=Linear System : Improper Node ,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

37、, 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

38、(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)

39、then printf(不稳定焦点);elif(s1)and(s7)then printf(中心);elif(s2)then printf(鞍点);else printf(无法断定);end if;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(低 )之间，找一条曲线使一质点沿该曲线无摩擦下

40、滑 的时间最短 。建立坐标系求曲线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)(2)( 12xgyxydtdxxgyxy

41、xyTa02)(2)( 1)(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);求解后给出的结果太复

42、杂求解后给出的结果太复杂,我们令我们令dxxuhxygxuhxyxuhxyTufa02)()(2)( )( (1)()()(bayycxyxy)(, 0)0(,)( 1)(2)2cos(1 (2)cos()cot(1)(),cot()( 22tctctcxytxydttctdtttcxydydx)2cos(1 ()cot()cos()sin(2)( )cos(1 (),sin(tcyttcx2.2 人体含铅量问题 铅是一种人体所必须的微量元素，但体内铅含量过多时就会引起铅中毒.铅是一种重金属元素，通过食物、饮料、空气进入人体，经过呼吸和消化系统后进入血液，再经过血液循环慢慢进入人体和骨头中.铅

43、可以经过人体的排泄系统、通出出汗、剪头发、剪指甲排出体外. 我们根据铅在人体内的变化情况将人体分为血液、组织、骨头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)根据前面的假设建立人体血液、组织和骨头中含铅量所满足的

44、微分方程组,并求其通解.(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时这些含铅量曲线的图形. 铅是一种重金属元素，在尘

45、土、空中飘浮物、汽车尾气、油漆、一些彩釉陶瓷制品以及塑料所含的增塑剂中含量较高，主要通过口腔-消化道进入生物体，对其神经系统和血液系统有较强的毒性，会破坏其正常功能，其中以神经毒性为主。 众所周知，人的思想和行为受神经系统控制，一旦神经系统遭到破坏，其后果是可想而知的。儿童正处于神经系统发育完善期，是铅污染的最大受害者。儿童好奇心强，喜欢用手去触摸各种物品，而且一般有“手-口接触（吮吸）”频繁的行为特点，这就给铅进入体内创造了条件。由于儿童正处于生长发育阶段，许多器官的功能都不完善，大脑的发育尚未完成，“血-脑屏障（一种保护性组织，有防止血液中物质随意进入脑细胞的功能）”不能阻止铅离子（2+）

46、进入大脑。虽然，环境因素造成的铅危害一般不会达到工业性中毒的程度，但低量、长期的铅接触会造成铅在体内积累。虽然血铅含量超标并不等于铅中毒，但研究表明，过量的铅接触会造成神经元能量代谢降低，使儿童，尤其是学龄前儿童的神经发育指数达不到应有的标准，造成记忆力和思维能力下降，影响智力的发展这种损害是不可逆的。 儿童体内的铅水平可分从血、骨、齿、尿、发检测到，儿童血铅水平分为5级。 级：血铅值低于10微克分升，身体处于相对安全状态； 级：血铅值为10微克19微克分升，属于轻度铅中毒。影响造血、神经传导和认知能力，儿童易出现头昏、烦躁、注意力涣散、多动、厌食、腹胀、轻度贫血； 级：血铅值达到20微克44

47、微克分升，为中度铅中毒。引起缺钙、缺锌、缺铁、免疫力低下，运动不协调，视力和听力受损，学习困难、智力下降，生长发育迟缓，贫血、腹绞痛、食癖、反应迟钝等； 级：血铅值达到45微克69微克分升，就是重度铅中毒。可出现性格改变、易激怒、攻击性行为、运动失调、贫血、腹绞痛、高血压、心律失常和痴呆等； 级：当血铅值大于70微克分升时，为极重度铅中毒，可导致脏器损害、铅性脑病、瘫痪、昏迷等。 铅中毒会引起神经、消化、循环系统紊乱，表现为贫血、腹痛、高血压等。血铅高组的儿童的总智商、操作智商、语言智商分别比低血铅组落后分、分和分，而每升血液中的铅浓度上升微克，儿童的身高将降低厘米。 低剂量的铅接触可以对人体

48、的红细胞、肾脏、免疫系统、骨髓和中枢神经的功能产生不良影响，而所有这些影响发生前都可能没有明显的临床症状。铅中毒会影响婴幼儿最初站立、行走和说话的年龄，也可能引起孩子注意力涣散、记忆力减退、理解力降低及学习困难等。 铅通过各种方式进入人体后，可使人的智力下降，学习、工作成绩低落；蓄积到一定程度时会使人出现精神障碍、噩梦、失眠、头痛等慢性中毒症状；严重者还可有乏力、食欲不振、恶心、腹胀、腹痛或腹泻等。图简化图模型, 0) 0(, 0) 0(, 0) 0()()()()()()()()()()()()(321313131321202121231321213121011xxxtxatxadttdxt

49、xaatxadttdxbtxatxatxaaadttdx模型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.000035 = 0.000035from blood to bone and backa a0101 = 0.0

50、21 = 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 :=matrix(3,1,L, 0, 0);#定义方程

51、定义方程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 := e

52、valm(-(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);#得到解的表达式得到解的表

53、达式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:=sim

54、plify(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(lambd

55、a1*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:=simp

56、lify(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, xx

57、2, xx3, t=365.1460, color=red,blue,green, thickness=2); plot2:=%:#将两个图画在一个图上display(plot1,plot2);3.振动与摆的动画演示振动与摆的动画演示o 研究多种摆和质点振动问题研究多种摆和质点振动问题o 建立问题的微分方程模型建立问题的微分方程模型o 求出模型的解析解或数值解求出模型的解析解或数值解o 应用应用Maple 给出模型求解过程给出模型求解过程o 应用应用Maple 给出模型动画演示给出模型动画演示3.1 一个弹簧质量系统一个弹簧质量系统问题问题 一个弹性系数为一个弹性系数为k的弹簧左端固定，右端连

58、的弹簧左端固定，右端连接一个质量为接一个质量为m的滑块，放在一个光滑的水平的滑块，放在一个光滑的水平面上，给出一定的初始位移后让其运动，求该面上，给出一定的初始位移后让其运动，求该弹簧质量系统的运动规律。弹簧质量系统的运动规律。问题的数学模型为： +k6初始条件为：x(0)=8, x(0)=0m=k=1,求解此问题的Maple语句为：dsolve(diff(x(t),t,t)+x(t)=6, x(0)=(8),D(x)(0)=0,x(t);Maple给出的解为22d xmkxdt ( )x t62( )cos t 这里建立一个这样的模型，两个物体之间连这里建立一个这样的模型，两个物体之间连有弹

59、簧，两个物体两端又分别连接有弹簧，两个物体两端又分别连接1个弹簧，个弹簧，两个弹簧的另一端又分别固定在墙上。两个弹簧的另一端又分别固定在墙上。 将两个物体分别拉伸一段距离后放开物体，将两个物体分别拉伸一段距离后放开物体，弹簧弹簧 -质量系统就会做往复运动，观察它们的质量系统就会做往复运动，观察它们的运动情况，列出它们的运动方程。运动情况，列出它们的运动方程。2t2( )x1 t2( )x1 t( )x2 t2t2( )x2 t2( )x2 t( )x1 t( )x2 t52( )cos t12()cos3 t( )x1 t52( )cos t12()cos3 tMaple语句为k1:=1;k2:=1;k3:=1;m1:=1; m2:=1; m3:=1;x10:=2;x20:=3;s10:=0;s20:=0;deq1:=m1*diff(x1