打靶法含Matlab程序

1、西京学数学软件实验任务书课程名称数学软件实验班级数0901学号0912020107姓名李亚强实验课题微分方程组边值问题数值算法（打靶法，有限差分法）实验目的熟悉微分方程组边值问题数值算法（打靶法，有限差分法）实验要求运用Matlab/C/C+/Java/Maple/Mathematica等其中一种语言完成实验内容微分方程组边值问题数值算法（打靶法，有限差分法）成绩教师动方向控制减速的推力，主要的控制量只有一个减速推力，减速还会消耗燃料让登月器的质量减小。所以在极坐标下系统的状态就是x=质量m，角度theta，高度r，角速度omega，径速度v这五个量，输入就是减速力F。先列微分方程，dx/dt

2、=f(x)+B*F，其中x是5*1的列向量，质量dm/dt=-F/2940，剩下几个翻下极坐标的手册。把这个动力学模型放到matlab里就能求解了，微分方程数值解用ode45。第一问F=0，让你求椭圆轨道非常容易。注意附件1里说15公里的时候速度是1.7km/s。算完以后验证一下对不对，对的话就是他了，不对的话说明这个椭圆轨道有进动，到时再说。(2) 算出轨道就能计算减速力了。这时候你随便给个常数减速力到方程里飞船八成都能降落，但不是最优解。想想整个过程，开始降落之前飞船总机械能就那么多，你需要对飞船做负功让机械能减到0。题目里写发动机喷出翔的相对速度是一定的，直觉告诉我飞船速度快的

3、时候多喷一些速度慢的时候少喷一些，可以提高做负功的效率。但是多喷也不能超过上限7500N，所以这就是一个带约束优化问题，matlab里边有专用的优化函数，用fmincon就好。找出最优解以后把过程画出来，看看F可不可以是那5个状态量的线性组合，如果是的话就非常happy，不是的话再说。三四阶段你可以扯点图像识别，什么二维复利叶分解找平坦区域，怎么一边下降一边根据自身状态调整路径之类的。五六阶段还真不知道说什么。一二阶段肯定是重点啦(3) 误差分析其实还挺难的。可能的误差来源是地球的引力，月亮绕地球向心加速度，太阳的引力(可能会很小)，对自身速度、角度的测量误差(比如你测出自身当前速度

4、100m/s但实际上是105m/s)，控制的时候F大小以及角度的误差(比如你想朝正前方向喷2000N但实际上偏了2度而且F=2010N之类)。上一问已经求出了最优控制策略和飞船路线，把这些扰动加进去以后算出新的路线减掉理想路线求偏差，然后随便用个卡尔曼滤波器把误差给校正All for Joy2014/9/13 11:14:38老师的思路，求大神解答给我一份呀实验二十七实验报告1、 实验名称：微分方程组边值问题数值算法（打靶法，有限差分法）。2、 实验目的：进一步熟悉微分方程组边值问题数值算法（打靶法，有限差分法）。3、 实验要求：运用Matlab/C/C+/Jav

5、a/Maple/Mathematica等其中一种语言完成程序设计。4、 实验原理：1 打靶法：对于线性边值问题 (1)假设是一个微分算子使：则可得到两个微分方程： ， (2) ， (3)方程（2），（3）是两个二阶初值问题.假设是问题（2）的解，是问题（3）的解，且，则线性边值问题（1）的解为: 。2 有限差分法：基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组 ，

6、 解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。5、 实验内容：%线性打靶法function k,X,Y,wucha,P=xxdb(dydx1,dydx2,a,b,alpha,beta,h)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

7、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;endu=Y1(:,1)for k=1:nk1=feval(dydx2,X(k),Y2(k,:) x2=X(k)+h/2;y2=Y

8、2(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;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

9、=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)的曲线','是初值问题1的数值解u(x)的曲线', '是初值问题2的数值解v(x)的曲线')title('用线性打靶法求线性边值问题的数值解的图形

10、9;)%有限差分法function k,A,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:n X=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-(h.2)*k2(k); A3(k,k+

11、1)= 1-h*k1(k)/2; k3(k)=feval(q3,X(k);endfor k=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=AB1;Y1=Y' y=alpha;Y;beta;for k=2:n+1wucha(k)=norm(y(k)-y(k-1); k=k+1;endX=X

12、(1:n+1); y=y(1:n+1,1); k=1:n+1; wucha=wucha(1:k,:); plot(X,y(:,1),'mp')xlabel('轴it x'); ylabel('轴it y'),legend('是边值问题的数值解y(x)的曲线')title('用有限差分法求线性边值问题的数值解的图形'),p=k',X',y,wucha'打靶法3.Matlab 源代码：创建M 文件：function ys=dbf(f,a,b,alfa,beta,h,eps) ff=(x,y)y(

13、2),f(y(1),y(2),x); xvalue=a:h:b;%x取值范围n=length(xvalue) s0=a-0.01;%选取适当的s的初值x0=alfa,s0;%迭代初值flag=0;%用于判断精度y0=rk4(ff,a,x0,h,a,b); if abs(y0(1,n)-beta)<=eps flag=1; y1=y0; else s1=s0+1; x0=alfa,s1; y1=rk4(ff,a,x0,h,a,b); if abs(y1(1,n)-beta)<=eps flag=1; end end if flag=1 while abs(y1(1,n)-beta)&

14、gt;eps s2=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; end end xvalue=a:h:b; yvalue=y1(1,:); ys=xvalue',yvalue' function x=rk4(f,t0,x0,h,a,b)%rung-kuta法求每个点的近似值（参考大作业一）t=a:h:b;%迭代区间m=length(t);%区间长度t(1)=t0; x(:,1)=x0;%迭代初值for i=1:m-1 L1=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); end 4.举例求二阶非线性方程的边值问题：在matlab 控制台中输入：f=(x,y,z)(x2+z*x2); x0l=0; x0u=2*exp(-1); alfa=0; beta=2;