第九講常微分方程模型2

1、6.1 人口增長(zhǎng)模型人口增長(zhǎng)模型英國(guó)人口學(xué)家malthus (1766-1834)模型假設(shè)模型假設(shè)人口自然增長(zhǎng)率人口自然增長(zhǎng)率 r 為常數(shù)為常數(shù)即單位時(shí)間內(nèi)人口的增長(zhǎng)量與當(dāng)時(shí)的人口呈正比即單位時(shí)間內(nèi)人口的增長(zhǎng)量與當(dāng)時(shí)的人口呈正比。模型建立模型建立1.指數(shù)增長(zhǎng)模型指數(shù)增長(zhǎng)模型人口以幾何級(jí)數(shù)增加！人口以幾何級(jí)數(shù)增加！模型分析模型分析0r( )x t 人口將人口將按指數(shù)規(guī)律無限增長(zhǎng)按指數(shù)規(guī)律無限增長(zhǎng)！ 0r 0( )x tx人口將人口將始終保持不變始終保持不變！ 0r ( )0 x t 人口將人口將按指數(shù)規(guī)律減少直至絕滅按指數(shù)規(guī)律減少直至絕滅！ 人口倍增時(shí)間人口倍增時(shí)間模型求解模型求解malthus

2、模型預(yù)測(cè)美國(guó)人口malthus模型預(yù)測(cè)美國(guó)人口malthus模型預(yù)測(cè)的優(yōu)缺點(diǎn)優(yōu)點(diǎn)優(yōu)點(diǎn)短期預(yù)報(bào)比較準(zhǔn)確缺點(diǎn)缺點(diǎn)不適合中長(zhǎng)期預(yù)報(bào)原因原因預(yù)報(bào)時(shí)假設(shè)人口增長(zhǎng)率 r 為常數(shù)。沒有考慮環(huán)境對(duì)人口增長(zhǎng)的制約作用。2.阻滯增長(zhǎng)模型阻滯增長(zhǎng)模型假設(shè)人口增長(zhǎng)率假設(shè)人口增長(zhǎng)率 r(t) 是是 t 時(shí)刻時(shí)刻人口人口 x(t) 的減函數(shù)的減函數(shù) ：其中，其中，xm 為考慮到受自然資源和環(huán)境條為考慮到受自然資源和環(huán)境條件限制所能容納的最大人口數(shù)量件限制所能容納的最大人口數(shù)量（稱（稱最大人口容量最大人口容量） 模型假設(shè)模型假設(shè)模型建立模型建立模型分析（定性分析）模型分析（定性分析）0mxx( )mx tx人口將人口將遞減

3、并趨向于遞減并趨向于xm！ 0mxx( )mx tx人口將人口將始終保持始終保持xm不變不變！ 00mxx( )mx tx人口將人口將遞增并趨向于遞增并趨向于xm！ 無論在哪種情況下，人口最終將趨向于最大人口容量！無論在哪種情況下，人口最終將趨向于最大人口容量！模型求解模型求解 xm/2 xm x dtdx t xm /2 x m x 2mxx人口增長(zhǎng)率達(dá)到最大值人口增長(zhǎng)率達(dá)到最大值阻滯增長(zhǎng)模型預(yù)測(cè)美國(guó)人口阻滯增長(zhǎng)模型預(yù)測(cè)美國(guó)人口阻滯增長(zhǎng)模型預(yù)測(cè)的優(yōu)缺點(diǎn)優(yōu)點(diǎn)優(yōu)點(diǎn) 中期預(yù)報(bào)比較準(zhǔn)確缺點(diǎn)缺點(diǎn) 理論上很好，實(shí)用性不強(qiáng)原因原因預(yù)報(bào)時(shí)假設(shè)固有人口增長(zhǎng)率 r 以及最大人口容量 xm 為定值。實(shí)際上這兩個(gè)參

4、數(shù)（特別是 xm ）很難確定，而且會(huì)隨著社會(huì)發(fā)展情況變化而變化。前面圖中曲線末端分叉就是由于這個(gè)原因。6.2 藥物在體內(nèi)的分布與排除藥物在體內(nèi)的分布與排除 藥物進(jìn)入機(jī)體形成藥物進(jìn)入機(jī)體形成血藥濃度血藥濃度( (單位體積血液的藥物量單位體積血液的藥物量) ) 血藥濃度需保持在一定范圍內(nèi)血藥濃度需保持在一定范圍內(nèi)給藥方案設(shè)計(jì)給藥方案設(shè)計(jì) 藥物在體內(nèi)吸收、分布和排除過程藥物在體內(nèi)吸收、分布和排除過程 藥物動(dòng)力學(xué)藥物動(dòng)力學(xué) 建立建立房室模型房室模型藥物動(dòng)力學(xué)的基本步驟藥物動(dòng)力學(xué)的基本步驟 房室房室機(jī)體的一部分，藥物在一個(gè)房室內(nèi)均勻機(jī)體的一部分，藥物在一個(gè)房室內(nèi)均勻分布分布( (血藥濃度為常數(shù)血藥濃度為

5、常數(shù)) )，在房室間按一定規(guī)律轉(zhuǎn)移，在房室間按一定規(guī)律轉(zhuǎn)移 本節(jié)討論本節(jié)討論二室模型二室模型中心室中心室( (心、肺、腎等心、肺、腎等) )和和周邊室周邊室( (四肢、肌肉等四肢、肌肉等) ) 中心室中心室周邊室周邊室給藥給藥排出排出( )f t111)(),(vtxtc222)(),(vtxtc12k21k13k1121131212( )( )x tk xk xkxf t 模型假設(shè)模型假設(shè) 中心室中心室(1)和周邊室和周邊室(2), ,容積不變?nèi)莘e不變 藥物在房室間轉(zhuǎn)移速率及向體外排藥物在房室間轉(zhuǎn)移速率及向體外排除除 速率，與該室血藥濃度成正比速率，與該室血藥濃度成正比 藥物從體外進(jìn)入中心室

6、，在二室間藥物從體外進(jìn)入中心室，在二室間 相互轉(zhuǎn)移相互轉(zhuǎn)移, ,從中心室排出體外從中心室排出體外模型建立模型建立2 , 1)()(ivtctxiii容積濃度藥量f給藥速率2211122)(xkxktxttttebeatcebeatc222111)()(211213121211121212122( )( )()( )vf tctkkckcvvvctk ckcv 2 , 1),()(itcvtxiii線性常系數(shù)線性常系數(shù)非齊次方程非齊次方程對(duì)應(yīng)齊次對(duì)應(yīng)齊次方程通解方程通解模型建立模型建立1210000( ),( ),( )dftccv幾種常見的給藥方式幾種常見的給藥方式1. .快速靜脈注射快速靜脈

7、注射t=0 瞬時(shí)瞬時(shí)注射劑量注射劑量d 的的藥物進(jìn)入中心室藥物進(jìn)入中心室, ,血藥血藥濃度立即為濃度立即為 d/v12211122121022112113121)()()()(ckckvvtcvtfckvvckktc1321132112kkkkk給藥速率給藥速率 f (t) 和初始條件和初始條件0lim ( )tc t11113112222211321121311213212121221200( ),( ),()(),ttttkcta eb ettk vkkcta eb ettkk vvkkvkkaabbk vk v 12000( ),( )( )f tk cc2. .恒速靜脈滴注恒速靜脈滴注

8、211213121211121212122( )( )()( )vf tc tkkck cvvvctk ck cv t t時(shí)時(shí), c1(t)和和 c2(t)按指數(shù)規(guī)律衰減趨于零按指數(shù)規(guī)律衰減趨于零藥物以恒定速率藥物以恒定速率k 進(jìn)入中心室進(jìn)入中心室0tt 010fk x)(0tx吸收室中心室00 1000( )()xtkxxd tktteebeaetc01)(1010( )k tx tde01001001( )( )k tftk x tdk e3. .口服或肌肉注射口服或肌肉注射相當(dāng)于藥物相當(dāng)于藥物( 劑量劑量d)先進(jìn)入吸收室，吸收后再進(jìn)入中心室先進(jìn)入吸收室，吸收后再進(jìn)入中心室吸收室藥量吸收室

9、藥量x0(t)211213121211121212122( )( )()( )vf tc tkkck cvvvctk ck cv ebacc,0)0(, 0)0(2110lim( )tc tttbeaetctc)()(11參數(shù)估計(jì)參數(shù)估計(jì)各種給藥方式下的各種給藥方式下的 c1(t), c2(t) 取決于參數(shù)取決于參數(shù)k12, k21, k13, v1,v2以以快速靜脈注射快速靜脈注射為例為例 , ,在在ti(i=1,2,n)測(cè)得測(cè)得c1(ti)121211( )()()()ttdc tkekev 充分大設(shè)t ,由由較大的較大的 用用最小二乘法最小二乘法定定a a, , 1,()iitct由由較

10、小的較小的 用用最小二乘法最小二乘法定定b, , 1,( )iit c t2111()( )()ttd kcteaev 211312kkk11110( )dcabv13110( )dk vct dt10,tc 1321132112kkkkk11131abdk v 111 311()abkba 1321kk參數(shù)估計(jì)法一參數(shù)估計(jì)法一進(jìn)入中心室的藥物全部排除進(jìn)入中心室的藥物全部排除參數(shù) 估計(jì)法二n% 構(gòu)造非線性擬合函數(shù)構(gòu)造非線性擬合函數(shù)ntwoexps.mnfunction e=twoexps(a,x,y)nx=x(:);y=y(:);ny=a(1)*exp(-a(3)*x)+a(2)*exp(-a

11、(4)*x);ne=sum(y-y).2) na0=100 1 1 ;noptions=optimset(fminsearch);noptions.tolx=0.01;noptions.display=off;na=fminsearch(ps,a0,options,x,y) a= 112.2378 0.1823 2.1773 6.3 傳染病模型傳染病模型問題問題 描述傳染病的傳播過程描述傳染病的傳播過程 分析受感染人數(shù)的變化規(guī)律分析受感染人數(shù)的變化規(guī)律 預(yù)報(bào)傳染病高潮到來的時(shí)刻預(yù)報(bào)傳染病高潮到來的時(shí)刻 預(yù)防傳染病蔓延的手段預(yù)防傳染病蔓延的手段 按照傳播過程的一般規(guī)律，按照傳播過程的一般規(guī)律，用

12、機(jī)理分析方法建立模型用機(jī)理分析方法建立模型 已感染人數(shù)已感染人數(shù) (病人病人) i(t) 每個(gè)病人每天有效接觸每個(gè)病人每天有效接觸(足以使人致病足以使人致病)人數(shù)為人數(shù)為 模型模型1 1假設(shè)假設(shè)ttititti)()()(若有效接觸的是病人，若有效接觸的是病人，則不能使病人數(shù)增加則不能使病人數(shù)增加必須區(qū)分已感染者必須區(qū)分已感染者(病人病人)和未感染者和未感染者(健康人健康人)建模建模0)0(iiidtdiitteiti0)(?sidtdi1)()(tits模型模型2 2區(qū)分已感染者區(qū)分已感染者(病人病人)和未感染者和未感染者(健康人健康人)假設(shè)假設(shè)1）總?cè)藬?shù)）總?cè)藬?shù)n不變，病人和健康不變，病人

13、和健康 人的人的 比例分別為比例分別為)(),(tsti 2）每個(gè)病人每天有效接觸人數(shù)）每個(gè)病人每天有效接觸人數(shù)為為 , 且且使接觸的健康人致病使接觸的健康人致病建模建模ttnitstittin)()()()(0)0()1(iiiidtdi 日日接觸率接觸率si 模型模型teiti1111)(00)0()1(iiiidtdi模型模型21/2tmii010t11ln01itmtm傳染病高潮到來時(shí)刻傳染病高潮到來時(shí)刻 (日接觸率日接觸率) tm 1itlogistic 模型病人可以治愈！病人可以治愈！?t=tm, di/dt 最大最大模型模型3傳染病無免疫性傳染病無免疫性病人治愈成病人治愈成為健康

14、人，健康人可再次被感染為健康人，健康人可再次被感染增加假設(shè)增加假設(shè)sis 模型模型3）病人每天治愈的比例為）病人每天治愈的比例為 日日治愈率治愈率ttnittitnstittin)()()()()(建模建模/ 日接觸率日接觸率1/ 感染期感染期 一個(gè)感染期內(nèi)一個(gè)感染期內(nèi)每個(gè)病人的每個(gè)病人的有效接觸人數(shù)，稱為有效接觸人數(shù)，稱為接觸數(shù)接觸數(shù)。0)0()1(iiiiidtdi1,01,11)(i)11 (iidtdi模型模型3i0i0接觸數(shù)接觸數(shù) =1 閾值閾值/1)(ti形曲線增長(zhǎng)按sti )(感染期內(nèi)感染期內(nèi)有效接觸感染的有效接觸感染的健康者人數(shù)不超過病人數(shù)健康者人數(shù)不超過病人數(shù)小01i1-1/

15、 i0iiidtdi)1 (模型模型2(si模型模型)如何看作模型如何看作模型3(sis模型模型)的特例的特例idi/dt01 10ti 11-1/ i0t 1di/dt 1/ i(t)先升后降至先升后降至0p2: s01/ i(t)單調(diào)降至單調(diào)降至01/ 閾值閾值p3p4p2s0ssss00lnln模型模型4sir模型模型預(yù)防傳染病蔓延的手段預(yù)防傳染病蔓延的手段 (日接觸率日接觸率) 衛(wèi)生水平衛(wèi)生水平 (日日治愈率治愈率) 醫(yī)療水平醫(yī)療水平 傳染病不蔓延的條件傳染病不蔓延的條件s01/ 的估計(jì)的估計(jì)0ln1000sssis0i忽略 降低降低 s0提高提高 r0 1000ris 提高閾值提高閾

16、值 1/ 降低降低 (= / ) , 群體免疫群體免疫sir模型模型被傳染人數(shù)的估計(jì)法一被傳染人數(shù)的估計(jì)法一 0ln1000sssis記被傳染人數(shù)比例記被傳染人數(shù)比例ssx00)211 (200sxsx0)1ln(10sxx)1(200ssx2xxs0i0s/1p10ssi0 0, s0 1 小小, s0 1提高閾值提高閾值 1/ 降低降低被傳染人數(shù)比例被傳染人數(shù)比例 xs0 - 1/ = 被傳染人數(shù)的估計(jì)法二被傳染人數(shù)的估計(jì)法二nx=fzero(x-1.2*log(x/0.96)-0.99,0.5)n x=0.86510ln1000sssis000.96,0.03,1.2si6.4多種群生態(tài)

17、數(shù)學(xué)模型多種群生態(tài)數(shù)學(xué)模型 意大利生物學(xué)家dancona曾致力于魚類種群相互制約關(guān)系的研究，他從第一次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚類捕獲量百分比的資料中，發(fā)現(xiàn)鯊魚等的比例有明顯增加（見下表），而供其捕食的食用魚的百分比卻明顯下降.顯然戰(zhàn)爭(zhēng)使捕魚量下降，食用魚增加，鯊魚等也隨之增加，但為何鯊魚的比例大幅增加呢？ 他無法解釋這個(gè)現(xiàn)象，于是求助于其岳父，著名的意大利數(shù)學(xué)家v.volterra，希望建立一個(gè)食餌捕食系統(tǒng)的數(shù)學(xué)模型，定性或定量地回答這個(gè)問題.年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分

18、比27.316.015.914.819.72基基本本假假設(shè)設(shè)：（1） 食餌由于捕食者的存在使增長(zhǎng)率降低，假設(shè)降低的程度與捕食者數(shù)量成正比； （2）捕食者由于食餌為它提供食物的作用使其死亡率降低或使之增長(zhǎng)，假定增長(zhǎng) 的程度與食餌數(shù)量成正比。 該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系，沒有考慮食餌和捕食者自身的阻滯作用，是volterra提出的最簡(jiǎn)單的模型.0, , ,dxaxbxydta b c ddycydxydt 模型（一）模型（一）不考慮捕獲定理定理 volterra微分方程組對(duì)應(yīng)初值問題0)0(, 0)0(yx的解)(),(tytx是周期函數(shù)，且解的周期平均值為

19、001d,1dttcxx tttdayy tttb首先，建立m-文件shier.m如下： function dx=shier(t,x) dx=zeros(2,1); dx(1)=x(1)*(1-0.1*x(2); dx(2)=x(2)*(-0.5+0.02*x(1);其次，建立主程序shark.m如下： t,x=ode45(shier,0 15,25 2); plot(t,x(:,1),-,t,x(:,2),*) plot(x(:,1),x(:,2)to matlab(shark)相圖),(21xx為：0510150102030405060708090100020406080100051015

20、202530數(shù)值解如下圖：)(1tx為實(shí)線，)(2tx為“*”線.求解結(jié)果： 左圖反映了x1（t）與x2（t）的關(guān)系。 可以猜測(cè)： x1（t）與x2（t）都是周期函數(shù)。模型（二）模型（二） 考慮人工捕獲 設(shè)表示捕獲能力的系數(shù)為e，相當(dāng)于食餌的自然增長(zhǎng)率由a降為a-e，捕食者的死亡率由c增為 c+e() ()dxx aebydtdyycedxdt10201,0.1,0.5,0.02,25,2abcdxx仍取設(shè)戰(zhàn)前捕獲能力系數(shù)e=0.3, 戰(zhàn)爭(zhēng)中降為e=0.1, 則戰(zhàn)前與戰(zhàn)爭(zhēng)中的模型分別為:2)0(,25)0()02. 08 . 0()1 . 07 . 0(21122211xxxxdtdxxxdt

21、dx2)0(,25)0()02. 06 . 0()1 . 09 . 0(21122211xxxxdtdxxxdtdxvolterra原理原理模型求解:1、分別用m-文件shier1.m和shier2.m定義上述兩個(gè)方程2、建立主程序shark1.m, 求解兩個(gè)方程，并畫出兩種情況下鯊魚數(shù)在魚類總數(shù)中所占比例 x2(t)/x1(t)+x2(t)to matlab(shark1)05101500.10.20.30.40.50.60.70.8 實(shí)線為戰(zhàn)前的鯊魚比例，“*”線為戰(zhàn)爭(zhēng)中的鯊魚比例結(jié)論：戰(zhàn)爭(zhēng)中鯊魚的比例比戰(zhàn)前高！結(jié)論：戰(zhàn)爭(zhēng)中鯊魚的比例比戰(zhàn)前高！function y=shier(t,x)r=

22、1;d=0.5;a=0.1;b=0.02;y=diag(r-a*x(2),-d+b*x(1)*x;shier.mts=0:0.1:35;x0=25,2;t,x=ode45(shier,ts,x0);t,x,plot(t,x),grid,gtext(x1(t),gtext(x2(t),pause,plot(x(:,1),x(:,2),grid,xlabel(x1),ylabel(x2)shiyan42注：ts 中終值(=15)和步長(zhǎng)=(0.1)的確定 t x1 x2 0 25.0000 2.0000 0.1000 27.0818 2.0041 0.2000 29.3344 2.0170 0.30

23、00 31.7689 2.0394 0.8000 46.9360 2.3503 0.9000 50.6072 2.4683 1.0000 54.5301 2.6106 10.3000 17.8063 2.0898 10.4000 19.2755 2.0627 10.5000 20.8708 2.0422 10.6000 22.6016 2.0287 10.7000 24.4779 2.0228 10.8000 26.5102 2.0249 10.9000 28.7093 2.0355 11.0000 31.0844 2.0558 14.7000 9.9982 25.9788 14.8000 8

24、.5540 25.1736 14.9000 7.3762 24.3300 15.0000 6.4158 23.46450510150102030405060708090100 x1(t)x2(t)0102030405060708090100051015202530 x1x2計(jì)算結(jié)果（數(shù)值，圖形）x(t),y(t)是周期函數(shù)，相圖是周期函數(shù)，相圖(x,y)是封閉曲線；是封閉曲線；axyrxxayrtx)()(bxydyybxdty)()(00) 0(,) 0(yyxx2,25,02. 0, 1 . 0, 5 . 0, 100yxbadr觀察，猜測(cè)觀察，猜測(cè)x(t),y(t)的周期約為的周期約為1

25、0.7;xmax=99.3, xmin=2.0, ymax=28.4, ymin=2.0. 用數(shù)值積分可算出用數(shù)值積分可算出x(t)一周期的平均值為一周期的平均值為25, y(t)一周期的平均值為一周期的平均值為10.6.5其它生態(tài)數(shù)學(xué)模型其它生態(tài)數(shù)學(xué)模型存在一大類生態(tài)模型源于對(duì)volterra模型的改造0 , , ,dxaxbxydta b c ddycydxydt 模型模型1考慮食餌種群與外界有遷入或遷出 g.r.walsh（1978） 0 外界有食餌遷入0 外界有食餌遷出 也可以表示人工干預(yù)，如投放或捕獲模型討論模型討論20-, , , ,-dxaax bxyxdtkk a b c dd

26、ycydxydt模型模型2考慮食餌種群內(nèi)部存在生存競(jìng)爭(zhēng) g.bojadziev0k表示當(dāng)沒有捕食者存在時(shí)食餌種群的環(huán)境容納量220 , , , , ,dxaxbxykxdtk l a b c ddycydxylydt模型模型3考慮食餌和捕食者種群內(nèi)部都存在生存競(jìng)爭(zhēng) 張錦炎張錦炎（1979） 0akk表示當(dāng)沒有捕食者存在時(shí)食餌種群的環(huán)境容納量20, , ,/dxaxbxydta b c ddycydyxdt模型模型4考慮雙方內(nèi)部都存在生存競(jìng)爭(zhēng)，且捕食者另有食物來源 e.c.pielou 平衡點(diǎn)平衡點(diǎn)：,adaebcb1112221211( )()( )()xyx tr xnnxyy tr ynn

27、 模型模型兩種群競(jìng)爭(zhēng)模型兩種群競(jìng)爭(zhēng)模型競(jìng)爭(zhēng)排斥原理競(jìng)爭(zhēng)排斥原理（competition exclution law） 多個(gè)種群依靠同一個(gè)生存資源而生活，如果生活在同一個(gè)地理空間，獵取相同食物或營(yíng)養(yǎng)物。在有限的相同生存資源條件下，如果存在競(jìng)爭(zhēng)關(guān)系，它們必然相互排斥，展開激烈的生存競(jìng)爭(zhēng)。結(jié)局是競(jìng)爭(zhēng)力較弱的種群滅絕，競(jìng)爭(zhēng)力最強(qiáng)的種群達(dá)到其環(huán)境容納量。 模型討論模型討論121 種群 y 最終將被滅絕,種群 x 最終趨于最大容量 211 種群 x 最終將被滅絕,種群 y 最終趨于最大容量 121,1存在過正平衡點(diǎn)的一條分界線，將第一象限分成種群 x 和種群 y 的兩個(gè)吸引域。121,1種群 x 和種群

28、y 最終達(dá)到穩(wěn)定的正平衡態(tài)競(jìng)爭(zhēng)排斥原理競(jìng)爭(zhēng)排斥原理 是針對(duì)前三種情形得出的結(jié)論， 第四種情況極為罕見。n% m 函數(shù)nfunction dy=cwf1 (t,y)n dy=zeros(2,1);n dy(1)=0.1*y(1)*(1-0.001*y(1)-0.0008*y(2);n dy(2)=0.2*y(2)*(1-0.0012*y(1)-0.001*y(2);n 主程序t,x=ode45(cwf1,0 200,200 200);t,y=ode45(cwf1,0 200,500 200); t,z=ode45(cwf1,0 200,1200 500);plot(x(:,1),x(:,2),

29、(y(:,1),y(:,2),(z(:,1),z(:,2)200300400500600700800900100011001200050100150200250300350400450500 11n 1112221211 ( )()( )()xyx tr xnnxyy tr ynn 模型模型兩種群互惠模型兩種群互惠模型研究多個(gè)種群之間相互依賴、共生現(xiàn)象。模型討論模型討論模型有三個(gè)平衡點(diǎn)，分別為模型有三個(gè)平衡點(diǎn)，分別為)1) 1(,1)1 (),0 ,(),0 , 0(211221113121nnpnpp1212:1,1,1a 兩種群最終達(dá)到穩(wěn)定平衡態(tài)兩種群共生p3 為正平衡點(diǎn)p3 穩(wěn)定1212

30、:1,1,1b 兩種群最終達(dá)不到穩(wěn)定平衡態(tài)p3 不穩(wěn)定不穩(wěn)定%m程序nfunction dy=cwf2 (t,y)n dy=zeros(2,1);n dy(1)=0.1*y(1)*(1-0.001*y(1)+0.0005*y(2);n dy(2)=0.2*y(2)*(-1+0.0015*y(1)-0.001*y(2);nt,x=ode45(cwf2,0,100,1600,2800);nt,y=ode45(cwf2,0,100,1200,2500);nt,z=ode45(cwf2,0,100,2500,2200);nplot(x(:,1),x(:,2),y(:,1),y(:,2),z(:,1),

31、z(:,2)ntext(2000,2600, 2000,2600,sigma_11,sigma_1sigma_21)1200140016001800200022002400260012001400160018002000220024002600280011,1216.6常微分方程的數(shù)值解及實(shí)驗(yàn)常微分方程的數(shù)值解及實(shí)驗(yàn)（一）常微分方程數(shù)值解的定義（一）常微分方程數(shù)值解的定義 在生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實(shí)際上對(duì)初值問題，一般是要求得到解在若干個(gè)點(diǎn)上滿足規(guī)定精確度的近似值，或者得到一個(gè)滿足精確度要求的便于計(jì)算的表達(dá)式。因此，研究常微分方程的數(shù)值解法是十分必要的因

32、此，研究常微分方程的數(shù)值解法是十分必要的。的相應(yīng)近似值求出準(zhǔn)確值，值處，即對(duì)的若干離散的開始其數(shù)值解是指由初始點(diǎn)，：對(duì)常微分方程nnnyyxyxyxxxxxy,y )(,),(),y(x x )y(xy)f(x,y 2121210000歐拉公式001i)y(xy)f(x,y , 1, 2 , 1 , 0 , xynihxi解微分方程：可用以下離散化方法求設(shè)1、用差商代替導(dǎo)數(shù)、用差商代替導(dǎo)數(shù) 若步長(zhǎng)h較小，則有hxyhxyxy)()()( 故有公式：1-n,0,1,2,i )(),(001xyyyxhfyyiiii此即歐拉法（歐拉法（向前歐拉法向前歐拉法）。2、使用數(shù)值積分、使用數(shù)值積分對(duì)方程y

33、=f(x,y), 兩邊由xi到xi+1積分，并利用梯形公式，有：)(,()(,(2)(,()()(11111iiiiiixxiixyxfxyxfxxdttytfxyxyii實(shí)際應(yīng)用時(shí)，與歐拉公式結(jié)合使用：, 2 , 1 , 0 ),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyykiiiiikiiiii的計(jì)算。然后繼續(xù)下一步，取時(shí)，當(dāng)滿足，對(duì)于已給的精確度）（ y y 2i111i)(1)1(1kikikiyyy此即改進(jìn)的歐拉法改進(jìn)的歐拉法。故有公式：)(),(),(200111xyyyxfyxfhyyiiiiii例1 求解初值問題0 . 1)0(8 . 102yxxy

34、y)(y21y0.1)y0.2(xyy)0.2x(1yycp1ipiiciip龍格龍格庫特方法庫特方法 n考慮微分中值定理 10),( )()(1hxyhxyxyiii)(,()()(1hxyhxhfxyxyiiii)(,(*hxyhxfkii2階龍格庫特公式階龍格庫特公式n 1,0),(),()(12122111hkyhxfkyxfkkkhyyiiiiiin利用泰勒展式得到 )()(,()(,()(,(212hoxyxfhkxyxhfxyxfkiiyiixii211223( ) ()( )( ( , ( )( , ( ) ( , ( )( )iiixiiiiyiiyy xhy xh f x y xf x y xf x y xoh1,21, 12214階龍格庫特公式階龍格庫特公式n ),(),(),(),()(3625143423122311121443322111hkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkkhyyiiiiiiiiii),()21,21()21,21(),()22(61342312