Numerical Solution using MATLAB

To evaluate the Torricelli solution for the Big Gulp dribble, open an m-file and put in the right hand side of the differential equation:

       function hdot=biggulp(t,h)
       %     right hand side of Torricelli model for a matlab numerical
       %     solution of Big-Gulp Emptying.

       %     H is the initial height of the fluid (cm)
       %     a is the area of  the small hole (cm squared)
       %     g is the constant accelerationof gravity (cm/sec-sec) 
       %     rb is the X-section radius  at the hole (cm)
       %     rT is the X-section radius  at the initial top of fluid (cm)
       %     
       %     A is the current cross section  at top of fluid
       %     r is current radius at top of fluid
       %     h is the current height above hole of fluid (dependent variable)
       %     t is the time

       rb=3.5; rT=5;  H=16; g=981; 
       a=pi*(.2)^2;
       
       h=max(h,0);          %  To avoid negative heights
                            %  Next calculate denominator, d/dh [h A(h)]
       denom=pi*(rb^2 + 4*rb*(rT-rb)/H*h + 3*((rT-rb)/H*h)^2);
                            %  return value of RHS as hdot
       hdot=-a*sqrt(2*g*h)/denom;

Now save this as a file named biggulp.m and go back to the command window. The following command will run the mlb ODE solver:

       [t,h] = ode45('biggulp', [0,100], 16);

In this command 16 is the initial height and the time interval for the solution runs from 0 to 100 seconds. Plot this solution; does it look right? Can you plot the numerical and analytic solution on the same axes?

To actually find the time at which the Big Gulp emptied from the numerical solution we need to find the first index at which the vector h is zero. For this we can use the find command, which returns all indices of a vector for which some test is true. Try

       find(h<.001)

This command should return all the indices for which an element of the vector h has value smaller than .001. To find the index at which this first happens we want the smallest index, or

       emptyindex = min(find(h<.001))

Now we know what index corresponds to the element in the vector for h at which the element is smaller than .001. To see what the time was for this index, use

       t(emptyindex)

which should tell you the time at which the numerical solution got to be (pretty close to) zero.

	EXERCISE 7: Adapt the MATLAB function biggulp.m to solve the Torricelli model when the container has constant cross section. Pick some parameter values and find when the numerical solution predicts that the bucket will empty. How does this compare with the analytic prediction?

	EXERCISE 8: Upside-Down Big Gulp Suppose the Big Gulp is emptied through the same size hole, but in the top as opposed to the bottom. Should it empty more rapidly or more slowly than the right-side-up Big Gulp? Make a prediction in your head, using only physical principles, and then use MATLAB to solve the upside down Gulp problem, compare right side up and up side down solutions, and see if your intuitive expectations were correct.