1、常微分方程及數(shù)學(xué)建模, 考慮一階常微分方程的初值問(wèn)題 /* Initial-Value Problem */:,只要 f (x, y) 在a, b R1 上連續(xù)，且關(guān)于 y 滿(mǎn)足 Lipschitz 條件，即存在與 x, y 無(wú)關(guān)的常數(shù) L 使 對(duì)任意定義在 a, b 上的 y1(x) 和 y2(x) 都成立，則上述IVP存在唯一解。,要計(jì)算出解函數(shù) y(x) 在一系列節(jié)點(diǎn) a = x0 x1 xn= b 處的近似值,節(jié)點(diǎn)間距 為步長(zhǎng)，通常采用等距節(jié)點(diǎn)，即取 hi = h (常數(shù))。,一、常微分方程數(shù)值解法,1 歐拉方法 /* Eulers Method */, 歐拉公式：, 歐拉法的局部截?cái)?/p>

2、誤差：,歐拉法具有 1 階精度。,Ri 的主項(xiàng) /* leading term */,亦稱(chēng)為歐拉折線(xiàn)法 /* Eulers polygonal arc method*/,1 Eulers Method, 歐拉公式的改進(jìn)：, 隱式歐拉法 /* implicit Euler method */,由于未知數(shù) yi+1 同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱(chēng)為隱式 /* implicit */ 歐拉公式，而前者稱(chēng)為顯式 /* explicit */ 歐拉公式。,一般先用顯式計(jì)算一個(gè)初值，再迭代求解。, 隱式歐拉法的局部截?cái)嗾`差：,即隱式歐拉公式具有 1 階精度。,Hey! Isnt the lead

3、ing term of the local truncation error of Eulers method ? Seems that we can make a good use of it ,1 Eulers Method, 梯形公式 /* trapezoid formula */, 顯、隱式兩種算法的平均,注：的確有局部截?cái)嗾`差 ， 即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計(jì)算時(shí)不得不用到迭代法，其迭代收斂性與歐拉公式相似。, 中點(diǎn)歐拉公式 /* midpoint formula */,假設(shè) ，則可以導(dǎo)出 即中點(diǎn)公式具有 2 階精度。,需要2個(gè)初值 y0

4、和 y1來(lái)啟動(dòng)遞推 過(guò)程，這樣的算法稱(chēng)為雙步法 /* double-step method */，而前面的三種算法都是單步法 /* single-step method */。,1 Eulers Method,簡(jiǎn)單,精度低,穩(wěn)定性最好,精度低, 計(jì)算量大,精度提高,計(jì)算量大,精度提高, 顯式,多一個(gè)初值, 可能影響精度,Cant you give me a formula with all the advantages yet without any of the disadvantages?,Do you think it possible?,Well, call me greedy,OK,

5、 lets make it possible., 改進(jìn)歐拉法 /* modified Eulers method */,注：此法亦稱(chēng)為預(yù)測(cè)-校正法 /* predictor-corrector method */。可以證明該算法具有 2 階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過(guò)程簡(jiǎn)單。后面將看到，它的穩(wěn)定性高于顯式歐拉法。,1 Eulers Method,2 龍格 - 庫(kù)塔法 /* Runge-Kutta Method */,建立高精度的單步遞推格式。,單步遞推法的基本思想是從 ( xi , yi ) 點(diǎn)出發(fā)，以某一斜率沿直線(xiàn)達(dá)到 ( xi+1 , yi+1 ) 點(diǎn)。歐拉

6、法及其各種變形所能達(dá)到的最高精度為2階。,斜率 一定取K1 K2 的平均值嗎？,步長(zhǎng)一定是一個(gè)h 嗎？,2 Runge-Kutta Method,首先希望能確定系數(shù) 1、2、p，使得到的算法格式有2階精度，即在 的前提假設(shè)下，使得,Step 1: 將 K2 在 ( xi , yi ) 點(diǎn)作 Taylor 展開(kāi),Step 2: 將 K2 代入第1式，得到,2 Runge-Kutta Method,Step 3: 將 yi+1 與 y( xi+1 ) 在 xi 點(diǎn)的泰勒展開(kāi)作比較,要求 ，則必須有：,這里有 個(gè)未知數(shù)， 個(gè)方程。,3,2,存在無(wú)窮多個(gè)解。所有滿(mǎn)足上式的格式統(tǒng)稱(chēng)為2階龍格 - 庫(kù)塔格

7、式。,注意到， 就是改進(jìn)的歐拉法。,Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？,其中i ( i = 1, , m )，i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 ) 均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。,2 Runge-Kutta Method, 最常用為四級(jí)4階經(jīng)典龍格-庫(kù)塔法 /* Classical Runge-Kutta Method */ ：,2 Runge-Kutta Method, 由于龍格-庫(kù)塔法的導(dǎo)出基于泰勒展開(kāi)，故精度主要受解函數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好采用低階算法而將步長(zhǎng)h 取小。,龍格庫(kù)塔法:

8、基于龍格庫(kù)塔法，MATLAB提供了求常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為： t,y=ode23(fname,ts,y0, options) t,y=ode45(fname,ts,y0,options) 前者是二階和三階龍格-庫(kù)塔公式,后者是四階和五階龍格-庫(kù)塔公式. 其中fname是待解的方程寫(xiě)成的M文件名，該函數(shù)文件必須返回一個(gè)列向量。ts=t0,t1, ,tf, t0和tf是自變量的初值和終值,并且t0,t1,tf指定在時(shí)刻t0,t1,tf輸出結(jié)果;t和y分別給出時(shí)間向量和相應(yīng)的狀態(tài)向量。Options用于設(shè)定誤差限,缺省時(shí)的相對(duì)誤差為10-3,絕對(duì)誤差為10-6,程序如下: optio

9、ns=odeset(reltol, rt, abstol, at) 這里,rt和at分別是設(shè)定的相對(duì)誤差和絕對(duì)誤差.,二、常微分方程數(shù)值解的Matlab實(shí)現(xiàn),大家知道導(dǎo)數(shù)的定義,它的含義通常是改變率、速度等，隨具體問(wèn)題而定。,（1）當(dāng)y=f(x)是一條曲線(xiàn)時(shí)，y(x)是x點(diǎn)處切線(xiàn)的斜率。,（2）當(dāng)x=x(t)是一個(gè)隨時(shí)間改變的量， 是單位時(shí)間相應(yīng)的量。,（3）當(dāng)某一個(gè)量為常數(shù)時(shí)，問(wèn)題 就簡(jiǎn)單了。 （4）導(dǎo)數(shù)也是函數(shù)。,三、微分方程建模,方法及案例,1.放水過(guò)程,有一個(gè)柱狀水箱，水平截面積為常數(shù)A，,原來(lái)水高H，t=0時(shí)刻下面一個(gè)面積為B的門(mén)打開(kāi)開(kāi)始放水，,求之后水位與t的關(guān)系，何時(shí)水放光？,解

10、：,開(kāi)門(mén)后，水自然向外流，開(kāi)始快，后來(lái)越來(lái)越慢，,記t時(shí)刻水位為h(t)小門(mén)處流速為v(t)，由能量守恒定律， 門(mén)口處水的勢(shì)能要轉(zhuǎn)化成動(dòng)能，兩者相等。,則,在t，t+t之間，-(h(t+t)-h(t)A=Bv(t)t， 即體積的變化等于流出的水量。,令t0則有,是變量分離的常微分方程。,初始情況是t=0時(shí)，h(0)=H,這就是水位與時(shí)間的關(guān)系。,在h=0，即水放光時(shí),2.耐用消費(fèi)品的銷(xiāo)售新產(chǎn)品的銷(xiāo)售量,一種耐用消費(fèi)品進(jìn)入市場(chǎng)后，一般是開(kāi)始銷(xiāo)得慢，逐漸加快，,當(dāng)普及了之后，速度又逐步減小，Product Life Cycle產(chǎn)品生 命周期。有人認(rèn)為應(yīng)該是鐘型曲線(xiàn)，請(qǐng)建模分析一下PLC曲線(xiàn)。,新產(chǎn)

11、品銷(xiāo)售的規(guī)律對(duì)于制定計(jì)劃指揮生產(chǎn)，促進(jìn)銷(xiāo)售有指導(dǎo)意義。,問(wèn)題分析：,未購(gòu)買(mǎi)者的購(gòu)買(mǎi)常通過(guò)兩種渠道的宣傳：,1.廠(chǎng)家或銷(xiāo)售商的廣告宣傳，另外在商店看到商品的表演，總之是消費(fèi)者之外的信息傳播；,2.已購(gòu)買(mǎi)者對(duì)未購(gòu)買(mǎi)者的宣傳，這是消費(fèi)者內(nèi)部的信息交流。,耐用消費(fèi)品通常一個(gè)家庭有一個(gè)即可，比如洗衣機(jī)、電冰箱、微 波爐、熱水器、DVD所以銷(xiāo)售數(shù)量即為銷(xiāo)售人數(shù)，不重復(fù)購(gòu)買(mǎi),建立模型：,記潛在市場(chǎng)人數(shù)為K，n(t)為t時(shí)刻已購(gòu)買(mǎi)者的人數(shù)，,t到t+t之間n1為完全由消費(fèi)者外信息交流造成的增加量, n2為由消費(fèi)者內(nèi)部造成的消費(fèi)增加量。,n1與未購(gòu)買(mǎi)者成正比，即n1=a(K-n(t)t,n2與未購(gòu)買(mǎi)者成正比，也

12、與已購(gòu)買(mǎi)者成正比，,所以 n2=bn(t)(K-n(t)t,其中a,b為比例系數(shù)（常數(shù)）,總增量 n(t)=n1+ n2=a(K-n(t) t+ bn(t)(K-n(t) t,設(shè)n(0)=0，微分方程的解為：,假設(shè) 假設(shè),3.二氧化碳的吸收,空氣通過(guò)盛有CO2的吸收劑的圓柱形器皿，已知它吸收CO2的量與 CO2的濃度及吸收層的厚度成正比，今有含CO28%的空氣通過(guò)厚度 為10cm的吸收層后濃度為2%，求： （1）若吸收層變?yōu)?0cm厚，出口濃度是多少？ （2）要使出口濃度為1%，應(yīng)該設(shè)多厚的吸收層？,解：,記吸收層厚度為d，等分n份，每小層d/n厘米。入口濃度為8%，在每小層看吸收量，第一層后

13、被吸收量為：k8%d/n，含量變?yōu)椋?第二層吸收了,第二層后濃度,=,依此類(lèi)推，最后第n層后的濃度為,從而n即無(wú)限細(xì)分通過(guò)d厘米后出口濃度,這就是我們要求的表達(dá)式。,在通過(guò)10cm厚吸收層后濃度為2%：8%e-k*10=2%得到,（1）在通過(guò)30cm厚吸收層后濃度為8%e-30*ln2/5=8%/26=0.125%,（2）要使出口濃度為1%，8%e-d*ln2/5=1%，則d=15cm,練習(xí)： 1.用處理放水問(wèn)題和耐用消費(fèi)品銷(xiāo)售量的方法推出出口 濃度與吸收層厚度的關(guān)系模型。 2.推出線(xiàn)長(zhǎng)度為L(zhǎng)的單擺的周期計(jì)算公式。,引入變量t：,表示厚度的變化，,引入函數(shù)f(t)表示通過(guò)厚度t后的濃度：,在(

14、t,t+t)中f(t)的增量的相反數(shù)即為濃度的吸收量，它與 濃度和厚度成正比：,微分方程解決問(wèn)題時(shí)常用4種方法建立模型： 1.直接分析y與x之間的關(guān)系，建立模型； 2.無(wú)限細(xì)分利用趨近于e的特殊極限建立模型； 3.把變量和它的各階導(dǎo)數(shù)都看作不同的變量，綜合建模。比如運(yùn)動(dòng)學(xué)中t表示時(shí)間，x表示走過(guò)的路程， 表示速度， 表示加速度，各看成各量建模。 4.追擊問(wèn)題有特殊的一套建模方法，下面我們介紹幾個(gè)這種方法的建模例題。,在追擊問(wèn)題里，追擊者常瞄著逃跑者而追去，二者之間的坐標(biāo) 形成一個(gè)差向量(x2-x1,y2-y1)，它就是追擊曲線(xiàn)的切向量方向。,令y=y2-y1， x=x2-x1， 注意到：y=

15、y/x 利用這一點(diǎn)，很快就 建立起微分方程來(lái)。,導(dǎo)彈打敵艦,一艘導(dǎo)彈驅(qū)逐艦在距敵艦a時(shí)發(fā)射一枚自動(dòng)跟蹤的導(dǎo)彈，與此同時(shí) 敵艦以v0a速度向與兩船聯(lián)線(xiàn)垂直的方向逃走，導(dǎo)彈速度是5v0a， 求導(dǎo)彈追擊軌跡與擊中時(shí)間。,解,以?xún)纱?lián)線(xiàn)為橫軸，驅(qū)逐艦為原點(diǎn)建坐標(biāo)架。,引入單位a，將所有的距離及速度都除以a，方便計(jì)算。,O,記t=0導(dǎo)彈發(fā)射時(shí)刻，t時(shí)刻敵船在Q，坐標(biāo)為：,(1,v0t)，導(dǎo)彈位置P(x,y)，差向量(1-x,v0t-y),就是曲線(xiàn)的切向量，,模型里y(t)，x(t)都是t的函數(shù)，但是三個(gè) 變量不好處理，注意我們要求的是y(x)。,再建立一個(gè)y(t)，x(t)，t的關(guān)系：t時(shí)間里導(dǎo)彈已飛行

16、的距離是可求的。,消去t得到,這個(gè)微分方程用我們所學(xué)過(guò)的知識(shí)解不了， 只能用MatLab求解。,語(yǔ)句格式是y=dsolve(方程1，方程n， 初始條件，自變量) 其中 Y用Dy表示，y用D2y表示，有關(guān)情況 可參考嗎MatLab的書(shū)籍。,求出結(jié)果,x=1時(shí)y=5/24即為擊中位置，擊中時(shí)間5/(24v0),A 題: 中國(guó)人口增長(zhǎng)預(yù)測(cè),2007 高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽題目,中國(guó)是一個(gè)人口大國(guó), 人口問(wèn)題始終是制約我國(guó)發(fā)展的關(guān)鍵因素之一。根據(jù)已有數(shù)據(jù)，運(yùn)用數(shù)學(xué)建模的方法，對(duì)中國(guó)人口做出分析和預(yù)測(cè)是一個(gè)重要問(wèn)題。,機(jī)動(dòng) 目錄 上頁(yè) 下頁(yè) 返回 結(jié)束,近年來(lái)中國(guó)的人口發(fā)展出現(xiàn)了一些新的特點(diǎn)，例

17、如，老齡化進(jìn)程加速、出生人口性別比持續(xù)升高，以及鄉(xiāng)村人口城鎮(zhèn)化等因素，這些都影響著中國(guó)人口的增長(zhǎng)。2007 年初發(fā)布的國(guó)家人口發(fā)展戰(zhàn)略研究報(bào)告(附錄1) 還做出了進(jìn)一步的分析。 關(guān)于中國(guó)人口問(wèn)題已有多方面的研究，并積累了大量數(shù)據(jù)資料。附錄2就是從中國(guó)人口統(tǒng)計(jì)年鑒上收集到的部分?jǐn)?shù)據(jù)。,費(fèi)馬 目錄 上頁(yè) 下頁(yè) 返回 結(jié)束,試從中國(guó)的實(shí)際情況和人口增長(zhǎng)的上述特點(diǎn)出發(fā)，參考附錄2中的相關(guān)數(shù)據(jù)（也可以搜索相關(guān)文獻(xiàn)和補(bǔ)充新的數(shù)據(jù)），建立中國(guó)人口增長(zhǎng)的數(shù)學(xué)模型，并由此對(duì)中國(guó)人口增長(zhǎng)的中短期和長(zhǎng)期趨勢(shì)做出預(yù)測(cè)；特別要指出你們模型中的優(yōu)點(diǎn)與不足之處。,費(fèi)馬 目錄 上頁(yè) 下頁(yè) 返回 結(jié)束,費(fèi)馬 目錄 上頁(yè) 下頁(yè) 返回 結(jié)束,費(fèi)馬 目錄 上頁(yè) 下頁(yè) 返回 結(jié)束,費(fèi)馬 目錄 上頁(yè) 下頁(yè) 返回 結(jié)束,費(fèi)馬 目錄 上頁(yè) 下頁(yè) 返回 結(jié)束,費(fèi)馬 目錄 上頁(yè) 下頁(yè) 返回 結(jié)束,費(fèi)馬 目錄 上頁(yè) 下頁(yè) 返回