Please help to solve the error in evaluating a system of odes using ode45

3 次查看（过去 30 天）
Daniel
Daniel 2024-3-7
评论： Torsten 2024-3-10
Ran in:
The following code keeps on returning the following error.
Error using odearguments
@(T,Y)[ODE1(T,Y);ODE2(T,Y)] returns a vector of length 366, but the length of initial conditions vector is 2. The vector returned by @(T,Y)[ODE1(T,Y);ODE2(T,Y)] and the initial conditions vector must have the same number of elements.
Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
% Solar constant [MJ m-2 min-1]
Gsc = 0.08;
lat = 0; % Latitude
% Define simulation period
startdate = datetime(2025, 01, 01, 'Format', 'uuuu-MM-dd');
finishdate = datetime(2025, 12, 31, 'Format', 'uuuu-MM-dd');
tspans = [1 days(finishdate - startdate)+ 1];
% Initialize arrays to store MM, doy, and dd values
MM_array = zeros(1, days(finishdate - startdate) + 1);
doy_array = zeros(1, days(finishdate - startdate) + 1);
dd_array = 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 = datenum(current_date) - datenum(year(current_date), 1, 1) + 1;
% Store values in arrays
MM_array(idx) = MM;
doy_array(idx) = doy;
dd_array(idx) = dd;
end
% 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 * doy_array / 365 - 1.39);
% Calculate sunset hour angle (ws)
ws = acos(-tan(phi) * tan(delta));
% Calculate extraterrestrial radiation (Ra)
Ra = extraterrestrialRadiation(Gsc, dr, pi, phi, ws, delta);
% Display the result
disp(['Extraterrestrial Radiation (Ra): ' num2str(Ra(1)) ' MJ/m^2/day']); % Displaying only the first value
Extraterrestrial Radiation (Ra): 34.8742 MJ/m^2/day
% -------------------------------------
% Constants and parameters
As = 0.06;
as = 0.25;
bs = 0.5;
N = (2 / 15) .* acosd(-tan(lat) .* tan(rad2deg(delta)));
n = 11;
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;
T_a = 25; % Later read from file
T_ak = Absolute_temp1(T_a);
r = 0.5; % Change the value as necessary
C_c = 0.5; % Read from file
epslon_a = 0.937 * 10^-5 * T_ak^2 * (1 + 0.17 * C_c^2);
Rho_an = (1 - r) * epslon_a * sigma * T_ak^4;
% Calculate Water Surface Radiation (HO2_heat_Emm)
T_w = 25; % To be simulated first
epslon_w = 0.97;
T_wk = Absolute_temp2(T_w);
Rho_HO2_heat_Emm = epslon_w * sigma * T_wk^4;
% Calculate Evaporative Heat Loss (Evap_HLoss)
U2 = 1.3; % Read from file
RH = 0.75; % Read from file
bo = 368.61;
lambda = 311.02;
es = 25.37 * exp(17.62 - 5271 / T_wk);
ea = RH * 25.37 * exp(17.62 - 5271 / T_ak);
z = 1000; % Entered in the app
P = 760 / (10^(z / 19748.2));
T_wv = (T_wk / (1.0 - (0.378 * es / P)));
T_av = (T_ak / (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 * ((T_wk - T_ak) / (es - ea));
% Calculate interfacial heat transfer due to various processes that occur at the pond (Rho-net)
Rho_net = Rho_SWR + Rho_an - Rho_HO2_heat_Emm - Rho_Evap_HLoss + Rho_Heat_cond;
% Computation of water balance components
A = 100;
pond_depth = 1;
V = A * pond_depth;
dmin = 0.8;
dcurr = 0.6;
dmax = pond_depth;
Ir = 20;
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;
CWR = 20; % To be read from Aquacrop
A_Irr = 200;
Irr = A_Irr * CWR / 1000;
% Initial conditions for the ODEs
T_win = 25; % Initial water temperature
T_wout = 26;
Cpw = 5;
% Define your ODEs
ode1 = @(t, Y) Qi + PCP - Qe + Evap + Seep - Irr;
ode2 = @(t, Y) Qi * T_win / Y(1) - Qe * T_wout / Y(1) + Rho_net' / (Density_w * Cpw * pond_depth) - Y(2) / Y(1) * ode1(t, Y);
% Combine the ODEs into a single function
dYdt = @(t, Y) [ode1(t, Y); ode2(t, Y)];
% Set initial conditions
initial_conditions = [V; T_w]; % Ensure that this is a column vector
% Solve the system of ODEs using ode45
[t, Y] = ode45(dYdt, tspans, initial_conditions);
Error using odearguments
@(T,Y)[ODE1(T,Y);ODE2(T,Y)] returns a vector of length 366, but the length of initial conditions vector is 2. The vector returned by @(T,Y)[ODE1(T,Y);ODE2(T,Y)] and the initial conditions vector must have the same number of elements.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
% 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(t, V_solution);
title('Volume (V) vs Time');
xlabel('Time');
ylabel('Volume (V)');
subplot(2, 1, 2);
plot(t, T_w_solution);
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
function Absolute_1 = Absolute_temp1(T_a)
Absolute_1 = 273 + T_a;
end
function Absolute_2 = Absolute_temp2(T_w)
Absolute_2 = 273 + T_w;
end

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

采纳的回答

Torsten
Torsten 2024-3-7
编辑：Torsten 2024-3-7
Rho_net is an array of length 365, but it must be a single scalar value in ode_2.
If you want to evaluate the array at time t given by the integrator, use
interp1(time_in_Rho_net,Rho_net,t)
where the array "time_in_Rho_net" is of the same size as "Rho_net" and are the time points at which Rho_net is given.
  8 个评论
显示 6更早的评论
Daniel
Daniel 2024-3-10
编辑：Daniel 2024-3-10
However my challange is that for the integration of the following two ODEs, some of the variables are known at different time intervals of the period of integration. These variales include preipitation (PCP), Irrigation amount(Irr) and Rho_net. How do I provide/input these values to the intergation prrocess at the specified time intervals so that the integartion puts into consideration these knowns while approximating the unknowns (water temperature and pond volume) at the various time steps . ode1 = @(t, Y) Qi + PCP - Qe + Evap + Seep - Irr;
ode2 = @(t, Y) Qi * T_win / Y(1) - Qe * T_wout / Y(1) + Rho_net / (Density_w * Cpw * pond_depth) - Y(2) / Y(1) * ode1(t, Y);
Torsten
Torsten 2024-3-10
My answer to interpolate the variables for which you have time-dependent arrays of values to the time t where the integrator wants their respective values remains the same.
I'm sorry that you can't understand how to do it from the answer given in this link:
Example
"ODE with Time-Dependent Terms"
under
https://uk.mathworks.com/help/matlab/ref/ode45.html

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by