Solving system of ODE with dataset of variables
3 views (last 30 days)
Show older comments
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?
0 Comments
Accepted Answer
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,..)
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!