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

3 ビュー (過去 30 日間)
Daniel 2024 年 3 月 31 日
コメント済み: Daniel 2024 年 4 月 3 日
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);
% 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));
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;
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)
Ra = (24 * 60 / pi) * Gsc * dr .* (ws .* sin(phi) .* sin(delta) + cos(phi) .* cos(delta) .* sin(ws));
end
##### 4 件のコメント2 件の古いコメントを表示2 件の古いコメントを非表示
Torsten 2024 年 4 月 1 日
Don't use this loop over the elements of "tspans".
Take a look at the example
ODE with Time-Dependent Terms
under
to see how to proceed in your case.
Daniel 2024 年 4 月 3 日
Thank you @Torsten. I will check out the example

サインインしてコメントする。

### 採用された回答

Star Strider 2024 年 3 月 31 日
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.
##### 2 件のコメントなしを表示なしを非表示
Daniel 2024 年 4 月 1 日
Thank you for your enligtenment. the terms of the odes are changing at everytime step inside the loop.
Star Strider 2024 年 4 月 1 日
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
.

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

R2023a

### Community Treasure Hunt

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

Start Hunting!

Translated by