我正在写一个脚本

ode45

为了整合

卫星在火星附近双曲线轨道上的运动方程

.

我需要整合地球上的整个通道:从SOI半径开始(

576000km

)向行星前进,然后穿过大气层直到卫星到达

opposite

“大气边界(设置在

250km

从表面)。

当它接收到输入A时

tspan

比大约高

200000

秒(我需要大约

400000

秒),Matlab给出了以下信息:

'无法执行分配,因为左侧的大小为

4乘2,右边的尺寸是4乘5。

error happens in line 488 of Ode45

.

我搜索了一些类似的案例,但找不到任何东西,也不知道如何使用

条件断点

在一个复杂的函数上计算出一些东西,如ode45。我也尝试了不同的选择

odeset

“但没什么变化。我不知道错误可能在哪里。

这是脚本,其中我使用两个附加函数来获取一些参数:

Vinf=[2.7 4 6 8];

mu_m=42828;

R_msoi=.576e6;

[afas,efas,pfas]=parasint(Vinf);

an_vera0=zeros(1,length(Vinf));

Vr0=zeros(1,length(Vinf));

Vt0=zeros(1,length(Vinf));

for j=1:length(Vinf)

c_tstar0=(1/efas(j))*((pfas(j)/R_msoi)-1);

an_vera0(j)=-acos(c_tstar0);

s_tstar0=sin(an_vera0(j));

Vr0(j)=sqrt(mu_m/pfas(j))*efas(j)*s_tstar0;

Vt0(j)=sqrt(mu_m/pfas(j))*(1+efas(j)*c_tstar0);

end

ti=time(an_vera0,efas,afas,mu_m);

for n=1:length(ti)

I=(0:3600:ti(n));

options=odeset('Vectorized','on','RelTol',1e-8,'AbsTol',1e-9);

[t,Y]=ode45(@(t,y) eqMotoCur(t,y),I,y0,options);

end

这是“

parasint

“功能:

function [afas,efas,pfas]=parasint(Vinf)

mu_m=42828;

R_m=3396.2;

Rp0=R_m+400;

Rpf=R_m+50;

a0=zeros(1,length(Vinf));

e0=zeros(1,length(Vinf));

p0=zeros(1,length(Vinf));

afas=zeros(1,length(Vinf));

efas=zeros(1,length(Vinf));

pfas=zeros(1,length(Vinf));

for j=1:length(Vinf)

a0(j)=(-mu_m)/(Vinf(j)^2);

e0(j)=1-Rp0/a0(j);

p0(j)=a0(j)*(1-e0(j)^2);

c(1)=p0(j);

c(2)=Rpf*(1-e0(j)^2);

c(3)=Rpf*(1-e0(j)^2)-p0(j);

Efas=roots(c);

ind=(Efas>1 & isreal(Efas));

efas(j)=Efas(ind);

afas(j)=Rpf/(1-efas(j));

pfas(j)=afas(j)*(1-efas(j)^2);

end

这是“

time

“功能:

function [ti] = time(an_vera0,efas,afas,mu_m)

ti=zeros(1,length(an_vera0));

for j=1:length(an_vera0)

F=2*atanh(sqrt((efas(j)-1)/(efas(j)+1))*tan(an_vera0(j)/2));

T=sqrt((-afas(j))^3/mu_m)*(efas(j)*sinh(F) - F);

ti(j)=-2*T;

end

这是ODE45的输入功能:

function [dydt] = eqMotoCur(t,y)

R_m=3396.2;

mu_m=42828;

Cd=2.2;

S=7.065e-6;

m=3699.046;

wE=0.24117/R_m;

if abs(y(1))>R_m+250 && y(3)>0

dydt(1:4)=0;

else

dydt=zeros(4,1);

dydt(1)=y(3);

dydt(2)=y(4)/y(1);

dydt(3)=(-mu_m/(y(1)^2))+((y(4)^2)/y(1))-(.5*Cd)*(S/m)*...

(dens(y(1)-R_m))*(sqrt((y(3)^2)+(y(4)-(wE*y(1)))^2))*y(3);

dydt(4)=-(y(4)*(y(3)/y(1)))-(.5*Cd)*(S/m)*...

(dens(y(1)-R_m))*(sqrt((y(3)^2)+(y(4)-(wE*y(1)))^2))*...

(y(4)-(wE*y(1)));

end

end

这是“

dens

“通过大气插入数据的功能:

function [rho] = dens(z)

load Density.mat h d

d=d*10^9;

if any(h)==abs(z)

j=h==abs(z);

rho=d(j);

elseif abs(z)<250

c_tot=(find(h>=abs(z)));

c=c_tot(1);

H=-(h(c)-h(c-1))/(log(d(c)/d(c-1)));

rho=d(c-1)*exp(-(abs(z)-h(c-1))/H);

else

rho=0;

end

end

这是Matlab给出的错误:

无法执行分配,因为左侧的大小为

4×2,右侧尺寸为4×5。

中的错误

ode45 (line 488)

yout(:,idx) = yout_new

;

中的错误

MainCur (line 59)

[t,Y]=ode45(@(t,y) eqMotoCur(t,y),I,y0,options)

;

提前谢谢。