打靶法(含Matlab程序)

1、西京學(xué)數(shù)學(xué)軟件實(shí)驗(yàn)任務(wù)書(shū)課程名稱數(shù)學(xué)軟件實(shí)驗(yàn)班級(jí)數(shù) 0901學(xué)號(hào)0912020107姓名李亞強(qiáng)實(shí)驗(yàn)課題微分方程組邊值問(wèn)題數(shù)值算法（打靶法，有限差分法）實(shí)驗(yàn)?zāi)康氖煜の⒎址匠探M邊值問(wèn)題數(shù)值算法（打靶法，有限差分法）實(shí)驗(yàn)要求運(yùn)用 Matlab/C/C+/Java/Maple/Mathematica等其中一種語(yǔ)言完成實(shí)驗(yàn)內(nèi)容微分方程組邊值問(wèn)題數(shù)值算法（打靶法，有限差分法）成績(jī)教師動(dòng)方向控制減速的推力，主要的控制量只有一個(gè)減速推力，減速 還會(huì)消耗燃料讓登月器的質(zhì)量減小。所以在極坐標(biāo)下系統(tǒng)的狀態(tài) 就是x =質(zhì)量m ,角度theta，高度r,角速度omega，徑速度v這五個(gè)量，輸入就是減速力F。先列微分方程

2、，dx/dt=f(x)+B*F , 其中 x 是 5*1 的列向量， 質(zhì)量 dm/dt=-F/2940 ，剩下幾個(gè)翻下極 坐標(biāo)的手冊(cè)。 把這個(gè)動(dòng)力學(xué)模型放到 matlab 里就能求解了， 微分 方程數(shù)值解用 ode45 。第一問(wèn) F=0 ，讓你求橢圓軌道非常容易。 注意附件 1 里說(shuō) 15 公里的時(shí)候速度是 1.7km/s 。算完以后驗(yàn)證一 下對(duì)不對(duì)，對(duì)的話就是他了，不對(duì)的話說(shuō)明這個(gè)橢圓軌道有進(jìn)動(dòng)， 到時(shí)再說(shuō)。(2) 算出軌道就能計(jì)算減速力了。這時(shí)候你隨便給個(gè)常數(shù)減速力 到方程里飛船八成都能降落，但不是最優(yōu)解。想想整個(gè)過(guò)程，開(kāi) 始降落之前飛船總機(jī)械能就那么多，你需要對(duì)飛船做負(fù)功讓機(jī)械 能減到

3、0。題目里寫(xiě)發(fā)動(dòng)機(jī)噴出翔的相對(duì)速度是一定的， 直覺(jué)告訴 我飛船速度快的時(shí)候多噴一些速度慢的時(shí)候少噴一些，可以提高 做負(fù)功的效率。但是多噴也不能超過(guò)上限 7500N ，所以這就是一 個(gè)帶約束優(yōu)化問(wèn)題， matlab 里邊有專用的優(yōu)化函數(shù)， 用 fmincon就好。找出最優(yōu)解以后把過(guò)程畫(huà)出來(lái)，看看 F可不可以是那5個(gè) 狀態(tài)量的線性組合，如果是的話就非常 happy ，不是的話再說(shuō)。 三四階段你可以扯點(diǎn)圖像識(shí)別，什么二維復(fù)利葉分解找平坦區(qū)域， 怎么一邊下降一邊根據(jù)自身狀態(tài)調(diào)整路徑之類(lèi)的。 五六階段還真不知道說(shuō)什么。一二階段肯定是重點(diǎn)啦（3） 誤差分析其實(shí)還挺難的。可能的誤差來(lái)源是地球的引力，月 亮繞

4、地球向心加速度，太陽(yáng)的引力 （可能會(huì)很小 ），對(duì)自身速度、角 度的測(cè)量誤差 （比如你測(cè)出自身當(dāng)前速度 100m/s 但實(shí)際上是 105m/s） ，控制的時(shí)候 F 大小以及角度的誤差 （比如你想朝正前方 向噴 2000N 但實(shí)際上偏了 2 度而且 F=2010N 之類(lèi) ）。上一問(wèn)已 經(jīng)求出了最優(yōu)控制策略和飛船路線，把這些擾動(dòng)加進(jìn)去以后算出 新的路線減掉理想路線求偏差，然后隨便用個(gè)卡爾曼濾波器把誤 差給校正All for Joy2014/9/13 11:14:38老師的思路，求大神解答給我一份呀實(shí)驗(yàn)二十七實(shí)驗(yàn)報(bào)告一、實(shí)驗(yàn)名稱： 微分方程組邊值問(wèn)題數(shù)值算法(打靶法，有限差分 法)。二、實(shí)驗(yàn)?zāi)康模?進(jìn)

5、一步熟悉微分方程組邊值問(wèn)題數(shù)值算法 (打靶法， 有限差分法)。三、實(shí)驗(yàn)要求： 運(yùn)用 Matlab/C/C+/Java/Maple/Mathematica 等其中一種語(yǔ)言完成程序設(shè)計(jì)。四、實(shí)驗(yàn)原理：1打靶法：對(duì)于線性邊值問(wèn)題y p(x)y q(x)y f(x) x a,b (1)y(a) , y(b) (1)假設(shè) L 是一個(gè)微分算子使： Ly y p(x)y q(x)y 則可得到兩個(gè)微分方程：Ly1 f (x) ， y1(a)， y1(a) 0y1 p(x)y1 q(x)y1 f (x) ， y1(a)， y1(a) 0 (2)Ly2 0， y2(a) 0， y2(a) 1y2 p(x)y2 q

6、(x)y2 0， y2(a) 0， y2(a) 1 (3)方程(2), (3)是兩個(gè)二階初值問(wèn)題 假設(shè)yi是問(wèn)題(2)的解，y是問(wèn)題(3)的解，且y2(b) 0,則線性邊值問(wèn)題(1) 的解為：y(x) %(x) 警 y2(x)。y2(b)2 .有限差分法：基本思想是把連續(xù)的定解區(qū)域用有限個(gè)離散點(diǎn)構(gòu)成的網(wǎng)格來(lái)代替，這些離散點(diǎn)稱作網(wǎng)格的節(jié)點(diǎn)；把連續(xù)定解 區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù) 來(lái)近似；把原方程和定解條件中的微商用差商來(lái)近似，積分用積分和來(lái)近似，于是原微分方程和定解條件就近似地 代之以代數(shù)方程組，即有限差分方程組，解此方程組就可以得到原問(wèn)題在離散點(diǎn)上的近似解。然后再利用插

7、值方 法便可以從離散解得到定解問(wèn)題在整個(gè)區(qū)域上的近似解。五、實(shí)驗(yàn)內(nèi)容：%線性打靶法fun ctio n k,X,Y,wucha,P=xxdb(dydx1,dydx2,a,b,alpha,beta,h)n=fix(b-a)/h); X=zeros(n+1,1); CT仁alpha,。;Y=zeros( n+1,le ngth(CTI); Y1= zeros( n+1,le ngth(CTI);Y2=zeros( n+1,le ngth(CT1);X=a:h:b;Y1(1,:)= CT1;CT2=0,1;Y2(1,:)= CT2;for k=1:nk1=feval(dydx1,X(k),Y1(k,

8、:)x2=X(k)+h/2;y2=Y1(k,:)'+k1*h/2;k2=feval(dydx1,x2,y2);k3=feval(dydx1,x2,Y1(k,:)'+k2*h/2);k4=feval(dydx1, X(k)+h,Y1(k,:)'+k3*h);Y1(k+1,:)=Y1(k,:)+h*(k1'+2*k2'+2*k3'+k4')/6,k=k+1;endu=Y1(:,1)for k=1:nk1=feval(dydx2,X(k),Y2(k,:)x2=X(k)+h/2;y2=Y2(k,:)'+k1*h/2;k2=feval(dy

9、dx2,x2,y2);k3=feval(dydx2,x2,Y2(k,:)'+k2*h/2);k4=feval(dydx2, X(k)+h,Y2(k,:)'+k3*h);Y2(k+1,:)=Y2(k,:)+h*(k1'+2*k2'+2*k3'+k4')/6,k=k+1;endv=Y2(:,1)Y=u+(beta-u(n+1)*v/v(n+1)for k=2:n+1wucha(k)=norm(Y(k)-Y(k-1); k=k+1;endX=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:);P=k'

10、;,X',Y,wucha'plot(X,Y(:,1),'ro',X,Y1(:,1),'g*',X,Y2(:,1),'mp')xlabel('軸it x'); ylabel(' 軸it y')legend(' 是邊值問(wèn)題的數(shù)值解 y(x) 的曲線 ',' 是初值問(wèn)題 1 的數(shù)值解 u(x) 的曲線','是初值問(wèn)題2的數(shù)值解v(x)的曲線')title(' 用線性打靶法求線性邊值問(wèn)題的數(shù)值解的圖形 ')%有限差分法function k,A,

11、B1,X,Y,y,wucha,p=yxcf(q1,q2,q3,a,b,alpha,beta,h)n=fix(b-a)/h); X=zeros(n+1,1);Y=zeros(n+1,1); A1=zeros(n,n);A2=zeros(n,n); A3=zeros(n,n); A=zeros(n,n);B= zeros(n,1);for k=1:nX=a:h:b;k1(k)=feval(q1,X(k); A1(k+1,k)=1+h*k1(k)/2;k2(k)=feval(q2,X(k);A2(k,k)=-2-(h42)*k2(k);A3(k,k+1)= 1-h*k1(k)/2; k3(k)=fe

12、val(q3,X(k);endfor k=2:nB(k,1)=(h.A2)*k3(k);endB(1,1)=(h.A2)*k3(1)-(1+h*k1(1)/2)*alpha;B(n-1,1)=(h.A2)*k3(n-1)-(1+h*k1(n-1)/2)*beta;A=A1(1:n-1,1:n-1)+A2(1:n-1,1:n-1)+A3(1:n-1,1:n-1);B1=B(1:n-1,1);Y=AB1;Y1=Y' y=alpha;Y;beta;for k=2:n+1wucha(k)=norm(y(k)-y(k-1); k=k+1;endX=X(1:n+1); y=y(1:n+1,1);

13、k=1:n+1;y(x) 的曲線 ')wucha=wucha(1:k,:); plot(X,y(:,1),'mp')xlabel('軸it x'); ylabel('軸it y'),legend('是邊值問(wèn)題的數(shù)值解 title（' 用有限差分法求線性邊值問(wèn)題的數(shù)值解的圖形 '）, p=k',X',y,wucha'打靶法3.Matlab源代碼：創(chuàng)建 M 文件：function ys=dbf(f,a,b,alfa,beta,h,eps) ff=(x,y)y(2),f(y(1),y(2),x);

14、 xvalue=a:h:b;%x 取值范圍 n=length(xvalue) s0=a-0.01;% 選取適當(dāng)?shù)?s 的初值x0=alfa,s0;%迭代初值 flag=0;%用于判斷精度 y0=rk4(ff,a,x0,h,a,b);if abs(y0(1,n)-beta)<=epsflag=1;y1=y0;elses1=s0+1;x0=alfa,s1;y1=rk4(ff,a,x0,h,a,b);if abs(y1(1,n)-beta)<=eps flag=1;endendif flag=1while abs(y1(1,n)-beta)>epss2=s1-(y1(1,n)-be

15、ta)*(s1-s0)/(y1(1,n)-y0(1,n);x0=alfa,s2;y2=rk4(ff,a,x0,h,a,b);s0=s1;s1=s2;y0=y1;y1=y2;endend xvalue=a:h:b; yvalue=y1(1,:); ys=xvalue',yvalue' function x=rk4(f,t0,x0,h,a,b) %rung-kuta 法求每個(gè)點(diǎn)的近似值(參考大作業(yè)一) t=a:h:b;% 迭代區(qū)間 m=length(t); %區(qū)間長(zhǎng)度t(1)=t0;x(:,1)=x0;%迭代初值for i=1:m-1L1=f(t(i),x(:,i);L2=f(t(i)+h/2,x(:,i)'+(h/2)*L1);L3=f(t(i)+h/2,x(:,i)'+(h/2)*L2);L4=f(t(i)+h,x(:,i)'+h*L3); x(:,i+1)=x(:,i)'+(h/6)*(L1+2*L2+2*L3+L4); end4.舉例求二階非線性方程的邊值問(wèn)題： 在 matlab 控制臺(tái)中輸入： f=(x,y,z)(xA2+z*xA2);x0l=0;x0u=2*exp(-1);alfa=0;be