How to impose conditions on a system of ODEs and solve them using ode45?

Hi all,
I'm trying to model water flowing in and out of a cascade of 3 tanks. I am able to solve the system of ODE's using ode45, but I want to make sure that when the tank height reaches a certain minimum or maximum, that the height does not change. For some reason, in the function below, Matlab is unable to recognize the variables h(1), h(2) and h(3) for evaluation in the if statements, but ode45 seems to have no problem using them in the evaluation of dh(1), dh(2) and dh(3). Thanks in advance for your help.
function dh = threetanks(t,h)
dh = zeros(3,1);
Vdot0 = .00005;%m^3/s
Apipe = 3.14/4*.02^2; %m^2
Atank = 3.14/4*.1^2; %m^2
g = 9.81; %m/s^2
ff = 0.005;
L = .1; %m
D = .02; %m
F1 = Apipe*sqrt((2*g*(h(1)-h(2)))/(1+ff*L/D));
F2 = Apipe*sqrt((2*g*(h(2)-h(3)))/(1+ff*L/D));
F3 = Apipe*sqrt((2*g*(h(3)-0))/(1+ff*L/D));
dh(1) = (Vdot0 - F1 )/Atank;
dh(2) = (F1 - F2)/Atank;
dh(3) = (F2 - F3)/Atank;
if h(1) <= 1.25 && dh(1) < 0
h(1) = 1.25;
disp('Tank 1 hit bottom!')
end
if h(1) >= 1.75 && dh(1) > 0
h(1) = 1.75;
disp('Tank 1 overflowed!')
end
if h(2) <= .75 && dh(2) < 0
h(2) = 0.75;
disp('Tank 2 hit bottom!')
end
if h(2) >= 1.25 && dh(2) > 0
h(2) = 1.25;
disp('Tank 2 overflowed!')
end
if h(3) <= .25 && dh(3) < 0
h(3) = .25;
disp('Tank 3 hit bottom!')
end
if h(3) >= .75 && dh(3) > 0
h(3) = 0.75;
disp('Tank 3 overflowed!')
end
end

3 Comments

Never ever change the values for h in your function routine !!!
In the if-Statements, set dh to zero for the variable in question.
Best wishes
Torsten.
Torsten,
Thanks, I did that and it worked as I wanted it to. It seems obvious now that you told me. Just out of curiosity, why can't h be changed in the function routine? Does it have something to do with the way ode45 works?
It looks like you're maybe expecting h to be changed after you exit the function, which will not be the case since all input arguments are "pass by value", not "pass by reference". I hope you know what those basic computer science terms mean. Thus your code is deceptive and confusing to the reader. It's not that you can't do it - it will work - it's just that it's confusing to maintain and if you think h will stay changed after the function exits (which is what it looks like to us), then you're wrong about that.

Sign in to comment.

Answers (0)

Products

Tags

Asked:

BG
on 29 Jan 2015

Commented:

on 30 Jan 2015

Community Treasure Hunt

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

Start Hunting!