Solving system of ODE with dataset of variables

4 views (last 30 days)
I´m trying to solve a Temperature in a reservoir with a system of two Differential Equation, the code for the system is:
function dTsdt = Temp_EH(t,Ts,Qin,Qout,Tin,rad,Tair,Vol_hip,Vol_ep,H_ep,AreaThermo,eair,f,...
dens,Cp,sigma,A,RL,eee,c1,TransferCoef)
dTsdt = zeros(2,1);
dTsdt(1) = Qin*Tin/Vol_ep + rad/(dens*Cp*H_ep)+(sigma*(Tair+273)^4)*(A+0.031*sqrt(eair))*(1-RL)/(dens*Cp*H_ep)
- Qout*Ts(1)/Vol_ep - (eee*sigma*(Ts(1)+273)^4)/(dens*Cp*H_ep)
- c1*f*(Ts(1)-Tair)/(dens*Cp*H_ep)
- f*(4.596*exp(17.27*Ts(1)/(237.3+Ts(1)))-eair)/(dens*Cp*H_ep);
dTsdt(2) =((TransferCoef*AreaThermo)/Vol_hip)*(Ts(1)-Ts(2))
end
I have a dataset, for time, vol_tot_hm3, AS_km2, Qin_m3s1, Qout_m3s1 ,Tin_C, rad_Wm2, Tair_C, Uw_ms1 and Rhum (2922 values for each one)
for i=1:(length(time)-1)
%Import data
vol=vol_tot_hm3(i)*1000000; % a m3
AS=AS_km2(i)*1000000; % a m2
Qin=Qin_m3s1(i)*(3600*24); % a m3/d
Qout=Qout_m3s1(i)*(3600*24); % a m3/d
Tin=Tin_C(i); % °C
rad=rad_Wm2(i)/41867.280720117*86400; % cal/cm2d
Tair=Tair_C(i); % °C
Uw=Uw_ms1(i); % m/s
RH=Rhum(i); % /100
%Constants
dens = 1; % Densidad del agua
Cp =1; % Calor Especifico del agua
sigma = 11.7*(10^-8); % Cte de Stefan Boltzmann
A =0.6; % Coef atenuacion atmosferica
RL=0.03; % Coef de Reflexion
eee=0.97; % Emisividad de un cuerpo radiante(agua en este caso)
c1=0.47; % Coef de Bowen (Conduccion y conveccion)
TransferCoef = 0.034; %m/d %(vt)
%Variables
Vol_ep= vol/2; %m3
Vol_hip = vol - Vol_ep; %m3
Termocl = -7.548*log((vol/1000000)/2) + 39.773; % Altura en metros de epilimnio o prof de termoclina
% log(x) es ln(x) en MATLAB
H_ep = Vol_ep/AS;
AreaThermo = (-0.0002*Termocl^4 + 0.0097*Termocl^3 - 0.1741*Termocl^2 + 0.2523*Termocl + 14.298)*1000000; %m2
Eigenvalorh= (TransferCoef*AreaThermo)/Vol_hip; %dias (para Hipolimnio)
esat=4.596*exp(17.27*Tair/(237.3+Tair)); %Presion de vapor de saturacion de aire
eair=RH*esat ; %Presion de vapor del aire
% "es" corresponde a Presion de vapor de saturacion de agua
f=19+0.95*Uw^2; %Func de transferencia de Vel viento hacia el agua
%Solving ODE system
tspan=[0 2922];
Tsi=[13 8.5];
[t,Ts]=ode45(@Temp_EH,tspan,Tsi)
end
Before use "For", I test the function Temp_EH, and it generate error:
Not enough input arguments.
Error in Temp_EH (line 4)
dTsdt(1) = Qin*Tin/Vol_ep + rad/(dens*Cp*H_ep)+(sigma*(Tair+273)^4)*(A+0.031*sqrt(eair))*(1-RL)/(dens*Cp*H_ep)
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);
I don´t realize what the problem is with the Temp_EH function, and the solution for the ODE system for the 2922 values of the dataset is correct?

Accepted Answer

Torsten
Torsten on 29 Apr 2019
[t,Ts]=ode45(@(t,Ts)Temp_EH(t,Ts,time,Qin,Qout,Tin,rad,Tair,Vol_hip,Vol_ep,H_ep,AreaThermo,eair,f,dens,Cp,sigma,A,RL,eee,c1,TransferCoef),tspan,Tsi)
function dTsdt = Temp_EH(t,Ts,time,Qin,Qout,Tin,rad,Tair,Vol_hip,Vol_ep,H_ep,AreaThermo,eair,f,dens,Cp,sigma,A,RL,eee,c1,TransferCoef)
Qin_at_t = interp1(time,Qin,t);
Qout_at_t = interp1(time,Qout,t);
...
dTsdt = zeros(2,1);
dTsdt(1) = Qin_at_t ...
end
And use elementwise multiplication and division when working with arrays, e.g.
H_ep = Vol_ep./AS
instead of
H_ep = Vol_ep/AS
(same for AreaThermo,Eigenvalorh,..)
  1 Comment
Albert Mamani
Albert Mamani on 30 Apr 2019
Thanks a lot
My function works now:
function dTsdt = Temp_EH(t,Ts,time,Qin,Qout,Tin,rad,Tair,Vol_hip,Vol_ep,H_ep,AreaThermo,eair,f)
%Sistema de EDO para Temp epilimnio e Hipolimnio
Qin_at_t = interp1(time,Qin,t);
Qout_at_t = interp1(time,Qout,t);
Tin_at_t = interp1(time,Tin,t);
rad_at_t = interp1(time,rad,t);
Tair_at_t = interp1(time,Tair,t);
Vol_hip_at_t = interp1(time,Vol_hip,t);
Vol_ep_at_t = interp1(time,Vol_ep,t);
H_ep_at_t = interp1(time,H_ep,t);
AreaThermo_at_t = interp1(time,AreaThermo,t);
eair_at_t = interp1(time,eair,t);
f_at_t = interp1(time,f,t);
%Constants
dens = 1; % Densidad del agua
Cp =1; % Calor Especifico del agua
sigma = 11.7*(10^-8); % Cte de Stefan Boltzmann
A =0.6; % Coef atenuacion atmosferica
RL=0.03; % Coef de Reflexion
eee=0.97; % Emisividad de un cuerpo radiante(agua en este caso)
c1=0.47; % Coef de Bowen (Conduccion y conveccion)
TransferCoef = 0.034; %m/d %(vt)
dTsdt = zeros(2,1);
dTsdt(1) = Qin_at_t.*Tin_at_t./Vol_ep_at_t + rad_at_t./(dens.*Cp.*H_ep_at_t)+(sigma.*(Tair_at_t+273).^4)*(A+0.031.*sqrt(eair_at_t)).*(1-RL)./(dens.*Cp.*H_ep_at_t)...
-Qout_at_t.*Ts(1)./Vol_ep_at_t - (eee.*sigma.*(Ts(1)+273).^4)/(dens.*Cp.*H_ep_at_t)...
- c1.*f_at_t.*(Ts(1)-Tair_at_t)/(dens.*Cp.*H_ep_at_t) - f_at_t.*(4.596*exp(17.27.*Ts(1)/(237.3+Ts(1)))-eair_at_t)/(dens.*Cp.*H_ep_at_t);
dTsdt(2) =((TransferCoef.*AreaThermo_at_t)./Vol_hip_at_t).*(Ts(1)-Ts(2));
dTsdt =dTsdt(:);
end

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!