Why does my ode45 return a solution with the first values as zero despite I supplying non-zero initial conditions?

4 views (last 30 days)
Daniel
Daniel on 27 Mar 2024
Answered: Thomas Males on 16 May 2024
Ran in:
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
Extraterrestrial Radiation (Ra): 32589.928 kJ/m^2/day
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:');
Volume (V) solution:
disp(V_solution);
0 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982
disp('Water Temperature (T_w) solution:');
Water Temperature (T_w) solution:
disp(T_w_solution);
0 27.2105 27.2138 27.2173 27.2210 27.2250 27.2291 27.2334 27.2379 27.2426 27.2475 27.2525 27.2577 27.2630 27.2686 27.2742 27.2800 27.2859 27.2919 27.2980 27.3043 27.3106 27.3171 27.3236 27.3301 27.3368 27.3434 27.3502 27.3569 27.3637 27.3705 27.3773 27.3841 27.3908 27.3976 27.4043 27.4109 27.4176 27.4241 27.4306 27.4370 27.4433 27.4495 27.4555 27.4615 27.4673 27.4730 27.4786 27.4839 27.4892 27.4942 27.4991 27.5037 27.5082 27.5125 27.5165 27.5203 27.5239 27.5273 27.5304 27.5333 27.5359 27.5383 27.5404 27.5422 27.5438 27.5451 27.5461 27.5468 27.5472 27.5473 27.5472 27.5467 27.5460 27.5449 27.5435 27.5419 27.5399 27.5376 27.5351 27.5322 27.5290 27.5255 27.5218 27.5177 27.5133 27.5087 27.5037 27.4985 27.4930 27.4872 27.4811 27.4748 27.4682 27.4614 27.4543 27.4470 27.4394 27.4316 27.4236 27.4153 27.4069 27.3982 27.3893 27.3803 27.3711 27.3617 27.3522 27.3425 27.3326 27.3227 27.3126 27.3024 27.2920 27.2816 27.2711 27.2606 27.2499 27.2393 27.2285 27.2178 27.2070 27.1962 27.1854 27.1746 27.1638 27.1531 27.1424 27.1317 27.1211 27.1106 27.1001 27.0898 27.0795 27.0694 27.0593 27.0494 27.0397 27.0301 27.0206 27.0113 27.0022 26.9932 26.9844 26.9759 26.9675 26.9594 26.9515 26.9438 26.9363 26.9291 26.9222 26.9155 26.9090 26.9028 26.8969 26.8913 26.8859 26.8809 26.8761 26.8716 26.8675 26.8636 26.8600 26.8568 26.8538 26.8512 26.8489 26.8469 26.8452 26.8439 26.8429 26.8422 26.8418 26.8417 26.8420 26.8426 26.8435 26.8447 26.8462 26.8481
% 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

Sign in to answer this question.

Accepted Answer

Dyuman Joshi
Dyuman Joshi on 27 Mar 2024
Edited: Dyuman Joshi on 27 Mar 2024
Ran in:
"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
Extraterrestrial Radiation (Ra): 32589.928 kJ/m^2/day
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:');
Volume (V) solution:
disp(V_solution);
100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 100.7982 0
disp('Water Temperature (T_w) solution:');
Water Temperature (T_w) solution:
disp(T_w_solution);
27.2105 27.2138 27.2173 27.2210 27.2250 27.2291 27.2334 27.2379 27.2426 27.2475 27.2525 27.2577 27.2630 27.2686 27.2742 27.2800 27.2859 27.2919 27.2980 27.3043 27.3106 27.3171 27.3236 27.3301 27.3368 27.3434 27.3502 27.3569 27.3637 27.3705 27.3773 27.3841 27.3908 27.3976 27.4043 27.4109 27.4176 27.4241 27.4306 27.4370 27.4433 27.4495 27.4555 27.4615 27.4673 27.4730 27.4786 27.4839 27.4892 27.4942 27.4991 27.5037 27.5082 27.5125 27.5165 27.5203 27.5239 27.5273 27.5304 27.5333 27.5359 27.5383 27.5404 27.5422 27.5438 27.5451 27.5461 27.5468 27.5472 27.5473 27.5472 27.5467 27.5460 27.5449 27.5435 27.5419 27.5399 27.5376 27.5351 27.5322 27.5290 27.5255 27.5218 27.5177 27.5133 27.5087 27.5037 27.4985 27.4930 27.4872 27.4811 27.4748 27.4682 27.4614 27.4543 27.4470 27.4394 27.4316 27.4236 27.4153 27.4069 27.3982 27.3893 27.3803 27.3711 27.3617 27.3522 27.3425 27.3326 27.3227 27.3126 27.3024 27.2920 27.2816 27.2711 27.2606 27.2499 27.2393 27.2285 27.2178 27.2070 27.1962 27.1854 27.1746 27.1638 27.1531 27.1424 27.1317 27.1211 27.1106 27.1001 27.0898 27.0795 27.0694 27.0593 27.0494 27.0397 27.0301 27.0206 27.0113 27.0022 26.9932 26.9844 26.9759 26.9675 26.9594 26.9515 26.9438 26.9363 26.9291 26.9222 26.9155 26.9090 26.9028 26.8969 26.8913 26.8859 26.8809 26.8761 26.8716 26.8675 26.8636 26.8600 26.8568 26.8538 26.8512 26.8489 26.8469 26.8452 26.8439 26.8429 26.8422 26.8418 26.8417 26.8420 26.8426 26.8435 26.8447 26.8462 26.8481 0
% 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

More Answers (2)

Saurabh
Saurabh on 27 Mar 2024
Edited: Saurabh on 27 Mar 2024
Hi Daniel,
It seems like you were expecting to have a non-zero initial value, but it appears you received a zero initial value instead.
Upon very closely analyzing your code, I find that the value of ‘tspans’ is 181 and the size of matrix ‘Y’ is also 181.
In the looping part, the variable ‘t’ is ranging from 1 to 180. So, the maximum time matrix ‘Y’ is getting updated is only 180 times.
If you change this part of the code,
Y (t+1, :) = Y_integrated(end, :);
Initial_conditions10=Y (t+1, :);
Rho_net_values(t+1) = Rho_net;
With this:
Y (t, :) = Y_integrated(end, :);
Initial_conditions10=Y (t, :);
Rho_net_values(t) = Rho_net;
The last value of the vector will hold a zero value.
I hope this was helpful.
  2 Comments
Daniel
Daniel on 27 Mar 2024
Its true but my wish the code could execute in such a manner such that both the first and the last values are not zeros.
Saurabh
Saurabh on 28 Mar 2024
In that case, simply reduce the size of the matrix 'Y' by 1. This will exclude the possibility of zero in the code.

Sign in to comment.

Thomas Males
Thomas Males on 16 May 2024
Hi Daniel,
I would suggest looking into timetables for processing of timeseries data. This would likely allow you to get rid of the for loops and calulate the new variables with dot notation like this:
myTable.NewVar = mytable.OldVar1 .* myTable.Oldvar2 + Const1.
This also opens up a lot of funtionality through the retime function and simplifies plotting of multiple variables through the stacked plot function.
You've probably got everything working with your loop by now but it could be a consideration for the future.
Cheers,
Tom

Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

Tags

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!