打靶法(含matlab程序)1

《打靶法(含matlab程序)1》由会员分享，可在线阅读，更多相关《打靶法(含matlab程序)1（10页珍藏版）》请在金锄头文库上搜索。

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

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

5、组边值问题数值算法（打靶法，有限差分法） 。二、实验目的实验目的：进一步熟悉微分方程组边值问题数值算法（打靶法，有限差分法） 。三、实验要求：实验要求：运用 Matlab/C/C+/Java/Maple/Mathematica 等其中一种语言完成程序设计。四、实验原理实验原理：1打靶法：对于线性边值问题 )(,)(,)()()(byaybaxxfyxqyxpy(1)假设L是一个微分算子使：( )( )Lyyp x yq x y则可得到两个微分方程：)(1xfLy ，)(1ay，0)(1 ay)()()(111xfyxqyxpy ，)(1ay，0)(1 ay(2)02Ly，0)(2ay，1)(2

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

7、方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。五、五、实验内容：实验内容：%线性打靶法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; fo

8、r k=1:n k1=feval(dydx1,X(k),Y1(k,:) x2=X(k)+h/2;y2=Y1(k,:)+k1*h/2;- 4 -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=

9、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; plo

10、t(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,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

11、);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);- 5 -A2(k,k)=-2-(h.2)*k2(k); A3(k,k+1)= 1-h*k1(k)/2; k3(k)=feval(q3,X(k); end for k=2:n B(k,1)=(h.2)*k3(k); end B(1,1)=(h.2)*k3(1)-(1+h*k1(1)/2)*alpha; B(n-1,1)=(h.2)*k3(n-1

12、)-(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+1 wucha(k)=norm(y(k)-y(k-1); k=k+1; end X=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)的曲线)title(

13、用有限差分法求线性边值问题的数值解的图形), p=k,X,y,wucha;打靶法3.Matlab3.Matlab源代码：源代码：创建创建M M 文件：文件：functionfunctionys=dbf(f,a,b,alfa,beta,h,eps)ys=dbf(f,a,b,alfa,beta,h,eps)ff=(x,y)y(2),f(y(1),y(2),x);ff=(x,y)y(2),f(y(1),y(2),x);- 6 -xvalue=a:h:b;xvalue=a:h:b;%x%x取值范围取值范围n=length(xvalue)n=length(xvalue)s0=a-0.01;s0=a-0.

14、01;% %选取适当的选取适当的s s的初值的初值x0=alfa,s0;x0=alfa,s0;% %迭代初值迭代初值flag=0;flag=0;% %用于判断精度用于判断精度y0=rk4(ff,a,x0,h,a,b);y0=rk4(ff,a,x0,h,a,b);ifif abs(y0(1,n)-beta)epsabs(y1(1,n)-beta)epss2=s1-(y1(1,n)-beta)*(s1-s0)/(y1(1,n)-y0(1,n);s2=s1-(y1(1,n)-beta)*(s1-s0)/(y1(1,n)-y0(1,n);x0=alfa,s2;x0=alfa,s2;y2=rk4(ff,

15、a,x0,h,a,b);y2=rk4(ff,a,x0,h,a,b);s0=s1;s0=s1;s1=s2;s1=s2;y0=y1;y0=y1;y1=y2;y1=y2;endendendendxvalue=a:h:b;xvalue=a:h:b;yvalue=y1(1,:);yvalue=y1(1,:);ys=xvalue,yvalue;ys=xvalue,yvalue;functionfunctionx=rk4(f,t0,x0,h,a,b)x=rk4(f,t0,x0,h,a,b)- 8 -%rung-kuta%rung-kuta法求每个点的近似值（参考大作业一）法求每个点的近似值（参考大作业一）t

16、=a:h:b;t=a:h:b;% %迭代区间迭代区间m=length(t);m=length(t);% %区间长度区间长度t(1)=t0;t(1)=t0;x(:,1)=x0;x(:,1)=x0;% %迭代初值迭代初值forfor i=1:m-1i=1:m-1L1=f(t(i),x(:,i);L1=f(t(i),x(:,i);L2=f(t(i)+h/2,x(:,i)+(h/2)*L1);L2=f(t(i)+h/2,x(:,i)+(h/2)*L1);L3=f(t(i)+h/2,x(:,i)+(h/2)*L2);L3=f(t(i)+h/2,x(:,i)+(h/2)*L2);L4=f(t(i)+h,x(:,i)+h*L3);L4=f(