/home/coin/SVN-release/OS-2.0.1/Couenne/src/convex/operators/trigNewton.cpp

Go to the documentation of this file.
00001 /* $Id: trigNewton.cpp 141 2009-06-03 04:19:19Z pbelotti $ */
00002 /*
00003  * Name:    trigNewton.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: numerically find tangents to (co)sines
00006  *
00007  * (C) Carnegie-Mellon University, 2006. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include <math.h>
00012 #include <stdlib.h>
00013 #include <stdio.h>
00014 #include "CouenneTypes.hpp"
00015 
00016 #define MAX_ITER 1000
00017 #define COU_TRIG_TOLERANCE 1e-12
00018 
00019 CouNumber trigNewton (CouNumber a, CouNumber l, CouNumber u) {
00020 
00021   // Find a zero to the function
00022   //
00023   // F(x) = cos x - (sin x - sin a) / (x - a)
00024   // 
00025   // whose derivative is
00026   //
00027   // F'(x) = - sin x - cos x / (x-a) + (sin x - sin a) / (x - a)^2
00028 
00029   if (l>u) {
00030     register CouNumber swap = l;
00031     l = u;
00032     u = swap;
00033   }
00034 
00035   register CouNumber xk = 0.5 * (u+l);
00036 
00037   CouNumber sina  = sin (a),
00038             sinxk = sin (xk),
00039             cosxk = cos (xk),
00040             dy    = sinxk - sina,
00041             dx    = xk - a,
00042             dydx  = dy/dx,
00043             F     = cosxk - dydx;
00044 
00045   // Newton loop. Tolerance is set above
00046   for (register int k = MAX_ITER; (fabs (F) > COU_TRIG_TOLERANCE) && k--;) {
00047 
00048     CouNumber Fp = sinxk + (cosxk - dydx) / dx;
00049 
00050     xk += F/Fp;
00051 
00052     if      (xk < l) xk = l;
00053     else if (xk > u) xk = u;
00054 
00055     sinxk = sin (xk);
00056     cosxk = cos (xk);
00057     dy    = sinxk - sina;
00058     dx    = xk - a;
00059     dydx  = dy/dx;
00060     F     = cosxk - dydx;
00061   }
00062 
00063   return xk;
00064 }
00065 
00066 /*
00067 int main (int argc, char **argv) {
00068 
00069   CouNumber a = atof (argv [1]),
00070             l = atof (argv [2]),
00071             u = atof (argv [3]), r;
00072 
00073   for (register int i=100000; i--;)
00074     r = trigNewton (a, l, u);
00075 
00076   printf ("b0 = %.14f: slope %.15f, derivative %.15f\n", 
00077           r, (sin (r) - sin (a)) / (r-a), cos (r));
00078 
00079   return 0;
00080 }
00081 */

Generated on Thu Oct 8 03:02:56 2009 by  doxygen 1.4.7