How to fix assignment error using ode45 in MatLab (line 488 of ode45 function)

Question

I'm trying to write a script that uses ode45 in order to integrate the equations of motion of a satellite on an hyperbolic trajectory near Mars.

I need to integrate the entire passage around the planet: starting from SOI radius (576000km) and proceeding towards the planet, then through the atmosphere until the satellite reaches the "opposite" atmosphere boundary (set at 250km from the surface).

When it receives in input a tspan higher than about 200000 seconds (I need approximately 400000 seconds), Matlab gives me this message:

'Unable to perform assignment because the size of the left side is 4-by-2 and the size of the right side is 4-by-5.' It says that the error happens in line 488 of Ode45.

I searched for some similar cases but i couldn't find anything and i don't know how to use conditional break points to figure out something on a complex function such as ode45. I also tried different options using "odeset" but nothing changed. I have no clue where the error could be.

This is the script, in which i use two additional functions to obtain some of the parameters:

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

This is "parasint" function:

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

This is "time" function:

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

This is the input function for 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

This is "dens" function for interpolation of data through atmosphere:

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

This is the error Matlab gives me:

Unable to perform assignment because the size of the left side is 4-by-2 and the size of the right side is 4-by-5.

Error in ode45 (line 488) yout(:,idx) = yout_new;

Error in MainCur (line 59) [t,Y]=ode45(@(t,y) eqMotoCur(t,y),I,y0,options);

Thank you in advance.

Answer 1

• ode45(@(t,y) eqMotoCur(t,y),I,y0,options) returns a 5 by 4 matrix
• [t,Y]is a 2 by 4 matrix
• Thus dimension mismatch as 5 by 4 # 2 by 4
• Replace [t,Y] by [t,Y, ~,~, ~]
• Only keep the first two variables you want, and ignore the remaining by using ~

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

Make sure you define y0 first before running ode45

function dens(z) should be defined inside eqMotoCur(t,y) as you need it to perform some calculations inside 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

The main issue is the time t, is quite large if it's evaluated in second and that's where ode45 gets stuck.
Alternately convert second to hours(devide it by 3600), it will work. Also the result obtained is struct data type where x denotes t and y denotes y. let's say t = sol{1}.x and y = sol{1}.y
Since the array will be changing every iteration I use a cell to keep them separately

The main code is as follow

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

Answer 2

I figured out that the problem was, as you said, in the stopping condition that i manually wrote in EqMotoCur function, but considering my purpose I couldn't correct it in the way you suggested: the computation would have stopped in the moment when the radius would become zero. So, I deleted that condition and I wrote an event function to use as a stopping condition. My modification didn't solve the "division for a possible zero" issue, so I also switched the first two elements of dydt and y in EqMotoCur. In this way, the "computation radius" (y(2) or second row in sol{1}.y now) will never become neither zero nor lower than Mars radius and the hyperbolic trajectory will be simulated realistically using any value for tspan in seconds. This is the event function I used, if anyone is interested (it's not adapted to the cell form because i didn't have time to adapt it):

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

I could never have done it without your help and time anyway, thank you so much!