My Project
PolynomialUtils.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_POLYNOMIAL_UTILS_HH
28 #define OPM_POLYNOMIAL_UTILS_HH
29 
30 #include <cmath>
31 #include <algorithm>
32 
34 
35 namespace Opm {
50 template <class Scalar, class SolContainer>
51 unsigned invertLinearPolynomial(SolContainer& sol,
52  Scalar a,
53  Scalar b)
54 {
55  if (std::abs(scalarValue(a)) < 1e-30)
56  return 0;
57 
58  sol[0] = -b/a;
59  return 1;
60 }
61 
78 template <class Scalar, class SolContainer>
79 unsigned invertQuadraticPolynomial(SolContainer& sol,
80  Scalar a,
81  Scalar b,
82  Scalar c)
83 {
84  // check for a line
85  if (std::abs(scalarValue(a)) < 1e-30)
86  return invertLinearPolynomial(sol, b, c);
87 
88  // discriminant
89  Scalar Delta = b*b - 4*a*c;
90  if (Delta < 0)
91  return 0; // no real roots
92 
93  Delta = sqrt(Delta);
94  sol[0] = (- b + Delta)/(2*a);
95  sol[1] = (- b - Delta)/(2*a);
96 
97  // sort the result
98  if (sol[0] > sol[1])
99  std::swap(sol[0], sol[1]);
100  return 2; // two real roots
101 }
102 
104 template <class Scalar, class SolContainer>
105 void invertCubicPolynomialPostProcess_(SolContainer& sol,
106  int numSol,
107  Scalar a,
108  Scalar b,
109  Scalar c,
110  Scalar d)
111 {
112  // do one Newton iteration on the analytic solution if the
113  // precision is increased
114  for (int i = 0; i < numSol; ++i) {
115  Scalar x = sol[i];
116  Scalar fOld = d + x*(c + x*(b + x*a));
117 
118  Scalar fPrime = c + x*(2*b + x*3*a);
119  if (std::abs(scalarValue(fPrime)) < 1e-30)
120  continue;
121  x -= fOld/fPrime;
122 
123  Scalar fNew = d + x*(c + x*(b + x*a));
124  if (std::abs(scalarValue(fNew)) < std::abs(scalarValue(fOld)))
125  sol[i] = x;
126  }
127 }
129 
147 template <class Scalar, class SolContainer>
148 unsigned invertCubicPolynomial(SolContainer* sol,
149  Scalar a,
150  Scalar b,
151  Scalar c,
152  Scalar d)
153 {
154  // reduces to a quadratic polynomial
155  if (std::abs(scalarValue(a)) < 1e-30)
156  return invertQuadraticPolynomial(sol, b, c, d);
157 
158  // normalize the polynomial
159  b /= a;
160  c /= a;
161  d /= a;
162  a = 1;
163 
164  // get rid of the quadratic term by subsituting x = t - b/3
165  Scalar p = c - b*b/3;
166  Scalar q = d + (2*b*b*b - 9*b*c)/27;
167 
168  if (std::abs(scalarValue(p)) > 1e-30 && std::abs(scalarValue(q)) > 1e-30) {
169  // At this point
170  //
171  // t^3 + p*t + q = 0
172  //
173  // with p != 0 and q != 0 holds. Introducing the variables u and v
174  // with the properties
175  //
176  // u + v = t and 3*u*v + p = 0
177  //
178  // leads to
179  //
180  // u^3 + v^3 + q = 0 .
181  //
182  // multiplying both sides with u^3 and taking advantage of the
183  // fact that u*v = -p/3 leads to
184  //
185  // u^6 + q*u^3 - p^3/27 = 0
186  //
187  // Now, substituting u^3 = w yields
188  //
189  // w^2 + q*w - p^3/27 = 0
190  //
191  // This is a quadratic equation with the solutions
192  //
193  // w = -q/2 + sqrt(q^2/4 + p^3/27) and
194  // w = -q/2 - sqrt(q^2/4 + p^3/27)
195  //
196  // Since w is equivalent to u^3 it is sufficient to only look at
197  // one of the two cases. Then, there are still 2 cases: positive
198  // and negative discriminant.
199  Scalar wDisc = q*q/4 + p*p*p/27;
200  if (wDisc >= 0) { // the positive discriminant case:
201  // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
202  Scalar u = - q/2 + sqrt(wDisc);
203  if (u < 0) u = - pow(-u, 1.0/3);
204  else u = pow(u, 1.0/3);
205 
206  // at this point, u != 0 since p^3 = 0 is necessary in order
207  // for u = 0 to hold, so
208  sol[0] = u - p/(3*u) - b/3;
209  // does not produce a division by zero. the remaining two
210  // roots of u are rotated by +- 2/3*pi in the complex plane
211  // and thus not considered here
212  invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d);
213  return 1;
214  }
215  else { // the negative discriminant case:
216  Scalar uCubedRe = - q/2;
217  Scalar uCubedIm = sqrt(-wDisc);
218  // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
219  Scalar uAbs = pow(sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3);
220  Scalar phi = atan2(uCubedIm, uCubedRe)/3;
221 
222  // calculate the length and the angle of the primitive root
223 
224  // with the definitions from above it follows that
225  //
226  // x = u - p/(3*u) - b/3
227  //
228  // where x and u are complex numbers. Rewritten in polar form
229  // this is equivalent to
230  //
231  // x = |u|*e^(i*phi) - p*e^(-i*phi)/(3*|u|) - b/3 .
232  //
233  // Factoring out the e^ terms and subtracting the additional
234  // terms, yields
235  //
236  // x = (e^(i*phi) + e^(-i*phi))*(|u| - p/(3*|u|)) - y - b/3
237  //
238  // with
239  //
240  // y = - |u|*e^(-i*phi) + p*e^(i*phi)/(3*|u|) .
241  //
242  // The crucial observation is the fact that y is the conjugate
243  // of - x + b/3. This means that after taking advantage of the
244  // relation
245  //
246  // e^(i*phi) + e^(-i*phi) = 2*cos(phi)
247  //
248  // the equation
249  //
250  // x = 2*cos(phi)*(|u| - p / (3*|u|)) - conj(x) - 2*b/3
251  //
252  // holds. Since |u|, p, b and cos(phi) are real numbers, it
253  // follows that Im(x) = - Im(x) and thus Im(x) = 0. This
254  // implies
255  //
256  // Re(x) = x = cos(phi)*(|u| - p / (3*|u|)) - b/3 .
257  //
258  // Considering the fact that u is a cubic root, we have three
259  // values for phi which differ by 2/3*pi. This allows to
260  // calculate the three real roots of the polynomial:
261  for (int i = 0; i < 3; ++i) {
262  sol[i] = cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
263  phi += 2*M_PI/3;
264  }
265 
266  // post process the obtained solution to increase numerical
267  // precision
268  invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);
269 
270  // sort the result
271  std::sort(sol, sol + 3);
272 
273  return 3;
274  }
275  }
276  // Handle some (pretty unlikely) special cases required to avoid
277  // divisions by zero in the code above...
278  else if (std::abs(scalarValue(p)) < 1e-30 && std::abs(scalarValue(q)) < 1e-30) {
279  // t^3 = 0, i.e. triple root at t = 0
280  sol[0] = sol[1] = sol[2] = 0.0 - b/3;
281  return 3;
282  }
283  else if (std::abs(scalarValue(p)) < 1e-30 && std::abs(scalarValue(q)) > 1e-30) {
284  // t^3 + q = 0,
285  //
286  // i. e. single real root at t=curt(q)
287  Scalar t;
288  if (-q > 0) t = pow(-q, 1./3);
289  else t = - pow(q, 1./3);
290  sol[0] = t - b/3;
291 
292  return 1;
293  }
294 
295  assert(std::abs(scalarValue(p)) > 1e-30 && std::abs(scalarValue(q)) > 1e-30);
296 
297  // t^3 + p*t = 0 = t*(t^2 + p),
298  //
299  // i. e. roots at t = 0, t^2 + p = 0
300  if (p > 0) {
301  sol[0] = 0.0 - b/3;
302  return 1; // only a single real root at t=0
303  }
304 
305  // two additional real roots at t = sqrt(-p) and t = -sqrt(-p)
306  sol[0] = -sqrt(-p) - b/3;;
307  sol[1] = 0.0 - b/3;
308  sol[2] = sqrt(-p) - b/3;
309 
310  return 3;
311 }
312 
330 template <class Scalar, class SolContainer>
331 unsigned cubicRoots(SolContainer* sol,
332  Scalar a,
333  Scalar b,
334  Scalar c,
335  Scalar d)
336 {
337  // reduces to a quadratic polynomial
338  if (std::abs(scalarValue(a)) < 1e-30)
339  return invertQuadraticPolynomial(sol, b, c, d);
340 
341  // We need to reduce the cubic equation to its "depressed cubic" form (however strange that sounds)
342  // Depressed cubic form: t^3 + p*t + q, where x = t - b/3*a is the transform we use when we have
343  // roots for t. p and q are defined below.
344  // Formula for p and q:
345  Scalar p = (3.0 * a * c - b * b) / (3.0 * a * a);
346  Scalar q = (2.0 * b * b * b - 9.0 * a * b * c + 27.0 * d * a * a) / (27.0 * a * a * a);
347 
348  // Check if we have three or one real root by looking at the discriminant, and solve accordingly with
349  // correct formula
350  Scalar discr = 4.0 * p * p * p + 27.0 * q * q;
351  if (discr < 0.0) {
352  // Find three real roots of a depressed cubic, using the trigonometric method
353  // Help calculation
354  Scalar theta = (1.0 / 3.0) * acos( ((3.0 * q) / (2.0 * p)) * sqrt(-3.0 / p) );
355 
356  // Calculate the three roots
357  sol[0] = 2.0 * sqrt(-p / 3.0) * cos( theta ) - b / (3.0 * a);
358  sol[1] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((2.0 * M_PI) / 3.0) ) - b / (3.0 * a);
359  sol[2] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((4.0 * M_PI) / 3.0) ) - b / (3.0 * a);
360 
361  // Sort in ascending order
362  std::sort(sol, sol + 3);
363 
364  // Return confirmation of three roots
365 
366  return 3;
367  }
368  else if (discr > 0.0) {
369 
370  // Find one real root of a depressed cubic using hyperbolic method. Different solutions depending on
371  // sign of p
372  Scalar t = 0;
373  if (p < 0) {
374  // Help calculation
375  Scalar theta = (1.0 / 3.0) * acosh( ((-3.0 * abs(q)) / (2.0 * p)) * sqrt(-3.0 / p) );
376 
377  // Root
378  t = ( (-2.0 * abs(q)) / q ) * sqrt(-p / 3.0) * cosh(theta);
379  }
380  else if (p > 0) {
381  // Help calculation
382  Scalar theta = (1.0 / 3.0) * asinh( ((3.0 * q) / (2.0 * p)) * sqrt(3.0 / p) );
383  // Root
384  t = -2.0 * sqrt(p / 3.0) * sinh(theta);
385 
386  }
387  else {
388  std::runtime_error(" p = 0 in cubic root solver!");
389  }
390 
391  // Transform t to output solution
392  sol[0] = t - b / (3.0 * a);
393  return 1;
394 
395  }
396  else {
397  // The discriminant, 4*p^3 + 27*q^2 = 0, thus we have simple (real) roots
398  // If p = 0 then also q = 0, and t = 0 is a triple root
399  if (p == 0) {
400  sol[0] = sol[1] = sol[2] = 0.0 - b / (3.0 * a);
401  }
402  // If p != 0, the we have a simple root and a double root
403  else {
404  sol[0] = (3.0 * q / p) - b / (3.0 * a);
405  sol[1] = sol[2] = (-3.0 * q) / (2.0 * p) - b / (3.0 * a);
406  std::sort(sol, sol + 3);
407  }
408  return 3;
409  }
410 }
411 } // end Opm
412 
413 #endif
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
unsigned invertQuadraticPolynomial(SolContainer &sol, Scalar a, Scalar b, Scalar c)
Invert a quadratic polynomial analytically.
Definition: PolynomialUtils.hpp:79
unsigned invertLinearPolynomial(SolContainer &sol, Scalar a, Scalar b)
Invert a linear polynomial analytically.
Definition: PolynomialUtils.hpp:51
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:331
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148