one of the terms "Rho_net" on the RHS of an set is varying with time over the tspan. How do I implement the set of the ODEs such that the term updates at every time step?

I have been trying the code as below but it returns the error "dYdt = @(t, Y) [ode1(t, Y); ode2(t, Y)] must return a column vector"

clc; clear; close all;

% Solar constant [MJ m-2 min-1]
Gsc = 0.08;
lat = 0; % Latitude

% Define simulation period
startdate = datetime(2025, 01, 01, 'Format', 'uuuu-MM-dd');
finishdate = datetime(2025, 12, 31, 'Format', 'uuuu-MM-dd');
tspans = [1, days(finishdate - startdate) + 1];

% Initialize arrays to store MM, doy, and dd values
MM_array = zeros(1, days(finishdate - startdate) + 1);
doy_array = zeros(1, days(finishdate - startdate) + 1);
dd_array = zeros(1, days(finishdate - startdate) + 1);

% Convert date to datetime object
for idx = 1:days(finishdate - startdate) + 1
    current_date = startdate + days(idx - 1);
    MM = month(current_date);
    dd = day(current_date);
    doy = days(current_date - datetime(year(current_date), 1, 1)) + 1;

    % Store values in arrays
    MM_array(idx) = MM;
    doy_array(idx) = doy;
    dd_array(idx) = dd;
end

% Convert latitude to radians
phi = deg2rad(lat);

% Calculate inverse relative distance Earth-Sun (dr)
dr = 1 + 0.033 * cos(2 * pi * doy_array / 365);

% Calculate solar declination (delta)
delta = 0.409 * sin(2 * pi * doy_array / 365 - 1.39);

% Calculate sunset hour angle (ws)
ws = acos(-tan(phi) * tan(delta));

% Calculate extraterrestrial radiation (Ra)
Ra = 1000*extraterrestrialRadiation(Gsc, dr, pi, phi, ws, delta);

% Display the result
disp(['Extraterrestrial Radiation (Ra): ' num2str(Ra(1)) ' MJ/m^2/day']); % Displaying only the first value

Extraterrestrial Radiation (Ra): 34874.1716 MJ/m^2/day

% Constants and parameters
As = 0.06;
as = 0.25;
bs = 0.5;
N = (2 / 15) .* acosd(-tan(lat) .* tan(rad2deg(delta)));
n = 11;
Rs = (as + (bs .* (n ./ N))) .* Ra;

% Calculate net shortwave solar radiation (and photoperiod/24(phi_Light))
pi_light = N / 24;
Rho_SWR = Rs * (1 - As);

% Calculate rho_an (the difference between the incident and reflected components of long-wave radiation)
sigma = 4.896 * 10^(-6);
T_a = 25 * ones(size(dd_array)); % Later read from file/ it should be a data series
T_ak = Absolute_temp1(T_a);
r = 0.5; % Change the value as necessary
C_c = 0.9; % Read from file/enter
epslon_a = 0.937 * 10^-5 * T_ak.^2 * (1 + 0.17 * C_c^2);
Rho_an = (1 - r) * epslon_a * sigma .* T_ak.^4;

% Calculate Water Surface Radiation (HO2_heat_Emm)
T_w = 18; % To be simulated first
epslon_w = 0.97;
T_wk = Absolute_temp2(T_w);
Rho_HO2_heat_Emm = epslon_w * sigma * T_wk.^4;

% Calculate Evaporative Heat Loss (Evap_HLoss)
U2 = 1.3; % Read from file
RH = 0.75; % Read from file
bo = 368.61;
lambda = 311.02;
es = 25.37 * exp(17.62 - 5271 ./ T_wk);
ea = RH * 25.37 * exp(17.62 - 5271 ./ T_ak);
z = 1000; % Entered in the app
P = 760 / (10^(z / 19748.2));
T_wv = (T_wk / (1.0 - (0.378 * es / P)));
T_av = (T_ak / (1.0 - (0.378 * ea / P)));
Rho_Evap_HLoss = (es - ea) * (lambda * ((T_wv - T_av))^(1 / 3) + bo * U2);

% Calculate Conductive Heat Loss or Gain (Heat_cond)
Rho_Heat_cond = Rho_Evap_HLoss * 0.61 * 10^-3 * P * ((T_wk - T_ak) / (es - ea));

Rho_net_data = Rho_SWR + Rho_an - Rho_HO2_heat_Emm - Rho_Evap_HLoss + Rho_Heat_cond;
% Define Rho_net as an anonymous function of time
Rho_net = @(t) interp1(tspans(1):tspans(end), Rho_net_data, t, 'linear', 'extrap');
figure
plot(tspans(1):tspans(end), Rho_net_data(tspans(1):tspans(end)), '.', 'markersize', 3), grid on, 

Warning: Imaginary parts of complex X and/or Y arguments ignored.

title('\rho_{net}'), xlabel('Day number'), xlim(tspans)

% Computation of water balance components
A = 100;
pond_depth = 1;
V = A * pond_depth;
dmin = 0.8;
dcurr = 0.6;
dmax = pond_depth;
Ir = 12.923;

if Ir == 0
    Qi = A * (dmin - dcurr) / tspans(2);
else
    Qi = Ir / 100 * V;
end
Er = 10;
if Er == 0
    Qe = A * (dcurr - dmax) / tspans(2);
else
    Qe = Er / 100 * V;
end

pd = 10; % To be read from a file
PCP = A * pd / 1000;
Density_w = 1000;
Latent_vapw = 2260;
Evap = A * Rho_Evap_HLoss / (Density_w * Latent_vapw);
Sr = 5;
Seep = A * Sr / 1000;
CWR = 20; % To be read from Aquacrop
A_Irr = 200;
Irr = A_Irr * CWR / 1000;

% Initial conditions for the ODEs
T_win = 18; % Initial water temperature
T_wout = 18;
Cpw = 5;

% Define your ODEs
ode1 = @(t, Y) Qi + PCP - Qe + Evap + Seep - Irr;
ode2 = @(t, Y) Qi * T_win / Y(1) - Qe * T_wout / Y(1) + Rho_net(t) / (Density_w * Cpw * pond_depth) - Y(2) / Y(1) * ode1(t, Y);

% Combine the ODEs into a single function
dYdt = @(t, Y) [ode1(t, Y); ode2(t, Y)];

% Set initial conditions
initial_conditions = [V; T_w]; % Ensure that this is a column vector

% Solve the system of ODEs using ode45
[t, Y] = ode45(dYdt, tspans, initial_conditions);

Error using odearguments
@(T,Y)[ODE1(T,Y);ODE2(T,Y)] must return a column vector.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% Extract results
V_solution = Y(:, 1);
T_w_solution = Y(:, 2);

% Display or use the solutions as needed
disp('Volume (V) solution:');
disp(V_solution);
disp('Water Temperature (T_w) solution:');
disp(T_w_solution);

% Plot the results
figure;
subplot(2, 1, 1);
plot(t, V_solution); grid on
title('Volume (V) vs Time');
xlabel('Time');
ylabel('Volume (V)');

subplot(2, 1, 2);
plot(t, T_w_solution); grid on
title('Water Temperature (T_w) vs Time');
xlabel('Time');
ylabel('Water Temperature (T_w)');

function Ra = extraterrestrialRadiation(Gsc, dr, pi, phi, ws, delta)
    % Calculate extraterrestrial radiation (Ra)
    Ra = (24 * 60 / pi) * Gsc * dr .* (ws .* sin(phi) .* sin(delta) + cos(phi) .* cos(delta) .* sin(ws));
end

function Absolute_1 = Absolute_temp1(T_a)
    Absolute_1 = 273 + T_a;
end

function Absolute_2 = Absolute_temp2(T_w)
    Absolute_2 = 273 + T_w;
end

Answer 1

[Torsten]

If you insert the line

size(dYdt(tspans(1),initial_conditions))

before the call to ode45, you will see that your function dYdt still returns a matrix of size 2x365 instead of a vector of size 2x1.

So here seem to be other variables in the function definition for which the same interpolation procedure has to be applied as you did for Rho_net. The one I could identify is Evap.

If I were you, I wouldn't try to supply dYdt as a function handle, but compute the time derivatives in an extra function. Things will become much easier this way.

[Daniel]

Thank you very much @Torsten for looking into my code and for your suggestion. Actually i have realized as per your suggestion that I should do the same for Evap as i did for Rho_net.

Answer 2

[Sam Chak]

Hi @Daniel

The majority of your code remains unchanged. I simply relocated your ODEs into a function named 'dYdt' and incorporated the interpolation code, similar to what you learned from the example involving time-dependent terms in the ode45 documentation. The remaining modifications involve passing the additional parameters.

As mentioned previously in your other thread, the simulated results obtained by interpolating the day-dependent Rho_net should show a close resemblance to the system response when using the average value of Rho_net, because it varies around throughout the entire year.

% Solar constant [MJ m-2 min-1]
Gsc         = 0.08;
lat         = 0; % Latitude

% Define simulation period
startdate   = datetime(2025, 01, 01, 'Format', 'uuuu-MM-dd');
finishdate  = datetime(2025, 12, 31, 'Format', 'uuuu-MM-dd');
tspans      = [1, days(finishdate - startdate) + 1]

tspans = 1×2

1 365

% Initialize arrays to store MM, doy, and dd values
MM_array    = zeros(1, days(finishdate - startdate) + 1);
doy_array   = zeros(1, days(finishdate - startdate) + 1);
dd_array    = zeros(1, days(finishdate - startdate) + 1);

% Convert date to datetime object
for idx = 1:days(finishdate - startdate)+1
    current_date    = startdate + days(idx - 1);
    MM              = month(current_date);
    dd              = day(current_date);
    doy             = datenum(current_date) - datenum(year(current_date), 1, 1) + 1;
    
    % Store values in arrays
    MM_array(idx)   = MM;
    doy_array(idx)  = doy;
    dd_array(idx)   = dd;
end

% Convert latitude to radians
phi = deg2rad(lat);

% Calculate inverse relative distance Earth-Sun (dr)
dr = 1 + 0.033 * cos(2 * pi * doy_array / 365);

% Calculate solar declination (delta)
delta = 0.409 * sin(2 * pi * doy_array / 365 - 1.39);

% Calculate sunset hour angle (ws)
ws = acos(-tan(phi) * tan(delta));

% Calculate extraterrestrial radiation (Ra)
Ra = extraterrestrialRadiation(Gsc, dr, pi, phi, ws, delta);

% Display the result
disp(['Extraterrestrial Radiation (Ra): ' num2str(Ra(1)) ' MJ/m^2/day']); % Displaying only the first value

Extraterrestrial Radiation (Ra): 34.8742 MJ/m^2/day

% Constants and parameters
As          = 0.06;
as          = 0.25;
bs          = 0.5;
N           = (2 / 15) .* acosd(-tan(lat) .* tan(rad2deg(delta)));
n           = 11;
Rs          = (as + (bs .* (n ./ N))) .* Ra;

% Calculate net shortwave solar radiation (and photoperiod/24(phi_Light))
pi_light    = N / 24;
Rho_SWR     = Rs * (1 - As);

% Calculate rho_an (the difference between the incident and reflected components of long-wave radiation)
sigma       = 4.896 * 10^-6;
T_a         = 25; % Later read from file
T_ak        = Absolute_temp1(T_a);
r           = 0.5; % Change the value as necessary
C_c         = 0.5; % Read from file
epslon_a    = 0.937 * 10^-5 * T_ak^2 * (1 + 0.17 * C_c^2);
Rho_an      = (1 - r) * epslon_a * sigma * T_ak^4;

% Calculate Water Surface Radiation (HO2_heat_Emm)
T_w         = 25; % To be simulated first
epslon_w    = 0.97;
T_wk        = Absolute_temp2(T_w);
Rho_HO2_heat_Emm = epslon_w * sigma * T_wk^4;

% Calculate Evaporative Heat Loss (Evap_HLoss)
U2          = 1.3; % Read from file
RH          = 0.75; % Read from file
bo          = 368.61;
lambda      = 311.02;
es          = 25.37 * exp(17.62 - 5271 / T_wk);
ea          = RH * 25.37 * exp(17.62 - 5271 / T_ak);
z           = 1000; % Entered in the app
P           = 760 / (10^(z / 19748.2));
T_wv        = (T_wk / (1.0 - (0.378 * es / P)));
T_av        = (T_ak / (1.0 - (0.378 * ea / P)));
Rho_Evap_HLoss = (es - ea) * (lambda * ((T_wv - T_av))^(1 / 3) + bo * U2);

% Calculate Conductive Heat Loss or Gain (Heat_cond)
Rho_Heat_cond = Rho_Evap_HLoss * 0.61 * 10^-3 * P * ((T_wk - T_ak) / (es - ea));

% Calculate interfacial heat transfer due to various processes that occur at the pond (Rho-net)
Rho_net = Rho_SWR + Rho_an - Rho_HO2_heat_Emm - Rho_Evap_HLoss + Rho_Heat_cond;

figure
plot(tspans(1):tspans(end), Rho_net, '.', 'markersize', 3), grid on, 
title('\rho_{net}'), xlabel('Day numder'), xlim(tspans)

% Computation of water balance components
A           = 100;
pond_depth  = 1;
V           = A * pond_depth;
dmin        = 0.8;
dcurr       = 0.6;
dmax        = pond_depth;
Ir          = 20;

if Ir == 0
    Qi = A * (dmin - dcurr) / tspans(2);
else
    Qi = Ir / 100 * V;
end
Er = 10;
if Er == 0
    Qe = A * (dcurr - dmax) / tspans(2);
else
    Qe = Er / 100 * V;
end

pd          = 10; % To be read from a file
PCP         = A * pd / 1000;
Density_w   = 1000;
Latent_vapw = 2260;
Evap        = A * Rho_Evap_HLoss / (Density_w * Latent_vapw);
Sr          = 5;
Seep        = A * Sr / 1000;
CWR         = 20; % To be read from Aquacrop
A_Irr       = 200;
Irr         = A_Irr * CWR / 1000;

% Initial conditions for the ODEs
T_win       = 25; % Initial water temperature
T_wout      = 26;
Cpw         = 5;

% Set initial conditions
initial_conditions = [V; T_w];  % Ensure that this is a column vector

% Solve the system of ODEs using ode45
[t, Y]      = ode45(@(t, Y) dYdt(t, Y, Qi, PCP, Qe, Evap, Seep, Irr, T_win, T_wout, Rho_net, Density_w, Cpw, pond_depth), tspans, initial_conditions);

% Extract results
V_solution  = Y(:, 1);
T_w_solution= Y(:, 2);

% Display or use the solutions as needed
disp('Volume (V) solution:');

Volume (V) solution:

disp(V_solution);

   1.0e+03 *

    0.1000
    0.1021
    0.1042
    0.1063
    0.1084
    0.1189
    0.1295
    0.1400
    0.1505
    0.1755
    0.2005
    0.2255
    0.2505
    0.2996
    0.3487
    0.3978
    0.4468
    0.5170
    0.5871
    0.6573
    0.7274
    0.7975
    0.8677
    0.9378
    1.0079
    1.0781
    1.1482
    1.2184
    1.2885
    1.3586
    1.4288
    1.4989
    1.5690
    1.6392
    1.7093
    1.7795
    1.8496
    1.9197
    1.9899
    2.0600
    2.1302
    2.2003
    2.2704
    2.3406
    2.4107
    2.4808
    2.5510
    2.6211
    2.6913
    2.7448
    2.7984
    2.8520
    2.9055

disp('Water Temperature (T_w) solution:');

Water Temperature (T_w) solution:

disp(T_w_solution);

   25.0000
   23.7557
   22.5336
   21.3325
   20.1512
   14.5029
    9.2112
    4.1946
   -0.6059
  -11.3766
  -21.5055
  -31.2118
  -40.6296
  -58.5450
  -75.9564
  -93.0565
 -109.9559
 -133.8571
 -157.5613
 -181.1334
 -204.6136
 -228.0244
 -251.3831
 -274.7018
 -297.9891
 -321.2511
 -344.4926
 -367.7171
 -390.9275
 -414.1257
 -437.3135
 -460.4920
 -483.6624
 -506.8256
 -529.9821
 -553.1327
 -576.2780
 -599.4184
 -622.5545
 -645.6868
 -668.8158
 -691.9420
 -715.0656
 -738.1873
 -761.3072
 -784.4256
 -807.5427
 -830.6586
 -853.7734
 -871.4261
 -889.0781
 -906.7295
 -924.3801

% Plot the results
figure;
subplot(2, 1, 1);
plot(t, V_solution); grid on
title('Volume (V) vs Time');
xlabel('Time');
ylabel('Volume (V)');

subplot(2, 1, 2);
plot(t, T_w_solution); grid on
title('Water Temperature (T_w) vs Time');
xlabel('Time');
ylabel('Water Temperature (T_w)');

%% Local functions
% function 1
function dYdt = dYdt(t, Y, Qi, PCP, Qe, Evap, Seep, Irr, T_win, T_wout, Rho_net, Density_w, Cpw, pond_depth);
    dYdt    = zeros(2, 1);
    
    % perform interpolations
    Rho_net = interp1(1:365, Rho_net, t);   % time-dependent Rho_net for 365 days
    
    % state equations
    ode1    = Qi + PCP - Qe + Evap + Seep - Irr;
    ode2    = Qi * T_win / Y(1) - Qe * T_wout / Y(1) + Rho_net / (Density_w * Cpw * pond_depth) - Y(2) / Y(1) * ode1;
    
    % full dynamics
    dYdt(1) = ode1;
    dYdt(2) = ode2;
end

% function 2
function Ra = extraterrestrialRadiation(Gsc, dr, pi, phi, ws, delta)
    % Calculate extraterrestrial radiation (Ra)
    Ra = (24 * 60 / pi) * Gsc * dr .* (ws .* sin(phi) .* sin(delta) + cos(phi) .* cos(delta) .* sin(ws));
end

% function 3
function Absolute_1 = Absolute_temp1(T_a)
    Absolute_1 = 273 + T_a;
end

% function 4
function Absolute_2 = Absolute_temp2(T_w)
    Absolute_2 = 273 + T_w;
end

[Torsten]

It's physically impossible that "Evap" remains constant over the year.

[Daniel]

Thank you very much @Sam Chak for your correction and suggestions for improving my code. It's also true as to what @Torsten has mentioned about "Evap" being not constant. I have i addition tried to make Evap time varying by introducing time varying T_a although i have not included the actual values but just an array of thesame value for all time steps. i hope to load the actual data for PCP, T_a, and Irr since they are all time varying as i procees further. Than

[Sam Chak]

I presume that 'Evap' stands for the evaporation rate or something similar. If my assumption is correct, it appears that the model did not capture the dynamics of the evaporation process, resulting in unrealistic results. Nevertheless, it is a positive step forward to have the code now working for time-dependent terms such as 'Rho_net'. If the 'Evap' data is discrete and measured per day throughout the whole year, it is indeed possible to set up interpolation using the same method described earlier.

[Daniel]

I totally agree with you @Sam Chak. i have already updated Evap and will be updating other time dependent variables to make my simulations more realistick. In case of any difficulties i will be sougting for further assistance in this community. Thank You once more.