using ode45 to solve 3 ODEs with solar radiation

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.

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.

this is my code and when i run it i get this error.

Thanks in advance!

Error using odearguments (line 95)
SOL returns a vector of length 17, 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)

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
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

[Jan]

Omit useless lines as:

costhetaz;
costheta_w;
costheta_g;

The line "hour=t1:1:t2;" can be removed also, because "hour" is not used anymore.

Replace

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

by the nicer and faster:

totalrad_w = sum(rad_w .* Aw);
totalrad_g = sum(rad_g .* Ag);

While 5.67*10^-8 is an expensive power operation and a mutliplication, 5.67e-8 is a cheap constant.

[Hasan Abbas]

thanks Jan, I already slove it but i have another problem if you could help me

i have a vector of (1*n) as input for the ode45, how can i do that so the ode45 takes one value at each time step?

[Jan]

You can't. You cannot know in advance, which time steps are chosen and the integrator can and does reject steps, when the step size controller detects an exceeded limit.

Please explain, which problem you want to solve. Providing different parameters for each time step sound like you try to integrate a non-smooth function. This would be a bad idea, because ODE45 requires a smooth function.

[Hasan Abbas]

hello Jan

it is the same problem mentioned above, i solved the pervious error, i found that the terms "Qabssolar" and "Qtranssolar" were vectors of length (1*n) where n is the same as time span, so i searched for a solution and wrote this code and it worked well

Qabs_t = linspace(t1*3600,t2*3600,h);
Qtrans_t = linspace(t1*3600,t2*3600,h);
Qabssolar=abswall*totalrad_w;
Qtranssolar=trans*absmass*totalrad_g;
Qabssolar = interp1(Qabs_t, Qabssolar, t); 
Qtranssolar = interp1(Qtrans_t, Qtranssolar, t); 

i guess it`s right, because the result were as i expected. do you think it`s right or wrong?

this is my full code

t1=8;
t2=15;
h=((t2-t1)*3600)/60;
tspan=linspace(t1*3600,t2*3600,h);
IC=[32,33,35];
[time,temp]=ode45(@sol,tspan,IC);

function tdot =sol(t,T)
t1=8;
t2=15;
h=((t2-t1)*3600)/60;
n=((t2-t1)*3600)/h;
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=46;
cpwall=600;
cpmass=1450;
mwall=900;
mmass=50;
emiss=0.9;
Tsky=22.6;
abswall=0.19;
absmass=1;
trans=0.8;
Aw=[2.479 %roof
    0.864  %windsheild
    0.787 %rear window
    1.533+0.224+0.266 %right side
    1.533+0.224+0.266]; %left side
Ag=[0.864 %windsheild
    0.787 %rear window
    0.224+0.266 %right windows
    0.224+0.266];%left windows
Awall=sum(Aw,'all');
AHwall=Aw(1,:);
lat=33.3152;
long=44.3661;
gmt=3;
dn=212; % 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;
mint=0;
tiltang_w=[0
    60
    45
    90
    90];
zs_w=[0
    0
    180
    90
    270];
tiltang_g=[60
    45
    90
    90];
zs_g=[0
    180
    90
    270];
costhetaz=[];
costheta_w=[];
costheta_g=[];

for hour=t1:1/n:t2-1/60
    lot=hour+mint/60;%local time of the city
    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);
    
    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);
    
    costhetaz=[costhetaz;costhetaz_i];
    costheta_w=[costheta_w;costheta_wi];
    costheta_g=[costheta_g;costheta_gi];
    
end

idir=a.*exp(-b./costhetaz); % direct solar radiation
idif=c.*idir;
i_tot=idir+idif;

costheta_w(costheta_w<0)=0;
costheta_g(costheta_g<0)=0;


rad_w=zeros(5,h);
for i=1:h
    for j=1:5
        rad_w(j,i)=rad_w(j,i)+ i_tot(i,1).*costheta_w(j,1);
    end
end

rad_g=zeros(4,h);
for i=1:h
    for j=1:4
        rad_g(j,i)=rad_g(j,i)+ i_tot(i,1).*costheta_g(j,1);
    end
end

totalrad_w=zeros(5,h);
totalrad_g=zeros(5,h);
for i=1:h
    for j=1:1:5
        totalrad_w(j,i)=totalrad_w(j,i)+rad_w(j,i).*Aw(j,1);
    end
end
totalrad_w=sum(totalrad_w,1);

for i=1:h
    for j=1:1:4
        totalrad_g(j,i)=totalrad_g(j,i)+rad_g(j,i).*Ag(j,1);
    end
end
totalrad_g=sum(totalrad_g,1);

Qabs_t = linspace(t1*3600,t2*3600,h);
Qtrans_t = linspace(t1*3600,t2*3600,h);
Qabssolar=abswall*totalrad_w;
Qtranssolar=trans*absmass*totalrad_g;
Qabssolar = interp1(Qabs_t, Qabssolar, t); 
Qtranssolar = interp1(Qtrans_t, Qtranssolar, t); 

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.67e-8*Amass*((T(3)+273)^4-(T(2)+273)^4));
Qout_rad=(5.67e-8*emiss*AHwall*((T(2)+273)^4-(Tsky+273)^4));

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

[Jan]

I do not have information about the purpose of your code. Therefore I cannot estimate, if the code is working as wanted or not. At least the variables Tac and mac are not used anywhere.

The expression "(T(3)+273)^4" seems to be very sensitive. What about using more digits like 273.15?

Qabssolar = interp1(Qabs_t, Qabssolar, t); 
Qtranssolar = interp1(Qtrans_t, Qtranssolar, t);

Linear interpolations are not smooth. The stepsize control of ODE45 requires a smooth function. This can cause serious troubles, e.g. results, which are dominated by rounding errors. A stopping integration is possble also. Running a numerical integrator with not matching inputs is not a reliable method from a scientific point of view. If this is a part of a PhD thesis, I would reject it, if I am your professor.

The code does not allow to check the stability of the solution by varying the initial conditions and parameters. Did you validate the applied model?

Answer 1

[Mischa Kim]

Hi Hasan,

based on the error message I suggest looking into tdot at the bottom of your ode function. Apparently, this is a 17x1 vector and by browsing through your code it seems that Qabssolar and Qtranssolar are also vectors, not scalars.

I suggest you set a breakpoint next to the last end command in your code above and run the code. The execution of the code will halt when the breakpoint is reached and you can hover over the variables to see value and size. Once you have confirmed that Qabssolar and Qtranssolar are indeed the cause of the issue you can work your way back up and isolate the source of the problem.

[Hasan Abbas]

Thanks for your answer, i checked it and it is indeed a vector. but i have another problem, i will appreciateit if you can help me. I want the values of rad_w and rad_g to change by taking values from t1 to t2 and i want the ode45 to evaluate the 'sol' function at each time span corresponding to that time. so i have time span from t1 to t2 and i have for loop from t1 to t2.

Also, in that for loop i get vector of 8*1 for costhetaz but i get a vector of 40*1 for costheta_w and a vector of 32*1 for costheta_g because i have vectors of 5*1 for the walls and 4*1 for the glass, how can i make it all simultaneous so at each time span it take 1 value of costheaz and 5 values of costheta_w and 4 values of costheat_g?

if i have done this wrong, please explain it to me

thanks again