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

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');

Answer 1

[Cris LaPierre]

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.

[JoeKay]

Thank you for you response. Whilst I can add values for f2, f3 and f41that were missing how do I find the initial values x0?

[Cris LaPierre]

Your initial values are likely 0. You need to add in the constant rates of f2 and f3 that I pointed out to the corresponding ponds, and then include the rates corresponding to the stream that connects ponds 4 & 1.

[Torsten]

You can't find the initial conditions in your assignment.

But if you add the external sources f2 and f3, you can assume x0 = [0 0 0 0].

[JoeKay]

I have add the items you mentioned but I still get a zero answer.

% Define parameters
f2 = 1; % kg/h
f3 = 4; % kg/h
f12 = 0.05; % kg/h
f23 = 0.02; % kg/h
f24 = 0.03; % kg/h
f34 = 0.02; % kg/h
f41 = 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);

% 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');

[Cris LaPierre]

You have defined additional variables, though the value for f41 is wrong and V4 is missing. However, you have not updated ode_system, so they are still not incorporated into your simulation.

The items outlined in green are not included in ode_system.

% 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);
    ];

[Steven Lord]

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]

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');