打靶法含Matlab程序Word下载.docx
《打靶法含Matlab程序Word下载.docx》由会员分享，可在线阅读，更多相关《打靶法含Matlab程序Word下载.docx（8页珍藏版）》请在冰豆网上搜索。
这时候你随便给个常数减速力到方程里飞船八成都能降落，但不是最优解。
想想整个过程,开始降落之前飞船总机械能就那么多,你需要对飞船做负功让机械能减到0。
题目里写发动机喷出翔的相对速度是一定的，直觉告诉我飞船速度快的时候多喷一些速度慢的时候少喷一些,可以提高做负功的效率。
但是多喷也不能超过上限7500N,所以这就是一个带约束优化问题，matlab里边有专用的优化函数,用ｆmincｏn就好。
找出最优解以后把过程画出来，看看F可不可以是那5个状态量的线性组合，如果是的话就非常haｐpｙ，不是的话再说。
ﻫ三四阶段你可以扯点图像识别,什么二维复利叶分解找平坦区域，怎么一边下降一边根据自身状态调整路径之类的。
五六阶段还真不知道说什么。
一二阶段肯定是重点啦ﻫ（3）
误差分析其实还挺难的。
可能的误差来源是地球的引力，月亮绕地球向心加速度,太阳的引力（可能会很小）,对自身速度、角度的测量误差（比如你测出自身当前速度１00m/ｓ但实际上是１０5ｍ/s）,控制的时候F大小以及角度的误差（比如你想朝正前方向喷200０Ｎ但实际上偏了2度而且F=201０N之类）。
上一问已经求出了最优控制策略和飞船路线，把这些扰动加进去以后算出新的路线减掉理想路线求偏差,然后随便用个卡尔曼滤波器把误差给校正ﻫＡll
for
Joy
2０１４/9/13
１1:
14:
３8ﻫﻫ老师的思路,求大神解答给我一份呀
实验二十七实验报告
1、实验名称:
微分方程组边值问题数值算法（打靶法,有限差分法）。
2、实验目的：
进一步熟悉微分方程组边值问题数值算法（打靶法，有限差分法）。
3、实验要求：
运用Matlaｂ/Ｃ／Ｃ++／Java／Maplｅ／Ｍatｈematicａ等其中一种语言完成程序设计。
4、实验原理:
1．打靶法：
对于线性边值问题
（１）
假设是一个微分算子使:
则可得到两个微分方程：
，,
，， （２）
,
，,（3）
方程（２）,（3）是两个二阶初值问题.假设是问题
（2）的解,是问题（3）的解,且，则线性边值问题（１）的解为:
。
2．有限差分法:
基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点；
把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；
把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解。
然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。
5、实验内容：
%线性打靶法
ｆuncｔｉon[k,X,Y,wucha,P]＝xｘdb（dydx1,ｄｙｄx2,a，b,alｐha，beta，h）
n=ｆｉx（（b-a）／ｈ）;
Ｘ＝zeros（n+1，1）;
CT１=［alpｈa,0];
Y=zｅros（n＋1，length（CT1））;
Ｙ1=zeros（n＋1，lenｇtｈ（CＴ１））；
Y2=zeｒｏs（ｎ＋１,lenｇtｈ（CT1））；
Ｘ=a:
ｈ：
b；
Y1（１，：
）＝ＣT1;
CT2=[0,1];
Y2（１,:
）=ＣT２;
fｏrk=１:
n
k1=fｅvaｌ（dydx1，Ｘ（k）,Y1（k,：
））
x2=X（k）+h/2;
y2=Y1（k，:
）＇+k1*h／2;
k2＝feval（dｙdx1，x2,y2）;
k3=fｅｖal（ｄydｘ1,x2,Ｙ1（ｋ,:
）'
+k2*h/２）;
ｋ4=fevａｌ（dydx１,X（k）＋h,Ｙ1（k,:
+k3＊ｈ）；
Y1（ｋ+1,：
）=Ｙ1（k,:
）+h*（k1'
+2＊ｋ2＇＋2*ｋ３'
+k4'
）/6,ｋ=k+１;
end
ｕ=Y1（：
1）
ｆork=１:
k1=ｆｅvａl（ｄydｘ2,X（k）,Ｙ2（k，:
ｘ２=X（ｋ）+ｈ/2;
y2=Ｙ２（ｋ，:
+k１＊ｈ/２;
k2=feｖal（dydx２，x２，y2）；
k3=feｖal（dｙdx２,x2,Y２（k,:
+ｋ２＊h/２）;
k4＝ｆeval（ｄydx2, X（ｋ）+ｈ,Y2（k,：
+ｋ3*ｈ）;
Ｙ2（k+1,:
）=Y2（k，：
）＋h＊（k１'
＋２*k２'
+2*k３'
+k４＇）/６,k=k+1;
eｎd
v=Y２（：
Y=u+（beta-ｕ（n＋1））＊v／v（n＋1）
ｆork=2:
n+１
wucha（k）＝norm（Ｙ（k）-Ｙ（k-1））;
ｋ=k+１；
enｄ
X=X（1：
n+１）；
Y=Y（1:
n＋1,:
）;
ｋ=1:
ｎ+1;
ｗucha=wucha（１：
k,:
）；
P=[k＇,X＇,Ｙ，ｗucｈａ＇];
ploｔ（Ｘ,Y（：
1）,'
ro'
X,Y１（：
1），＇g*'
，Ｘ,Y2（:
mｐ'
）
xlａbｅｌ（'
轴\ｉtｘ'
ylabel（'
轴＼ity'
legend（'
是边值问题的数值解y（x）的曲线'
＇是初值问题1的数值解u（x）的曲线'
'
是初值问题2的数值解v（x）的曲线'
ｔitle（'
用线性打靶法求线性边值问题的数值解的图形'
%有限差分法
ｆunction[ｋ,Ａ,B１,X，Y,y，ｗucｈa,p]=yxcf（q1,ｑ2,ｑ3,a，b,ａlｐha,beta,h）
n＝fｉx（（b－a）/h）;
X＝ｚeｒｏs（n+1,1）;
Y=zerｏｓ（n+1,１）;
A1=ｚeｒos（ｎ,n）;
A2=ｚeｒoｓ（n,n）;
A3＝ｚeros（n,n）;
Ａ=zｅros（n,n）;
B＝ ｚeｒｏs（n，1）;
fork=１:
ｎ
ｈ:
b;
k１（ｋ）=fｅval（q1,X（k））；
A１（k+1,k）＝1+ｈ*ｋ1（ｋ）/２;
k2（k）=fｅvａｌ（q2,X（k））；
A２（k,k）=－2-（h.^2）*k2（ｋ）;
A３（k,k+1）=１-ｈ＊k1（k）/2;
ｋ3（k）=fevaｌ（ｑ３，Ｘ（k））；
for ｋ＝２:
Ｂ（k,1）=（ｈ.^２）*k3（k）;
B（１，1）＝（ｈ.^２）＊k3
（1）－（１＋ｈ＊k１
（1）/2）*aｌpha；
B（n-1,1）=（ｈ.^2）*k3（n-1）-（1+h＊k1（ｎ-1）/2）*bｅta;
A=Ａ1（１:
n－1,１:
ｎ-1）+A２（1：
n-１，１:
n-1）+A3（１:
n-1,１：
n-1）；
B1＝B（1：
ｎ－1,１）;
Y=A\B1;
Y1=Y'
;
y=[aｌpha;
Y;
bｅta］;
n＋1
wｕcha（k）=nｏrｍ（ｙ（k）-ｙ（k-1））;
k=k＋1;
Ｘ=X（1:
ｎ+１）;
ｙ=y（1：
n＋1,1）;
k=１:
ｎ+1；
wucha＝ｗucha（1:
plot（X，y（:
1），＇mｐ'
xlaｂel（'
轴\itx'
ｙlabel（'
轴\ｉt y'
）,lｅgend（'
是边值问题的数值解ｙ（x）的曲线'
用有限差分法求线性边值问题的数值解的图形'
）,
p=[k'
，Ｘ'
y，ｗｕcha'
]；
打靶法
３.Ｍatlａb
源代码：
创建M文件：
ｆｕｎction
ys=dbf（f,ａ,b，alfa,ｂeta,h，eｐｓ）
fｆ=＠（x,y）[y（２），f（y
（1）,y
（2）,x）];
xvalue=ａ:
h:
%ｘ取值范围
n=lｅngth（xｖalｕe）
s0=a-0.01;
％选取适当的ｓ的初值
x0=［ａｌfa,s０］;
%
迭代初值
flａg=0；
%用于判断精度
y0＝rk4（fｆ，a,ｘ０,h,a,b）;
ifabs（y０（1,n）-ｂeta）<
＝eps
flａg=1；
y1=y0;
else
s1=s0+１;
ｘ0=［ａlfa,ｓ１]；
y1=rk4（ff，ａ,x０，h,a，b）;
ifａｂs（ｙ1（1,ｎ）-beta）<
=eps
flag＝1;
enｄ
end
iffｌag～＝１
whiｌe abs（y1（1,n）－bｅｔa）>
ｅps
s2=s1-（ｙ1（1,n）－beta）*（ｓ1－s0）/（y1（1,n）－y0（1,n））;
x0=[aｌfa,ｓ2];
y2=rk4（fｆ,ａ,x０，ｈ，a,ｂ）;
s0=s1;
ｓ１＝s２;
y０＝y１;
y1=y2;
xvalｕe=a:
h：
ｙvａlue=y１（１,:
ys=［xvａｌue'
yvaluｅ'
];
fｕｎcｔion
x=rk4（f,t０，ｘ0,h，a,b）
%rｕng－ｋutａ
法求每个点的近似值（参考大作业一）
t＝a:
%迭代区间ｍ=ｌenｇｔｈ（t）;
%区间长度
t（１）=t0;
x（:
１）=ｘ0；
%迭代初值
ｆoｒi=1：
ｍ-1
L1＝f（ｔ（ｉ）,x（：
ｉ））；
L2=ｆ（t（ｉ）+ｈ/2，ｘ（：
ｉ）＇＋（h/２）＊L1）；
L３＝f（ｔ（ｉ）＋h／2,x（:
ｉ）'
+（h／２）＊L2）;
L4＝f（t（i）+h，x（:
i）'
＋h*L３）;
x（：
，i＋１）=ｘ（：
+（h/6）*（L１+２＊L2＋2*L3+L4）；
4.举例求二阶非线性方程的边值问题:
在ｍatlab控制台中输入:
f=@（x,ｙ,z）（ｘ^２+ｚ*x＾２）;
x0l＝0;
x０u=2＊exp（-１）;
alfa=0;
beｔa=2;
h=0.0１
ｄbｆ（f,ｘ0l,x0ｕ，ｙ0ｌ，ｙ0ｕ,ｈ,1e-6）;
>
y=anｓ（：
，2）;
x=ans（:
1）;
plot（x,ｙ,'
-r'
）
＞
结果:
再输入：
m=0:
０.０1:
2;
＞>
n＝m.*eｘp（－1/2*m）;
pｌot（n,m）
plｏｔ（ｘ,ｙ,'
-ｒ'
n,ｍ,'
－b'
升级会员 