Solving system of ODE with dataset of variables
4 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Torsten
2019-4-29
[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,..)
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!