共回答了19个问题采纳率：89.5%

function dx=appollo(t,x)

mu=1/82.45;

mustar=1-mu;

r1=sqrt((x(1)+mu)^2+x(3)^2);

r2=sqrt((x(1)-mustar)^2+x(3)^2);

dx=[x(2)

2*x(4)+x(1)-mustar*(x(1)+mu)/r1^3-mu*(x(1)-mustar)/r2^3

x(4)

-2*x(2)+x(3)-mustar*x(3)/r1^3-mu*x(3)/r2^3];

--------------------------------------------------------------------------------------

x0=[1.2;0;0;-1.04935751];%x0(i)对应与xi的初值

options=odeset('reltol',1e-8);%该命令的另一种写法是options=odeset;options.reltol=1e-8;

tic

[t,y]=ode45(@appollo,[0,20],x0,options);%t是时间点,y的第i列对应xi的值,t和y的行数相同

toc

plot(y(:,1),y(:,3))%绘制x1和x3,也就是x和y的图形

title('Appollo卫星运动轨迹')

xlabel('X')

ylabel('Y')

1年前

10

!using["XSLSF"];                //使用命名空间XSLSF

//数组xArray存放x的值；ti为当前有效值的个数；tmax为ti对应的时间；tmin为起始时间。

xt(t:k:xArray,ti,tmax,tmin)=

{

k=(t-tmin)/(tmax-tmin)*ti-1,

if[k<0, k=0], if[k>=ti, k=ti-1],

xArray.getrai[k]

};

//数组yArray存放y的值；ti为当前有效值的个数；tmax为ti对应的时间；tmin为起始时间。

yt(t:k:yArray,ti,tmax,tmin)=

{

k=(t-tmin)/(tmax-tmin)*ti-1,

if[k<0, k=0], if[k>=ti, k=ti-1],

yArray.getrai[k]

};

//函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数f。

//t为自变量，x,y为函数值，dx,dy为右端函数值(即微分方程的值)。

f(t,x,y,dx,dy)=

{

dx=-(x^3)-y,

dy=x-y^3

};

//用连分式法对微分方程组进行积分，获得理论值。

//t1,t2为积分的起点和终点。

//h,s为自动变量。

//模块变量：hf为函数f的句柄，要预先获得该句柄；Array为工作数组；step为积分步长；eps为积分精度。

积分(t1,t2:h,s:hf,Array,step,eps)=

{

s=ceil[(t2-t1)/step],        //计算积分步数

h=(t2-t1)/s,                 //重新计算积分步长

{   pbs1[hf,t1,Array,h,eps], //对微分方程组积分一步

t1=t1+h                  //积分步长增加

}.until[abs(t1-t2)

};

数据(:i,h,t1,t2,x,y,static,free:hf,Array,step,eps,max,xArray,yArray,ti,tmax,tmin)= //微分方程组积分获得数据

{

if[free,delete(Array),delete(xArray),delete(yArray),return(0)],

hf=HFor("f"),                               //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄

step=0.01,eps=1e-6,                         //积分步长step要合适，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1

Array=new[rtoi(real_s),rtoi(2*15)],         //申请工作数组

max=1001,

xArray=new[rtoi(real_s),rtoi(max)],         //申请工作数组

yArray=new[rtoi(real_s),rtoi(max)],         //申请工作数组

Array.setra(0,1,0.5),                         //设置积分初值，通过模块变量Array传递，Array是一个数组

xArray.setra(0,1),                          //设置xArray的第一个值

yArray.setra(0,0.5),                          //设置yArray的第一个值

ti=1, h=step*3, tmin=0, tmax=0, t1=0, t2=h,

i=1,(i

积分(&t1,t2),                             //从t1积分到t2

Array.getra(0,&x,&y),                     //从数组Array获得t2时的积分值

xArray.setra(i,x), yArray.setra(i,y),     //将积分值保存到数组

ti=i+1, tmax=t1, t2=t2+h,

i++

}

};

//绘制函数图形

(::hxt,hyt)= hxt=HFor("xt"), hyt=HFor("yt");  //获得函数xt和yt的句柄，保存在模块变量中

gleDrawScene[HFor("Scene")],stop();      //设置场景绘制函数后退出

Scene(::hxt,hyt,tmax,tmin)=

{

glClear[],                           //清除屏幕以及深度缓存

glLoadIdentity[],                    //重置视图

glTranslated[0,0,-20],               //移动坐标，向屏幕里移动10个单元

glColor3d[0,0,1],                    //设置颜色

fgPlot[hxt,tmin,tmax,FG_Y,-1,2],   //绘制一元函数图象

glColor3d[1,0,0],                    //设置颜色

fgPlot[hyt,tmin,tmax,FG_Y,-1,2]    //绘制一元函数图象

};