Below I tried to solve Population Growth problem, I provide you the mathematical solution. Also, I create MATLAB code, however i feel like that i repeat this code to each part of this code as you can see. How could I improve this code?

Assume a continuous growth model where the function that gives the rate of change is defined as . Provide the function , which represents the population size.
If the growth of a population follows the logistic law but there is also a term for population reduction (e.g., from a predator while the members of the population are the prey), the equation that expresses the population size is of the form
where the term is given here to be and it represents a term for population reduction. The constants are positive. Normalize the equation by setting and .Find the steady solutions of the problem (real and positive) and then implement a program in MATLAB to numerically solve the dimensionless form of the equation. Make the graphical representation of the numerical solution for , with initial conditions , and respectively. What do you observe?

So we have

We are going to find the function :
To find , we integrate :

where is the integration constant, which can be determined if an initial condition is given, such as .
Now we are going to normalize the equation:
Let and . Then:

We substitute into the original equation:

Then we are going to calculate the Stationary Points:
Set :
]

For :

% Define the normalized logistic equation with reduction term
% K/a = 35, ar/b = 0.4
Ka = 35;
ar_b = 0.4;

% Define the function for the normalized logistic equation
f = @(tau, N) ar_b * N * (1 – N/Ka) – N / (1 + N);

% Time span for the simulation
tspan = [0 50]; % Adjust if needed to observe long-term behavior

% Initial conditions
N0 = [40, 25, 2, 0.05];

% Prepare the figure for plotting
figure(1);
hold on;

% Colors for different trajectories
colors = [‘r’, ‘g’, ‘b’, ‘k’]; % Red, Green, Blue, Black

% Solve the equation for each initial condition and plot
for i = 1:length(N0)
[T, N] = ode45(f, tspan, N0(i));
plot(T, N, ‘Color’, colors(i), ‘LineWidth’, 2);
legendInfo{i} = [‘N_0 = ‘ num2str(N0(i))]; % Create legend entry
end

% Add graph details
title(‘Numerical solutions of the normalized logistic equation’);
xlabel(‘Normalized time τ’);
ylabel(‘Normalized population size N’);
legend(legendInfo, ‘Location’, ‘northeastoutside’);
grid on;
hold off;

% Equilibrium analysis
syms N
eqn = ar_b * N * (1 – N / Ka) – N / (1 + N) == 0;
equilibria = double(solve(eqn, N));

% Stability analysis
f = @(N) ar_b * N * (1 – N / Ka) – N / (1 + N);
J = diff(f(N), N);
eigenvalues = double(subs(J, N, equilibria));

disp(‘Equilibrium Points and their Stability:’);
for i = 1:length(equilibria)
if real(eigenvalues(i)) < 0
stability = ‘Stable’;
else
stability = ‘Unstable’;
end
fprintf(‘N = %f: Eigenvalue = %f, %sn’, equilibria(i), eigenvalues(i), stability);
end

% Bifurcation analysis setup
ar_b_values = linspace(0.1, 1, 100);
bifurcation_results = zeros(length(ar_b_values), 10); % Adjust the columns based on expected number of roots

% Compute equilibrium points for varying ar/b
for i = 1:length(ar_b_values)
eqn = ar_b_values(i) * N * (1 – N / Ka) – N / (1 + N) == 0;
solutions = double(solve(eqn, N));
bifurcation_results(i, 1:length(solutions)) = solutions’;
end

% Plot bifurcation diagram
figure;
plot(ar_b_values, bifurcation_results, ‘b*’);
title(‘Bifurcation Diagram’);
xlabel(‘ar/b’);
ylabel(‘Equilibrium Points’);

% Set up the ODE function for phase plane analysis
f = @(N) ar_b * 0.4 * N .* (1 – N / Ka) – N ./ (1 + N);

% Create a grid of N values and compute their time derivatives
N_values = linspace(0, 50, 400);
dN_dtau = f(N_values);

% Plot phase plane
figure;
plot(N_values, dN_dtau);
title(‘Phase Plane Analysis’);
xlabel(‘Population Size N’);
ylabel(‘dN/dτ’);
grid on;
hold on;

% Mark zero crossings (equilibrium points)
zero_crossings = N_values(abs(dN_dtau) < 1e-3);
for i = 1:length(zero_crossings)
plot(zero_crossings(i), 0, ‘ro’);
end
hold off;Below I tried to solve Population Growth problem, I provide you the mathematical solution. Also, I create MATLAB code, however i feel like that i repeat this code to each part of this code as you can see. How could I improve this code?

Assume a continuous growth model where the function that gives the rate of change is defined as . Provide the function , which represents the population size.
If the growth of a population follows the logistic law but there is also a term for population reduction (e.g., from a predator while the members of the population are the prey), the equation that expresses the population size is of the form
where the term is given here to be and it represents a term for population reduction. The constants are positive. Normalize the equation by setting and .Find the steady solutions of the problem (real and positive) and then implement a program in MATLAB to numerically solve the dimensionless form of the equation. Make the graphical representation of the numerical solution for , with initial conditions , and respectively. What do you observe?

So we have

We are going to find the function :
To find , we integrate :

where is the integration constant, which can be determined if an initial condition is given, such as .
Now we are going to normalize the equation:
Let and . Then:

We substitute into the original equation:

Then we are going to calculate the Stationary Points:
Set :
]

For :

% Define the normalized logistic equation with reduction term
% K/a = 35, ar/b = 0.4
Ka = 35;
ar_b = 0.4;

% Define the function for the normalized logistic equation
f = @(tau, N) ar_b * N * (1 – N/Ka) – N / (1 + N);

% Time span for the simulation
tspan = [0 50]; % Adjust if needed to observe long-term behavior

% Initial conditions
N0 = [40, 25, 2, 0.05];

% Prepare the figure for plotting
figure(1);
hold on;

% Colors for different trajectories
colors = [‘r’, ‘g’, ‘b’, ‘k’]; % Red, Green, Blue, Black

% Solve the equation for each initial condition and plot
for i = 1:length(N0)
[T, N] = ode45(f, tspan, N0(i));
plot(T, N, ‘Color’, colors(i), ‘LineWidth’, 2);
legendInfo{i} = [‘N_0 = ‘ num2str(N0(i))]; % Create legend entry
end

% Add graph details
title(‘Numerical solutions of the normalized logistic equation’);
xlabel(‘Normalized time τ’);
ylabel(‘Normalized population size N’);
legend(legendInfo, ‘Location’, ‘northeastoutside’);
grid on;
hold off;

% Equilibrium analysis
syms N
eqn = ar_b * N * (1 – N / Ka) – N / (1 + N) == 0;
equilibria = double(solve(eqn, N));

% Stability analysis
f = @(N) ar_b * N * (1 – N / Ka) – N / (1 + N);
J = diff(f(N), N);
eigenvalues = double(subs(J, N, equilibria));

disp(‘Equilibrium Points and their Stability:’);
for i = 1:length(equilibria)
if real(eigenvalues(i)) < 0
stability = ‘Stable’;
else
stability = ‘Unstable’;
end
fprintf(‘N = %f: Eigenvalue = %f, %sn’, equilibria(i), eigenvalues(i), stability);
end

% Bifurcation analysis setup
ar_b_values = linspace(0.1, 1, 100);
bifurcation_results = zeros(length(ar_b_values), 10); % Adjust the columns based on expected number of roots

% Compute equilibrium points for varying ar/b
for i = 1:length(ar_b_values)
eqn = ar_b_values(i) * N * (1 – N / Ka) – N / (1 + N) == 0;
solutions = double(solve(eqn, N));
bifurcation_results(i, 1:length(solutions)) = solutions’;
end

% Plot bifurcation diagram
figure;
plot(ar_b_values, bifurcation_results, ‘b*’);
title(‘Bifurcation Diagram’);
xlabel(‘ar/b’);
ylabel(‘Equilibrium Points’);

% Set up the ODE function for phase plane analysis
f = @(N) ar_b * 0.4 * N .* (1 – N / Ka) – N ./ (1 + N);

% Create a grid of N values and compute their time derivatives
N_values = linspace(0, 50, 400);
dN_dtau = f(N_values);

% Plot phase plane
figure;
plot(N_values, dN_dtau);
title(‘Phase Plane Analysis’);
xlabel(‘Population Size N’);
ylabel(‘dN/dτ’);
grid on;
hold on;

% Mark zero crossings (equilibrium points)
zero_crossings = N_values(abs(dN_dtau) < 1e-3);
for i = 1:length(zero_crossings)
plot(zero_crossings(i), 0, ‘ro’);
end
hold off; Below I tried to solve Population Growth problem, I provide you the mathematical solution. Also, I create MATLAB code, however i feel like that i repeat this code to each part of this code as you can see. How could I improve this code?

Assume a continuous growth model where the function that gives the rate of change is defined as . Provide the function , which represents the population size.
If the growth of a population follows the logistic law but there is also a term for population reduction (e.g., from a predator while the members of the population are the prey), the equation that expresses the population size is of the form
where the term is given here to be and it represents a term for population reduction. The constants are positive. Normalize the equation by setting and .Find the steady solutions of the problem (real and positive) and then implement a program in MATLAB to numerically solve the dimensionless form of the equation. Make the graphical representation of the numerical solution for , with initial conditions , and respectively. What do you observe?

So we have

We are going to find the function :
To find , we integrate :

where is the integration constant, which can be determined if an initial condition is given, such as .
Now we are going to normalize the equation:
Let and . Then:

We substitute into the original equation:

Then we are going to calculate the Stationary Points:
Set :
]

For :

% Define the normalized logistic equation with reduction term
% K/a = 35, ar/b = 0.4
Ka = 35;
ar_b = 0.4;

% Define the function for the normalized logistic equation
f = @(tau, N) ar_b * N * (1 – N/Ka) – N / (1 + N);

% Time span for the simulation
tspan = [0 50]; % Adjust if needed to observe long-term behavior

% Initial conditions
N0 = [40, 25, 2, 0.05];

% Prepare the figure for plotting
figure(1);
hold on;

% Colors for different trajectories
colors = [‘r’, ‘g’, ‘b’, ‘k’]; % Red, Green, Blue, Black

% Solve the equation for each initial condition and plot
for i = 1:length(N0)
[T, N] = ode45(f, tspan, N0(i));
plot(T, N, ‘Color’, colors(i), ‘LineWidth’, 2);
legendInfo{i} = [‘N_0 = ‘ num2str(N0(i))]; % Create legend entry
end

% Add graph details
title(‘Numerical solutions of the normalized logistic equation’);
xlabel(‘Normalized time τ’);
ylabel(‘Normalized population size N’);
legend(legendInfo, ‘Location’, ‘northeastoutside’);
grid on;
hold off;

% Equilibrium analysis
syms N
eqn = ar_b * N * (1 – N / Ka) – N / (1 + N) == 0;
equilibria = double(solve(eqn, N));

% Stability analysis
f = @(N) ar_b * N * (1 – N / Ka) – N / (1 + N);
J = diff(f(N), N);
eigenvalues = double(subs(J, N, equilibria));

disp(‘Equilibrium Points and their Stability:’);
for i = 1:length(equilibria)
if real(eigenvalues(i)) < 0
stability = ‘Stable’;
else
stability = ‘Unstable’;
end
fprintf(‘N = %f: Eigenvalue = %f, %sn’, equilibria(i), eigenvalues(i), stability);
end

% Bifurcation analysis setup
ar_b_values = linspace(0.1, 1, 100);
bifurcation_results = zeros(length(ar_b_values), 10); % Adjust the columns based on expected number of roots

% Compute equilibrium points for varying ar/b
for i = 1:length(ar_b_values)
eqn = ar_b_values(i) * N * (1 – N / Ka) – N / (1 + N) == 0;
solutions = double(solve(eqn, N));
bifurcation_results(i, 1:length(solutions)) = solutions’;
end

% Plot bifurcation diagram
figure;
plot(ar_b_values, bifurcation_results, ‘b*’);
title(‘Bifurcation Diagram’);
xlabel(‘ar/b’);
ylabel(‘Equilibrium Points’);

% Set up the ODE function for phase plane analysis
f = @(N) ar_b * 0.4 * N .* (1 – N / Ka) – N ./ (1 + N);

% Create a grid of N values and compute their time derivatives
N_values = linspace(0, 50, 400);
dN_dtau = f(N_values);

% Plot phase plane
figure;
plot(N_values, dN_dtau);
title(‘Phase Plane Analysis’);
xlabel(‘Population Size N’);
ylabel(‘dN/dτ’);
grid on;
hold on;

% Mark zero crossings (equilibrium points)
zero_crossings = N_values(abs(dN_dtau) < 1e-3);
for i = 1:length(zero_crossings)
plot(zero_crossings(i), 0, ‘ro’);
end
hold off; differential equations, equation, plotting, plot, ode, homework MATLAB Answers — New Questions

​

I am trying to use dtm for solving SIR model, Although my code is run but I think the DTM part is wrong. I need help for DTM part
here is my code
function sir_model_dtm()

% Parameters
alpha = 0.3; % Infection rate
beta = 0.1; % Recovery rate
N = 1000; % Total population
I0 = 1; % Initial number of infected individuals
R0 = 0; % Initial number of recovered individuals
S0 = N – I0 – R0; % Initial number of susceptible individuals

% Time parameters
T = 100; % Total simulation time
dt = 1; % Time step

% Initialize arrays to store results
S = zeros(1, T);
I = zeros(1, T);
R = zeros(1, T);

% Initial conditions
S(1) = S0;
I(1) = I0;
R(1) = R0;

% Simulate the SIR model using DTM
for t = 2:T
% Compute new values
S(t) = S(t-1) – alpha*S(t-1)*I(t-1)/N * dt;
I(t) = I(t-1) + (alpha*S(t-1)*I(t-1)/N – beta*I(t-1)) * dt;
R(t) = R(t-1) + beta*I(t-1) * dt;
end

% Plot the results
t = 0:dt:T-dt;
plot(t, S, ‘b’, t, I, ‘r’, t, R, ‘g’, ‘LineWidth’, 2);
legend(‘Susceptible’, ‘Infectious’, ‘Recovered’);
xlabel(‘Time’);
ylabel(‘Number of individuals’);
title(‘SIR Model’);

endI am trying to use dtm for solving SIR model, Although my code is run but I think the DTM part is wrong. I need help for DTM part
here is my code
function sir_model_dtm()

% Parameters
alpha = 0.3; % Infection rate
beta = 0.1; % Recovery rate
N = 1000; % Total population
I0 = 1; % Initial number of infected individuals
R0 = 0; % Initial number of recovered individuals
S0 = N – I0 – R0; % Initial number of susceptible individuals

% Time parameters
T = 100; % Total simulation time
dt = 1; % Time step

% Initialize arrays to store results
S = zeros(1, T);
I = zeros(1, T);
R = zeros(1, T);

% Initial conditions
S(1) = S0;
I(1) = I0;
R(1) = R0;

% Simulate the SIR model using DTM
for t = 2:T
% Compute new values
S(t) = S(t-1) – alpha*S(t-1)*I(t-1)/N * dt;
I(t) = I(t-1) + (alpha*S(t-1)*I(t-1)/N – beta*I(t-1)) * dt;
R(t) = R(t-1) + beta*I(t-1) * dt;
end

% Plot the results
t = 0:dt:T-dt;
plot(t, S, ‘b’, t, I, ‘r’, t, R, ‘g’, ‘LineWidth’, 2);
legend(‘Susceptible’, ‘Infectious’, ‘Recovered’);
xlabel(‘Time’);
ylabel(‘Number of individuals’);
title(‘SIR Model’);

end I am trying to use dtm for solving SIR model, Although my code is run but I think the DTM part is wrong. I need help for DTM part
here is my code
function sir_model_dtm()

% Parameters
alpha = 0.3; % Infection rate
beta = 0.1; % Recovery rate
N = 1000; % Total population
I0 = 1; % Initial number of infected individuals
R0 = 0; % Initial number of recovered individuals
S0 = N – I0 – R0; % Initial number of susceptible individuals

% Time parameters
T = 100; % Total simulation time
dt = 1; % Time step

% Initialize arrays to store results
S = zeros(1, T);
I = zeros(1, T);
R = zeros(1, T);

% Initial conditions
S(1) = S0;
I(1) = I0;
R(1) = R0;

% Simulate the SIR model using DTM
for t = 2:T
% Compute new values
S(t) = S(t-1) – alpha*S(t-1)*I(t-1)/N * dt;
I(t) = I(t-1) + (alpha*S(t-1)*I(t-1)/N – beta*I(t-1)) * dt;
R(t) = R(t-1) + beta*I(t-1) * dt;
end

% Plot the results
t = 0:dt:T-dt;
plot(t, S, ‘b’, t, I, ‘r’, t, R, ‘g’, ‘LineWidth’, 2);
legend(‘Susceptible’, ‘Infectious’, ‘Recovered’);
xlabel(‘Time’);
ylabel(‘Number of individuals’);
title(‘SIR Model’);

end sir, matlab, dtm, differential equations, transform method, nonlinear, covid19 MATLAB Answers — New Questions

​

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

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(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));
endIn 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

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(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 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

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(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 ode45 MATLAB Answers — New Questions

​

Hello:
I am trying to apply the concept and sample code of the tutorial "3D brain tumor segmentation using deep learning" found on the below link –
https://in.mathworks.com/help/deeplearning/examples/segment-3d-brain-tumor-using-deep-learning.html
to a custom CT dataset in order to segment the lymph nodes. In my use case also there are 2 classes i.e background and lymph nodes. The current dataset has less data points (70 in total) and the split after preprocessing is 56 for train, 10 for validation and 4 for test.

I have modified the preprocessing code to accept the custom dataset and am cropping the ROI (i.e. the cubiod region around the abdomen which has the ground truth for the lumph nodes from the full CT Torso Volume)
I have visualized the patches that are extracted by the randomPatchExtractorDatastore and found that there is a class imbalance i.e. 4 out of 5 times the extracted patch has a Lymph node portion and 1 out of 5 times – it has the background portion. So, this may be causing an issue is what I am suspecting. (The dice accuracy suggest the same – as mentioed at the later part of this text)
In my case the volume image are 1 channel where as in the tutorial the volume image are 4 channel, so I have changed the input_layer to [64 64 64 1]
The code is working and there are no errors. The training process is also working but the accuracy starts at 100% and the loss starts at approx 1. So, this is unlike the graph shown in the tutorial where the accuracy gradually increase and reaches about 72 %. To train 1 EPOCH it takes like 75 mins on RTX 2080Ti so I tried training for 10 epochs as I was suspicious about the class imbalance of the randomly extracted patches.
I used PatchPerImage = 32 and the initial learning rate at 0.001 for the 1st 5 EPOCHs. The validation freq is 10.

Training Graph is as below –

The segmentation test shows that the model is just identifing the background and not learning to segment the lymph node. The Quantify Segmentation Accuracy section of the code suggest the same thing that the lymph nodes are not identified by the newtork.
Average Dice Score of the background across 4 test volumes = 0.99685 and
Average Dice Score of the lymph nodes across 4 test volumes = 0

Box plot of the Dice Score is as below –

I need some help and guidance on the following –

So, what can be the reason for this. Is class imbalance the reason or lack of training or few datapoint ?

How can i ensure that the patches in the randomPatchExtractorDatastore are equally balance accross both the classes i.e. Background and Lymph Nodes ?
Also, should i train the network for more epochs (like 50 or 100) before evaluating its efficacy ?
Is the limited data point the problem. I only have 175 total datapoints in the comprehensive dataset. So, should i use a higher patchPerImage value like 64 (instead of 16 as used in the tutorial)

Thanks for your guidance.Hello:
I am trying to apply the concept and sample code of the tutorial "3D brain tumor segmentation using deep learning" found on the below link –
https://in.mathworks.com/help/deeplearning/examples/segment-3d-brain-tumor-using-deep-learning.html
to a custom CT dataset in order to segment the lymph nodes. In my use case also there are 2 classes i.e background and lymph nodes. The current dataset has less data points (70 in total) and the split after preprocessing is 56 for train, 10 for validation and 4 for test.

I have modified the preprocessing code to accept the custom dataset and am cropping the ROI (i.e. the cubiod region around the abdomen which has the ground truth for the lumph nodes from the full CT Torso Volume)
I have visualized the patches that are extracted by the randomPatchExtractorDatastore and found that there is a class imbalance i.e. 4 out of 5 times the extracted patch has a Lymph node portion and 1 out of 5 times – it has the background portion. So, this may be causing an issue is what I am suspecting. (The dice accuracy suggest the same – as mentioed at the later part of this text)
In my case the volume image are 1 channel where as in the tutorial the volume image are 4 channel, so I have changed the input_layer to [64 64 64 1]
The code is working and there are no errors. The training process is also working but the accuracy starts at 100% and the loss starts at approx 1. So, this is unlike the graph shown in the tutorial where the accuracy gradually increase and reaches about 72 %. To train 1 EPOCH it takes like 75 mins on RTX 2080Ti so I tried training for 10 epochs as I was suspicious about the class imbalance of the randomly extracted patches.
I used PatchPerImage = 32 and the initial learning rate at 0.001 for the 1st 5 EPOCHs. The validation freq is 10.

Training Graph is as below –

The segmentation test shows that the model is just identifing the background and not learning to segment the lymph node. The Quantify Segmentation Accuracy section of the code suggest the same thing that the lymph nodes are not identified by the newtork.
Average Dice Score of the background across 4 test volumes = 0.99685 and
Average Dice Score of the lymph nodes across 4 test volumes = 0

Box plot of the Dice Score is as below –

I need some help and guidance on the following –

So, what can be the reason for this. Is class imbalance the reason or lack of training or few datapoint ?

How can i ensure that the patches in the randomPatchExtractorDatastore are equally balance accross both the classes i.e. Background and Lymph Nodes ?
Also, should i train the network for more epochs (like 50 or 100) before evaluating its efficacy ?
Is the limited data point the problem. I only have 175 total datapoints in the comprehensive dataset. So, should i use a higher patchPerImage value like 64 (instead of 16 as used in the tutorial)

Thanks for your guidance. Hello:
I am trying to apply the concept and sample code of the tutorial "3D brain tumor segmentation using deep learning" found on the below link –
https://in.mathworks.com/help/deeplearning/examples/segment-3d-brain-tumor-using-deep-learning.html
to a custom CT dataset in order to segment the lymph nodes. In my use case also there are 2 classes i.e background and lymph nodes. The current dataset has less data points (70 in total) and the split after preprocessing is 56 for train, 10 for validation and 4 for test.

I have modified the preprocessing code to accept the custom dataset and am cropping the ROI (i.e. the cubiod region around the abdomen which has the ground truth for the lumph nodes from the full CT Torso Volume)
I have visualized the patches that are extracted by the randomPatchExtractorDatastore and found that there is a class imbalance i.e. 4 out of 5 times the extracted patch has a Lymph node portion and 1 out of 5 times – it has the background portion. So, this may be causing an issue is what I am suspecting. (The dice accuracy suggest the same – as mentioed at the later part of this text)
In my case the volume image are 1 channel where as in the tutorial the volume image are 4 channel, so I have changed the input_layer to [64 64 64 1]
The code is working and there are no errors. The training process is also working but the accuracy starts at 100% and the loss starts at approx 1. So, this is unlike the graph shown in the tutorial where the accuracy gradually increase and reaches about 72 %. To train 1 EPOCH it takes like 75 mins on RTX 2080Ti so I tried training for 10 epochs as I was suspicious about the class imbalance of the randomly extracted patches.
I used PatchPerImage = 32 and the initial learning rate at 0.001 for the 1st 5 EPOCHs. The validation freq is 10.

Training Graph is as below –

The segmentation test shows that the model is just identifing the background and not learning to segment the lymph node. The Quantify Segmentation Accuracy section of the code suggest the same thing that the lymph nodes are not identified by the newtork.
Average Dice Score of the background across 4 test volumes = 0.99685 and
Average Dice Score of the lymph nodes across 4 test volumes = 0

Box plot of the Dice Score is as below –

I need some help and guidance on the following –

So, what can be the reason for this. Is class imbalance the reason or lack of training or few datapoint ?

How can i ensure that the patches in the randomPatchExtractorDatastore are equally balance accross both the classes i.e. Background and Lymph Nodes ?
Also, should i train the network for more epochs (like 50 or 100) before evaluating its efficacy ?
Is the limited data point the problem. I only have 175 total datapoints in the comprehensive dataset. So, should i use a higher patchPerImage value like 64 (instead of 16 as used in the tutorial)

Thanks for your guidance. 3d, unet, semantic segmentation, deep learning, custom dataset, own dataset, class imbalance, randompatchextractordatastore, random patch extractor datastore MATLAB Answers — New Questions

​