Differential Equations with Inputs (Tank Draining Problem)

23 views (last 30 days)
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)

Accepted Answer

Bjorn Gustavsson
Bjorn Gustavsson on 3 Dec 2022
Edited: Bjorn Gustavsson on 3 Dec 2022
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
  2 Comments
QUINLAN
QUINLAN on 3 Dec 2022
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
Torsten on 3 Dec 2022
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.

Sign in to comment.

More Answers (1)

Sam Chak
Sam Chak on 8 Dec 2022
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

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!