Как исправить ошибку присваивания, используя ode45 в MatLab (строка 488 функции ode45) - PullRequest
Новая
       51

Как исправить ошибку присваивания, используя ode45 в MatLab (строка 488 функции ode45)

1 голос
/ 20 мая 2019

Я пытаюсь написать скрипт, который использует ode45, чтобы интегрировать уравнения движения спутника по гиперболической траектории около Марса .

Мне нужно объединить весь проход вокруг планеты: начиная с радиуса SOI (576000km) и далее по направлению к планете, затем через атмосферу, пока спутник не достигнет границы атмосферы "opposite" (установлен на 250km с поверхности).

Когда он получает на входе 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

Это ошибка, которую мне дает Матлаб:

Невозможно выполнить назначение, потому что размер левой стороны 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);

Заранее спасибо.

Ответы [ 2 ]

1 голос
/ 20 мая 2019
  • ode45(@(t,y) eqMotoCur(t,y),I,y0,options) возвращает матрицу 5 by 4
  • [t,Y] - матрица 2 by 4
  • Таким образом, несоответствие размеров как 5 by 4 # 2 by 4
  • Заменить [t,Y] на [t,Y, ~,~, ~]
  • Оставьте только первые две нужные вам переменные и игнорируйте оставшиеся используя ~
[t,Y, ~,~, ~]=ode45(@(t,y) eqMotoCur(t,y),I,y0,options);

Убедитесь, что вы определили y0 сначала перед запуском ode45

функция dens(z) должна быть определена внутри eqMotoCur(t,y), поскольку вам нужно выполнить некоторые вычисления внутри eqMotoCur(t,y)

function [dydt] = eqMotoCur(t,y)
    %% define dens here 
    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
    %% Now start eqMotoCur
    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) ||y(1)==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

Основной проблемой является время t, довольно большое, если оно оценивается в секунда и вот где ode45 застревает.
Альтернативное преобразование секунд в часов (делим на 3600) , это будет работать. Также результат полученный тип данных структуры, где x обозначает t, а y обозначает y. Давайте скажем t = sol{1}.x и y = sol{1}.y
, так как массив будет изменяя каждую итерацию, я использую cell, чтобы хранить их отдельно

Основной код следующий: 

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);
ti = ti/3600;
sol = cell(1,4);
for n=1:length(ti)
   y0=[.576e6; an_vera0(n); Vr0(n) ;Vt0(n)];


    I=0:1:ti(n);

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

sol{n} = ode45(@eqMotoCur,I, y0,options);

end
0 голосов
/ 21 мая 2019

Я понял, что проблема, как вы сказали, была в условии остановки, которое я вручную написал в функции EqMotoCur, но, учитывая мою цель, я не смог исправить ее так, как вы предлагали: вычисление остановилось бы в тот момент, когда радиус станет нулевым. Итак, я удалил это условие и написал функцию события для использования в качестве условия остановки. Моя модификация не решила проблему «деления для возможного нуля», поэтому я также переключил первые два элемента dydt и y в EqMotoCur. Таким образом, «радиус вычисления» (y(2) или вторая строка в sol{1}.y сейчас) никогда не станет ни нулем, ни ниже радиуса Марса, и гиперболическая траектория будет реалистично смоделирована с использованием любого значения для tspan в секундах. Это функция события, которую я использовал, если кому-то интересно (она не адаптирована к форме cell, потому что у меня не было времени ее адаптировать): 

function [value, isterminal, direction] = atmexit(t, y)

    R_m=3396.2;
    value=((y(2)>=R_m+250) && (y(3)>0))-0.5;
    isterminal=1;
    direction=0;

end

Я бы никогда не сделал это без твоей помощи и времени, спасибо большое!

...