source: branches/Couenne/Couenne/src/convex/operators/powNewton.cpp @ 506

Last change on this file since 506 was 506, checked in by pbelotti, 13 years ago

avoid cuts with huge coefficients. Added powNewton to compute deepest cut for powers

File size: 1.7 KB
Line 
1/*
2 * Name:    powNewton.cpp
3 * Author:  Pietro Belotti
4 * Purpose: numerically find tangents to power functions
5 *
6 * This file is licensed under the Common Public License (CPL)
7 */
8
9#include <math.h>
10#include <CouenneTypes.h>
11
12#define MAX_ITER 1000
13#define COU_POW_TOLERANCE 1e-12
14
15CouNumber powNewton (CouNumber xc, CouNumber yc, 
16                     unary_function f, 
17                     unary_function fp, 
18                     unary_function fpp) {
19
20  // Find a zero to the function
21  //
22  // F(x) = x - xc + f'(x) (f(x) - yc)
23  //
24  // where f(x) is either x^k, exp(x), or log(x).
25  // The derivative of F(x) is
26  //
27  // F'(x) = 1 + f''(x) (f(x) - yc) + (f'(x))^2
28  //
29  // Follow usual update:
30  //
31  // x(k+1) = x(k) - f(x(k))/f'(x(k))
32
33  register CouNumber xk = xc;
34
35  CouNumber fk  = f (xk) - yc,
36            fpk = fp (xk),
37            F   = fpk * fk,
38            Fp  = 1 + fpp (xk) * fk + fpk * fpk;
39
40  // Newton loop. Tolerance is set above
41  for (register int k = MAX_ITER; (fabs (F) > COU_POW_TOLERANCE) && k--;) {
42
43    xk -= F / Fp;
44
45    fk  = f (xk) - yc;
46    fpk = fp (xk);
47    F   = xk - xc + fpk * fk;
48    Fp  = 1 + fpp (xk) * fk + fpk * fpk;
49  }
50
51  return xk;
52}
53
54/*
55CouNumber expon = 4;
56
57inline CouNumber f (CouNumber x)
58{return pow (x, expon);}
59
60inline CouNumber fp (CouNumber x)
61{return expon * pow (x, expon-1);}
62
63inline CouNumber fpp (CouNumber x)
64{return expon * (expon-1) * pow (x, expon-2);}
65
66int main (int argc, char **argv) {
67
68  CouNumber r,
69    xc = atof (argv [2]),
70    yc = atof (argv [3]);
71
72  expon = atof (argv [1]);
73
74  for (register int i=10000; i--;)
75    r = powNewton (xc, yc, f, fp, fpp);
76
77  printf ("xc = %.14f: xk = %.15f, slope %.15f -- %.15f ==> %.15f\n",
78          xc, r, fp (r), (yc - f (r)) / (xc - r), fp (r) * (yc - f (r)) / (xc - r));
79
80  return 0;
81}
82*/
Note: See TracBrowser for help on using the repository browser.