What am I missing here? I have included all the parameters into my script but I still get zero as my answer. See script below

1 次查看（过去 30 天）
JoeKay
JoeKay 2024-5-1
评论： Cris LaPierre 2024-5-15
Ran in:
Consider four ponds connected by streams. Initially, the ponds were pristine.
At time t=0, after environmental accident, the second and third ponds are getting polluted with the sources, releasing pollutant substance with the following rates: f2=1 kg/h and f3=4 kg/h. Pollution spreads via the connecting streams to the other ponds, as shown in the Figure.
Formulate associated differential equations in terms of x1(t), x2(t), x3(t) and x4(t) - time functions, describing pollutant amounts in each lake. Solve problem numerically, using ode45 MATLAB procedure, and DETERMINE the amount of pollutant (in kg) in lake #2 towards the end of 48 hours of contamination. For a specific numerical example, take f12/V1=f41/V4=0.05; f23/V2=f34/V3=0.02; f24/V2=0.03 [all in 1/h].
Figure-2: Four ponds #1, #2, #3 and #4 of volumes V1, V2, V3, V4 [all in m^3], connected by streams with specific flow rates f12, f13, f23, f34, f41 [m^3/h].
Here is the script I came up with but when I run it I get zero as my anser. What am I missing?
% Define parameters
f12 = 0.05; % kg/h
f23 = 0.02; % kg/h
f24 = 0.03; % kg/h
f34 = 0.02; % kg/h
V1 = 1; % m^3
V2 = 1; % m^3
V3 = 1; % m^3
% Define time span
tspan = [0 48]; % 48 hours
% Define initial conditions
x0 = [0; 0; 0; 0]; % Initial pollutant amounts in each pond
% Define the function for the system of ODEs
ode_system = @(t, x) [
-f12/V1 * x(1);
f12/V1 * x(1) - (f23/V2 + f24/V2) * x(2);
f23/V2 * x(2) - f34/V3 * x(3);
f24/V2 * x(2) + f34/V3 * x(3);
];
% Solve the system of ODEs
[t, x] = ode45(ode_system, tspan, x0);
% Extract pollutant amount in pond 2 at the end of 48 hours
x2_end = x(end, 2);
% Display result
fprintf('Pollutant amount in pond 2 at the end of 48 hours: %.4f kg\n', x2_end);
Pollutant amount in pond 2 at the end of 48 hours: 0.0000 kg
% Plot results
plot(t, x);
xlabel('Time (hours)');
ylabel('Pollutant amount (kg)');
legend('Pond 1', 'Pond 2', 'Pond 3', 'Pond 4');
title('Pollutant Amounts in Ponds Over Time');

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

回答（1 个）

Cris LaPierre
Cris LaPierre 2024-5-1
编辑：Cris LaPierre 2024-5-1
All things being equal, your answer will always be zero based on how ode_system is defined if x0 is a vector of zeros.
The main issue is that your system does not account for f2 and f3: "f2=1 kg/h and f3=4 kg/h". You need to add these fixed rates to your equations.
You also appear to have left one of the streams out of your system of equations. You need to model the entire system to get the expected resutls.
  7 个评论
显示 5更早的评论
Steven Lord
Steven Lord 2024-5-1
What is "t = 0" for your problem?
If the concentration is 0 for all the lakes, that means it's right before the pollution occurs and so you're solving the system without any pollution. In that case it makes sense for the concentration to always be 0, as the pollution doesn't occur.
So you probably want t = 0 to be right when the pollution hits the first lake, in which case one of the components of your x0 should be non-zero.
Cris LaPierre
Cris LaPierre 2024-5-15
Ran in:
Try this
% Define parameters
f12 = 0.05; % kg/h
f23 = 0.02; % kg/h
f24 = 0.03; % kg/h
f34 = 0.02; % kg/h
V1 = 1; % m^3
V2 = 1; % m^3
V3 = 1; % m^3
%% Define missing parameters
f2 = 1;
f3 = 4;
f41 = 0.05;
V4 = 1;
% Define time span
tspan = [0 48]; % 48 hours
% Define initial conditions
x0 = [0; 0; 0; 0]; % Initial pollutant amounts in each pond
%% Add missing components (f(2), f(3), f41/v4*x(4)
% Define the function for the system of ODEs
ode_system = @(t, x) [
-f12/V1 * x(1) + f41/V4 * x(4); % need to account for input from #4
f2 + f12/V1 * x(1) - (f23/V2 + f24/V2) * x(2); % need to account for input f2
f3 + f23/V2 * x(2) - f34/V3 * x(3); % need to account for input f3
f24/V2 * x(2) + f34/V3 * x(3) - f41/V4 * x(4); % need to account for output to #1
];
% Solve the system of ODEs
[t, x] = ode45(ode_system, tspan, x0);
% Extract pollutant amount in pond 2 at the end of 48 hours
x2_end = x(end, 2);
% Display result
fprintf('Pollutant amount in pond 2 at the end of 48 hours: %.4f kg\n', x2_end);
Pollutant amount in pond 2 at the end of 48 hours: 31.0563 kg
% Plot results
plot(t, x);
xlabel('Time (hours)');
ylabel('Pollutant amount (kg)');
legend('Pond 1', 'Pond 2', 'Pond 3', 'Pond 4');
title('Pollutant Amounts in Ponds Over Time');

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Scopes and Data Logging 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by