飞行器建模与仿真

飞行器建模与仿真作业

班级：飞设12

姓名：潘周周

学号：2110702015

作业一

一、题目

已知微分方程组

x

x 2

1=∙

x x 32=∙ x

x 4

3=∙

x

x 5

4=∙

x x 6

5

=∙

⎥⎥⎦

⎤⎢⎢⎣⎡------=∙

w x a w x a w x a w x a w x a x x a w x 5

6

54054304320320

211660

6 用四阶龙哥库塔法求解该方程组。 二、程序

1.编写rk4f.m 函数

function r=rk4f(x1,x2,x3,x4,x5,x6,xs,ts,tf,N) %四阶龙格库塔法计算微分方程组程序 %x1,x2,x3,x4,x5,x6分别是初值 %xs 是方程组中的x

%ts ，tf 代表代表左端点和右端点 %N 是迭代步数

w0=50; a1=3.86; a2=7.46; a3=9.13; a4=7.46; a5=3.86; a6=1;

aa=zeros(1,N+1); ab=zeros(1,N+1); ac=zeros(1,N+1); ad=zeros(1,N+1); ae=zeros(1,N+1); af=zeros(1,N+1);

aa(1)=x1;

ab(1)=x2;

ac(1)=x3;

xd(1)=x4;

xe(1)=x5;

af(1)=x6;

h=(tf-ts)/N;

t=ts:h:tf;

for i=1:1:N;

%ki1(i=1,2,...6)

k11=x2;

k21=x3;

k31=x4;

k41=x5;

k51=x6;

k61=w0^6/a6*(xs-x1-a1*x2/w0-a2*x3/w0^2-a3*x4/w0^3-a4*x5/w0^4-a5*x6/w0 ^5);

%ki2(i=1,2,...6)

k12=x2+h/2*k21;

k22=x3+h/2*k31;

k32=x4+h/2*k41;

k42=x5+h/2*k51;

k52=x6+h/2*k61;

k62=w0^6/a6*(xs-(x1+h/2*k11)-a1*k12/w0-a2*k22/w0^2-a3*k32/w0^3-a4*k42 /w0^4-a5*k52/w0^5);

%%ki3(i=1,2,...6)

k13=x2+h/2*k22;

k23=x3+h/2*k32;

k33=x4+h/2*k42;

k43=x5+h/2*k52;

k53=x6+h/2*k62;

k63=w0^6/a6*(xs-(x1+h/2*k12)-a1*k13/w0-a2*k23/w0^2-a3*k33/w0^3-a4*k43 /w0^4-a5*k53/w0^5);

%ki4(i=1,2,...6)

k14=x2+h*k23;

k24=x3+h*k33;

k34=x4+h*k43;

k44=x5+h*k53;

k54=x6+h*k63;

k64=w0^6/a6*(xs-(x1+h*k13)-a1*k14/w0-a2*k24/w0^2-a3*k34/w0^3-a4*k44/w 0^4-a5*k54/w0^5);

%计算下一步的值

x1=x1+h*(k11+2*k12+2*k13+k14)/6;

x2=x2+h*(k21+2*k22+2*k23+k24)/6;

x3=x3+h*(k31+2*k32+2*k33+k34)/6;

x4=x4+h*(k41+2*k42+2*k43+k44)/6;

x5=x5+h*(k51+2*k52+2*k53+k54)/6;

x6=x6+h*(k61+2*k62+2*k63+k64)/6;

%形成向量

aa(i+1)=x1;

ab(i+1)=x2;

ac(i+1)=x3;

xd(i+1)=x4;

xe(i+1)=x5;

af(i+1)=x6;

end

%画出图型

subplot(2,3,1)

plot(t,aa,'-')

subplot(2,3,2)

plot(t,ab,'-')

subplot(2,3,3)

plot(t,ac,'-')

subplot(2,3,4)

plot(t,ad,'-')

subplot(2,3,5)

plot(t,ae,'r')

subplot(2,3,6)

plot(t,af,'b')

2．编写主函数调用该函数

rk4f(0,0,0,0,0,0,1,0,0.5,100)

三、运行结果

作业2

一、 题目

弹道式再入轨迹仿真。已知太空舱质量为350Kg ，从近地点高度为200Km 的轨道返回。偏心率为0.2，轨道倾角为80度，近地点角距为265度，升焦点赤经100度。已知:

()8489182889.79)0(;10)0(;89967.65790-=︒-==δλkm r

︒===955734.99)0(;54681217.0)0(;/30433867.7589)0(A s m v φ

求轨迹基本参数 )(),(),(),(),(),(t A t t v t t t r φδλ。 二、用龙格库塔法编写程序求解此题

在此函数中要求考虑非球形假设以及国际标准大气。所以在编写此函数之前应先编写标准大气函数、阻力函数以及重力加速度函数。 （1）大气函数