Error using odearguments (line 95) SOL returns a vector of length 18, but the length of initial conditions vector is 3. The vector returned by SOL and the initial conditions vector must have the same number of elements.

I am trying to solve 3 ODEs by using ode45. i dont know if im doing this right or not, I appreciate any small help.

this is my code.

my main problem is how to include the transient effect of solar radiation. I tried to solve it with constant solar radiation and it worked fine.

t1=8; %start time
t2=15; %end time 
tspan=[0,(t2-t1)*3600];
IC=[38,40,37]; %initial temperatures 

[time,temp]=ode45(@sol,tspan,IC);
function tdot =sol(t,T)
t1=8;
t2=15;
vol=4;
mair=vol*1.225;
cp=1005;
Qh=0;%220;
Tac=15;
mac=0.09;
Uinwall=4;
Uoutwall=4.5;
Uinmass=12;
Amass=3;
Tamb=45;
cpwall=600;
cpmass=1450;
mwall=900;
mmass=50;
emiss=0.9;
Tsky=22.6;
abswall=0.19;
absmass=1;
trans=0.8;
A=[2.479,0 %roof
    0.864,0.864  %windsheild
    0.787,0.787 %rear window
    1.533+0.224+0.266,0.224+0.266 %right side
    1.533+0.224+0.266,0.224+0.266]; %left side

Aw=A(:,1); %area of walls 
Ag=A(:,2); % area of glass
Awall=sum(Aw,'all');
AHwall=Aw(1,:);
lat=33.3152; %latitude  
long=44.3661; % longitude 
gmt=3;
mint=0;
tiltang_w=[0 % tilt angle of the walls
    60
    45
    90
    90];
zs_w=[0 % zenith angle of the walls
    0
    180
    90
    270];
tiltang_g=[60  % tilt angle of the glass
    45
    90
    90];
zs_g=[0  % zenith angle of the glass
    180
    90
    270];
costhetaz=[]; % cosine of the zenith angle
costheta_w=[]; %incident angle of solar radiation on the walls
costheta_g=[]; %incident angle of solar radiation on the glass
for hour=t1:1:t2
    lot=hour+mint/60;%local time of the city
    dn=190; % day of the year
    a=1160+75*sind((dn-275)*360/365); %constants of the equation
    b=0.174+0.035*sind((dn-100)*360/365); %constants of the equation
    c=0.095+0.04*sind((dn-100*360/365)); %constants of the equation
    delta=23.45*sind((dn+284)*360/365); % Declination angle
    bn=(dn-1)*360/365;
    eot=2.2918*(0.0075+0.1686*cosd(bn)-...
        3.2077*sind(bn)-1.4615*cosd(2*bn)-...
        4.089*sind(2*bn));  % equation of time
    lst=15*gmt; % local standard time
    st=(lot+(eot/60)+(long-lst)/15); % solar time
    omega=15*(12-st); % hour angle
    costhetaz_i=cosd(delta).*cosd(lat).*cosd(omega)+...
        sind(delta).*sind(lat); % cosine of the zenith angle
        
    costheta_wi=sind(delta).*sind(lat).*cosd(tiltang_w)-...
        sind(delta).*cosd(lat).*sind(tiltang_w).*cosd(zs_w)+...
        cosd(delta).*cosd(lat).*cosd(tiltang_w).*cosd(omega)+...
        cosd(delta).*sind(lat).*sind(tiltang_w).*cosd(zs_w).*cosd(omega)+...
        cosd(delta).*sind(tiltang_w).*sind(zs_w).*sind(omega);  %incident angle of solar radiation on the walls
    
    costheta_gi=sind(delta).*sind(lat).*cosd(tiltang_g)-...
        sind(delta).*cosd(lat).*sind(tiltang_g).*cosd(zs_g)+...
        cosd(delta).*cosd(lat).*cosd(tiltang_g).*cosd(omega)+...
        cosd(delta).*sind(lat).*sind(tiltang_g).*cosd(zs_g).*cosd(omega)+...
        cosd(delta).*sind(tiltang_g).*sind(zs_g).*sind(omega);%incident angle of solar radiation on the glass 
    
    costhetaz=[costhetaz;costhetaz_i];
    costheta_w=[costheta_w;costheta_wi];
    costheta_g=[costheta_g;costheta_gi];
    
end
hour=t1:1:t2;
costhetaz;
costheta_w;
costheta_g;
minLength = min(length(costhetaz), length(costheta_w));  % i wrote this line becasue i didnt get the same number of values for costheatz,costheta_w and costheta_g 
costhetaz = costhetaz(1:minLength);
costheta_w = costheta_w(1:minLength);
costheta_g = costheta_g(1:minLength);
[idir]=a.*exp(-b./costhetaz); % direct solar radiation
costheta_w(costheta_w<0)=0; % this line to neglect the negative values of costheta_w and costheta_g
costheta_g(costheta_g<0)=0;
rad_w=idir.*costheta_w; % solar radiation on the walls 
rad_g=idir.*costheta_g; % solar radiation on the glass
totalrad_w=0;
totalrad_g=0;
for i=1:1:5
    totalrad_w=totalrad_w+rad_w.*Aw(i,1);
end
for i=1:1:5
    totalrad_g=totalrad_g+rad_g.*Ag(i,1);
end

Q_ac=0;%(mac*cp*(Tac-T(1)));
Qin_w=(Uinwall*Awall*(T(1)-T(2)));
Qin_m=(Uinmass*Amass*(T(1)-T(3)));
Qout_w=(Uoutwall*Awall*(Tamb-T(2)));
Qin_rad=(5.67*10^-8*Amass*((T(3)+273)^4-(T(2)+273)^4));
Qout_rad=-(5.67*10^-8*emiss*AHwall*((T(2)+273)^4-(Tsky+273)^4));
Qabssolar=abswall*totalrad_w;
Qtranssolar=trans*absmass*totalrad_g;
tdot=[(Qh+Q_ac-Qin_w-Qin_m)/(mair*cp)
    (Qin_w+Qout_w+Qin_rad+Qabssolar-Qout_rad)/(mwall*cpwall)
    (Qin_m+Qtranssolar-Qin_rad)/(mmass*cpmass)];
end

Answer 1

[Alan Stevens]

Your calculations of rad_w and rad_g produce vectors of size 8x1 that lead to tdot being larger than 3x1. Should these two (i.e. rad_w and rad_g) be the sum of the values in the vectors?

[Hasan Abbas]

rad_w and rad_g are the solar radiation on each side (i.e the walls and the glass, which are Aw and Ag) , then i need to sum rad_w and rad_g on all sides at each time step (in this loop : hour=t1:1:t2), because the solar radiation is chainging at each hour which makes the temperature to rise

[Alan Stevens]

So, my question was: instead of

rad_w=idir.*costheta_w; % solar radiation on the walls 
rad_g=idir.*costheta_g; % solar radiation on the glass

should you have

rad_w=sum(idir.*costheta_w); % solar radiation on the walls 
rad_g=sum(idir.*costheta_g); % solar radiation on the glass

[Hasan Abbas]

No, that will not work (i tried it), because it will sum them at all time steps.

Can i make the whole solar radiation change with the time span?

[Alan Stevens]

Currently, each time ode45 calls sol with a value of t, you use a range of times (t1 to t2) within sol. So, for example, if ode45 is considering a time prior to t1, sol still does a calculation from t1 to t2. Surely, this can't be correct. Instead of looping from t1 to t2 within sol, shouldn't you have

if t>=t1 & t<=t2
    ...% all the calcualtions within the existing loop
end

If so, you would also need to set appropriate values parameters at times that aren't between t1 and t2.

[Hasan Abbas]

excuse me, can you explain more?

[Alan Stevens]

Each time ode45 calls function sol it does so for a single time, t. Therefore, one might expect that all the calculations within sol should apply for that time only. However, you have written sol such that times between t1 and t2 are used within it to calculate some parameters (such as lot, st, omega etc.). But, you can't have the cumulative value of, say, lot, accumulated over times between t1 and t2, applying for a specific time, t (certainly not if t precedes t2).

[Hasan Abbas]

So, what i understand is that i cant make loop within sol between t1 and t2 because ode45 calls the function only once. so how can i make rad_w and rad_g changes with every step of the time span?

thanks for all your replies!

[Alan Stevens]

As I noted before, where you currently have

for hour=t1:1:t2

replace this with

if t>=t1 & t<=t2

then work out lot etc for this specific time using the value of t.

[Hasan Abbas]

well, I have this error now

Unrecognized function or variable 'a'.

Error in test>sol (line 205)
[idir]=a.*exp(-b./costhetaz); % direct solar radiation

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in test (line 51)
[time,temp]=ode45(@sol,tspan,IC);

but when I replaced the 'a' and 'b' with their numerical values another error has appeared

Error using odearguments (line 95)
SOL returns a vector of length 1, but the length of initial conditions vector is 3. The vector returned by SOL and
the initial conditions vector must have the same number of elements.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in test (line 51)
[time,temp]=ode45(@sol,tspan,IC);