September 25, 2014 --- Class 10 --- Details of secant and false position methods,
adapt.c, Oscillatory Motion, Non-conservation of
Energy in the Euler Algorithm, Sources of Numerical Analysis Code,
Final Project
Activities:
Details of secant and false position methods
In each method, we use two points on the curve
of the function we want to find the root of. We connect these
points by a straight line, and find the root of the line. That
x value is used as our next guess of the root. It is easy to
show using simlar triangles that if the two points are (x1,y1)
and (x2,y2). the value of x on the straight line for which y=0
is given by
x2 - y2* (x1-x2)/(y1-y2)
In ~sg/chap3, you will find code in secant.c and
false_position.c for testing these methods.
They implement the algorithms and print the
results after each step. You may vary the function f(x).
These are set up as examples of the algorithm, not as general
purpose root finding subroutines.
We tried to find the root of f(x) = tan(x) - 6 using the
secant method. It didn't work for one set of initial guesses,
but did for another. Look at the difference between secant.c
and secant2.c. The code secant does not converge, but
secant2, with different initial conditions does converge:
1.325900e+00 -1.998603e+00
1.374589e+00 -9.689273e-01
1.420405e+00 5.991223e-01
1.402900e+00 -1.000287e-01
1.405404e+00 -8.996072e-03
1.405652e+00 1.485253e-04
1.405648e+00 -2.169805e-07
1.405648e+00 -5.230483e-12
1.405648e+00 3.552714e-15
1.405648e+00 3.552714e-15
Notice how quickly this converges near the end. False position
does not converge as quickly:
8.878937e-01 -4.770704e+00
1.207087e+00 -3.372870e+00
1.328629e+00 -1.951659e+00
1.375649e+00 -9.408807e-01
1.393945e+00 -4.046209e-01
1.401079e+00 -1.645238e-01
1.403864e+00 -6.530093e-02
1.404951e+00 -2.566518e-02
1.405376e+00 -1.004793e-02
1.405541e+00 -3.927745e-03
1.405606e+00 -1.534439e-03
1.405631e+00 -5.993137e-04
1.405641e+00 -2.340556e-04
1.405645e+00 -9.140465e-05
1.405647e+00 -3.569534e-05
1.405647e+00 -1.393967e-05
1.405648e+00 -5.443676e-06
1.405648e+00 -2.125846e-06
1.405648e+00 -8.301780e-07
1.405648e+00 -3.241982e-07
1.405648e+00 -1.266047e-07
Note that in each set of output the variables are double
precision, but only a limited number of digits are printed out.
Putting it all together: adapt.c
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 as we saw above.
We will talk more about root finding later in the course.
There is a routine zbrent in Numerical Recipes that is
quite sophisticated.
The secant method is just fine for the problem solved in
adapt.c because the function is very close to a straight
line.
Gnu Scientific Library
We briefly looked at documentation in the Gnu Scientific Library
(GSL) for root finding.
http://www.gnu.org/software/gsl/manual/ is the top level link.
We looked, in particular at Sec. 33.8 Root Bracketing Algorithms.
The three routines gsl_root_fsolver_bisection, gsl_root_fsolver_falsepos,
and gsl_root_fsolver_brent are described. We used the first
two algorithms in class.
Chapter 4---Oscillatory Motion
We used the code ~sg/chap4/sho2.c and executable sho2 to
calculate the motion of a simple harmonic oscillator.
This code is based on the Euler method. We found that energy
is not conserved. In fact, we found that the energy grows
exponentially with time, no matter how small the
integration step size. sho_ec.c implements the Euler-Cromer
algorithm. We found that the energy is much better
conserved in this case. Here is a simple script
~sg/chap4/x_vs_dt.csh
to run the Euler-Cromer algorithm and plot the position at
t=2 sec as a function of step size.
#!/bin/csh
# script to plot position at t= 2 s vs dt
foreach dt ( .1 .05 .01 .005 .001 .0005 )
echo -n $dt " "
sho_ec <