打靶法(含Matlab程序)

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

2、速度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) 算出轨道就能计算减速力了。这时候你随便给个常数减速力到方程里飞船八成都能降落，但不是最优解。想想整个过程，开始降落之前飞船总机械能就那么多，你需要对飞

3、船做负功让机械能减到 0。题目里写发动机喷出翔的相对速度是一定的， 直觉告诉我飞船速度快的时候多喷一些速度慢的时候少喷一些，可以提高做负功的效率。但是多喷也不能超过上限 7500N ，所以这就是一个带约束优化问题， matlab 里边有专用的优化函数，用 fmincon- 2 -就好。找出最优解以后把过程画出来，看看F 可不可以是那 5 个状态量的线性组合，如果是的话就非常happy ，不是的话再说。三四阶段你可以扯点图像识别，什么二维复利叶分解找平坦区域，怎么一边下降一边根据自身状态调整路径之类的。五六阶段还真不知道说什么。一二阶段肯定是重点啦(3) 误差分析其实还挺难的。可能的误差来源是地

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

5、目的：进一步熟悉微分方程组边值问题数值算法（打靶法，有限差分法）。三、实验要求： 运用 Matlab/C/C+/Java/Maple/Mathematica等其中一种语言完成程序设计。四、实验原理 ：1打靶法：对于线性边值问题y p( x) yq( x) yf ( x)x a,b(1)y(a),y(b)假设 L 是一个微分算子使： Lyy p( x) yq(x) y则可得到两个微分方程：Ly1 f (x) ， y1 (a)， y1 (a) 0y1p( x) y1q(x) y1f (x) ， y1 (a)， y1 (a)0(2)Ly 20 ， y2 (a)0 ， y2 ( a) 1y2p( x)

6、 y2q( x) y20 ， y2 (a)0 ， y2 ( a)1(3)方程（ 2 ），（3）是两个二阶初值问题.假设 y1 是问题（ 2 ）- 4 -的解， y2 是问题（ 3 ）的解，且 y2 (b)0 ，则线性边值问题（ 1）的解为 : y( x) y1 ( x)y1 (b) y2 ( x) 。y2 (b)2 有限差分法：基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似； 把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限

7、差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。五、实验内容：% 线性打靶法functionk,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;- 5 -Y1(1,:)= CT1;CT2=0,1;Y2(1,:)= CT2;for

8、 k=1:nk1=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=

9、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- 6 -v=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,:

10、);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) 的曲线 ',' 是初值问题1 的数值解 u(x) 的曲线 ' , ' 是初值问题 2 的数值解 v(x) 的曲线 ')title

11、( '用线性打靶法求线性边值问题的数值解的图形')% 有限差分法functionk,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:nX=a:h:b;- 7 -k1(k)=feval(q1,X(k); A1(k+1,k)=1+h*k1(k)/2;k2(k)=feval(q2,X(k

12、);A2(k,k)=-2-(h.2)*k2(k);A3(k,k+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)

13、=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( '轴it x' ); ylabel( ' 轴 it y' ),legend('是边值问题的数值解y(x) 的曲线 ')- 8 -title( '用有限差分法求线性边值问题的数值解的图形'),p=k',X',y,wucha'打靶法3.Matlab源代码：创建 M文件：funct

14、ionys=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;% 选取适当的 s的初值x0=alfa,s0;- 9 -%迭代初值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)<=epsflag=1;end

15、endif flag=1while abs(y1(1,n)-beta)>eps s2=s1-(y1(1,n)-beta)*(s1-s0)/(y1(1,n)-y0(1,n);-10-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法求每个点的近似值（参考大作业一）t=a:h:b;% 迭代区间 m=length(t);% 区间长度t(1)=t0;-11-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. 举例求二阶非线性方程的边值问题：在matlab 控制台中输入：f=(x,y,z)(x2+z*x2);x0l=0;x0u=2*exp(