數(shù)學(xué)實(shí)驗(yàn)之微分方程

數(shù)學(xué)實(shí)驗(yàn)之微分方程第一頁，共四十七頁，編輯于2023年，星期六

MATLAB軟件的實(shí)現(xiàn)引例：增長與傳播模型

實(shí)驗(yàn)內(nèi)容

微分方程模型與方法

解析解

數(shù)值解主要內(nèi)容第二頁，共四十七頁，編輯于2023年，星期六

1)掌握微分方程求解的三種解法：解析法、數(shù)值解法以及圖形表示解的方法;2)掌握求微分方程數(shù)值解的歐拉方法,了解龍格---庫塔方法的思想;3）學(xué)會(huì)使用MATLAB軟件求解析解、數(shù)值解和圖形解；

4）通過范例學(xué)習(xí)怎樣建立微分方程模型和分析問題的思想；實(shí)驗(yàn)?zāi)康姆祷氐谌?，共四十七頁，編輯?023年，星期六首先研究某些物種增長的一階微分方程.其中x(t)表示一種給定物種在時(shí)刻t的總數(shù);r(t,x)表示該物種的凈增長率;進(jìn)一步假定r(t,x)是常數(shù),于是得到簡(jiǎn)化模型：解：引例:群體增長模型第四頁，共四十七頁，編輯于2023年，星期六t=0:12;x=2*exp(0.4*t);plot(t,x,'r:o')思考:該模型是否符合物種的自然增長規(guī)律?第五頁，共四十七頁，編輯于2023年，星期六

該模型用于預(yù)測(cè)地球上的人口總數(shù),例如1961年地球上的人口總數(shù)為3,060,000,000，又假定人口總數(shù)以2%的速度增長(凈增長率)，故檢驗(yàn)：1700~1961年期間，實(shí)際數(shù)字與理論數(shù)值吻合！預(yù)測(cè)：p(2510)=2,000,000(億)；

p(2670)=36,000,000(億)；星球的總表面積約18,600,000億ft2，80%被水覆蓋，到2510年，每人占地面積只有9.3

ft2,到2670年，每人占地面積只有0.5167ft2。第六頁，共四十七頁，編輯于2023年，星期六分析以上模型,并修改模型.一般地，b<<r，若x不是很大，則可以忽略-bx2。相反,當(dāng)x很大時(shí),-bx2不能忽略,它將減慢群體迅速的增長率.解出結(jié)果：競(jìng)爭(zhēng)項(xiàng)第七頁，共四十七頁，編輯于2023年，星期六分析解

稱該曲線為邏輯曲線或S形曲線。符合很多實(shí)際情況。

物種發(fā)展到一定時(shí)期,會(huì)受到其它因素的制約,增長速度漸趨緩慢.第八頁，共四十七頁，編輯于2023年，星期六引例：新技術(shù)的傳播模型

一項(xiàng)新技術(shù)如何在同類企業(yè)中推廣?傳播速度怎樣？哪些因素影響傳播速度？第九頁，共四十七頁，編輯于2023年，星期六1）記p(t)表示t時(shí)刻采用了該技術(shù)的企業(yè)數(shù)，并假定它連續(xù)可微；2）該項(xiàng)技術(shù)具有成效,值得傳播；模型假設(shè)3）當(dāng)t=0時(shí)，有一個(gè)企業(yè)已經(jīng)采用該技術(shù)，共有N個(gè)同類企業(yè)數(shù)；4）在[t,t+△t]內(nèi)，△p(t)與p(t)(N-p(t))△t成正比；第十頁，共四十七頁，編輯于2023年，星期六由條件4），得到模型建立新技術(shù)傳播模型

某行業(yè)內(nèi)采納一項(xiàng)新技術(shù)的企業(yè)數(shù)必然依賴于該技術(shù)的效益L和所需要的投資資金s，一般情況下,一個(gè)企業(yè)的日傳播率r=f(L,s).但這里假定r為常數(shù).第十一頁，共四十七頁，編輯于2023年，星期六模型求解計(jì)算結(jié)果：邏輯曲線①P(t)>0，并且t→+∞時(shí)，p(t)→N;②對(duì)任意的t，p(t)<N;③當(dāng)p=N/2時(shí)，dp/dt達(dá)到最大；結(jié)果分析第十二頁，共四十七頁，編輯于2023年，星期六模型檢驗(yàn)

該模型需要進(jìn)行實(shí)證分析，檢驗(yàn)是否符合實(shí)際情況，檢驗(yàn)?zāi)Ｐ图僭O(shè)是否合理。

例如,考察過去的內(nèi)燃機(jī)火車頭、集中化的交通控制、車廂減速器等三種技術(shù)從一個(gè)企業(yè)推廣到另一些企業(yè)的速率問題，如下圖所示。第十三頁，共四十七頁，編輯于2023年，星期六理論曲線N=25t0=1925實(shí)測(cè)數(shù)據(jù)p0=1r1=4%r2=3.1%r3=1.5%第十四頁，共四十七頁，編輯于2023年，星期六模型修正

隨著通訊能力的提高和大眾媒介的普及，廣告將對(duì)這項(xiàng)新技術(shù)的傳播起到一定的作用。因此，在上述模型中增加一個(gè)因素——廣告效應(yīng)的作用。即修改模型如下：計(jì)算結(jié)果：仍然是邏輯曲線，與實(shí)際情況應(yīng)該是更加吻合。第十五頁，共四十七頁，編輯于2023年，星期六1、上述模型哪些假設(shè)需要修改？

2、流行性感冒的傳播模型；思考:第十六頁，共四十七頁，編輯于2023年，星期六定義：含有導(dǎo)數(shù)的方程稱為微分方程。如上例，

f(x,V(x),V’(x))=0微分方程模型1、微分方程的一般形式：

F(x,y,y’,…,y(n))=0

隱式或

y(n)=f(x,y,y’,…,y(n-1))

顯式n階第十七頁，共四十七頁，編輯于2023年，星期六特殊情形：2、一階微分方程組的一般形式：1階初始條件：y(x0)=y0微分方程模型第十八頁，共四十七頁，編輯于2023年，星期六滿足一定的初始條件，稱為偏微分方程。微分方程模型第十九頁，共四十七頁，編輯于2023年，星期六微分方程求解方法簡(jiǎn)介③圖形解xyo①簡(jiǎn)單的微分方程②復(fù)雜、大型的微分方程返回①

解析解

y=f(x)②

數(shù)值解

(xi,yi)歐拉方法改進(jìn)歐拉方法

梯形法龍格-庫塔法第二十頁，共四十七頁，編輯于2023年，星期六MATLAB軟件實(shí)現(xiàn)解析解dsolve('eqn1','eqn2',…,'c1',…,'var')微分方程組初值條件自變量名注意：①y'Dy，y''D2y②自變量名可以省略，默認(rèn)變量名‘t’。第二十一頁，共四十七頁，編輯于2023年，星期六例①輸入：y=dsolve('Dy=1+y^2')y1=dsolve('Dy=1+y^2','y(0)=1','x')輸出：y=tan(t-C1)（通解，一簇曲線）

y1=tan(x+1/4*pi)（特解，一條曲線）MATLAB軟件實(shí)現(xiàn)第二十二頁，共四十七頁，編輯于2023年，星期六例②常系數(shù)的二階微分方程y=dsolve('D2y-2*Dy-3*y=0','x')y=dsolve('D2y-2*Dy-3*y=0','y(0)=1,Dy(0)=0','x')輸入：x=dsolve('D2x-(1-x^2)*Dx+x=0','x(0)=3,Dx(0)=0')上述兩例的計(jì)算結(jié)果怎樣？由此得出什么結(jié)論？例③非常系數(shù)的二階微分方程例③無解析表達(dá)式！第二十三頁，共四十七頁，編輯于2023年，星期六x=dsolve('(Dx)^2+x^2=1','x(0)=0')例④非線性微分方程x=[sin(t)][-sin(t)]若欲求解的某個(gè)數(shù)值解，如何求解？t=pi/2;eval(x)MATLAB軟件實(shí)現(xiàn)第二十四頁，共四十七頁，編輯于2023年，星期六輸入：[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=0,y(0)=1')例④輸出：（li3.m）MATLAB軟件實(shí)現(xiàn)返回第二十五頁，共四十七頁，編輯于2023年，星期六數(shù)值解1、歐拉法2、龍格—庫塔法數(shù)值求解思想：(變量離散化)

引入自變量點(diǎn)列{xn}→{yn}，在x0x1x2…xn…上求y(xn)的近似值yn.通常取等步長h，即xn=x0+n×h，或

xn

=xn-1+h，（n=1,2,…)。研究常微分方程的數(shù)值解法是十分必要的。解:y=y(x)第二十六頁，共四十七頁，編輯于2023年，星期六1)向前歐拉公式:（y’=f(x,y)）

y(xn+1)

y(xn)+hf(xn,y(xn))(迭代式)

yn+1

yn+hf(xn,yn)(近似式)

特點(diǎn)：f(x,y)取值于區(qū)間[xn,xn+1]的左端點(diǎn).1、歐拉方法在小區(qū)間[xn,

xn+1]上用差商代替微商(近似),第二十七頁，共四十七頁，編輯于2023年，星期六

yn+1

yn

+hf(xn+1,yn

+1)特點(diǎn)：①f(x,y)取值于區(qū)間[xn,xn+1]的右端點(diǎn).②非線性方程,稱‘隱式公式’。yn+1

=

yn+hf(xn,yn)2)向后歐拉公式方法：迭代（y’=f(x,y)）x=[];y=[];x(1)=x0;y(1)=y0;forn=1:kx(n+1)=x(n)+n*h;y(n+1)=

y(n)+h*f(x(n),y(n));(向前）end,x,y,1、歐拉方法第二十八頁，共四十七頁，編輯于2023年，星期六例1

觀察向前歐拉、向后歐拉算法計(jì)算情況。與精確解進(jìn)行比較。誤差有多大？解：1）解析解：y=x+e-xy=dsolve('Dy=-y+x+1','y(0)=1','x')1、歐拉方法第二十九頁，共四十七頁，編輯于2023年，星期六2）向前歐拉法：

yn+1=yn+h(-yn+xn+1)=(1-h)yn+hxn+h3）向后歐拉法：

yn+1=yn+h(-yn+1+xn+1+1)

轉(zhuǎn)化yn+1=(yn+hxn+1+h)/(1+h)y’=f(x,y)=-y+x+1;1、歐拉方法第三十頁，共四十七頁，編輯于2023年，星期六x1(1)=0;y1(1)=1;y2(1)=1;h=0.1;%(died.m)fork=1:10x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)+h*x1(k)+h;y2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);endx1,y1,y2,%（y1——向前歐拉解，y2——向后歐拉解）x=0:0.1:1;y=x+exp(-x)%（解析解）plot(x,y,x1,y1,'k:',x1,y2,'r--')1、歐拉方法第三十一頁，共四十七頁，編輯于2023年，星期六x精確解向前歐拉向后歐拉01110.11.004811.00910.21.01871.011.02640.31.04081.0291.05130.41.07031.05611.08300.51.10651.09051.12090.61.14881.13141.16450.71.19661.17831.21320.81.24931.23051.26650.91.30661.28741.324111.36791.34871.3855（1）步長h=0.1的數(shù)值解比較表計(jì)算結(jié)果第三十二頁，共四十七頁，編輯于2023年，星期六（2）步長h=0.01的數(shù)值解比較表x精確解向前歐拉向后歐拉01110.11.00481.00441.00530.21.01871.01791.01950.31.04081.03971.04190.41.07031.06901.07170.51.10651.10501.10800.61.14881.14721.15040.71.19661.19481.19830.81.24931.24751.25110.91.30661.30471.308411.36791.36601.3697結(jié)論：顯然迭代步長h的選取對(duì)精度有影響。第三十三頁，共四十七頁，編輯于2023年，星期六圖形顯示

有什么方法可以使精度提高？向后歐拉法返回第三十四頁，共四十七頁，編輯于2023年，星期六

對(duì)方程y’=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：2、使用數(shù)值積分即梯形法：第三十五頁，共四十七頁，編輯于2023年，星期六梯形公式

改進(jìn)歐拉公式y(tǒng)n+hf(xn,yn)返回第三十六頁，共四十七頁，編輯于2023年，星期六

以例1為例，用改進(jìn)歐拉公式編程計(jì)算，再與精確解的比較。yn+1=yn+(h/2)*[(-yn+xn+1)+(-yn+1+xn+1+1)]=yn+(h/2)*[(-yn+xn+1)-(yn+h*(-yn+xn+1))+xn+1+1]=yn+(h/2)*[(1-h)*xn+xn+1+2-h+(h-2)*yn]died1.m使用改進(jìn)歐拉公式的例第三十七頁，共四十七頁，編輯于2023年，星期六x精確解向前歐拉向后歐拉改進(jìn)歐拉011110.11.004811.00911.0050.21.01871.011.02641.0190.31.04081.0291.05131.04120.41.07031.05611.08301.07080.51.10651.09051.12091.10710.61.14881.13141.16451.14940.71.19661.17831.21321.19720.81.24931.23051.26651.25000.91.30661.28741.32411.307211.36791.34871.38551.3685步長

h=0.1的數(shù)值解比較表使用改進(jìn)歐拉公式的例第三十八頁，共四十七頁，編輯于2023年，星期六3、使用泰勒公式

以此方法為基礎(chǔ)，有龍格-庫塔法、線性多步法等方法。4、數(shù)值公式的精度

當(dāng)一個(gè)數(shù)值公式的截?cái)嗾`差可表示為O（hk+1）時(shí)（k為正整數(shù)，h為步長），稱它是一個(gè)k階公式。

k越大，則數(shù)值公式的精度越高。

歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。龍格-庫塔法有二階公式和四階公式。線性多步法有四階阿達(dá)姆斯外插公式和內(nèi)插公式。返回第三十九頁，共四十七頁，編輯于2023年，星期六[t，x]=solver（’f’,ts,x0）ode23

ode45由待解方程寫成的m-函數(shù)文件ts=[t0，tf]，t0、tf為自變量的初值和終值函數(shù)的初值ode23：組合的2/3階龍格-庫塔算法ode45：運(yùn)用組合的4/5階龍格-庫塔算法自變量值函數(shù)值Matlab軟件計(jì)算數(shù)值解第四十頁，共四十七頁，編輯于2023年，星期六1）首先建立M-文件（weif.m）

functionf=weif(x,y)f=-y+x+1;2）求解：[x，y]=ode23(‘weif’,[0,1],1)3)作圖形:plot(x,y,‘r’);4)與精確解進(jìn)行比較

holdonezplot(‘x+exp(-x)’,[0,1])例1

y’=-y+x+1，y(0)=1標(biāo)準(zhǔn)形式：y’=f(x,y)范例第四十一頁，共四十七頁，編輯于2023年，星期六1、在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-函數(shù)文件中的待解方程組應(yīng)以x的分量形式寫成.2、使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.注意:第四十二頁，共四十七頁，編輯于2023年，星期六注意1：1、建立M文件函數(shù)

functionxdot=fun(t,x,y)xdot=[f1(t,x(t),y(t))；

f2(t,x(t),y(t))]；2、數(shù)值計(jì)算（執(zhí)行以下命令）

[t,x,y]=ode23(‘fun',[t0,tf],[x0,y0])注意：執(zhí)行命令不能寫在M函數(shù)文件中。xd(1)=f1(t,x(t),y(t));xd(2)=f2(t,x(t),y(t));xdot=xd’;%