Question about using ODE45 to solve a differential equation using column vectors as inputs

Hello i am trying to setup an ode45 to solve this differential equation ρaCaVdTin/dT= Qshort − Qconv,cond − Qinfilt where Qshort Qconv and Qinflit are all equations describing the loss or gain of heat in a greenhouse.

What i am trying to achieve is by giving a column vector variable of external temperature wind speed and light intensity i want the ode45 to give me the internal temperature for the corresponding value of external temperature wind speed etc ie i want Tin(1) to come from ode45 which will use the first time value from my tspan the first external temperature etc.

The way i ve coded it by inputting a breakpoint inside the function i ve noticed that from the first execution of my code the internal temperature vector is already a 297x1 column vector same length as my tspan.going past the breakpoint of dT results in an error that the initial conditions vector must have the same number of elements as the vector returned. The question is, is there a way to achieve my intended result which is as i ve stated above for ode45 to use one value of my vectors to get one result? Sorry for my bad english it isn't my first language and thanks in advance.

the error i am getting which is saying that the length of the initial conditions vector needs to be the same length as the returning vector of elements.

I have included the txts from where i am getting the vectors of wind speed external temperature and solar radiation intensity.

What i am trying to achieve (and probably this isnt the right way to achieve ) is for ode45 to use for each of the 297 time variables the corresponding external temperature wind speed etc to create a 297x1 internal temperature model which is Tsol.I ve tried using a steady variable for wind speed external temperature etc and seen that this code works and creates a 297x1 internal temperature model the problem occurs when adding these column vectors as inputs. Is there another way to achieve the intended result?

below i am presenting the code

L=4;
W=3.7;
H=3;
volup=1*(1*4)/2; %ogkos trigwnikhs orofhs
voldown=W*L*2;  %ogkos katw or8ogwniou parellhlepipedou
V=voldown+volup; %Sunolikos ogkos
taxythtaaera = importdata('taxythta aera.txt');
exwterikh=importdata('metablhtes.txt');
ilios=importdata('hlios.txt');




S=L*W;

Vw=taxythtaaera*0.27777777777778;
Tac=15; %8ermokrasia se kelsioy
Ta=exwterikh+273.15; %se kelvin kai 15 se kelsiou kserw gw
Ta(265)=[];
Ta(143)=[];
Ta(53)=[];
Pa=1.137; %air density inside
Ca=1005;   %specific heat of air
sta=Pa*Ca*V;
ac=0.1;     %cover absorbtivity
tc=0.85;    %cover transmittance
I=ilios;    %kanonika pinakas olikh iliakh aktinobolia w/m^2
Sc=H*W; 
Tsky= 0.0552.*(Ta.^1.5);
Tskyc=Tsky-300;
R=0.75;        %air changes per hour APO FUSIKO AERISMO
Qshort=ac*tc*S*I;
stat=Qshort/sta;
stat2=S/sta;
stat3=5.2*(R/((Sc*L))^1/2);
stat4=R/V;
ho=2.8+1.2.*Vw;
ho(132)=[];
ho(53)=[];
ho(274)=[];
stat5=((ho*S)/sta)*(1-tc);
Kc=0.028;    %cover K w*m^-1*k^-1
Lc=0.03;    %cover thickness
%dT=Qshort/pa*ca*V-U*(S/Pa*Ca*V)(T-Tout)-(Pa*Ca*R/Pa*Ca*V)*(T-Ti)/3600-(ho*S*(1-tc)/(Pa*Ca*V)(T-Tsky)
t0=0;
To=286.15;
tend=35640;
n=297;  %giati allazei toso to apotelesma analoga me ta deigmata? sta 300 fainetai arketa sta8ero
tspan=linspace(t0,tend,n);
[tsol,Tsol]=ode45(@(t,T) funcODE(t,T,ho,Ta,I),tspan,To);
plot(tsol,Tsol-273.15)


function dT= funcODE(t,T,ho,Ta,I)



Tout=Ta;
ho=ho;
I=I;
L=4;
W=3.7;
H=3;
volup=1*(1*4)/2; %ogkos trigwnikhs orofhs
voldown=W*L*2;  %ogkos katw or8ogwniou parellhlepipedou
V=voldown+volup; %Sunolikos ogkos

S=L*W;
i=1;

Pa=1.137; %air density inside
Ca=1005;   %specific heat of air
sta=Pa*Ca*V;
ac=0.1;     %cover absorbtivity
tc=0.85;    %cover transmittance  
Sc=H*W; 
Tsky= 0.0552.*(Tout.^1.5);
R=0.75;        %air changes per hour APO FUSIKO AERISMO
Qshort=ac*tc*S*I;
stat=Qshort/sta;
stat2=S/sta;
stat3=5.2*(R/((Sc*L)^1/2));
stat4=R/V; 
stat5=((ho*S)/sta)*(1-tc);
Kc=0.028;    %cover K w*m^-1*k^-1
Lc=0.03;    %cover thickness
U=(1.0714+(1/ho)+(1/(1.52*(T-Tout).^1/3+stat3))).^-1;
dT=stat -(((R/(V*3600))*(T-Tout(i)))+(((ho*S*(1-tc))/sta)*(T-Tsky(i)))+(S*(T-Tout(i))*U(i)));
i=i+1;

Answer 1

[Torsten]

Use "interp1" to interpolate to the time t for which the solver requires the scalar values of ho, Ta and I.

Here is an example:

"ODE with time-dependent terms" under

https://de.mathworks.com/help/matlab/ref/ode45.html

[giorgos dimou]

thank you for the answer i will look into it and hopefuly see some results in the meantime can you be a bit more specific on how to implement this as i couldnt really understand a lot from the tutorial and some searching on this forum and stackoverflow on related topics?Thanks in advance

[Torsten]

Use

function main

L=4;
W=3.7;
H=3;
volup=1*(1*4)/2; %ogkos trigwnikhs orofhs
voldown=W*L*2;  %ogkos katw or8ogwniou parellhlepipedou
V=voldown+volup; %Sunolikos ogkos
taxythtaaera = importdata('taxythta aera.txt');
exwterikh=importdata('metablhtes.txt');
ilios=importdata('hlios.txt');

S=L*W;
Vw=taxythtaaera*0.27777777777778;
Tac=15; %8ermokrasia se kelsioy
Ta=exwterikh+273.15; %se kelvin kai 15 se kelsiou kserw gw
Ta(265)=[];
Ta(143)=[];
Ta(53)=[];
Pa=1.137; %air density inside
Ca=1005;   %specific heat of air
sta=Pa*Ca*V;
ac=0.1;     %cover absorbtivity
tc=0.85;    %cover transmittance
I=ilios;    %kanonika pinakas olikh iliakh aktinobolia w/m^2
Sc=H*W; 
Tsky= 0.0552.*(Ta.^1.5);
Tskyc=Tsky-300;
R=0.75;        %air changes per hour APO FUSIKO AERISMO
Qshort=ac*tc*S*I;
stat=Qshort/sta;
stat2=S/sta;
stat3=5.2*(R/((Sc*L))^1/2);
stat4=R/V;
ho=2.8+1.2.*Vw;
ho(132)=[];
ho(53)=[];
ho(274)=[];
stat5=((ho*S)/sta)*(1-tc);
Kc=0.028;    %cover K w*m^-1*k^-1
Lc=0.03;    %cover thickness
%dT=Qshort/pa*ca*V-U*(S/Pa*Ca*V)(T-Tout)-(Pa*Ca*R/Pa*Ca*V)*(T-Ti)/3600-(ho*S*(1-tc)/(Pa*Ca*V)(T-Tsky)
t0=0;
To=286.15;
tend=35640;
n=297;  %giati allazei toso to apotelesma analoga me ta deigmata? sta 300 fainetai arketa sta8ero
tspan=linspace(t0,tend,n);
[tsol,Tsol]=ode45(@(t,T) funcODE(t,T,tspan,ho,Ta,I),tspan,To);
plot(tsol,Tsol-273.15)
end

function dT= funcODE(t,T,tspan,ho,Ta,I)

Tout=interp1(tspan,Ta,t);
ho=interp1(tspan,ho,t);
I=interp1(tspan,I,t);
L=4;
W=3.7;
H=3;
volup=1*(1*4)/2; %ogkos trigwnikhs orofhs
voldown=W*L*2;  %ogkos katw or8ogwniou parellhlepipedou
V=voldown+volup; %Sunolikos ogkos
S=L*W;
Pa=1.137; %air density inside
Ca=1005;   %specific heat of air
sta=Pa*Ca*V;
ac=0.1;     %cover absorbtivity
tc=0.85;    %cover transmittance  
Sc=H*W; 
Tsky= 0.0552.*(Tout.^1.5);
R=0.75;        %air changes per hour APO FUSIKO AERISMO
Qshort=ac*tc*S*I;
stat=Qshort/sta;
stat2=S/sta;
stat3=5.2*(R/((Sc*L)^1/2));
stat4=R/V; 
stat5=((ho*S)/sta)*(1-tc);
Kc=0.028;    %cover K w*m^-1*k^-1
Lc=0.03;    %cover thickness
U=(1.0714+(1/ho)+(1/(1.52*(T-Tout).^1/3+stat3))).^-1;
dT=stat -(((R/(V*3600))*(T-Tout))+(((ho*S*(1-tc))/sta)*(T-Tsky))+(S*(T-Tout)*U));
end

[giorgos dimou]

Thank you so much this seems to be absolutely correct! i dont know if you have the answer to this as it is kinda off topic to the whole programming but do you have any idea what could be causing the spikes? there is a direct correlation between the number of samples i am taking and them but it is totaly randomised like adding samples might make it more stable adding a bit more and it makes it even more unstable... thank you for your time i wish i could give you more than the credit for the answer anyways have a great day and a happy new year!

[Torsten]

tspan=linspace(t0,tend,n);
options=odeset('RelTol',1e-8,'AbsTol',1e-8);
[tsol,Tsol]=ode45(@(t,T) funcODE(t,T,tspan,ho,Ta,I),tspan,To,options);
plot(tsol,Tsol-273.15)