Task

Try to run the code with different values of the time step \Delta t . How do you expect the approximate solution will change? Try with the values \Delta t = 0.01, 0.1, 1, 5, 10, 20 , while keeping the other parameters fixed.

The figure shows the solution (water level in the tank) obtained using different values of the time step \Delta t .

water level in the tank

Figure 1: Water level in the tank for the considered values of \Delta t .

When \Delta t gets smaller, the graphs approach the analytic solution. For a large value of \Delta t our approximation wildly oscillates up and down, showing a completely different behaviour from the real solution! Choosing \Delta t small enough is thus very important to get reliable and accurate results.

We recall that we obtained our approximation by substituting the time derivative \dfrac{dh}{dt} with a finite ratio:

\dfrac{dh}{dt}(\bar t) \approx \dfrac{h(\bar t + \Delta t) - h(\bar t)}{\Delta t}.

This approximation is justified by the definition of derivative:

\dfrac{dh}{dt}(\bar t) = \lim_{\Delta t \to 0} \dfrac{h(\bar t + \Delta t) - h(\bar t)}{\Delta t}.

So, as \Delta t becomes smaller, the ratio becomes a better approximation of the derivative and we expect that our approximation gets closer to the true solution.

Now that you have all the information you need to complete the task, use at least one of the two workspaces from previous lessons (Python or MATLAB) and write the code that gives the correct solution.

To experiment with different values of \Delta t , we use a for loop. Then, we save the simulated data in a list and plot them to compare the results. The full Python code for the solution is the following:


import matplotlib.pyplot as plt  # import package for plot

# Physical parameters
Q_in = 1.0      # inflow rate
A_t = 1.0       # tank cross-sectional area (m^2)
A_out = 0.1     # outlet area (m^2)
g = 9.81        # gravity
h0 = 1.0        # initial water level

# Simulation time parameters
T = 100.0       # total simulation time

# outflow rate Q_out
def Q_out(h):
    return A_out * (2 * g * h) ** 0.5

plt.figure(figsize=(10, 5))

# loop over different values for the time step
for dt in [0.01, 0.1, 1, 5, 10, 20]:
    N = int(T / dt)  # number of time steps

    # lists to store the solution
    times = [0.0] * (N + 1)
    h = [0.0] * (N + 1)
    h[0] = h0

    for i in range(N):
        dh = dt * (Q_in - Q_out(h[i]))
        h[i + 1] = h[i] + dh
        times[i + 1] = (i + 1) * dt

    plt.plot(times, h, label="dt = " + str(dt))

plt.xlabel("Time (s)")
plt.ylabel("Water Level (m)")
plt.title("Water Level Over Time for Different dt")
plt.legend()
plt.grid()
plt.show()

To experiment with different values of \Delta t , we use a for loop. Then, we can compare the results by reporting the graphs in the same figure. The full Matlab code for the solution is the following:


% Function to define the inflow rate function
Q_in = @(t) 1;

% Function to calculate the outflow rate based on the water level and outlet area
Q_out = @(outlet_area, water_level) sqrt(2 * 9.81 * max(water_level, 0)) * outlet_area;

A_out = 0.1;  % Outlet area (m^2)
A_t   = 1.0;  % Tank cross-sectional area (m^2)
h0    = 1.0;  % Initial water level

% Simulation time parameters
T = 100;      % Total simulation time (seconds)

% Create figure for plotting
figure();
grid on;
hold on;

% Loop over different time step sizes
for dt = [0.01, 0.1, 1, 5, 10, 20]
    N = round(T / dt);  % Number of time steps

    % Create arrays to store data for plotting
    times = zeros(1, N + 1);
    h     = zeros(1, N + 1);
    h(1)  = h0;
    t     = 0;

    % Simulate the water level over time
    for i = 1:N
        t = t + dt;
        h(i + 1) = h(i) + (Q_in(t) - Q_out(A_out, h(i))) * dt;
        times(i + 1) = t;
    end

    plot(times, h, '-');
end

legend('dt = 0.01', 'dt = 0.1', 'dt = 1', 'dt = 5', 'dt = 10', 'dt = 20');