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

2 次查看（过去 30 天）
Daniel
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);
% 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
  4 个评论
显示 2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
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 个评论
Star Strider
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
.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by