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

Go to the documentation of this file.
00001 /* $Id: powNewton.cpp 141 2009-06-03 04:19:19Z pbelotti $ */
00002 /*
00003  * Name:    powNewton.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: numerically find tangents to power functions
00006  *
00007  * (C) Carnegie-Mellon University, 2006-07.
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include <math.h>
00012 
00013 #include "CouenneTypes.hpp"
00014 
00015 #define MAX_ITER 10
00016 #define COU_POW_TOLERANCE 1e-12
00017 
00018 //#define DEBUG_POWNEW
00019 
00020 #ifndef DEBUG_POWNEW
00021 #include "funtriplets.hpp"
00022 #else
00023 #include <stdlib.h>
00024 #include <stdio.h>
00025 #endif
00026 
00027 CouNumber powNewton (CouNumber xc, CouNumber yc, 
00028                      unary_function f, 
00029                      unary_function fp, 
00030                      unary_function fpp) {
00031 
00032   // Find a zero to the function
00033   //
00034   // F(x) = x - xc + f'(x) (f(x) - yc)
00035   //
00036   // where f(x) is either x^k, exp(x), or log(x).
00037   // The derivative of F(x) is
00038   //
00039   // F'(x) = 1 + f''(x) (f(x) - yc) + (f'(x))^2
00040   //
00041   // Apply usual update:
00042   //
00043   // x(k+1) = x(k) - F(x(k))/F'(x(k))
00044 
00045   register CouNumber xk = xc;
00046 
00047   CouNumber fk  = f (xk) - yc,
00048             fpk = fp (xk),
00049             F   = fpk * fk,
00050             Fp  = 1 + fpp (xk) * fk + fpk * fpk;
00051 
00052   // Newton loop. Tolerance is set above
00053   for (int k = MAX_ITER; k--;) {
00054 
00055     xk -= F / Fp;
00056 
00057     fk  = f (xk) - yc;
00058     fpk = fp (xk);
00059     F   = xk - xc + fpk * fk;
00060 
00061     //    printf ("xk = %g; F = %g, fk = %g, fpk = %g\n", xk, F, fk, fpk);
00062 
00063     if (fabs (F) < COU_POW_TOLERANCE) break;
00064     Fp  = 1 + fpp (xk) * fk + fpk * fpk;
00065   }
00066 
00067   return xk;
00068 }
00069 
00070 #ifndef DEBUG_POWNEW
00071 
00073 CouNumber powNewton (CouNumber xc, CouNumber yc, funtriplet *tri) {
00074 
00075   // Find a zero to the function
00076   //
00077   // F(x) = x - xc + f'(x) (f(x) - yc)
00078   //
00079   // where f(x) is either x^k, exp(x), or log(x).
00080   // The derivative of F(x) is
00081   //
00082   // F'(x) = 1 + f''(x) (f(x) - yc) + (f'(x))^2
00083   //
00084   // Apply usual update:
00085   //
00086   // x(k+1) = x(k) - f(x(k))/f'(x(k))
00087 
00088   register CouNumber xk = xc;
00089 
00090   CouNumber fk  = tri -> F (xk) - yc,
00091             fpk = tri -> Fp (xk),
00092             F   = fpk * fk,
00093             Fp  = 1 + tri -> Fpp (xk) * fk + fpk * fpk;
00094 
00095   // Newton loop. Tolerance is set above
00096   for (int k = MAX_ITER; k--;) {
00097 
00098     xk -= F / Fp;
00099 
00100     fk  = tri -> F (xk) - yc;
00101     fpk = tri -> Fp (xk);
00102     F   = xk - xc + fpk * fk;
00103     if (fabs (F) < COU_POW_TOLERANCE) break;
00104     Fp  = 1 + tri -> Fpp (xk) * fk + fpk * fpk;
00105   }
00106 
00107   return xk;
00108 }
00109 #else
00110 
00112 inline CouNumber inv (register CouNumber arg) 
00113 {return 1.0 / arg;}
00114 
00115 
00117 inline CouNumber oppInvSqr (register CouNumber x) 
00118 {return (- inv (x*x));}
00119 
00120 
00122 inline CouNumber inv_dblprime (register CouNumber x) 
00123 {return (2 * inv (x*x*x));}
00124 
00125 
00126 int main (int argc, char **argv) {
00127 
00128   CouNumber r, 
00129     xc = atof (argv [2]),
00130     yc = atof (argv [3]);
00131 
00132   unary_function 
00133     f   = log,
00134     fp  = inv,
00135     fpp = oppInvSqr;
00136 
00137   //expon = atof (argv [1]);
00138 
00139   for (register int i=1; i--;)
00140     r = powNewton (xc, yc, f, fp, fpp);
00141 
00142   printf ("xc = %.14f: xk = %.15f, slope %.15f -- %.15f ==> [%.15f = -1?]\n", 
00143           xc, r, fp (r), 
00144                    (yc - f (r)) / (xc - r), 
00145           fp (r) * (yc - f (r)) / (xc - r));
00146 
00147   return 0;
00148 }
00149 
00150 #endif

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