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

15 | CouNumber 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 | /* |
---|

55 | CouNumber expon = 4; |
---|

56 | |
---|

57 | inline CouNumber f (CouNumber x) |
---|

58 | {return pow (x, expon);} |
---|

59 | |
---|

60 | inline CouNumber fp (CouNumber x) |
---|

61 | {return expon * pow (x, expon-1);} |
---|

62 | |
---|

63 | inline CouNumber fpp (CouNumber x) |
---|

64 | {return expon * (expon-1) * pow (x, expon-2);} |
---|

65 | |
---|

66 | int 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 | */ |
---|