常微分方程數(shù)值解與matlab.ppt

實(shí)驗(yàn)4-常微分方程數(shù)值解,1.求解常微分方程數(shù)值方法介紹（1）一階微分方程求方程(1)的數(shù)值解,就是計(jì)算(精確)解在一系列離散點(diǎn)的近似值.通常取相等的步長(zhǎng)h,于是xn=x0+nh(n=1,2,).(a)歐拉方法基本思想是在小區(qū)間xn,xn+1上用差商代替方程(1)左端的導(dǎo)數(shù)而方程右端函數(shù)f(x,y(x)中的x取xn,xn+1上得某一點(diǎn),公式為（2）,實(shí)驗(yàn)4-常微分方程數(shù)值解,(b)Runge-Kutta方法基本思想是用小區(qū)間xn,xn+1上的若干個(gè)點(diǎn)的導(dǎo)數(shù)的線性組合代替方程(2)右端的,一般形式為(3)滿足并使(3)的局部截?cái)嗾`差-L級(jí)p階Runge-Kutta公式,實(shí)驗(yàn)4-常微分方程數(shù)值解,(2)常微分方程組和高階方程的數(shù)值方法歐拉方法和Runge-Kutta方法可直接推廣到求常微分方程組,如對(duì)歐拉公式為Runge-Kutta公式有類似的形式.對(duì)高階方程(5)需先降階化為一階常微分方程組,降階方法不唯一.簡(jiǎn)單、常用的方法是令y1=y,將(5)化為,實(shí)驗(yàn)4-常微分方程數(shù)值解,2.Runge-Kutta方法的MatLab實(shí)現(xiàn)對(duì)微分方程(組)的初值問(wèn)題Runge-Kutta方法用MatLab命令實(shí)現(xiàn):t,x=ode23(f,ts,x0,options)%用3級(jí)2階Runge-Kutta公式t,x=ode45(f,ts,x0,options)%用5級(jí)4階Runge-Kutta公式命令的輸入f是待解方程寫成的函數(shù)M文件:functiondx=f(t,x)Dx=f1;f2;fn;,實(shí)驗(yàn)4-常微分方程數(shù)值解,2.Runge-Kutta方法的MatLab實(shí)現(xiàn)舉例：仿真模擬著名的Lorenz系統(tǒng)混沌圖其中，先建立一個(gè)函數(shù)M文件functionxdot=lorenz(t,x)sigma=10;r=28;row=8/3;xdot=-sigma*x(1)+sigma*x(2);(r-x(3)*x(1)-x(2);x(1)*x(2)-row*x(3);,實(shí)驗(yàn)4-常微分方程數(shù)值解,2.Runge-Kutta方法的MatLab實(shí)現(xiàn)畫出Lorenz系統(tǒng)圖clearall;clf;options=odeset(RelTol,1e-5,AbsTol,1e-5);tspan=0,100;x0=1,2,3;t,x=ode45(lorenz,tspan,x0,options);l=length(x(:,1);a=1;b=l;figure(1)plot3(x(a:b,3),x(a:b,1),x(a:b,2),b);gridon;%畫出三維相圖xlabel(z);ylabel(x);zlabel(y);figure(2)subplot(311);plot(t,x(a:b,1);%畫三分量演化圖subplot(312);plot(t,x(a:b,2)subplot(313);plot(t,x(a:b,3),實(shí)驗(yàn)4-常微分方程數(shù)值解,2.Runge-Kutta方法的MatLab實(shí)現(xiàn)作業(yè)報(bào)告：著名的Duffing系統(tǒng)（描述彈簧系統(tǒng)性質(zhì)）其中類似的，分別畫出F=1，2，3，4，6等時(shí)的相圖翻閱一些參考書，你能得到一些什么結(jié)論？,實(shí)驗(yàn)4-常微分方程數(shù)值解,3.實(shí)例問(wèn)題緝私艇追擊走私船海上邊防緝私艇發(fā)現(xiàn)距d公里處有一走私船正以勻速a沿直線行駛,緝私艇立即以最大勻速度v追趕,在雷達(dá)的引導(dǎo)下,緝私艇的方向始終指向走私船.問(wèn)緝私艇何時(shí)追趕上走私船?并求出緝私艇追趕的路線.,S,（1）建立模型,走私船初始位置在點(diǎn)(S0,0),行駛方向?yàn)閤軸正方向,緝私艇的初始位置在點(diǎn)(0,M0),在時(shí)刻t：走私船的位置到達(dá)點(diǎn):(S0+at,0)緝私艇到達(dá)點(diǎn)M(x,y),S,（2）模型求解（a）求解析解,令：,令：,（2）模型求解,(a)求解析解,當(dāng)y=0時(shí):,走私船a=0.4千米/秒，分別取v=0.6,0.8,1.0千米/秒時(shí)，緝私艇追趕路線的圖形。,clearall;clf;a=0.4;v=0.60.81.0;%取不同的速度r=0.4./v;t=20*r./(a*(1-r.2)%追上的時(shí)間fori=1:3y=20:-0.01:0;x(:,i)=-0.5*(-40*r(i)+20(-r(i)*(r(i)-1)*y.(1+r(i)+20r(i)*(r(i)+1)*y.(1-r(i)/(1-r(i)2);plot(x(:,i),y);axis(030020);holdonend,追趕時(shí)間分別為：T=60.0000,33.3333,23.8095（秒）,2),當(dāng),時(shí)，,緝私艇不可能追趕上走私船。,3）,，,，,當(dāng),時(shí)，,，,緝私艇不可能追趕上走私船。,（b）用MATLAB軟件求解析解,MATLAB軟件5.3以上版本提供的解常微分方程解析解的指令是Dsolve,完整的調(diào)用格式是:dsolve(eqn1,eqn2,.)其中eqn1,eqn2,.是輸入宗量,包括三部分:微分方程、初始條件、指定變量,若不指定變量,則默認(rèn)小寫字母t為獨(dú)立變量.書P-69,微分方程的書寫格式規(guī)定:當(dāng)y是因變量時(shí),用“Dny”表示y的n階導(dǎo)數(shù).,例求微分方程,的通解。,dsolve(Dy=x+x*y,x),Ans=-1+exp(1/2*x2)*C1,dsolve(Dx=1/2*(y/20)r-(20/y)r),x(20)=0,y),Ans=1/2*20(-r)*y(r+1)/(r+1)+1/2*20r/(r-1)*y*(1/y)r-20*r/(r2-1),（c）用MATLAB軟件防真,當(dāng)建立動(dòng)態(tài)系統(tǒng)的微分方程模型很困難時(shí),我們可以用計(jì)算機(jī)仿真法對(duì)系統(tǒng)進(jìn)行分析研究.所謂計(jì)算機(jī)仿真就是利用計(jì)算機(jī)對(duì)實(shí)際動(dòng)態(tài)系統(tǒng)的結(jié)構(gòu)和行為進(jìn)行編程、模擬和計(jì)算,以此來(lái)預(yù)測(cè)系統(tǒng)的行為效果.,追趕方向可用方向余弦表示為：%兩點(diǎn)形成的向量的方向余弦,時(shí)間步長(zhǎng)為,，,則在時(shí)刻,時(shí)：,仿真算法：,第一步：設(shè)置時(shí)間步長(zhǎng),速度a,v及初始距離d，,第二步：,計(jì)算動(dòng)點(diǎn)緝私艇D在時(shí)刻,時(shí)的坐標(biāo),，,計(jì)算走私船R在時(shí)刻,時(shí)的坐標(biāo),，,第三步：計(jì)算緝私艇與走私船這兩個(gè)動(dòng)點(diǎn)之間的距離：,根據(jù)事先給定的距離，判斷緝私艇是否已經(jīng)追上了走私船，從而判斷退出循環(huán)還是讓時(shí)間產(chǎn)生一個(gè)步長(zhǎng)，返回到第二步繼續(xù)進(jìn)入下一次循環(huán)；,第四步：當(dāng)從上述循環(huán)退出后，由點(diǎn)列,和,可分別繪,制成兩條曲線即為緝私艇和走私船走過(guò)的軌跡曲線。,緝私艇初始位置,，,走私船初始位置,追擊問(wèn)題的數(shù)值模擬(P-66)clear;clf;d=120;v=90;a=80;s0=8;%給出初始條件T=10;dt=0.001;%選取時(shí)間區(qū)間T(可以偏大一點(diǎn)),時(shí)間微元dtt=0:dt:T;%離散時(shí)間表tn=length(t);%離散時(shí)間表t長(zhǎng)度x(1)=0;y(1)=d;s(1)=s0;%初始位置、初始距離fori=1:nx(i+1)=x(i)+v*dt*(s0+a*t(i)-x(i)/sqrt(s0+a*t(i)-x(i)2+y(i)2);y(i+1)=y(i)+v*dt*(-y(i)/sqrt(s0+a*t(i)-x(i)2+y(i)2);%遞推算式、d=sqrt(s0+a*t(i)-x(i)2+y(i)2);%t(i)時(shí)刻的距離ifd10的方程便可認(rèn)為是剛性方程,實(shí)際問(wèn)題中可出現(xiàn)s達(dá)的情況.剛性是問(wèn)題本身的性質(zhì),與解法無(wú)關(guān).但正是由于這種性質(zhì),用數(shù)值方法求解時(shí)需要計(jì)算到最慢瞬態(tài)解衰減成可忽略的小量為止,使得積分區(qū)間很長(zhǎng),而為保證計(jì)算的穩(wěn)定性,當(dāng)最快瞬態(tài)解的很大時(shí),又必須使步長(zhǎng)充分小,這就出現(xiàn)了在大區(qū)間上用小步長(zhǎng)計(jì)算的困難情況.,實(shí)驗(yàn)4-常微分方程數(shù)值解,4.剛性現(xiàn)象與剛性方程MatLab求解Matlab中求解常微分方程的命令ode23,ode45.由于其步長(zhǎng)是按穩(wěn)定性要求和指定的精度加以調(diào)整的,所以用它們解剛性微分方程時(shí)步長(zhǎng)會(huì)自動(dòng)變小,對(duì)于大的區(qū)間會(huì)導(dǎo)致計(jì)算時(shí)間很長(zhǎng).Matlab中有專門求解剛性方程的命令ode23s,ode15s,用法與ode23,ode45相同.例解方程特征根,剛性比.解析解為,實(shí)驗(yàn)4-常微分方程數(shù)值解,4.剛性現(xiàn)象與剛性方程MatLab求解functiondx=stiff1(t,x)dx=x(1)+2*x(2);-(106+1)*x(1)-(106+2)*x(2);t=0:0.1:1;%t=0:0.1:10;x1=(106/4+1)*exp(-t)-exp(-106*t);x2=-(106/4+1)*exp(-t)+(106+1)/2*exp(-106*t);A=t;x1;x2x0=106/4,106/4-1/2;t,x=ode23s(stiff1,t,x0);%很快出結(jié)果B=t,x%t,y=ode23(stiff1,t,x0);%幾十秒出結(jié)果%C=t,y%要計(jì)算0,10才能保證精度,ode23薛要很長(zhǎng)很長(zhǎng)時(shí)間.,謝謝同學(xué)們！,