A Simple Finite-Difference Approach

Let us, instead, seek a numerical solution by replacing the derivatives with finite differences. In a way, this is a step backwards in mathematics, since we are going to replace an equation from calculus with one which is purely algebraic. Recall from your first course in calculus that you crossed the border between algebra and calculus by taking a limit to define a derivative. That is, the derivative df/dx was defined as:

So, we might approximate dy/dt in our equation as:

where is `small.' Rearranging our differential equation slightly gives:

Again the time-dependent major radius of the cone at the fluid level is given by:

If we are willing to approximate dy/dt with the finite differences , then

is the equation we need to solve on the computer. Notice that Equation (21) is non-linear; that is, the right-hand side contains the solution y and the square of the solution y² . Note also that just as the vessel is about to empty, the radical will go to zero if we are not careful to require that , since as .

Now there is clearly some ambiguity about what value of the fluid height y we should use in the right-hand side of Equation (21). Two approaches are obvious.

• We could evaluate completely in terms of y(t), which is always known from the last time step. This includes evaluating the fluid surface radius R(y) using the values computed at the previous time step.

• We could iterate at each time level. That is, we could first solve for using the value of y from the previous time level. Then we could repeat the calculation several times without advancing the time level and presumably get better and better agreement between on the left-hand side and the value of y in the terms on the right-hand side.

This second approach is the one we should take. No matter what programming language one might choose to use, a program to obtain numerical solutions to Equation (21) should take roughly the following form:

1. Specify physical parameters and g in a consistent set of units.
2. Specify a time step size .
3. Set initial values for the dependent variable y(t) and R(y).
4. Solve the finite-difference Equation (21) for , using values of y(t) in the right-hand side.
5. Repeat the previous step using the values of just computed to evaluate the right-hand side.
6. Iterate by repeating this last step until the value of stops changing significantly.
7. Advance the solution another time step by beginning with Step 3.
8. Stop and write out the solution when is sufficiently close to zero, or whatever the final value is.

Below is a simple program written in Fortran which solves our equation in the manner outlined above.

	program cone
	parameter (imax=1500) ! Set the maximum dimensions.
	dimension y(imax),t(imax),r(imax),dydt(imax)
	open(1,file='cone1.plt',form='formatted')! Open a file for plotting.
c	Below we compute 'machine epsilon,' the smallest number your
c	computer can add to one and get something different from one.
c	This is very useful for testing for 'machine zero', etc.
	eps=1.
100	eps=eps/2.
	if(eps +1. .gt. 1.)goto 100
	print *,'Machine epsilon: ',eps

c	Read constant parameters.
c	Here they have been set at values sufficient for evaluating
c	the behavior of the system.
 	print *, 'Enter R0,Y0,ro,re, and g in that order.'
	R0=1.  ! major  radius of the cone
	Y0=1.  ! initial fluid level
	ro=.1  ! minor radius of the cone
	re=.01 ! radius of the drain 
	g=9.81 ! acceleration of gravity
c	read(5,*)R0,Y0,ro,re,g ! un-comment this line to change the
c	fixed values.
	dr=R0-ro!  This is just a convenient definition to shorten the
c	Fortran statements below.
c	Read computational parameters.  Here we need to choose a time
c	step size, dt.  You will notice that the smaller dt, the more
c	nearly you can compute y=0, but the price for small dt is 
c	a huge data file, since many time steps will be required to empty
c	the vessel.
	print *, 'Enter dt and itmax in that order.'
	read(5,*)dt,itmax
c
c	Initialize. Here we set values at t=0.
	y(1)=Y0
	yold=Y0
	t(1)=0.
	r(1)=R0
	i=1
c	Begin time stepping.
1	i=i+1 !  This statement begins a time step.
	y(i)=y(i-1)!  This gives a non-zero value of y to get the computation
c	started at each new time level.
	if(i .gt. imax)goto 2!  Stop if we've exceeded the variable dimensions.
	t(i)=dt+t(i-1)! increment time
	do iter=1,itmax!  begin the iteration loop.  
c	Usually only 3 or 4 iterations are required at each time step.
 	r(i)=ro+y(i)*dr/Y0	! compute the radius of the fluid surface
c	The following statements just evaluate parts of Equation (21)
	term1=sqrt(2.*g*y(i))*re**2*Y0**2
	ratio=re**4/r(i)**4
	term2=sqrt(1.-ratio)
	term3=ro**2*Y0**2+2.*dr*ro*y(i)*Y0+ dr**2*y(i)**2
	y(i)=y(i-1)-dt*term1/(term2*term3)!  This is Equation (21)
c	Below is just a safety net to catch the solution should we
c	compute a non-physical negative fluid level as we work in the
c	neighborhood of y=0.
	if(y(i) .lt. 0.)then
		y(i)=0.
	 	r(i)=ro+y(i)*dr/Y0	
		goto 2
		endif
c	The statement below just lets you monitor the progress.
	print *,'r(i),y(i),i ',r(i),y(i),i
c	Below tests to see if we are very close to y=0.
	if(y(i) .lt. sqrt(eps).and. y(i-1).lt.sqrt(eps))goto 2
c	Below we test to see how closely two successive values of
c	y in the iteration agree.
	if(yold .gt. 0.)resid=abs((yold-y(i))/yold)
	yold=y(i)! yold is the value from the previous iteration
	if(resid .lt. sqrt(eps))goto 4!  if we have converged, take another
c	time step
	enddo
4	goto 1 ! start another time step after the maximum number of
c	iterations regardless of whether the solution has converged
2	n=i
c	Below we compute the derivative of the fluid level.  Clearly
c	dy/dt is the velocity of the fluid surface.  Due to the
c	singular behavior of Equation (21) as y-->0, dy/dt becomes
c	infinite, so we don't try to compute the last couple of values.
 	do i=2,n-1
	dydt(i)=(y(i+1)-y(i-1))/dt
	dmax=max(abs(dydt(i)),dmax)
	enddo
	dydt(1)=(y(2)-y(1))/dt
c	Now we write out the solution for plotting.
 	do i=1,n-1
	dydt(i)=dydt(i)/dmax
	write(1,*)t(i),y(i),r(i),-dydt(i)
	enddo
	stop
	end

A plot of the solution for the cone using this code is available here. Note that as the fluid level approaches zero, the rate at which it drops (dy/dt) increases rapidly, since the remain volume is decreasing very rapidly.