1 from Numeric import sqrt
2
4 """Cubic interpolation using f(a), f(b), g(a), and g(b).
5
6 Equations
7 ~~~~~~~~~
8
9 f(a) = a'a**3 + b'a**2 + c'a + d'
10 f(b) = a'b**3 + b'b**2 + c'b + d'
11 g(a) = 3a'a**2 + 2b'a + c'
12 g(b) = 3a'b**2 + 2b'b + c'
13
14
15 Interpolation
16 ~~~~~~~~~~~~~
17
18 The extrema are the roots of the quadratic equation:
19
20 3a'*alpha**2 + 2b'*alpha + c' = 0
21
22 The cubic interpolant is given by the formula:
23
24 g(b) + beta2 - beta1
25 ac = b - (b - a) . ---------------------
26 g(b) - g(a) + 2*beta2
27
28 where:
29 f(a) - f(b)
30 beta1 = g(a) + g(b) - 3 . -----------
31 a - b
32
33 if a < b:
34 beta2 = sqrt(beta1**2 - g(a).g(b))
35 else:
36 beta2 = -sqrt(beta1**2 - g(a).g(b))
37 """
38
39
40 beta1 = ga + gb - 3.0*(fa - fb)/(a - b)
41 if a < b:
42 beta2 = sqrt(beta1**2 - ga*gb)
43 else:
44 beta2 = -sqrt(beta1**2 - ga*gb)
45 alpha = b - (b - a)*(gb + beta2 - beta1)/(gb - ga + 2.0*beta2)
46 return alpha
47
48
49 -def cubic_ext(a, b, fa, fb, ga, gb, full_output=0):
50 """Cubic Extrapolation using f(a), f(b), g(a), and g(b).
51
52 Extrapolation
53 ~~~~~~~~~~~~~
54
55 The extrema are the roots of the quadratic equation:
56
57 3a'*alpha**2 + 2b'*alpha + c' = 0
58
59 The cubic extrapolant is given by the formula:
60
61 g(b) + beta2 - beta1
62 ac = b - (b - a) . ---------------------
63 g(b) - g(a) + 2*beta2
64
65 where:
66
67 f(a) - f(b)
68 beta1 = g(a) + g(b) - 3 . -----------
69 a - b
70
71 if a < b:
72 beta2 = sqrt(max(0.0, beta1**2 - g(a).g(b)))
73 else:
74 beta2 = -sqrt(max(0.0, beta1**2 - g(a).g(b)))
75
76 """
77
78
79 beta1 = ga + gb - 3.0*(fa - fb)/(a - b)
80 if a < b:
81 beta2 = sqrt(max(0.0, beta1**2 - ga*gb))
82 else:
83 beta2 = -sqrt(max(0.0, beta1**2 - ga*gb))
84
85 alpha = b - (b - a)*(gb + beta2 - beta1)/(gb - ga + 2.0*beta2)
86
87 if full_output:
88 return alpha, beta1, beta2
89 else:
90 return alpha
91
92
94 """Quadratic interpolation using f(a), f(b), and g(a).
95
96 The extremum of the quadratic is given by:
97
98 1 g(a)
99 aq = a + - . -------------------------
100 2 f(a) - f(b) - (a - b)g(a)
101 """
102
103 denom = fa - fb - (a - b)*ga
104 if denom == 0.0:
105 return 1e99
106 else:
107 return a + 0.5 * ga / denom
108
109
111 """Quadratic interpolation using g(a) and g(b).
112
113 The extremum of the quadratic is given by:
114
115 bg(a) - ag(b)
116 aq = -------------
117 g(a) - g(b)
118 """
119
120 denom = ga - gb
121 if denom == 0.0:
122 return 1e99
123 else:
124 return (b*ga - a*gb) / denom
125