September 28, 2004 --- Class 9 --- A Body Falling in a Position Dependent Potential (continued), Root finding Activities: Discussion of an object falling in a position dependent potential In view of the length of time since we last discussed this problem, we reviewed what we previously discussed. We wish to solve problem 3.4 of CSM. We drop a ball from high above the surface of the earth (neglecting air resistance!). From how great a height must we drop it so that it's velocity when it hits the ground is 1% different from what it would be assuming constant acceleration. There are several potential sources of error: 1) We know from the coffee cooling problem that when we do numerical integration there can be a step-size error. So we must be careful about this. 2) When we integrate the equations in time, we won't necessarily stop integrating exactly when the object reaches height 0. We may have to interpolate to find the velocity when the height is zero. 3) We also must find what height results in a 1% difference. Do we guess initial heights? How do we accurately find the height that corresponds to 1%? Step size error We have three approaches to the step size error issue. One possibility is to vary the step size and see how the answer varies. The accuracy on the velocity at height=0 must be much greater than 1%. The accuracy will depend on height and step size. We can also try to find a better method than Euler or Euler-Cromer. These are each first order methods. The Euler-Richardson method is one order higher. Try running this code with different step sizes. Note the accuracy of the energy and the position. (Interestingly, the difference is negligible for the velocity.) Another scheme involves using the Euler-Cromer algorithm with different step sizes to extrapolate to zero step size. This was described on the board and is coded in adapt_step.c. This code will increase or decrease the step size to achieve a desired level of accuracy during the integration. We examined the code in adapt_step.c to see how it increases or decreases the step size depending on how accurately our two estimates of the velocity agree with each other. We then ran the code with different initial step sizes and found that the step size could increase or decrease depending up whether we picked too large or too small an initial step size. Extrapolation back to height = 0 Since the integration goes to the next print time and we stop when the height is negative, we need to find the velocity when the height is zero. It is easy to do this neglecting the acceleration. However, it the print time is large, this may become a dominant source of error. We examined this by testing different print times. A second code, to extrapolate back to height = 0, using uniform acceleration was prepared by me the night before class. The code is in ~sg/adapt_step_nonlin_interp.c. We compiled the code and ran it to verify that it improves the calculation of the velocity at height = 0. Finding the ratio to the naive velocity The problem is stated in terms of the ratio of the velocity in the position dependent potential as compared with constant acceleration g. Of course, we can find the velocity with constant acceleration, so the code currently prints that out on a line by itself after the velocity from integration with the position dependent acceleration is printed. We may use AWK to get the ratio of velocities. We use the awk file ~/sg/vel_ratio.a. It contains two lines: /^V/{vel=-$4; rec=NR} NR==rec+1 && NR >1 { ratio =vel/$1; print ratio} We next ran the adapt_step program with several guesses for the initial height and found that as we increase the height the ratio of velocity in a position dependent potential to velocity assuming constant acceleration g is decreasing. It is a bit boring to keep guessing values that will get us to 0.99 as the ratio, so we realized that we are really trying to find the root of an equation: ratio( height) = 0.99, and we need to find the value of height that solves this equation. Root finding We devoted most of the rest of the class to a discussion of root finding. In the program adapt.c, the secant method is used to solve the problem of what height must the object be dropped from. The secant method is easy to implement because only the function, not its derviative is needed. Another advantage of this method is that it is not necessary to know that the root lies in a particular range. However, the secant method may be unstable. Some of the other methods require knowing that the root lies in a particular interval. One knows there is a root because the function changes sign. The bisection method evaluates the function at the middle of the interval. From the sign of the function at the midpoint, it can be determined which half of the interval contains the root. The function is then evaluated in the middle of that half interval, etc. This method doesn't have any stability problems, but it is slow to converge since the error only decreases by a factor of 2 each time the interval is cut in half. After all, in this method we are only using the sign of the function at the endpoints of the current interval. In the method of false position, we approximate the function by a straight line between the endpoints of the interval and use the root of the straight line to guess where the root is. We evaluate the function at the guessed root and replace the endpoint with the same sign of the function with the guess. We then use the two current endpoints in the same procedure. Convergence can be slow. The method of false position is a lot like the secant method except that in the secant method we always keep the two most recent points. When we do that, we might no longer be in the original interval that brackets a root. (In fact, we don't have to start off bracketing the root with the secant method.) This is what leads to the potential instability. In our next class, we will look further into these root finding algorithms and we will look at the adapt.c program.