when using ode45 with the for loop, should I put the initial conditions inside or outside the for loop?

Where should I placde the initial conditions for the ODE when using ode45 and the for loop to solve my system of ODEs? I have noticed i am getting different solutions when i either place inside or outside the for for loop. More so, is it okay to state the next set of initial conditions as I have stated in my code when using the for loop(i.e. Initial_conditions10=Y(t+1,:); ) such that the next time step used the solution of the previous time step as the initial conditions? Kindly assist.

clc; clear; close all;

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

% Define simulation period
startdate = datetime(2023, 01, 01, 'Format', 'uuuu-MM-dd');
finishdate = datetime(2023, 12, 31, '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

Initial_T_w = 27.5;
A = 100;
pond_depth = 1;
V = A * pond_depth;
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
n_data = 11*ones(length(tspans),1); %to be provided as a data series
pd_data = 5*ones(length(tspans),1); % precipitation depth; To be read from a 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);
    pd = pd_data(t);
    n = n_data(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 = (24/pi)*ws;
    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)
    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;%A constant with the units368.61 kJ m^2 d^-1·'mmHg^-1(m s^-1)^-1
    lambda  = 311.02;%A constant with the units kJ m^2 d^-1mmHg^-1K^-1/3
    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
    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
    
    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(1,1:2) =Initial_conditions10;%this prevents the first row entry from being zero
    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(V_solution); grid on
title('Volume (V) vs Time');
xlabel('Time');
ylabel('Volume (V)');
subplot(2, 1, 2);
plot(T_w_solution); grid on
title('Water Temperature (T_w) vs Time');
xlabel('Time');
ylabel('Water Temperature (T_w)');

%% Local function in the script
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

[Torsten]

dVdt = @(t, Y) Qi + PCP - Qe + Evap + Seep - Irr;

So dVdt is a constant ? Then

V = V0 + (Qi + PCP - Qe + Evap + Seep - Irr)*(t-t0)

No need to integrate an ODE for V.

[Daniel]

Thank you for the suggestion. However, dVdt isn't s constant as all the terms are time dependend and evap too varies with temperature.

[Torsten]

Don't use this loop over the elements of "tspans".

Take a look at the example

ODE with Time-Dependent Terms

under

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

to see how to proceed in your case.

[Daniel]

Thank you @Torsten. I will check out the example

Answer 1

[Star Strider]

It is difficult to understand what your code does. If you are changing the parameters of the differential equations inside the loop (and producing step changes as the result), the end result of the previous integration should be the initial conditions for the next integration step.

[Daniel]

Thank you for your enligtenment. the terms of the odes are changing at everytime step inside the loop.

[Star Strider]

My pleasure!

I cannot run anything here today, however this is the sort of approach I was thinking of —

odefcn = @(t,y,k) [y(1); exp(y(1) + (-1)^k)];
ic = [0 1];

for k = 0:5;
    [t,y] = ode45(@(t,y)odefcn(t,y,k), [0 10]+10*k, ic);
    ic = y(end,:);
    tc{k+1,:} = t;
    yc{k+1,:} = y;
end

tv = cell2mat(tc);
ym = cell2mat(yc);

figure
plot(tv, ym)
grid

.