Solving system of ODE with dataset of variables

Question

Albert Mamani 2019-4-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/459062-solving-system-of-ode-with-dataset-of-variables

评论： Albert Mamani 2019-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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2019-4-29

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/459062-solving-system-of-ode-with-dataset-of-variables#answer_372729

[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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Albert Mamani 2019-4-30

在 MATLAB Online 中打开

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