1 /* $Id$
2 *
3 * Name: powNewton.cpp
4 * Author: Pietro Belotti
5 * Purpose: numerically find tangents to power functions
6 *
7 * (C) Carnegie-Mellon University, 2006-07.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include <math.h>
12
13 #include "CouenneFunTriplets.hpp"
14
15 #define MAX_ITER 10
16 #define COU_POW_TOLERANCE 1e-12
17
18 //#define DEBUG_POWNEW
19
20 #ifndef DEBUG_POWNEW
21 #include "CouenneTypes.hpp"
22 #else
23 #include <stdlib.h>
24 #include <stdio.h>
25 #endif
26
27 namespace Couenne {
28
powNewton(CouNumber xc,CouNumber yc,unary_function f,unary_function fp,unary_function fpp)29 CouNumber powNewton (CouNumber xc, CouNumber yc,
30 unary_function f,
31 unary_function fp,
32 unary_function fpp) {
33
34 // Find a zero to the function
35 //
36 // F(x) = x - xc + f'(x) (f(x) - yc)
37 //
38 // where f(x) is either x^k, exp(x), or log(x).
39 // The derivative of F(x) is
40 //
41 // F'(x) = 1 + f''(x) (f(x) - yc) + (f'(x))^2
42 //
43 // Apply usual update:
44 //
45 // x(k+1) = x(k) - F(x(k))/F'(x(k))
46
47 register CouNumber xk = xc;
48
49 CouNumber fk = f (xk) - yc,
50 fpk = fp (xk),
51 F = fpk * fk,
52 Fp = 1 + fpp (xk) * fk + fpk * fpk;
53
54 // Newton loop. Tolerance is set above
55 for (int k = MAX_ITER; k--;) {
56
57 xk -= F / Fp;
58
59 fk = f (xk) - yc;
60 fpk = fp (xk);
61 F = xk - xc + fpk * fk;
62
63 // printf ("xk = %g; F = %g, fk = %g, fpk = %g\n", xk, F, fk, fpk);
64
65 if (fabs (F) < COU_POW_TOLERANCE) break;
66 Fp = 1 + fpp (xk) * fk + fpk * fpk;
67 }
68
69 return xk;
70 }
71
72 #ifndef DEBUG_POWNEW
73
74 ///
powNewton(CouNumber xc,CouNumber yc,funtriplet * tri)75 CouNumber powNewton (CouNumber xc, CouNumber yc, funtriplet *tri) {
76
77 // Find a zero to the function
78 //
79 // F(x) = x - xc + f'(x) (f(x) - yc)
80 //
81 // where f(x) is either x^k, exp(x), or log(x).
82 // The derivative of F(x) is
83 //
84 // F'(x) = 1 + f''(x) (f(x) - yc) + (f'(x))^2
85 //
86 // Apply usual update:
87 //
88 // x(k+1) = x(k) - f(x(k))/f'(x(k))
89
90 register CouNumber xk = xc;
91
92 CouNumber fk = tri -> F (xk) - yc,
93 fpk = tri -> Fp (xk),
94 F = fpk * fk,
95 Fp = 1 + tri -> Fpp (xk) * fk + fpk * fpk;
96
97 // Newton loop. Tolerance is set above
98 for (int k = MAX_ITER; k--;) {
99
100 xk -= F / Fp;
101
102 fk = tri -> F (xk) - yc;
103 fpk = tri -> Fp (xk);
104 F = xk - xc + fpk * fk;
105 if (fabs (F) < COU_POW_TOLERANCE) break;
106 Fp = 1 + tri -> Fpp (xk) * fk + fpk * fpk;
107 }
108
109 return xk;
110 }
111 #else
112
113 /// the operator itself
inv(register CouNumber arg)114 inline CouNumber inv (register CouNumber arg)
115 {return 1.0 / arg;}
116
117
118 /// derivative of inv (x)
oppInvSqr(register CouNumber x)119 inline CouNumber oppInvSqr (register CouNumber x)
120 {return (- inv (x*x));}
121
122
123 /// inv_dblprime, second derivative of inv (x)
inv_dblprime(register CouNumber x)124 inline CouNumber inv_dblprime (register CouNumber x)
125 {return (2 * inv (x*x*x));}
126
127
main(int argc,char ** argv)128 int main (int argc, char **argv) {
129
130 CouNumber r,
131 xc = atof (argv [2]),
132 yc = atof (argv [3]);
133
134 unary_function
135 f = log,
136 fp = inv,
137 fpp = oppInvSqr;
138
139 //expon = atof (argv [1]);
140
141 for (register int i=1; i--;)
142 r = powNewton (xc, yc, f, fp, fpp);
143
144 printf ("xc = %.14f: xk = %.15f, slope %.15f -- %.15f ==> [%.15f = -1?]\n",
145 xc, r, fp (r),
146 (yc - f (r)) / (xc - r),
147 fp (r) * (yc - f (r)) / (xc - r));
148
149 return 0;
150 }
151
152 #endif
153
154 }
155