數(shù)學(xué)建模

1、Maple 簡介,周義倉 西安交通大學(xué)理學(xué)院應(yīng)用數(shù)學(xué)系Tel. 82668741（O）中2-2210Email: zhouyc@mail.xjtu.edu.cn,主要內(nèi)容和參考文獻1.1 主要內(nèi)容參考資料Maple 軟件包介紹 引言，數(shù)值計算 ，代數(shù)運算 ， 圖形 ，微積分，線 性 代 數(shù) ， 微分方程， Maple 編程 Maple 應(yīng)用舉例 最速降線問題，人

2、體含鉛量,1.2 參考資料Maple與微分方程（復(fù)?。恢芰x倉、靳禎、秦軍林，常微分方程及其應(yīng)用，科學(xué)出版社，2003年6月。精英科技，孫非，Maple6實例教程，中國電力出版社，2001。李強，Maple8基礎(chǔ)教程，中國水利水電出版社，2004David Barrow等， Solving Differential Equations with Maple V, Brooks/Cole Publishing Company,

3、 1998.http://202.112.84.70/resource/maple.htmhttp://202.204.240.30/maple/tour/index.htm,Maple9教程“四部曲” Maple 9 Getting Started Guide , 834KB

4、 http://bolide.lamost.org/sheaweb/files/m9GettingStartedGuide.pdfMaple 9 Learning Guide, 4.49MB http://bolide.lamost.org/sheaweb/files/m9LearningGuide.pdf Maple 9 Introductory Programming Guide,

5、7.08MB http://bolide.lamost.org/sheaweb/files/m9IntroductoryProgrammingGuide.pdf 4. Maple 9 Advanced Programming Guide , 3.63MB http://bolide.lamost.

6、org/sheaweb/files/m9AdvancedProgrammingGuide.pdfMaple應(yīng)用http://www.maplesoft.com/applications/index.aspx,Maple 軟件包介紹2.1 引言 Maple是Waterloo大學(xué)開發(fā)的集數(shù)值計算、符號運算和圖形顯示于一體的軟件包。它可以把人們從煩瑣的代數(shù)運算中解脫出來，使得我們將主要精力集中與分析和解決問題的思路與

7、方法中去。數(shù)學(xué)運算分為數(shù)值運算、符號運算、邏輯運算等。 數(shù)值運算： 如函數(shù)求值、方程求根、矩陣的特征值、特征向量等，數(shù)值計算是一項簡單煩瑣的工作，計算機的出現(xiàn)解決了數(shù)值計算的困難，使得大規(guī)模、高速度的計算可以實現(xiàn)，使得數(shù)學(xué)在天氣預(yù)報、油藏模擬、航天等領(lǐng)域發(fā)揮了巨大的作用。,符號運算（代數(shù)運算）：這是一種更加智能化的計算，所處理的符號可以代表各種數(shù)、代數(shù)式、數(shù)學(xué)結(jié)構(gòu)（集合），代數(shù)運算是在一些代數(shù)規(guī)則下所進行

8、的數(shù)學(xué)處理，主要是尋求一個簡單完美的公式解答。 用計算機進行符號運算是數(shù)學(xué)與計算機領(lǐng)域中一個新的發(fā)展方向，是人們用計算機代替大腦進行數(shù)值運算后用計算機代替大腦進行代數(shù)運算的成功實踐，是數(shù)學(xué)走向機械化的重要步驟，這也使得計算機更加智能化。 目前數(shù)學(xué)軟件包很多，如 MathematicaMapleMatlabMathCADSASLindo……,Maple 的歷史1980年W

9、terloo大學(xué)符號計算研究小組成立；1985年：Maple3.3商業(yè)版發(fā)行；1992年：Windows操作系統(tǒng)下的Maple V Release2面世；1994年，Maple V Release3發(fā)行；1996年，Maple V Release4發(fā)行； 1997年，Maple V Release5發(fā)行； 1998年，Maple V Release5.2,5.2發(fā)行； 2000年，Maple 6發(fā)行； 2001年，Map

10、le 7發(fā)行； 200?年，Maple 8發(fā)行； 200?年，Maple 9發(fā)行； 200?年，Maple 9.5發(fā)行；,Maple提供了一整套問題解決的數(shù)學(xué)環(huán)境 它支持廣泛、大量的數(shù)學(xué)操作如：數(shù)值分析、符號代數(shù)、圖 形圖像。 Maple是一個廣泛應(yīng)用的符號計算軟件，它擁有強大的功能，主要表現(xiàn)在集符號運算、數(shù)值計算、可視化和程序設(shè)計于一體，這些功能是通過Maple提供的線性代數(shù)程序包、微積分方程程序包、統(tǒng)計程序包、偏微分方

11、程程序包以及畫圖程序包等實現(xiàn)的。 在和別的數(shù)學(xué)軟件相比，Maple在處理微分方程中問題的功能特別好。,在線幫助 Maple 允許通過如下方式獲得信息 ： 在由建菜單選擇 Help on XXXX,或 ctrl+F1得到鼠標(biāo)指針位置文 字的幫助 。 在 help菜單選擇 Introduction使用幫助導(dǎo)航器, 獲得層級結(jié)構(gòu)導(dǎo)航，單擊感興趣的標(biāo)題直到出現(xiàn)相應(yīng)幫助。 在help菜單選擇Topic Search 進行主題搜索。

12、在工作表中可使用?命令獲得明確主題的幫助頁。例如 ?plot 打開plot命令 的幫助頁。 在 help菜單選擇Full Text Search可以進行全文關(guān)鍵字檢索。 在 help菜單選擇History可以方便的回訪已訪問 的幫助頁。,2.2 數(shù)值計算 Maple可視為功能強大的計算器。請輸入下面的指令：32*12^13;200!;ifactor(%);expand(%); (2^30)/(3^20)*sqrt(3

13、);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 );,2.3 代數(shù)運算 展開

14、、分解、化簡表達式 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);

15、eval(expr2 , x=1 ); top := expr2; bottom := expand((3*x+5)*(2*x-1)); answer := normal( top/bottom );,表達式變形 my_expr := (a*x^2+b)/(x*(-3*x^2-x+4)); my_expr := (a*x^2+b)/(x*(-3*x^2-x+4)); 定義函數(shù) f := x -> x^2+1/2 ;

16、f(2); f(a+b); 使用 unapply命令將表達式轉(zhuǎn)化為函數(shù)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 );,定義如下方程組;eqn1 := a+2*b+3*c+4

17、*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來表示其他未知數(shù) a， b， c， d得到一組解 solve( {eqn1, eqn2, eqn3, eqn4}, {a, b, c, d} ); 解不等式ineq := x+y+4/(x+y) < 10: solve(

18、ineq, {x} );,2.4.圖 形 可以對顯式、隱式 、參數(shù)型函數(shù)及數(shù)據(jù)集作圖。 首先 用 with命令載入兩個作圖工具包, 列出可用的函數(shù)。with(plots); with(plottools); 2D作圖舉例 plot( tan(x), x=-2*Pi..2*Pi, y=-4..4); 隱函數(shù)作圖 implicitplot( { x^2+y^2=1, y=exp(x) }, x=-Pi..Pi

19、, y=-Pi..Pi,scaling=CONSTRAINED ); 線性不等數(shù)組的圖解 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),

20、optionsexcluded=(color=yellow) );,極坐標(biāo)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

21、^2-y^2), x=-2..2, y=-2..2, axes=BOXED,title=`A Surface Plot`); 球坐標(biāo)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,c

22、oords=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*P

23、i),t=0..1],n=1..360);,2.5 微積分f := x -> x*sin(a*x) + b*x^2; #定義函數(shù)diff( f(x), x ); #求導(dǎo) int(f(x), x); #積分int(f(x), x=1..2); #積分expr := (2*x+3)/(7*x+5); #定義表達式Limit( expr, x=infinity ); #顯示求表達式的極限limit( e

24、xpr, x=infinity ); #給出表達式的極限expr := sin(4*x)*cos(x): #定義表達式approx1 := series( expr, x=0 ); #展開為級數(shù)poly1 := convert( approx1, polynom ); #轉(zhuǎn)化為多項式 plot( {expr, poly1}, x=-1..1, y=-2..2, title = cat( convert(expr, s

25、tring)," vs. Series Approximation" ) ); #比較展開式與原函數(shù)的圖像,2.6 線 性 代 數(shù) 使用with(linalg)命令載入線 性代數(shù)工具包 restart; with(linalg); A:=matrix( 3, 3,[1/2, -1/3, 2, -5, 14/3, 9, 0, 11, -5/6] );#定義矩陣det(A); #計算行列式值 inver

26、se(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] ); #Van

27、dermonde矩陣 A := matrix( 3, 3 ); print(A);,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}

28、,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:=

29、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);,第21頁用逐次迭代法求下列初始值問題的解: 解. 該問題等價的積分方程為: 利用Maple去進行這些重復(fù)性的迭代:y0:=1;y1:=1+int(x^2+y0^2,x=0

30、..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);,第25頁DEtools[phaseportrait] #畫向量場及積分曲線；([diff(y(x),x)=-y(x)],y(x), #定義微分方程y’=-y;x=-2..2,

31、 #指出x的范圍[[y(-2)=2],[y(-2)=1],[y(-2)=-2]], #給出3個初始值dirgrid=[17,17], #定義網(wǎng)格點密度arrows=LINE, #定義線段類型axes=NORMAL);

32、 #定義坐標(biāo)系類型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);,第26頁DEtools[phaseportrait] #畫向量場及積分曲線；([diff(y(x),x)= x^2-y(x)],y(x),

33、 #定義微分方程x=-2..2, #指出x的范圍[[y(-2)=1.3],[y(-2)=-2]], #給出3個初始值dirgrid=[33,33], #定義網(wǎng)格點密度arrows=LINE,

34、 #定義線段類型axes=NORMAL); #定義坐標(biāo)系類型,第88-91: 常微分方程的數(shù)值解原理：函數(shù)的近似展開在此基礎(chǔ)上可以得到不同精度的算法,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)

35、->2*(x-y)+2*(y-x)*(1+(y-x)^2);for 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)+h^2*f2(x||n,y||n)/2;print (x||(n+1),y||(n+1),z||(n+1));od;,＃整理數(shù)據(jù)data1:=[[x0,y0],

36、[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]);plo

37、t({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*x^2)/(4-2*x)+1/(sqr

38、t(4-2*x));od;＃整理數(shù)據(jù)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]);,2.8 Maple 編程Maple提供了簡單

39、有效的編程語言.2.8.1 條件語句 if if 條件 then 語句組 fiif 條件 then 語句組 else 語句組 fiif 條件 then 語句組 elif 條件then 語句組 fiif 條件 then 語句組 elif 條件 then 語句組else 語句組 fix:=10:if x>2 then x^2 fi; if x0 then y:=1 else y:=0 fi;,2.8.2 循環(huán)語句fo

40、r,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 < 18 dos := s+ i;od;,2.8.3 控制流程語句 breakt:=1:while t>0 dot

41、:=t^2+1;if t>5 then break fi; od;2.8.4過程 Maple的絕大多數(shù)函數(shù)是以過程的形式存在，一個Maple語言的過程定義如下:lc := proc( s, u, t, v ) description "form a linear combination of the arguments"; s * u + t * vend proc;

42、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);,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:=

43、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("穩(wěn)定結(jié)點");elif((s1)and(s3)and(s8))the

44、n printf("不穩(wěn)定結(jié)點");elif((s1)and(s4)and(s6))then printf("穩(wěn)定的臨界或退化結(jié)點");elif((s1)and(s4)and(s8))then printf("不穩(wěn)定的臨界或退化結(jié)點");elif((s1)and(s5)and(s6))then printf("穩(wěn)定焦點"); elif((s1)and

45、(s5)and(s8))then printf("不穩(wěn)定焦點");elif((s1)and(s7))then printf("中心");elif(s2)then printf("鞍點");else printf("無法斷定");end if;end;,3.Maple 應(yīng)用舉例3.1最速降線問題 在數(shù)學(xué)史上最著名的問題之一就是最速降線問題

46、：在給定 點 P(高 )、Q(低 )之間，找一條曲線使一質(zhì)點沿該曲線無摩擦下滑 的時間最短 。建立坐標(biāo)系求曲線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+dy^2)/sqrt(2*g*y),x=0..a);Time=0.638

47、87restart;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,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

48、);with(student);middlesum(sqrt(1+dy^2)/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.經(jīng)過推導(dǎo)后得到方程利用 dsolve(y(x)*(diff(y(x),x)^2+1)=c,,y(x));求解后給出的結(jié)果太復(fù)雜,我們令,

49、3.2 人體含鉛量（第260頁 19題） 鉛是一種人體所必須的微量元素，但體內(nèi)鉛含量過多時就會引起鉛中毒.鉛是一種重金屬元素，通過食物、飲料、空氣進入人體，經(jīng)過呼吸和消化系統(tǒng)后進入血液，再經(jīng)過血液循環(huán)慢慢進入人體和骨頭中.鉛可以經(jīng)過人體的排泄系統(tǒng)、通出出汗、剪頭發(fā)、剪指甲排出體外. 我們根據(jù)鉛在人體內(nèi)的變化情況將人體分為血液、組織、骨頭3個倉室，鉛在這三個倉室中的轉(zhuǎn)化關(guān)系如圖所示.x1(t)表示t時刻血液中的含鉛量

50、，x2(t)表示t時刻組織中的含鉛量，x3(t)表示t時刻骨頭中的含鉛量.假設(shè)在單位時間內(nèi)從環(huán)境經(jīng)過消化、吸收系統(tǒng)進入血液的鉛為L，從血液進入組織、骨頭的鉛分別為a31*x1(t)和a21*x1(t)，從組織、骨頭再進入血液的鉛分別為a13*x2(t)、a12*x2(t)，血液和骨頭向外界排出的鉛分別為a01*x1(t)和a02*x3(t).. (1)根據(jù)前面的假設(shè)建立人體血液、組織和骨頭中含鉛量所滿足的微分方程組,并求其通解

51、.(2)取時間單位為天，對一個志愿者在一種環(huán)境中的生活情況進行測定得a21=0.011，a12=0.012, a31=0.0039，a13=0.000035,a01=0.021, a02=0.016,L=49.3. 假設(shè)該志愿者開始時體內(nèi)的含鉛量為0，求他體內(nèi)血液、組織和骨頭中含鉛量隨時間變化的關(guān)系，畫出這些解曲線的圖，并求當(dāng)x→+∞時這些函數(shù)的極限.(3)假設(shè)該志愿者在此環(huán)境中生活了365天后搬到了一個無鉛的環(huán)境中去(不再有外界

52、的鉛進入體內(nèi)，即L=0)，再討論他體內(nèi)血液、組織和骨頭中含鉛量隨時間變化的關(guān)系，并畫出0≤t≤1460時這些含鉛量曲線的圖形.,鉛是一種重金屬元素，在塵土、空中飄浮物、汽車尾氣、油漆、一些彩釉陶瓷制品以及塑料所含的增塑劑中含量較高，主要通過口腔---消化道進入生物體，對其神經(jīng)系統(tǒng)和血液系統(tǒng)有較強的毒性，會破壞其正常功能，其中以神經(jīng)毒性為主。 眾所周知，人的思想和行為受神經(jīng)系統(tǒng)控制，一旦神經(jīng)系統(tǒng)遭到破壞，其后果是可想而知的。兒童

53、正處于神經(jīng)系統(tǒng)發(fā)育完善期，是鉛污染的最大受害者。兒童好奇心強，喜歡用手去觸摸各種物品，而且一般有“手--口接觸（吮吸）”頻繁的行為特點，這就給鉛進入體內(nèi)創(chuàng)造了條件。由于兒童正處于生長發(fā)育階段，許多器官的功能都不完善，大腦的發(fā)育尚未完成，“血---腦屏障（一種保護性組織，有防止血液中物質(zhì)隨意進入腦細胞的功能）”不能阻止鉛離子（Ｐｂ2+）進入大腦。雖然，環(huán)境因素造成的鉛危害一般不會達到工業(yè)性中毒的程度，但低量、長期的鉛接觸會造成鉛在體內(nèi)積累

54、。雖然血鉛含量超標(biāo)并不等于鉛中毒，但研究表明，過量的鉛接觸會造成神經(jīng)元能量代謝降低，使兒童，尤其是學(xué)齡前兒童的神經(jīng)發(fā)育指數(shù)達不到應(yīng)有的標(biāo)準(zhǔn)，造成記憶力和思維能力下降，影響智力的發(fā)展——這種損害是不可逆的。,兒童體內(nèi)的鉛水平可分從血、骨、齒、尿、發(fā)檢測到，兒童血鉛水平分為5級。 Ⅰ級：血鉛值低于10微克／分升，身體處于相對安全狀態(tài)； Ⅱ級：血鉛值為10微克－19微克／分升，屬于輕度鉛中毒。影響造血、神經(jīng)傳導(dǎo)和認知能力，兒童易出現(xiàn)頭

55、昏、煩躁、注意力渙散、多動、厭食、腹脹、輕度貧血； Ⅲ級：血鉛值達到20微克－44微克／分升，為中度鉛中毒。引起缺鈣、缺鋅、缺鐵、免疫力低下，運動不協(xié)調(diào)，視力和聽力受損，學(xué)習(xí)困難、智力下降，生長發(fā)育遲緩，貧血、腹絞痛、食癖、反應(yīng)遲鈍等； Ⅳ級：血鉛值達到45微克－69微克／分升，就是重度鉛中毒?？沙霈F(xiàn)性格改變、易激怒、攻擊性行為、運動失調(diào)、貧血、腹絞痛、高血壓、心律失常和癡呆等； Ⅴ級：當(dāng)血鉛值大于70微克／分升時，為極重度

56、鉛中毒，可導(dǎo)致臟器損害、鉛性腦病、癱瘓、昏迷等。,鉛中毒會引起神經(jīng)、消化、循環(huán)系統(tǒng)紊亂，表現(xiàn)為貧血、腹痛、高血壓等。血鉛高組的兒童的總智商、操作智商、語言智商分別比低血鉛組落后１４分、１４分和１３分，而每升血液中的鉛濃度上升１００微克，兒童的身高將降低１．３厘米。 低劑量的鉛接觸可以對人體的紅細胞、腎臟、免疫系統(tǒng)、骨髓和中樞神經(jīng)的功能產(chǎn)生不良影響，而所有這些影響發(fā)生前都可能沒有明顯的臨床癥狀。鉛中毒會影響嬰幼兒最初站立、行走和說

57、話的年齡，也可能引起孩子注意力渙散、記憶力減退、理解力降低及學(xué)習(xí)困難等。 鉛通過各種方式進入人體后，可使人的智力下降，學(xué)習(xí)、工作成績低落；蓄積到一定程度時會使人出現(xiàn)精神障礙、噩夢、失眠、頭痛等慢性中毒癥狀；嚴重者還可有乏力、食欲不振、惡心、腹脹、腹痛或腹瀉等。,圖,簡化圖,模型,,模型x' = Ax + b,,數(shù)據(jù),時間分為兩段: [0,365]和[365,1460]#定義常數(shù)a01:= 0.021; a02:=

58、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+

59、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)});,結(jié)果太復(fù)雜,利用迭加原理、特

60、征值和特征向量法#非齊次方程的特解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.0

61、0746814];v2:=[-.5975248926, -.8018660787, .1178839056];v3:=[-.8256380196, .5640574390, .7307156328e-1];#特征向量組成矩陣 P:=augment(v1,v2,v3);,#得到解的表達式x1:=c[1]*P[1,1]*exp(lambda[1]*t)+c[2]*P[1,2]*exp(lambda[2]*t)+c[3]*P[1,3

62、]*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]

63、*t)+xe[3,1];#在解的表達式用t=0代入x10:=simplify(subs(t=0,x1));x20:=simplify(subs(t=0,x2));x30:=simplify(subs(t=0,x3));#求解常數(shù)solve({x10=0,x20=0,x30=0},{c[1],c[2],c[3]});assign(%);,#作圖plot([x1, x2, x3], t=0..365, color=[red,

64、blue,green], thickness=2); plot1:=%:,#重新顯示解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

65、(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);,# 計算解在t=

66、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));#求解常數(shù)so

67、lve({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);,#作圖plot([xx1, xx2, xx3], t=365..1460