數(shù)學(xué)實(shí)驗(yàn)第三講

數(shù)學(xué)實(shí)驗(yàn)之三微分方程解的形式①

解析解

y=f(x)②

數(shù)值解

(xi,yi)③圖形解xyo①簡單的微分方程。②復(fù)雜、大型的微分方程。一階微分方程:獲取解析解的方法歸類:①分離變量法；如dy/dx=x*y；②齊次方程的變換法；如dy/dx=f(y/x)③線性方程的常數(shù)變易法或公式法.

……解析解MATLAB軟件實(shí)現(xiàn)解析解dsolve('eqn1','eqn2',…,'c1',…,'var1',…)微分方程組初值條件變量組注意：①y'Dy，y''D2y②自變量名可以省略，默認(rèn)變量名‘t’。例①輸入：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)（特解，一條曲線）例②常系數(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')上述兩例的計算結(jié)果怎樣？由此得出什么結(jié)論？例③非常系數(shù)的二階微分方程例③無解析表達(dá)式！x=dsolve('(Dx)^2+x^2=1','x(0)=0')例④非線性微分方程x=[sin(t)][-sin(t)]若欲求解的某個數(shù)值解，如何求解？t=pi/2;eval(x)輸入：[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）數(shù)值解1、歐拉法2、龍格—庫塔法數(shù)值求解思想：(變量離散化)

引入自變量點(diǎn)列{xn}→{yn}，在x0

x1

x2

…

xn

…上求y(xn)的近似值yn.通常取等步長h，即xn=x0+n×h，或

xn

=xn-1+h，（n=1,2,…)。研究常微分方程的數(shù)值解法是十分必要的。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]上用差商代替微商(近似),

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例1

觀察向前歐拉、向后歐拉算法計算情況。與精確解進(jìn)行比較。誤差有多大？解：1）解析解：y=x+e-xy=dsolve('Dy=-y+x+1','y(0)=1','x')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;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--')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ù)值解比較表結(jié)果（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顯然迭代步長h的選取對精度有影響。圖形顯示

有什么方法可以使精度提高？向后歐拉法

對方程y’=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：使用數(shù)值積分即梯形法：梯形公式

改進(jìn)的歐拉公式y(tǒng)n+hf(xn,yn)

以例1為例，用改進(jìn)歐拉公式編程計算，再與精確解的比較。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.mx精確解向前歐拉向后歐拉改進(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ù)值解比較表結(jié)果使用泰勒公式

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

當(dāng)一個數(shù)值公式的截斷誤差可表示為O（hk+1）時（k為正整數(shù)，h為步長），稱它是一個k階公式。

k越大，則數(shù)值公式的精度越高。歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。龍格-庫塔法有二階公式和四階公式。線性多步法有四階阿達(dá)姆斯外插公式和內(nèi)插公式。[t，x]=solver（’f’,ts,x0,options）ode23

ode45ode113ode15sode23s由待解方程寫成的m-函數(shù)文件ts=[t0，tf]，t0、tf為自變量的初值和終值函數(shù)的初值ode23：組合的2/3階龍格-庫塔算法ode45：運(yùn)用組合的4/5階龍格-庫塔算法自變量值函數(shù)值用于設(shè)定誤差限(缺省時設(shè)定相對誤差10-3,絕對誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設(shè)定的相對誤差和絕對誤差.Matlab軟件計算數(shù)值解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)1、在解n個未知函數(shù)的方程組時，x0和x均為n維向量，m-函數(shù)文件中的待解方程組應(yīng)以x的分量形式寫成.2、使用Matlab軟件求數(shù)值解時，高階微分方程必須等價地變換成一階微分方程組.注意:注意1：1、建立M文件函數(shù)

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

f2(t,x(t),y(t))]；2、數(shù)值計算（執(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’;%列向量例如：令注意2：functionxdot=fun1(t,x,y)（fun1.m)xdot=[f(t,x(t),y(t))；x(t)]；[t,x,y]=ode23(‘fun1',[t0,tf],[x0,y0])M-文件函數(shù)如何寫呢？注意：y(t)是原方程的解。x(t)只是中間變量。如果方程形式是：z’’’=f(t,z,z’’)?例2Vanderpol方程：令y1=x(t),y2=x’(t);

該方程是否有解析解？（1）編寫M文件(文件名為vdpol.m):functionyp=vdpol(t,y);yp=[y(2);(1-y(1)^2)*y(2)-y(1)];（