WATERMODEL returns a vector of length 825, but the length of initial conditions vector is 3. The vector returned by WATERMODEL and the initial conditions vector must have the

2 次查看（过去 30 天）
Nick Austinos
Nick Austinos 2022-10-9
编辑： Torsten 2022-10-9
Ran in:
I am trying to solve the ODEs defined by equations:
dL1_dt=-fTr1-fDr1;
dL2_dt=fDr1-fTr2-fDr2;
dL3_dt=fDr2-fTr3-fDr3;
But I am getting an error, I cant figure out what is wrong. Parameters used in calculating the right hand sides of the equations change after every time step, depending on the data contained in "dummyData.csv". I attach the file
tspan=1:275; % [d] 275 is length of data in data table
%Initial conditions
L0=[54;80;144]; %[mm],
[t,L]=ode45(@waterModel,tspan,L0);
Error using odearguments
WATERMODEL returns a vector of length 825, but the length of initial conditions vector is 3. The vector returned by WATERMODEL and the initial conditions vector must have the same number of
elements.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
%% ODE FUNCTION
function dLdt=waterModel(t,Y)
L1=Y(1);
L2=Y(2);
L3=Y(3);
% Load file with variables temperature, ET etc,
%%
data=readtable("dummyData.csv");
% Use the data to calculate waterparamters
p=waterParameters(data.Temperature); % Call to the nested function for calculating water parameters
% p is a structure of parameters
%% Plant transpiration
WAI1=(L1-p.pwp(1))./(p.fc(1)-p.pwp(1));
WAI2=(L2-p.pwp(2))./(p.fc(2)-p.pwp(2));
WAI3=(L3-p.pwp(3))./(p.fc(3)-p.pwp(3));
if WAI1<p.WAIc
krt1=WAI1./p.WAIc;
else
krt1=1;
end
if WAI2<p.WAIc
krt2=WAI2./p.WAIc;
else
krt2=1;
end
if WAI3<p.WAIc
krt3=WAI3./p.WAIc;
else
krt3=1;
end
%% Transpiration from each soil layer 1-3[mm]
fTr1=krt1.*p.kra.*p.krf1.*data.ET;
fTr2=krt2.*p.kra.*p.krf2.*data.ET;
fTr3=krt3.*p.kra.*p.krf3.*data.ET;
%% Soil water Drainage from both layers
if L1>p.fc(1)
fDr1=L1-p.fc(1);
else
fDr1=0;
end
if L2>p.fc(2)
fDr2=L2-p.fc(2);
else
fDr2=0;
end
if L3>p.fc(3)
fDr3=L3-p.fc(3);
else
fDr3=0;
end
%% Total Water availability index
%WAI=(L1-p.pwp(1)+L2-p.pwp(2)+L3-p.pwp(3))./(p.fc(1)-p.pwp(1)+p.fc(2)-p.pwp(2)+p.fc(3)-p.pwp(3));
%WAI = max(WAI,0);
%% Controlled inputs
%firr = 0; % [mm d-1] Irrigation
%% Differential equations
dL1_dt=-fTr1-fDr1; % dL1_dt=fpe-fTr1-fEv-fDr1
dL2_dt=fDr1-fTr2-fDr2;
dL3_dt=fDr2-fTr3-fDr3;
dLdt=[dL1_dt;dL2_dt;dL3_dt];
function p=waterParameters(Tgrh)
% Function for calculating water parameters, takes temperature of
% greenhouse [Tgrh] and returns a structure of parametrs
p.WAIc=0.75; % [-] WAI critical, range (0.5-0.8)
p.theta_fc1=0.36; % [-] Field capacity of soil layer 1
p.theta_fc2=0.32; % [-] Field capacity of soil layer 2
p.theta_fc3=0.24; % [-] Field capacity of soil layer 3
p.theta_pwp1= 0.21; % [-] Permanent wilting point of soil layer 1
p.theta_pwp2=0.17; % [-] Permanent wilting point of soil layer 2
p.theta_pwp3=0.1; % [-] Permanent wilting point of soil layer 3
p.D1=150; % [mm] Depth of Soil layer 1
p.D2=250; % [mm] Depth of Soil layer 2
p.D3=600; % [mm] Depth of Soil layer 3
p.krf1=0.25; % [-] Rootfraction layer 1 (guess)
p.krf2=0.50; % [-] Rootfraction layer 2 (guess)
p.krf3=0.25; % [-] Rootfraction layer 3 (guess)
p.kra= 0.0408.* exp(0.19.*Tgrh); % [-] Root activity, varies with Tgrh (so its a vector of kra per time step)
% [mm] Field capacities of both soil layers
p.fc=[p.theta_fc1*p.D1, p.theta_fc2*p.D2, p.theta_fc3*p.D3];
% [mm] Permanent wilting points
p.pwp = [p.theta_pwp1*p.D1, p.theta_pwp2*p.D2, p.theta_pwp3*p.D3];
end
end

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

采纳的回答

Torsten
Torsten 2022-10-9
Ran in:
% Load file with variables temperature, ET etc,
%%
data=readtable("https://de.mathworks.com/matlabcentral/answers/uploaded_files/1150115/dummyData.csv");
tspan=data.TIME(1):data.TIME(end); % [d] 275 is length of data in data table
%Initial conditions
L0=[54;80;144]; %[mm],
[t,L]=ode45(@(t,y)waterModel(t,y,data),tspan,L0);
plot(t,L)
%% ODE FUNCTION
function dLdt=waterModel(t,Y,data)
L1=Y(1);
L2=Y(2);
L3=Y(3);
% Use the data to calculate waterparamters
p=waterParameters(data.TIME,data.Temperature,t); % Call to the nested function for calculating water parameters
% p is a structure of parameters
%% Plant transpiration
WAI1=(L1-p.pwp(1))./(p.fc(1)-p.pwp(1));
WAI2=(L2-p.pwp(2))./(p.fc(2)-p.pwp(2));
WAI3=(L3-p.pwp(3))./(p.fc(3)-p.pwp(3));
if WAI1<p.WAIc
krt1=WAI1./p.WAIc;
else
krt1=1;
end
if WAI2<p.WAIc
krt2=WAI2./p.WAIc;
else
krt2=1;
end
if WAI3<p.WAIc
krt3=WAI3./p.WAIc;
else
krt3=1;
end
%% Transpiration from each soil layer 1-3[mm]
fTr1=krt1.*p.kra.*p.krf1.*interp1(data.TIME,data.ET,t);
fTr2=krt2.*p.kra.*p.krf2.*interp1(data.TIME,data.ET,t);
fTr3=krt3.*p.kra.*p.krf3.*interp1(data.TIME,data.ET,t);
%% Soil water Drainage from both layers
if L1>p.fc(1)
fDr1=L1-p.fc(1);
else
fDr1=0;
end
if L2>p.fc(2)
fDr2=L2-p.fc(2);
else
fDr2=0;
end
if L3>p.fc(3)
fDr3=L3-p.fc(3);
else
fDr3=0;
end
%% Total Water availability index
%WAI=(L1-p.pwp(1)+L2-p.pwp(2)+L3-p.pwp(3))./(p.fc(1)-p.pwp(1)+p.fc(2)-p.pwp(2)+p.fc(3)-p.pwp(3));
%WAI = max(WAI,0);
%% Controlled inputs
%firr = 0; % [mm d-1] Irrigation
%% Differential equations
dL1_dt=-fTr1-fDr1; % dL1_dt=fpe-fTr1-fEv-fDr1
dL2_dt=fDr1-fTr2-fDr2;
dL3_dt=fDr2-fTr3-fDr3;
dLdt=[dL1_dt;dL2_dt;dL3_dt];
function p=waterParameters(Time,Tgrh,t)
% Function for calculating water parameters, takes temperature of
% greenhouse [Tgrh] and returns a structure of parametrs
p.WAIc=0.75; % [-] WAI critical, range (0.5-0.8)
p.theta_fc1=0.36; % [-] Field capacity of soil layer 1
p.theta_fc2=0.32; % [-] Field capacity of soil layer 2
p.theta_fc3=0.24; % [-] Field capacity of soil layer 3
p.theta_pwp1= 0.21; % [-] Permanent wilting point of soil layer 1
p.theta_pwp2=0.17; % [-] Permanent wilting point of soil layer 2
p.theta_pwp3=0.1; % [-] Permanent wilting point of soil layer 3
p.D1=150; % [mm] Depth of Soil layer 1
p.D2=250; % [mm] Depth of Soil layer 2
p.D3=600; % [mm] Depth of Soil layer 3
p.krf1=0.25; % [-] Rootfraction layer 1 (guess)
p.krf2=0.50; % [-] Rootfraction layer 2 (guess)
p.krf3=0.25; % [-] Rootfraction layer 3 (guess)
p.kra= 0.0408.* exp(0.19.*interp1(Time,Tgrh,t)); % [-] Root activity, varies with Tgrh (so its a vector of kra per time step)
% [mm] Field capacities of both soil layers
p.fc=[p.theta_fc1*p.D1, p.theta_fc2*p.D2, p.theta_fc3*p.D3];
% [mm] Permanent wilting points
p.pwp = [p.theta_pwp1*p.D1, p.theta_pwp2*p.D2, p.theta_pwp3*p.D3];
end
end
  2 个评论
Torsten
Torsten 2022-10-9
编辑：Torsten 2022-10-9
For any time t, you must calculate the suitable values from your data according to the time of integration.
You can't supply the complete column of the data vector.
And read the data once at the beginning of your code and pass them to the ODE function as done in the code above. It would be a waste of time to read them again every time the ODE function is called.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Agriculture 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by