"Why does my ode45 return a solution with the first values as zero despite I supplying non-zero initial conditions?"
Because you loop variable t is set from 1 to length(tspans)-1, but the indexing (for variables Y and Rho_net_values) and is done using t+1. And as you have preallocated the array using zeros(), the 1st value is not modified and remains as it is.
To resolve this issue change the index from t+1 to t -
% Solar constant [MJ m-2 min-1]
Gsc = 0.08;
lat = 0; % Latitude
% Define simulation period
startdate = datetime(2023, 01, 01, 'Format', 'uuuu-MM-dd');
finishdate = datetime(2023, 06, 30, 'Format', 'uuuu-MM-dd');
tspans = (1:1:days(finishdate - startdate) + 1);
% Initialize arrays to store MM, doy, and dd values
MM_array_s = zeros(1, days(finishdate - startdate) + 1);
doy_array_s = zeros(1, days(finishdate - startdate) + 1);
dd_array_s = 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_s(idx) = MM;
doy_array_s(idx) = doy;
dd_array_s(idx) = dd;
end
CWR = 20; % To be read from Aquacrop
A_Irr = 200;
Irr_data = A_Irr * CWR / 1000*ones(length(tspans),1);
T_a_data = 25 *ones(length(tspans),1); % Later read from file/ it should be a data series
U2_data = 1.3*ones(length(tspans),1); % Read from file
RH_data = 0.75*ones(length(tspans),1); % Read from file
Y= zeros(length(tspans), 2);
Rho_net_values =zeros(length(tspans), 1);
for t=1:length(tspans)-1
T_a =T_a_data(t);
U2 = U2_data(t);
RH = RH_data(t);
Irr = Irr_data(t);
doy_array = doy_array_s(t);
% 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 ./ 365)* doy_array) - 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
% Constants and parameters
As = 0.06;
as = 0.25;
bs = 0.5;
%N = (2 / 15) .* acos(-tand(lat) .* tand(rad2deg(delta)));
N = (24/pi)*ws;
n = 11; %to be provided as a data series
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);
r = 0.36; % Change the value as necessary
C_c = 0.5; % Read from file/enter
epslon_a = 0.937 * 10^-5 * (273+T_a).^2 * (1 + 0.17 * C_c^2);
Rho_an = (1 - r) * epslon_a * sigma .* (273+T_a).^4;
% Calculate Water Surface Radiation (HO2_heat_Emm)
Initial_T_w = 27;
T_w = Initial_T_w; % To be simulated first
epslon_w = 0.97;
Rho_HO2_heat_Emm = epslon_w * sigma * (273+T_w).^4;
% Calculate Evaporative Heat Loss (Evap_HLoss)
bo = 368.61;
lambda = 311.02;
es = 25.37 * exp(17.62 - 5271 ./ (273+T_w));
ea = RH .* 25.37 .* exp(17.62 - 5271 ./ (273+T_a));
z = 1000; % Entered in the app
P = 760 / (10^(z / 19748.2));
T_wv = ((273+T_w) / (1.0 - (0.378 * es / P)));
T_av = ((273+T_a) / (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 * (((273+T_w) - (273+T_a)) / (es - ea));
Rho_net= Rho_SWR + Rho_an - Rho_HO2_heat_Emm - Rho_Evap_HLoss + Rho_Heat_cond;
% Define Rho_net as an anonymous function of time
%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;
% Initial conditions for the ODEs
T_win = 27; % Initial water temperature
T_wout = 27;
Cpw = 4.18;
% Define your ODEs
dVdt = @(t, Y) Qi + PCP - Qe + Evap + Seep - Irr;
dTdt = @(t, Y) Qi * T_win / Y(1) - Qe * T_wout / Y(1) + Rho_net / (Density_w * Cpw * pond_depth) - Y(2) / Y(1) * dVdt(t,Y);
% Combine the ODEs into a single function
dYdt = @(t, Y) [dVdt(t, Y); dTdt(t, Y)];
% Set initial conditions
Initial_conditions10 = [V; Initial_T_w]; % Ensure that this is a column vector
% Solve the system of ODEs using ode45
options = odeset('AbsTol', 1e-6, 'RelTol', 1e-6);
[t_integrated10,Y_integrated]=ode45(@(t,Y) dYdt(t,Y),tspans(t:t+1),Initial_conditions10);
%%v
Y(t,:) = Y_integrated(end,:);
Initial_conditions10=Y(t+1,:);
%% v
Rho_net_values(t) = Rho_net;
end
disp(['Extraterrestrial Radiation (Ra): ' num2str(Ra(1)) ' kJ/m^2/day']); % Displaying only the first value
figure
plot(tspans(1):tspans(end), Rho_net_values(tspans(1):tspans(end)), '.', 'markersize', 3), grid on,
title('\rho_{net}')
% 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(tspans, V_solution); grid on
title('Volume (V) vs Time');
xlabel('Time');
ylabel('Volume (V)');
subplot(2, 1, 2);
plot(tspans, 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
