常微分方程的數(shù)值解

4 ——重要1. x2y+xy+(x2-n2)y=02x 22xy( )2

(2)=-

(n=0.5)y=sinxy1=yy2=y1=yn=0.5代入，則原轉(zhuǎn)化為：y=y1 2y 1y'2(12 x

)y4x2 1 2functiondy=fangcheng(x,y)dy=[y(2);-y(2)/x-(1-1/(4*x^2))*y(1)];xs=pi/2:0.1:100;y0=[2,-2/pi];[x,y]=ode45(@fangcheng,xs,y0);[x,y]plot(x,y(:,1)),gridy1(2)=2,y2(2)functiondy=fangcheng(x,y)dy=[y(2);-y(2)/x-(1-1/(4*x^2))*y(1)];xs=pi/2:0.1:100;y0=[2,-2/pi];[x,y]=ode45(@fangcheng,xs,y0);[x,y]plot(x,y(:,1)),grid20012y(x

y(x)x

ye

y(0)

y(0)

用法算出y(1)= （精確到4）,你用法是 ，調(diào)用的functiondy=fangcheng2(x,y)dy=[y(2);y(2)*sin(x)-y(1)*exp(x)];xs=0:0.01:1;functiondy=fangcheng2(x,y)dy=[y(2);y(2)*sin(x)-y(1)*exp(x)];xs=0:0.01:1;y0=[1,0];[x,y]=ode45(@fangcheng2,xs,y0);[x,y]20001. y(x

y(x)

x

y(0)

y(0)0,y(1)= ,你是 Matlab是 functiondy=fangcheng1(x,y)dy=[y(2);y(1)*sin(x)];xs=[0,1];y0=[1,0];[x,y]=ode45(@fangcheng1,xs,y0);functiondy=fangcheng1(x,y)dy=[y(2);y(1)*sin(x)];xs=[0,1];y0=[1,0];[x,y]=ode45(@fangcheng1,xs,y0);[x,y]題目11. x2y+xy+(x2-n2)y=02x 22xy( )2

(2)=-

()ynxy1=yy2=y1=yn=0.5代入則原程轉(zhuǎn)化：y=y1 2y'y2 1 )y2 x 4x2 1 2y1(2)=2,y2(2)=-向前歐拉公式：y(n+1)=y(n)+h*y(n)1 1 2y(n) 1y(n+1)=y(n)+h*[2

)y(n)]2 2 x(n)2

4x(n1)*x(n1) 1y)1

2(0)=-

)=

+n*h編程如下：clearallx=[pi/2:0.1:pi/2+100-0.05];y1=2;y2=-2/pi;fori=1:999y1(i+1)=y1(i)+0.1*y2(i);y2(i+1)=y2(i)-0.1*(y2(i)/x(i)+(1-0.25/x(i)^2)*y1(i));endplot(x,y1),grid3020100-10-20-30-400 20 40 60 80 100 120一樣2

100左右結(jié)束，中間的點(diǎn)也太一樣同：程序：clearallx=[pi/2:0.01*pi:30*pi];y1=2;fori=1:2950y1(i+1)=y1(i)+0.01*pi*y2(i);y2(i+1)=y2(i)-0.01*pi*(y2(i)/x(i)+(1-0.25/x(i)^2)*y1(i));endplot(x,y1),grid21.510.50051150 10 20 30 40 50 60 70 80 90 1001—庫(kù)塔y=y1 2y'y2 1 )y2 x 4x2 1 2y( )=2,y(1 2 2

)=-編寫(xiě)M文件：functiondx=suan(t,x) %suan函M表示組4級(jí)5—庫(kù)塔進(jìn)行編：clearts=[pi/2:0.1:100];%使—庫(kù)塔plot(t,x(:,1)),grid,title('庫(kù)塔'),pause2x%y=sin2xy=sin(t).*sqrt(2*pi./t)下：龍格庫(kù)塔方法21.510.50051150 10 20 30 40 50 60 70 80 90 100解21.510.50051150 10 20 30 40 50 60 70 80 90 100（4級(jí)5階2xy=sinx 隨x2x蕩x稍大的情況下x1階的度，45階5階的度一般情況下肯定用后者計(jì)算更接近真所還是應(yīng)該用—庫(kù)塔法好。題目2【題目 （課本習(xí)題第四章第 題）小型火箭初始重量為其中包括燃料?；鸺Q直向時(shí)燃料燃燒率為由此產(chǎn)生的推力設(shè)火箭升空氣阻力正于速度的平例系為高【模型假設(shè)】基于對(duì)題目的及需要作下假設(shè)：假設(shè)火箭飛行過(guò)中重力加速度的不變恒為 （由于火箭；；】：①盡到 個(gè)用fFFNkk=0.4kg/m=g初始s燒率t發(fā)射為 ；律：又由整理得有組據(jù) ②火箭燃料燃盡引擎關(guān)閉時(shí)火箭不再受到推力 的作用同理分析得到據(jù) 上式中為火箭體的質(zhì)量）【模型解】對(duì)模型進(jìn)行函數(shù)抽象題意對(duì)火箭飛行的兩個(gè)階段分別進(jìn)行函數(shù)抽象，并編寫(xiě)兩個(gè)函數(shù) 文件記火箭引擎工作階段： 即第一階段的度 即第一階的高度 火箭引擎關(guān)閉階段： 即第二階段的即第二階段的高度 得到火箭飛行的兩個(gè)階段的函數(shù) 文件如下：作業(yè)題 火箭引擎工作階段函數(shù)文件源程序火箭引擎工作階段： 即第一階段的度 即第一階段的高度設(shè)定各條件值平方的（）式列寫(xiě)微分方程： 度 、調(diào)整首先由已知可求出工時(shí)已知燃料總為 ，燃料燃燒率為 ，，為工時(shí)對(duì)于時(shí)達(dá)到最大為 用 編求解出對(duì)應(yīng)代碼略，得到下表：時(shí)間t與的數(shù)據(jù)關(guān)系（用于的終點(diǎn)試探）v0000000001001307020026170300393004005245050065636000026725760100257669602002486796030024023171000312471100214371200116371300018371.400-0.79771.500-1.77771.600-2.75879.800-120.86379.900-123.72180.000-126.671=71.3s00-6060s-71.3s。（在）值解中編寫(xiě)如下代碼：作業(yè)題腳本文件源程序設(shè)定第一階間終點(diǎn) 步長(zhǎng)設(shè)定第二階間終點(diǎn) 步第一階賦初值高度調(diào)用 計(jì)算第一階第一階高用第一階計(jì)算結(jié)果計(jì)算第二階初值調(diào)用 計(jì)算第二階第一階高將 組第一剪切否則與 合并第一階最后一會(huì)和第二階初始重復(fù)將第一階和第二階間高進(jìn)合并初始質(zhì)量；燃料燃燒率；；空氣阻力與平方比例系；重力加； 體質(zhì)量按式 求解第一階加

按式 求解第二階加速標(biāo)記標(biāo)題標(biāo)記標(biāo)題標(biāo)記標(biāo)題曲線用不同顏色以區(qū)分間 度 、標(biāo)記標(biāo)題引擎工引擎關(guān)閉例曲線用不同顏色以區(qū)分間 度 、標(biāo)記標(biāo)題引擎工引擎關(guān)閉例曲線用不同顏色以區(qū)分間 度 程 、標(biāo)記標(biāo)題引擎工引擎關(guān)閉例顯示、、、表為 文檔以】①由以上程序可得在以及、、如下面表、速和加速ht1(m)vt1(m/s)267at1(m/s2)09170ht1(m)vt1(m/s)0at1(m/s2)98000②由 得全過(guò)程隨變化情況如下面表 為節(jié)省篇幅只選取了部分以及圖 圖所示：第一階段：工作階段全過(guò)程第一階段：工作階段0000000000001305701000065130713086020002622617131140300058939301314104001048524513167050016386563131930600236078841321707003215920713240080042021053213263597001210964526698409195980012136350267075091859900121630632671660918二階段：火箭引擎關(guān)閉階段60000121897852672570917601001221602625766992792602001224133824867987101603001226578024023181938707001311348559099844708001311402649259830709001311447039439819710001311481529619811711001311506219809805712001311521110009801713001311526200209800全過(guò)程）全過(guò)程）加全過(guò)程）。得到下面：分段表示）度和加速情對(duì)比。【結(jié)果分析】數(shù)據(jù)表和曲線變化情整個(gè)過(guò)程中明顯分為兩個(gè)，對(duì)應(yīng)于引擎工作和引擎關(guān)閉兩個(gè)。見(jiàn)，從發(fā)射到燃料用盡，斷增大，但第一階增加速增加速明顯曲線斜率變化與速曲線變化是相一致。由行。由圖可見(jiàn)，加速度大小在第一階段逐漸減小，直至變?yōu)?。這是因?yàn)殡S可見(jiàn)，求解結(jié)果與理論分析是一致的。【求解解析表達(dá)式并與數(shù)值解對(duì)比】下面以火箭引擎關(guān)閉階段的模型為例求解其高度速度加速度值的解析解，并與 用龍格庫(kù)塔方法求出的數(shù)值解進(jìn)行對(duì)比。由之前的分析，火箭在第二階段滿足的微分方程組為加速度 。該方程組的第一個(gè)方程一個(gè)因變得由理論分析，時(shí)可以于 的，可式的解的解析式為第①問(wèn)結(jié)果所到的燃料燃盡時(shí)的速度值作為第二階段的值，即。，可得其與數(shù)計(jì)算數(shù)據(jù)結(jié)果對(duì)比如下面：數(shù)計(jì)算結(jié)果（5數(shù)計(jì)算結(jié)果（54龍格庫(kù)塔法ht(m)速vt(m/s)0at(m/s2)-9.8000htvt(m/s)at(m/s2)】。4tl裝0.08lbfs/ft現(xiàn)當(dāng)圓桶40ft/s時(shí)就會(huì)海底沖撞生破裂。建立決上述問(wèn)微分程模型。兩種法微分程并回答誰(shuí)贏了官司。問(wèn)分根牛頓律以上描述問(wèn)建立個(gè)關(guān)S(t)v微分程組利龍格庫(kù)塔法能給微分程組以利法模型假設(shè)】，v量GR8微分程組dsdtvdvGFRGFk*vdt m mS(0)v(0)0-。（一）求【算設(shè)計(jì)】有兩個(gè)思路：思路一是直接上面關(guān)于下沉距離S和下沉速度v的。思路二是-關(guān)于速度v 和時(shí)間t的一階dvdt

GFRm

GFkvm

再辛普森公式速度求積S；分別求SvSvS=300ft時(shí)速sst，從而判斷官司的勝負(fù)。functionwaste的函M文件m=527.436*0.4536;圓桶所受重力圓桶所受浮力%%clearall%終點(diǎn)試探%初始默認(rèn)誤差限5級(jí)4階%t,S(t),v(t)plot(t,x(:,2)),grid,title('v(t)-t線xlabel('t');ylabel('v(t)');figure; %在新窗口中作圖0%按照輸S(t),v(t)gtext('S(t)'),gtext('v(t)')holdon,Smax=300*0.3048;plot([0,25],[Smax,Smax],'r-'),holdoff%S=Smax的圖形以便比對(duì)do,'d%v=Vmax【計(jì)算結(jié)果】B：

v(t)-t曲線00 t1 速度時(shí)間像SmaxS(t)SmaxS(t)Vmaxv00 5 2functiondx=waste1(t,x,G,F,m,b)dx=[(G-F-b*x)/m];%列寫(xiě)微分方程clearall;%ts=0:h:2000;%%%%b=0.08*0.4536*9.8/0.3048; %比例系數(shù)opt=odeset;[t,x]=ode45(@waste1,ts,x0,opt,G,F,m,b);%調(diào)用5級(jí)4階龍格-庫(kù)塔方程%v(t)曲線xlabel('t');ylabel('v(t)');%用辛普森公式對(duì)速度積分求下沉深%間25s沉到海底fori=0:2:10*Ty=x(1:(i+1));end

k=length(y);%i=[0:2:10*T];figure;Smax=300*0.3048;Vmax=40.*0.3048;plot(x(i+1),z4'),grid;holdon,plot([0,25],[Smax,Smax],'r:'),grid %Smaxholdon,plot([Vmax,Vmax],[0,350],'m:'),grid,holdoff %legend('S(v)','Smax','Vmax');xlabel('v');ylabel('S(v)');figure;; %Stgrid;xlabel('t');ylabel('S(t)');【輸出結(jié)果】速度-時(shí)間圖像

00 t3S(v)S(v)S00 5 v4底會(huì)破裂。間S

00 t5【算設(shè)計(jì)】對(duì)于上面的組dsdtvdvGFRGFk*vdt m mS(0)0,v(0)0GF

527.436*0.4536*9.8470.327*0.4536*9.8令a= m

= 527.436*0.4536

=1.06111k 0.08*0.4536*9.8/0.3048b=m= 527.436*0.4536

=4.87676*10-3則原組等價(jià)于

dv

*）v(0)=0s(t)一階線性常(*)式的特征為

tvt)dt0得特征根為

+b=0*)A**9A0-s(t)

tvt)dt0-0【程序】clearallt=[0:0.1:25];Smax=300*0.3048;Vmax=40*0.3048;v=-217.59*exp(-4.87676*t*0.001)+217.59;s=44618*exp(-4.87676*t*0.001)+217.59*t-44618;plot(t,v,'b-',t,s,'g-'),gridlegend('速度-時(shí)間曲線','下沉深度-時(shí)間曲線')holdon,plot([0,25],[Smax,Smax],'r:'),gridholdon,plot([Vmax,Vmax],[0,350],'m:'),grid,holdoffgtext('Smax'),gtext('Vmax')【計(jì)算結(jié)果】圖速度速度時(shí)間曲線下沉深度時(shí)間曲線SmaxVmax00 5 圖51時(shí)間曲線看到開(kāi)時(shí)圓桶速度隨著時(shí)間逐漸增大加速度曲線t（。v(t)-t40ft/s(12.2m/s)所用的時(shí)間。結(jié)合S=300ft或判斷速度達(dá)到25。從圖中可見(jiàn)，圓桶還沒(méi)沉到海會(huì)與海底碰撞破裂。顯然是工程師們贏了官司！【結(jié)果的實(shí)際意義】構(gòu)造數(shù)學(xué)模型可以分析和模擬實(shí)際生題?？梢阅M圓桶下沉的速度和深度與下沉?xí)r間的上實(shí)際生的作。題目5為 圖 速與中的速度 為 。小的數(shù)學(xué)模型，其析。設(shè) ，作圖，并與析較。若速 為 結(jié)果將何。【題分析】

圖4.6度分到平與垂直向，利用數(shù)值和析列寫(xiě)微分程，即能到小程。【模型假】假速恒定，上無(wú)遇到障礙物阻礙素。依題意，v1,v2，為頭的向始終目Bv2。vv1v2BxyA0Bθx2y2x2y2則

yx2x2y2x2y2dxx2y2x2y2dyx2y2

1）(2)x(0)=0,y(0)=-d1÷ 去 得1x2y2k 1x2y2v1

v2x令u= ，ydx

d(u*

y*duu則dy dy dyu21ku21yu1

y

d( )k]2 d yy yx )k

d( )k]2 d y式即航析此外-庫(kù)塔列微程組就能得到問(wèn)題數(shù)值【程序】%以x(1)代表x，以x(2)代表yfunctiondx=boat(t,x,v1,v2)%建立名為boat的函M文件s=sqrt(x(1)^2+x(2)^2);dx=[v1-v2*x(1)/s;-v2*x(2)/s];%以向量形式表示clearallts=0:0.1:150;小船初始位置默認(rèn)誤差限%v1=1m/s[t,x]=ode15s(@boat,ts,x0,opt,1,2); %龍格-庫(kù)塔方%出t,x(t)plot(t,x),grid,title('v1=1m/s時(shí)的x(t),y(t)曲線')%按照輸出作x(t),y(t)gtext('x(t)'),gtext('y(t)')figure;plot(x(:,1),x(:,2)),grid,title('v1=1m/s時(shí)的航線')%v1=0[t,x]=ode15s(@boat,ts,x0,opt,0,2);[t,x]figure;plot(t,x),grid,title('v1=0m/s時(shí)的x(t),y(t)曲線')figure;plot(x(:,1),x(:,2)),grid,title('v1=0時(shí)的航線')%v1=0.5m/s[t,x]=ode15s(@boat,ts,x0,opt,0.5,2);[t,x]figure;plot(t,x),grid,title('v1=0.5m/s時(shí)的x(t),y(t)曲線')figure;plot(x(:,1),x(:,2)),grid,title('v1=0.5m/s時(shí)的航線')%v1=1.5m/s[t,x]=ode15s(@boat,ts,x0,opt,1.5,2);[t,x]figure;plot(t,x),grid,title('v1=1.5m/s時(shí)的x(t),y(t)曲線')figure;plot(x(:,1),x(:,2)),grid,title('v1=1.5m/s時(shí)的航線')%v1=2m/s[t,x]=ode15s(@boat,ts,x0,opt,2,2);[t,x]figure;plot(t,x),grid,title('v1=2m/sx(t),y(t)')figure;plot(x(:,1),x(:,2)),grid,title('v1=2m/s')v1=1m/s時(shí)y))y)：v1=1m/sx(t),y(t)線x(t)y(t)x(t)y(t)00 y=0x=0y=0By)000y0了B以渡河需66.65s。任意刻位置及行v1=1m/s線00 2 4 6 8 0x(t),y(t)0

vx(t),y(t)線0 5 0y。v線0-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 10。在在靜中度不變前提此以最短間到達(dá)對(duì)。所以渡最短間50s。（3）v1=0.5m/s輸出x(t),y(t)像所示：v1=0.5m/sx(t),y(t)線0-1000 syB。所以渡河所需間。任意刻的位置及航行曲線如下：v1=0.5m/s00 1 2 3 4 5 6 7 8 9 （4）v1=1.5m/sx(t),y(t)圖像如下所示：v1=1.5m/sx(t),y(t)曲線0-1000 ,1yB。所以渡河所需時(shí)間114.25s。任意時(shí)刻的位置及航行曲線如下：v1=1.5m/s時(shí)的航線00 5 （5）v1=2m/s，的x(t),y(t)圖像如下所示：v1=2m/s時(shí)的x(t),y(t)曲線0-1000 sBx0。當(dāng)∴B任意刻位置及航行曲線00 5 二用解法求解在前面我們已經(jīng)得了航線方程y yx ( )k

d( )k]2 d y面利用解法來(lái)求航線v1=1m/sv1=2m/s例?！境绦颉客ㄟ^(guò)試驗(yàn)發(fā)現(xiàn)本題ode45()求解則費(fèi)甚久所以以推測(cè)題微方程一剛性方程ode15s()求解。clearalld=100;k=0.5;y=-d:0.1:0;x=y/2.*((-y/d).^k-(-d./y).^k);plot(x,y,'r-'),gridxlabel('x'),ylabel('y'),title('v1=1m/s'),pause;%--ts=0:0.05:150;%-holdon,plot(x(:,1),x(:,2),’g:’),grid,holdoff%legend(‘,’)【計(jì)算結(jié)果】0-10

y -1000 2 4 6 8 xsk和s0-10

y 0 5 xv1v1v1Bxv≥在靜候無(wú)BK=1是否正分界上想正靜必須題目84-8．兩種群相互競(jìng)爭(zhēng)模型如下：rxx

y)t)

y

x y)1 n 1n1

2 2n n1 2y)分別為甲乙兩種群量rr為它們固有率nn為它1 2 1 2