Solving system of ODE with dataset of variables
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Albert Mamani
le 29 Avr 2019
Commenté : Albert Mamani
le 30 Avr 2019
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 commentaires
Réponse acceptée
Torsten
le 29 Avr 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,..)
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!