Why does my ode45 return a solution with the first values as zero despite I supplying non-zero initial conditions?
In the following code to simulate Pond volume and water temperature changes the code returns the solution (Y) with the first values of pondvolume and water temperature as 0,0. What could be the problem and ho can I solve it?
% 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);
Y(t+1,:) = Y_integrated(end,:);
Initial_conditions10=Y(t+1,:);
Rho_net_values(t+1) = 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));
endIn the following code to simulate Pond volume and water temperature changes the code returns the solution (Y) with the first values of pondvolume and water temperature as 0,0. What could be the problem and ho can I solve it?
% 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);
Y(t+1,:) = Y_integrated(end,:);
Initial_conditions10=Y(t+1,:);
Rho_net_values(t+1) = 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 In the following code to simulate Pond volume and water temperature changes the code returns the solution (Y) with the first values of pondvolume and water temperature as 0,0. What could be the problem and ho can I solve it?
% 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);
Y(t+1,:) = Y_integrated(end,:);
Initial_conditions10=Y(t+1,:);
Rho_net_values(t+1) = 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 ode45 MATLAB Answers — New Questions