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.

Adam · Answer 1 · 2019-05-20T17:59:12.663

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

edited May 20 '19 at 17:59

answered May 20 '19 at 03:33

Adam

2,726
1
9
22

I tried but it didn't work. It seems that the error happens during intermediate steps regardless of the outputs of `ode45`. When i use a `tspan` lower than `200000` seconds with your modification, it says `Output argument "varargout{4}" (and maybe others) not assigned during call to "ode45".` though. I truly don't understand the relation between `tspan` and the assignment error I have. – Giamarcus15 May 20 '19 at 11:35
Alright, I can't see ```y0``` in your code. Upload ```y0``` I will check it for you – Adam May 20 '19 at 11:38
You're right, I didn't notice. It's in the last `for` cycle of the main script. `y0=[.576e6;an_vera0(n);Vr0(n);Vt0(n)];` – Giamarcus15 May 20 '19 at 11:43
Also define function ```benz``` inside ```eqMotoCur```. I update my answer, check it. – Adam May 20 '19 at 12:14
I need also ```Density.mat h d``` to run the code for you – Adam May 20 '19 at 12:17
Okay, i modified `eqMotoCur` and defined function `dens` inside it, but it keeps giving me the same error. How can I give you `Density.mat h d`? – Giamarcus15 May 20 '19 at 13:19
Any link where I can download? – Adam May 20 '19 at 13:24
[Drive](https://drive.google.com/file/d/12s1tiHkolsO3F97atDE8mYI3H-MNd5ju/view?usp=sharing) – Giamarcus15 May 20 '19 at 13:31
Updated, check it – Adam May 20 '19 at 15:04
Okay, i tried it and it works now, thank you. But there's another problem: the function is executed as if time is still in seconds and gives me results that are equivalent of just a small portion of the hyperbolic trajectory. In fact the first row of `sol{1}.y` represents satellite distance from Mars center and its values go from `.576e6` to about `.575e6` now. – Giamarcus15 May 20 '19 at 15:24
I also tried to change the initial conditions in order to use a shorter `tspan` and it seems to show the error everytime `tspan` includes the moment when radius (first row of `sol{1}.y`) cross zero. Can it be caused by divisions for `y(1)` in `eqMotoCur`? – Giamarcus15 May 20 '19 at 16:04
```if abs(y(1))>R_m+250 && y(3)>0 dydt(1:4)=0;``` if ```y(1)== 0``` this condition will be false, I suggest you include it there, check the new update – Adam May 20 '19 at 17:56

Giamarcus15 · Answer 2 · 2019-05-20T22:06:06.763

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!

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

2 Answers2

The main code is as follow