數(shù)學(xué)實驗作業(yè)四

1、數(shù)學(xué)實驗作業(yè)四題目：P120. 1.d); 4); 5)日期:2003-3-23【實驗內(nèi)容】：用歐拉方法和龍格-庫塔方法求微分方程數(shù)值解,畫曲解的圖形，對結(jié)果進行比較分析。以下方程供選擇:a)y' = y +2x, y(0) =1,0 Ex E1,精確解 y =3ex2x 2.b)22y =x -y ,y(0) =0,或y(0) =1.c)x2y''+xy' + (x2n2)y=°,yC)=2,y'(2) = q(Bessel 方程，令n=1/2),精確解2 二 y=sinx x3d) y'' y cosx =0,y(0) =1

2、,y'(0) =0,哥級數(shù)解,1 22 49 6 55 8y = 1 - x '' x - x "x2!4!6!8!【模型建立】：原方程化為 y'' = -ycosx,y(0) fy'B)=。x(l) = y記兇2)=y',并用t代替x,則原方程化為:x(1)' =x(2)x(1)建=1x(2)' = x(1)cost .日x(2)|t包=0 .于是可以用龍格-庫塔法求解?！灸Ｐ颓蠼狻浚河肕atlab作龍格庫塔法求解：% chapter 4 -1.d%此函數(shù)是微分方程組function Xdot=ch41dfun

3、(t,x)Xdot=x(2),-cos*x'%ch41d.mfunction I=ch41d(a)x0=1,0't,x=ode45('ch41dfun',0,a,x0);y=x(:,1);plot(t,y,'r');hold on;gtext('y');y1=1-1/jiecheng(2)*t.A2+2/jiecheng(4)*t.A4-9/jiecheng(6)*t.A6+55/jiecheng(8)*t.A8;plot(t,y1,'b');gtext('y1');t,y,yihold off;運

4、行程序可以得到：(只保留圖像，數(shù)據(jù)省略)圖1>> ch41d(10)(數(shù)據(jù)省略，yl圖象也省略)圖25>> ch41d(15)>> ch41d(20)圖47【結(jié)果分析】由圖1得：y（龍格庫塔方法）和y1（級數(shù)近似解）在。到 大約1.5的區(qū)間內(nèi)是完全吻合的，從 x=1.5之后，兩條曲線 開始分離。之后y的變化趨勢可見圖2-圖4,呈幅度越來 越大的上下震蕩。對y取絕對值后再取對數(shù)，得到其圖象為：由此圖可見：|y|基本是隨x指數(shù)增加的，并上下振動，可 以猜測|y|接近于ex+osx的形式。四：求解“火箭發(fā)射”問題。【模型分析】從上圖可以看生，開始的幾個點基本呈線性

5、； 而以后的幾模型假設(shè):1) xwx0時，F(xiàn) 與 x 呈正比，F(xiàn)(x) = kx;模型建立:采用最小二乘法，目標函數(shù)為模型求解:為求生J的最小值，應(yīng)保證 編程語句如下：function k,a,b,c=number8(m) x = 1,2, 4, 7, 9,12,13,15,17;運行程序可以得到:x0 = 7.9342同時得到擬合曲線如下:結(jié)論:1.7000x x < 7,9342F =2-0,097僅2 3.2404x - 6.0417 x _ 7,9342五：一只小船渡過寬為 d的河流（見下圖），目標是起點 A正對著的另一岸 B點。已知河水流速 v1與船在靜水中的速度v2之比為ko

6、a)建立小船航線的方程，求其解析解。b)設(shè)d=100m,v1=1m/s,v2=2m/s,用數(shù)值解法求渡河所需時問、任意時刻小船的位置及航行曲線，作圖，并與解析解比較。c)若流速v1為0,0.5,1.5,2(m/s),結(jié)果將如何？【模型建立】：如圖，以B為原點，沿河岸向右為 x軸正向，垂直河岸 向下為y軸正向，建立坐標系。設(shè)在 t時刻，船在x方向上 的位移是x(t),在y方向上的位移是y(t),則：在t時刻，船在 X方向上的速度是X，在y方向上的速度是y'，將船的速 度V和靜水速度V1在x,y軸方向上分解，可得：Vx =V1 -V2sin FVy - -V 2cos 1又因為船頭始終指向

7、Bxtan =一點，所以 y,所以:Vx =V1 -V2Vy = -V2xx2 y2yx x'(t) =V1 -V2 x2 y2 y'(t) =N2= xy2這就是本題的微分方程。初始條件為:x(0)=0,y(0)=-d【模型求解】：1.解析解：x ;：；cosu令1y = Psine ,將直角坐標化成極坐標,1 P'cos - Psin 二二'=V1-V2cos1I i 'sin 二 、cosi 二'-V2sin1;則原微分方程化成化簡后為:'=V1cos -V2H'； '=V1sinu消去t,得:'V2二(-c

8、tg V1于是解得:一2sin日，其中V2R =V12.數(shù)值解法：根據(jù)上述微分方程，用龍格-庫塔方法求解:程序如下：%ch45fun.m%第四章第5題%此函數(shù)是微分方程組function Xdot=ch45fun(t,x,V1,V2)%x(1)代表x;x(2)代表y%加入限制條件，防止無限循環(huán)if (norm(x)>1e-5)Xdot=V1-V2*x(1)/sqrt(x(1)A2+x(2)A2),-V2*x(2)/sqrt(x(1)A2+x(2F2)'elseXdot=0,0'end%ch45.m%第四章第 5 題function I=ch45(V1,V2,d,a)x0=

9、0,-d'% 初始條件t,x=ode45(ch45fun,0,a,x0,V1,V2);plot(x(:,1),x(:,2),'r');hold on;%作圖t,x(:,1),x(:,2)%打印數(shù)據(jù)%下面是作出精確解的圖象seta=linspace(-pi/2,0,100);rou=d*(abs(tan(seta/2)(V2/V1)./sin(seta);xp=-rou.*cos(seta);yp=-rou.*sin(seta);plot(xp,yp,'b');hold off;需要說明：在ch45fun.m 中，加入了 (norm(x)>1e-5)

10、 的限制條件，以保證在船離B 點足夠近時中止運算。否則無法得出正確結(jié)果。運行結(jié)果：(中間數(shù)據(jù)省略)可知渡河時航線如圖，所用時間為： 66.7秒當(dāng)v1=0,0.5,1.5,2時，結(jié)果分別為（圖像）:v1=0 10 20 30 40-50 BO 70-80 901001-0.8-0 6 -a J-口N00 20.4 O S 0.81說明在靜水中，船沿直線到達B點。這與直觀經(jīng)驗相符合渡河時間為50秒；v1=0.5渡河時間為：53.3秒v1=1.517渡河時間為：114.3秒v1=2渡河時間：從圖上看由，到t=1000秒時，船已到達對岸,但是并沒有到達 B點，而是在B點下游50米處。由于船頭 指向B點

11、，即船頭指向逆流方向， 且船速（靜水）等于水速, 可知船將保持原地不動。也就是說，船永遠到達不了B點。所以渡河時間為無窮大。v1=2.5渡河時間：與v1=2時情況類似，船能到達對岸，但是是在 B點下游。由于船速（靜水）小于水速，船將被水沖得順流 而下，同樣永遠到達不了 B點?！窘Y(jié)果分析】1 . 以上圖中紅線是龍格庫塔方法的結(jié)果，藍線是解析解的結(jié)果。從圖中可以看出，兩種方法的結(jié)果基本上是相符的。但在接近 B 點時，解析解將無法得到正確解（因為會出現(xiàn)分母為零的情況，而Matlab 不會計算極限） 。當(dāng) V1=0 時，解析解也得不到正確解（同樣是因為出現(xiàn)了分母為零） 。2 .隨著V1的變化，船的航線也在變化。當(dāng)