## Solving system of ODE with dataset of variables

Albert Mamani

### Albert Mamani (view profile)

さんによって質問されました 2019 年 4 月 29 日

### Albert Mamani (view profile)

さんによって コメントされました 2019 年 4 月 30 日
Torsten

### Torsten (view profile)

さんの 回答が採用されました
I´m trying to solve a Temperature in a reservoir with a system of two Differential Equation, the code for the system is:
dens,Cp,sigma,A,RL,eee,c1,TransferCoef)
dTsdt = zeros(2,1);
- 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
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
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)
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 件のコメント

サインイン to comment.

## 1 件の回答

### Torsten (view profile)

2019 年 4 月 29 日
採用された回答

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
H_ep = Vol_ep/AS
(same for AreaThermo,Eigenvalorh,..)

Albert Mamani

### Albert Mamani (view profile)

2019 年 4 月 30 日
Thanks a lot
My function works now:
%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);
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
c1=0.47; % Coef de Bowen (Conduccion y conveccion)
TransferCoef = 0.034; %m/d %(vt)
dTsdt = zeros(2,1);