打靶法(含Matlab程序)

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

2、，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就好。找出最优解以后把过程画出来，看看 F可不可以是那5个 状态量的线性组合，如果是的话就非常 happy ，不是的话再说。 三四阶段你可以扯点图像识别，什么二维复利叶分解找平坦区域， 怎么一边下降一边根据自身状态调整路径之类的。 五六阶段还真不知道说什么。一二阶段肯定是重点啦（3） 误差分析其实还挺难的。可能的误差来源是地球的引力，月 亮绕

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

5、一步熟悉微分方程组边值问题数值算法 (打靶法， 有限差分法)。三、实验要求： 运用 Matlab/C/C+/Java/Maple/Mathematica 等其中一种语言完成程序设计。四、实验原理：1打靶法：对于线性边值问题y p(x)y q(x)y f(x) x a,b (1)y(a) , y(b) (1)假设 L 是一个微分算子使： Ly y p(x)y q(x)y 则可得到两个微分方程：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)是两个二阶初值问题 假设yi是问题(2)的解，y是问题(3)的解，且y2(b) 0,则线性边值问题(1) 的解为：y(x) %(x) 警 y2(x)。y2(b)2 .有限差分法：基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解 区域上的连续变量的函数用在网格上定义的离散变量函数 来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地 代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插

7、值方 法便可以从离散解得到定解问题在整个区域上的近似解。五、实验内容：%线性打靶法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(' 是边值问题的数值解 y(x) 的曲线 ',' 是初值问题 1 的数值解 u(x) 的曲线','是初值问题2的数值解v(x)的曲线')title(' 用线性打靶法求线性边值问题的数值解的图形 ')%有限差分法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('是边值问题的数值解 title（' 用有限差分法求线性边值问题的数值解的图形 '）, p=k',X',y,wucha'打靶法3.Matlab源代码：创建 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;% 选取适当的 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 法求每个点的近似值(参考大作业一) t=a:h:b;% 迭代区间 m=length(t); %区间长度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.举例求二阶非线性方程的边值问题： 在 matlab 控制台中输入： f=(x,y,z)(xA2+z*xA2);x0l=0;x0u=2*exp(-1);alfa=0;be