1-西京學(xué)數(shù)學(xué)軟件實(shí)驗(yàn)任務(wù)書課程名稱數(shù)學(xué)軟件實(shí)驗(yàn)班級數(shù)0901學(xué)號0912020107姓名李亞強(qiáng)實(shí)驗(yàn)課題微分方程組邊值問題數(shù)值算法（打靶法，有限差分法）實(shí)驗(yàn)?zāi)康氖煜の⒎址匠探M邊值問題數(shù)值算法（打靶法，有限差分法）實(shí)驗(yàn)要求運(yùn)用Matlab/C/C++/Java/Maple/Mathematica等其中一種語言完成實(shí)驗(yàn)內(nèi)容微分方程組邊值問題數(shù)值算法（打靶法，有限差分法）成績教師

動方向控制減速的推力，主要的控制量只有一個減速推力，減速還會消耗燃料讓登月器的質(zhì)量減小。所以在極坐標(biāo)下系統(tǒng)的狀態(tài)就是x‘=[質(zhì)量m，角度theta，高度r，角速度omega，徑速度v]這五個量，輸入就是減速力F。先列微分方程，dx/dt=f(x)+B*F，其中x是5*1的列向量，質(zhì)量dm/dt=-F/2940，剩下幾個翻下極坐標(biāo)的手冊。把這個動力學(xué)模型放到matlab里就能求解了，微分方程數(shù)值解用ode45。第一問F=0，讓你求橢圓軌道非常容易。注意附件1里說15公里的時候速度是1.7km/s。算完以后驗(yàn)證一下對不對，對的話就是他了，不對的話說明這個橢圓軌道有進(jìn)動，到時再說。

(2)

算出軌道就能計算減速力了。這時候你隨便給個常數(shù)減速力到方程里飛船八成都能降落，但不是最優(yōu)解。想想整個過程，開始降落之前飛船總機(jī)械能就那么多，你需要對飛船做負(fù)功讓機(jī)械能減到0。題目里寫發(fā)動機(jī)噴出翔的相對速度是一定的，直覺告訴我飛船速度快的時候多噴一些速度慢的時候少噴一些，可以提高做負(fù)功的效率。但是多噴也不能超過上限7500N，所以這就是一個帶約束優(yōu)化問題，matlab里邊有專用的優(yōu)化函數(shù)，用fmincon就好。找出最優(yōu)解以后把過程畫出來，看看F可不可以是那5個狀態(tài)量的線性組合，如果是的話就非常happy，不是的話再說。

三四階段你可以扯點(diǎn)圖像識別，什么二維復(fù)利葉分解找平坦區(qū)域，怎么一邊下降一邊根據(jù)自身狀態(tài)調(diào)整路徑之類的。

五六階段還真不知道說什么。一二階段肯定是重點(diǎn)啦

(3)

誤差分析其實(shí)還挺難的?？赡艿恼`差來源是地球的引力，月亮繞地球向心加速度，太陽的引力(可能會很小)，對自身速度、角度的測量誤差(比如你測出自身當(dāng)前速度100m/s但實(shí)際上是105m/s)，控制的時候F大小以及角度的誤差(比如你想朝正前方向噴2000N但實(shí)際上偏了2度而且F=2010N之類)。上一問已經(jīng)求出了最優(yōu)控制策略和飛船路線，把這些擾動加進(jìn)去以后算出新的路線減掉理想路線求偏差，然后隨便用個卡爾曼濾波器把誤差給校正

All

for

Joy

2014/9/13

11:14:38

老師的思路，求大神解答給我一份呀A2(k,k)=-2-(h.^2)*k2(k);A3(k,k+1)=1-h*k1(k)/2;k3(k)=feval(q3,X(k));endfork=2:nB(k,1)=(h.^2)*k3(k);endB(1,1)=(h.^2)*k3(1)-(1+h*k1(1)/2)*alpha;B(n-1,1)=(h.^2)*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=A\B1;Y1=Y';y=[alpha;Y;beta];fork=2:n+1wucha(k)=norm(y(k)-y(k-1));k=k+1;endX=X(1:n+1);y=y(1:n+1,1);k=1:n+1;wucha=wucha(1:k,:);plot(X,y(:,1),'mp')xlabel('軸\itx');ylabel('軸\ity'),legend('是邊值問題的數(shù)值解y(x)的曲線')title('用有限差分法求線性邊值問題的數(shù)值解的圖形'),p=[k',X',y,wucha'];打靶法3.Matlab源代碼：創(chuàng)建M文件：functionys=dbf(f,a,b,alfa,beta,h,eps)ff=@(x,y)[y(2),f(y(1),y(2),x)];xvalue=a:h:b;%x取值范圍n=length(xvalue)s0=a-0.01;%選取適當(dāng)?shù)膕的初值x0=[alfa,s0];%迭代初值flag=0;%用于判斷精度y0=rk4(ff,a,x0,h,a,b);ifabs(y0(1,n)-beta)<=epsflag=1;y1=y0;elses1=s0+1;x0=[alfa,s1];y1=rk4(ff,a,x0,h,a,b);ifabs(y1(1,n)-beta)<=epsflag=1;endendifflag~=1whileabs(y1(1,n)-beta)>epss2=s1-(y1(1,n)-beta)*(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;endendxvalue=a:h:b;yvalue=y1(1,:);ys=[xvalue',yvalue'];functionx=rk4(f,t0,x0,h,a,b)%rung-kuta法求每個點(diǎn)的近似值（參考大作業(yè)一）t=a:h:b;%迭代區(qū)間m=length(t);%區(qū)間長度t(1)=t0;x(:,1)=x0;%迭代初值fori=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.舉例求二階非線性方程的邊值問題：在matlab控制臺中輸入：f=@(x,y,z)(x^2+z*x^2);x0l=0;x0u=2*exp(-1);alfa=0;beta=2;h=0.01dbf(f,x0l,x0u,y0l,y0u,h,1e-