Differential Equations with Inputs (Tank Draining Problem)

All I am trying to do is let the user "customize" the parameters of the tank and the rate of water flowing out of it. I started with this Script/Function combo which works:

Script:

h0=str2num(input('Initial Height of Water\n', 's'));
tf=str2num(input('Time of Leak\n', 's'));
tspan=0:.05:tf
[t,h0]=ode45('Tank',tspan,h0)
plot(t,h0)
title('Height of Water vs Time Cylindrical Tank')
xlabel('time (s)')
ylabel('height (m)')

Function:

function dh = Tank(t,h0)
k = .1038
dh = -k*sqrt(h0);
end

BUT now I am trying to give the user some control and this combo is returning the error messages written below it:

Script:

cylD=str2num(input('Cylinder Diameter (m)\n', 's'));
holeD=str2num(input('Hole Diameter (m)\n', 's'));
h0=str2num(input('Initial Height of Water\n', 's'));
tf=str2num(input('Time of Leak\n', 's'));
g = 9.81
A = pi*(cylD/2)^2
hole = pi*(holeD/2)^2
tspan=0:.05:tf
[t,h0]=ode45(@TankIn,tspan,h0)
plot(t,h0)
title('Height of Water vs Time Cylindrical Tank')
xlabel('time (s)')
ylabel('height (m)')

Function:

function dh = TankIn(t,h0)
k = holeD/A*sqrt(2*g)
dh = -k*sqrt(h0);
end

Errors:

Error using odearguments

When the first argument to ode45 is a function handle, the tspan and y0 arguments must be supplied.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in Tank1 (line 12)

[t,h0]=ode45(@TankIn,tspan,h0)

Answer 1

[Bjorn Gustavsson]

What you should do is to make a small modification of your ODE-function - such that it takes an additional input argument where you can store all additional parameters that it needs. In this case a simple array will do the trick. Then you only need to collect these parameters in the script, and modify the ode45-call.

%% Script:
cylD=str2num(input('Cylinder Diameter (m)\n', 's'));
holeD=str2num(input('Hole Diameter (m)\n', 's'));
h0=str2num(input('Initial Height of Water\n', 's'));
tf=str2num(input('Time of Leak\n', 's'));
g = 9.81
A = pi*(cylD/2)^2
hole = pi*(holeD/2)^2
extraPars = [hole A g]; % Here we put the extra-parameters away, including 
                        % g so we can use this on the moon etc
tspan=0:.05:tf
[t,h0]=ode45(@(t,h) TankIn(t,h,extraPars),tspan,h0); % here we create an 
                                                     % "anonymous function" that 
                                                     % ony depends on t and
                                                     % h, that uses the
                                                     % current value of
                                                     % extraPars in TankIn
plot(t,h0)
title('Height of Water vs Time Cylindrical Tank')
xlabel('time (s)')
ylabel('height (m)')

%% Function:
function dh = TankIn(t,h0,extraPars) % here we have the extraPars as a 3rd input

holeD = extraPars(1); % then we only need to extract these parameters
A = extraPars(2);     % into the variables. I prefer to do it this way to keep
g = extraPars(3);     % different steps separate. 
k = holeD/A*sqrt(2*g);
dh = -k*sqrt(h0);
end

In my opinion it is cleaner to have all extra parameters in one additional array if possible. That way you dont need to modify the function-call in the ode45-calls you only need to modify the assignement to extraPars when you modify the TankIn-function.

HTH

[QUINLAN]

Super helpful thank you! Separate question but why am I getting a plot that goes into negative numbers? They shouldnt even be plotable on this axis I would have thought because they are complex numbers.

[Torsten]

h0 becomes complex, but integration can continue since ode45 also integrates complex-valued ODEs.

When plotting, the imaginary part of h is ignored.

Either detect the time when h becomes negative and only integrate up to this point or use the events detection facility of the ODE integrators. See the ballode example on how to choose this option.

Answer 2

[Sam Chak]

Hi @Quinlan Connell & @QUINLAN

You can also use the real() function to obtain only the Real part of the complex number.

% Parameters
cylD   = 1;                 % Cylinder Diameter (m)
holeD  = 0.2;               % Orifice (tech term for 'hole' in fluid flow) Diameter (m)
h0     = 50;                % Initial water height (m)
tf     = 300;               % Time of Leak
g      = 9.81;              % Earth's gravity
A      = pi*(cylD/2)^2;     % Tank cross-sectional area
hole   = pi*(holeD/2)^2;    % Area of Orifice 
Params = [hole A g]; 

tspan  = 0:0.05:tf;
[t, h] = ode45(@(t, h) TankIn(t, h, Params), tspan, h0);
plot(t, h), grid on
title('Water Level of the Cylindrical Tank')
xlabel('Time (s)')
ylabel('Water Level (m)')

% Water Tank
function dh = TankIn(t, h, Params)
    holeD =  Params(1); 
    A     =  Params(2);     
    g     =  Params(3);     
    k     =  holeD/A*sqrt(2*g);
    dh    = -k*real(sqrt(h));
end