飞行器建模与仿真.docx

飞行器建模与仿真作业班级：飞设12姓名：潘周周学号：2110702015作业一一、题目已知微分方程组用四阶龙哥库塔法求解该方程组。二、程序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;fori=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/w0^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度。已知:求轨迹基本参数。用龙格库塔法编写程序求解此题在此函数中要求考虑非球形假设以及国际标准大气。所以在编写此函数之前应先编写标准大气函数、阻力函数以及重力加速度函数。大气函数function Y=atmosdenty(h,vel,CL)%本函数用来求den,Kn,ma%输入h,vel,Cl分别表示高度，速度，特征长度R=287; %空气气体常数g0=9.806; %海平面重力加速度T0=288.15; %海平面温度re=6378140; %地球半径gama=1.405; %海平面比热比b=2/re;