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

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

moved include files to make them doxygenable. Introduced three-way branching, with fixed intervals for now. Added check for small bound interval within all generateCuts()

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 100
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.