Please help to solve the error in evaluating a system of odes using ode45

The following code keeps on returning the following error.

Error using odearguments

@(T,Y)[ODE1(T,Y);ODE2(T,Y)] returns a vector of length 366, but the length of initial conditions vector is 2. The vector returned by @(T,Y)[ODE1(T,Y);ODE2(T,Y)] 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);

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

% 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 = 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;

% 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;

% 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' / (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)] returns a vector of length 366, but the length of initial conditions vector is 2. The vector returned by @(T,Y)[ODE1(T,Y);ODE2(T,Y)] and the initial conditions vector must have the same number of elements.

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);
title('Volume (V) vs Time');
xlabel('Time');
ylabel('Volume (V)');

subplot(2, 1, 2);
plot(t, T_w_solution);
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]

Rho_net is an array of length 365, but it must be a single scalar value in ode_2.

If you want to evaluate the array at time t given by the integrator, use

interp1(time_in_Rho_net,Rho_net,t)

where the array "time_in_Rho_net" is of the same size as "Rho_net" and are the time points at which Rho_net is given.

[Sam Chak]

In my opinion, the simulated results obtained by interpolating the day-dependent points should show a close resemblance to the system response when using the average value of .

% 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

% 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) + mean(Rho_net) / (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)];

% Solve the system of ODEs using ode45
[t, Y]      = ode45(dYdt, 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.5691
    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.6060
  -11.3768
  -21.5058
  -31.2123
  -40.6303
  -58.5463
  -75.9583
  -93.0592
 -109.9595
 -133.8622
 -157.5679
 -181.1414
 -204.6231
 -228.0350
 -251.3943
 -274.7130
 -297.9998
 -321.2607
 -344.5004
 -367.7227
 -390.9303
 -414.1256
 -437.3104
 -460.4861
 -483.6540
 -506.8150
 -529.9700
 -553.1198
 -576.2649
 -599.4058
 -622.5430
 -645.6768
 -668.8076
 -691.9357
 -715.0613
 -738.1846
 -761.3059
 -784.4253
 -807.5430
 -830.6591
 -853.7738
 -871.4259
 -889.0773
 -906.7280
 -924.3780

% 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 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

[Daniel]

Thank you very much for your response @Torsten. It is really helping.

[Daniel]

Kindly elaborate on what you meant by this statements@Torsten. Thank you.

"If you want to evaluate the array at time t given by the integrator, use

interp1(time_in_Rho_net,Rho_net,t)"

[Torsten]

Look at the example

"ODE with Time-Dependent Terms"

under

https://uk.mathworks.com/help/matlab/ref/ode45.html

[Daniel]

Thank you for your guidance.

[Torsten]

Did you understand ? You must insert the value for Rho_net in your equation that is valid at time t. You usually get this value by interpolating your array values to the time t that the solver asks for.

[Daniel]

However my challange is that for the integration of the following two ODEs, some of the variables are known at different time intervals of the period of integration. These variales include preipitation (PCP), Irrigation amount(Irr) and Rho_net. How do I provide/input these values to the intergation prrocess at the specified time intervals so that the integartion puts into consideration these knowns while approximating the unknowns (water temperature and pond volume) at the various time steps . ode1 = @(t, Y) Qi + PCP - Qe + Evap + Seep - Irr;

ode2 = @(t, Y) Qi * T_win / Y(1) - Qe * T_wout / Y(1) + Rho_net / (Density_w * Cpw * pond_depth) - Y(2) / Y(1) * ode1(t, Y);

[Torsten]

My answer to interpolate the variables for which you have time-dependent arrays of values to the time t where the integrator wants their respective values remains the same.

I'm sorry that you can't understand how to do it from the answer given in this link:

Example

"ODE with Time-Dependent Terms"

under

https://uk.mathworks.com/help/matlab/ref/ode45.html