《常微分方程數(shù)值解與matlab.ppt》由會(huì)員分享，可在線閱讀，更多相關(guān)《常微分方程數(shù)值解與matlab.ppt（25頁珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

實(shí)驗(yàn)4-常微分方程數(shù)值解,1.求解常微分方程數(shù)值方法介紹（1）一階微分方程求方程(1)的數(shù)值解,就是計(jì)算(精確)解在一系列離散點(diǎn)的近似值.通常取相等的步長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ì)微分方程(組)的初值問題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í)例問題緝私艇追擊走私船海上邊防緝私艇發(fā)現(xiàn)距d公里處有一走私船正以勻速a沿直線行駛,緝私艇立即以最大勻速度v追趕,在雷達(dá)的引導(dǎo)下,緝私艇的方向始終指向走私船.問緝私艇何時(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))+20^r(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*x^2)*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*20^r/(r-1)*y*(1/y)^r-20*r/(r^2-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ì)算,以此來預(yù)測(cè)系統(tǒng)的行為效果.,追趕方向可用方向余弦表示為：%兩點(diǎn)形成的向量的方向余弦,時(shí)間步長為,，,則在時(shí)刻,時(shí)：,仿真算法：,第一步：設(shè)置時(shí)間步長,,速度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è)步長，返回到第二步繼續(xù)進(jìn)入下一次循環(huán)；,第四步：當(dāng)從上述循環(huán)退出后，由點(diǎ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長度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í)際問題中可出現(xiàn)s達(dá)的情況.剛性是問題本身的性質(zhì),與解法無關(guān).但正是由于這種性質(zhì),用數(shù)值方法求解時(shí)需要計(jì)算到最慢瞬態(tài)解衰減成可忽略的小量為止,使得積分區(qū)間很長,而為保證計(jì)算的穩(wěn)定性,當(dāng)最快瞬態(tài)解的很大時(shí),又必須使步長充分小,這就出現(xiàn)了在大區(qū)間上用小步長計(jì)算的困難情況.,,,,,,,,,,實(shí)驗(yàn)4-常微分方程數(shù)值解,4.剛性現(xiàn)象與剛性方程MatLab求解Matlab中求解常微分方程的命令ode23,ode45.由于其步長是按穩(wěn)定性要求和指定的精度加以調(diào)整的,所以用它們解剛性微分方程時(shí)步長會(huì)自動(dòng)變小,對(duì)于大的區(qū)間會(huì)導(dǎo)致計(jì)算時(shí)間很長.Matlab中有專門求解剛性方程的命令ode23s,ode15s,用法與ode23,ode45相同.例解方程特征根,剛性比.解析解為,,,,,,,,,,,實(shí)驗(yàn)4-常微分方程數(shù)值解,4.剛性現(xiàn)象與剛性方程MatLab求解functiondx=stiff1(t,x)dx=[x(1)+2*x(2);-(10^6+1)*x(1)-(10^6+2)*x(2)];t=0:0.1:1;%t=0:0.1:10;x1=(10^6/4+1)*exp(-t)-exp(-10^6*t);x2=-(10^6/4+1)*exp(-t)+(10^6+1)/2*exp(-10^6*t);A=[t;x1;x2]x0=[10^6/4,10^6/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薛要很長很長時(shí)間.,,,,,,,,,,,謝謝同學(xué)們！,