打靶法(含Matlab程序).docVIP

西京学数学软件实验任务书 课程名称 数学软件实验 班级 数0901 学号 0912020107 姓名 李亚强 实验课题 微分方程组边值问题数值算法（打靶法，有限差分法） 实验目的 熟悉微分方程组边值问题数值算法（打靶法，有限差分法） 实验要求 运用Matlab/C/C++/Java/Maple/Mathematica等其中一种语言完成 实验内容 微分方程组边值问题数值算法（打靶法，有限差分法） 成绩 教师  动方向控制减速的推力，主要的控制量只有一个减速推力，减速还会消耗燃料让登月器的质量减小。所以在极坐标下系统的状态就是x‘=[质量m，角度theta，高度r，角速度omega，径速度v]这五个量，输入就是减速力F。先列微分方程，dx/dt=f(x)+B*F，其中x是5*1的列向量，质量dm/dt=-F/2940，剩下几个翻下极坐标的手册。把这个动力学模型放到matlab里就能求解了，微分方程数值解用ode45。第一问F=0，让你求椭圆轨道非常容易。注意附件1里说15公里的时候速度是1.7km/s。算完以后验证一下对不对，对的话就是他了，不对的话说明这个椭圆轨道有进动，到时再说。(2)?算出轨道就能计算减速力了。这时候你随便给个常数减速力到方程里飞船八成都能降落，但不是最优解。想想整个过程，开始降落之前飞船总机械能就那么多，你需要对飞船做负功让机械能减到0。题目里写发动机喷出翔的相对速度是一定的，直觉告诉我飞船速度快的时候多喷一些速度慢的时候少喷一些，可以提高做负功的效率。但是多喷也不能超过上限7500N，所以这就是一个带约束优化问题，matlab里边有专用的优化函数，用fmincon就好。找出最优解以后把过程画出来，看看F可不可以是那5个状态量的线性组合，如果是的话就非常happy，不是的话再说。三四阶段你可以扯点图像识别，什么二维复利叶分解找平坦区域，怎么一边下降一边根据自身状态调整路径之类的。五六阶段还真不知道说什么。一二阶段肯定是重点啦(3)?误差分析其实还挺难的。可能的误差来源是地球的引力，月亮绕地球向心加速度，太阳的引力(可能会很小)，对自身速度、角度的测量误差(比如你测出自身当前速度100m/s但实际上是105m/s)，控制的时候F大小以及角度的误差(比如你想朝正前方向喷2000N但实际上偏了2度而且F=2010N之类)。上一问已经求出了最优控制策略和飞船路线，把这些扰动加进去以后算出新的路线减掉理想路线求偏差，然后随便用个卡尔曼滤波器把误差给校正All?for?Joy2014/9/13?11:14:38老师的思路，求大神解答给我一份呀 n=fix((b-a)/h); X=zeros(n+1,1); CT1=[alpha,0]; Y=zeros(n+1,length(CT1)); Y1=zeros(n+1,length(CT1)); Y2=zeros(n+1,length(CT1)); X=a:h:b; Y1(1,:)= CT1; CT2=[0,1];Y2(1,:)= CT2; for k=1:n k1=feval(dydx1,X(k),Y1(k,:)) 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; end u=Y1(:,1) for k=1:n k1=feval(dydx2,X(k),Y2(k,:)) x2=X(k)+h/2;y2=Y2(k,:)+k1*h/2; k2=feval(dydx2,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; end v=Y2(:,1) Y=u+(beta-u(n+1))*v/v(n+1) for k=2:n+1 wucha(k)=norm(Y(k)-Y(k-1)); k=k+1; end X=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:); P=[k,X,Y,wucha]; plot(X,Y(:,1),ro,X,Y1(:,1),g*,X,Y2(:,1),mp) xlabel(轴\it x); ylabel(轴\it y) legend(是边值问题的数值解y(x)的曲线,是