The trapezoid algorithm works well for most integrals, but cannot deal with singularities. One way of avoiding the problem is to use rectangles that intersect the curve at the interval midpoints. The sample program shows the construction on a diagram. The program deals with a ball bouncing on an impenetrable floor at y=0. The potential barrier at the floor is infinite, and above the floor the potential energy is U(y) = mgy. If the maximum height reached by the ball is h, the total energy can be written E = mgh, and the time for the ball to go from y = 0 to y = h is
The region between y = h/8 and y = 7h/8 has a green line between the endpoints to form a trapezoid, and a red line through the midpoint to form a rectangle. For a function with upward curvature, trapezoids overestimate the area, and rectangles underestimate it, but both become better approximations as the number of intervals increases. The program uses the rectangle algorithm to calculate the time to go from 0 to h, and prints the results in the form used previously for the trapezoid algorithm. The value of h is 4.9 m, so we expect the time to equal 1s.
It is not surprising that the values with 20 and 40 intervals both underestimate the time, but it is surprising that they converge so slowly. The rectangle algorithm is an example of a 'bad' algorithm - the answer to why this is so, and how a 'good' algorithm can be found has its answer in the physics of the problem. In numerical form, the time integral is
There is no approximation involved in writing time as displacement divided by average velocity - the sum is exact at this point. The approximation enters when the average velocity is written as half the sum of the endpoint velocities. Earlier we defined a good approximation as one that is exact when the acceleration is constant. In the Feynman algorithm we used midpoint velocities, but the midpoint was the midpoint in time, not space. If the acceleration is constant, the velocity at the midpoint of a time interval is the average velocity over that interval, but the velocity at the midpoint of a space interval is the root-mean-square velocity for the interval.
The average velocity is half the sum of the endpoint velocities no matter what the interval provided that the acceleration is constant. This is the basis for the approximation in what we refer to as the 'vbar' algorithm. By our definition it is a good algorithm. It is exact for the bouncing ball. The yellow line in the figure is drawn at an elevation equal to two divided by the sum of the velocities at h/8 and 7h/8. The area under it is exactly equal to the area under the blue curve. The example in the previous section printed results for the vbar algorithm as well as the trapezoid algorithm, and the results were very accurate.
The next example uses the algorithm to solve a problem with no easy analytic solution: the period of a pendulum with a 90-degree amplitude. If arc length s is used as the displacement variable and sm is the amplitude,
There is one subtle point in the implementation of the algorithm. At the turning point the velocity is zero, and the function which evaluates it involves a square root. If 's' is incremented by 'ds' until it reaches 'sm', exact equality will not be found, and the argument of the square-root function will be a small positive or negative number. If it is negative, the program will not run. If it is positive, the error after the square root is taken is larger than might be expected. The program avoids the problem by setting the velocity to zero at the turning points.
This is also a problem that the Feynman algorithm solves with no difficulty. The program prints the Feynman result as well as the vbar result.
The technique developed in this section applies equally well in two dimensions where a second constant of the motion, angular momentum, enters the equation.