The OpenFOAM Foundation
cubicEqn.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
9  This file is part of OpenFOAM.
10
11  OpenFOAM is free software: you can redistribute it and/or modify it
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23
24 \*---------------------------------------------------------------------------*/
25
26 #include "linearEqn.H"
28 #include "cubicEqn.H"
29
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31
33 {
34  /*
35
36  This function solves a cubic equation of the following form:
37
38  a*x^3 + b*x^2 + c*x + d = 0
39  x^3 + B*x^2 + C*x + D = 0
40
41  The following two substitutions are used:
42
43  x = t - B/3
44  t = w - P/3/w
45
46  This reduces the problem to a quadratic in w^3.
47
48  w^6 + Q*w^3 - P = 0
49
50  Where Q and P are given in the code below.
51
52  The properties of the cubic can be related to the properties of this
53  quadratic in w^3. If it has a repeated root a zero, the cubic has a tripl
54  root. If it has a repeated root not at zero, the cubic has two real roots,
55  one repeated and one not. If it has two complex roots, the cubic has three
56  real roots. If it has two real roots, then the cubic has one real root and
57  two complex roots.
58
59  This is solved for the most numerically accurate value of w^3. See the
60  quadratic function for details on how to pick a value. This single value of
61  w^3 can yield up to three cube roots for w, which relate to the three
62  solutions for x.
63
64  Only a single root, or pair of conjugate roots, is directly evaluated; the
65  one, or ones with the lowest relative numerical error. Root identities are
66  then used to recover the remaining roots, possibly utilising a quadratic
67  and/or linear solution. This seems to be a good way of maintaining the
68  accuracy of roots at very different magnitudes.
69
70  */
71
72  const scalar a = this->a();
73  const scalar b = this->b();
74  const scalar c = this->c();
75  const scalar d = this->d();
76
77  if (a == 0)
78  {
79  return Roots<3>(quadraticEqn(b, c, d).roots(), roots::nan, 0);
80  }
81
82  // This is assumed not to over- or under-flow. If it does, all bets are off.
83  const scalar p = c*a - b*b/3;
84  const scalar q = b*b*b*(2.0/27.0) - b*c*a/3 + d*a*a;
85  const scalar disc = p*p*p/27 + q*q/4;
86
87  // How many roots of what types are available?
88  const bool oneReal = disc == 0 && p == 0;
89  const bool twoReal = disc == 0 && p != 0;
90  const bool threeReal = disc < 0;
91  //const bool oneRealTwoComplex = disc > 0;
92
93  static const scalar sqrt3 = sqrt(3.0);
94
95  scalar x;
96
97  if (oneReal)
98  {
99  const Roots<1> r = linearEqn(- a, b/3).roots();
100  return Roots<3>(r.type(0), r[0]);
101  }
102  else if (twoReal)
103  {
104  if (q*b > 0)
105  {
106  x = - 2*cbrt(q/2) - b/3;
107  }
108  else
109  {
110  x = cbrt(q/2) - b/3;
111  const Roots<1> r = linearEqn(- a, x).roots();
112  return Roots<3>(Roots<2>(r, r), linearEqn(x*x, a*d).roots());
113  }
114  }
115  else if (threeReal)
116  {
117  const scalar wCbRe = - q/2, wCbIm = sqrt(- disc);
118  const scalar wAbs = cbrt(hypot(wCbRe, wCbIm));
119  const scalar wArg = atan2(wCbIm, wCbRe)/3;
120  const scalar wRe = wAbs*cos(wArg), wIm = wAbs*sin(wArg);
121  if (b > 0)
122  {
123  x = - wRe - mag(wIm)*sqrt3 - b/3;
124  }
125  else
126  {
127  x = 2*wRe - b/3;
128  }
129  }
130  else // if (oneRealTwoComplex)
131  {
132  const scalar wCb = - q/2 - sign(q)*sqrt(disc);
133  const scalar w = cbrt(wCb);
134  const scalar t = w - p/(3*w);
135  if (p + t*b < 0)
136  {
137  x = t - b/3;
138  }
139  else
140  {
141  const scalar xRe = - t/2 - b/3, xIm = sqrt3/2*(w + p/3/w);
142  x = - a*a*d/(xRe*xRe + xIm*xIm);
143
144  // This form of solving for the remaining roots seems more stable
145  // for this case. This has been determined by trial and error.
146  return
147  Roots<3>
148  (
149  linearEqn(- a, x).roots(),
150  quadraticEqn(a*x, x*x + b*x, - a*d).roots()
151  );
152  }
153  }
154
155  return
156  Roots<3>
157  (
158  linearEqn(- a, x).roots(),
159  quadraticEqn(- x*x, c*x + a*d, d*x).roots()
160  );
161 }
162
163
164 // ************************************************************************* //
dimensionedScalar sign(const dimensionedScalar &ds)
Roots< 2 > roots() const
Get the roots.
scalar c() const
Definition: cubicEqnI.H:67
dimensionedScalar sqrt(const dimensionedScalar &ds)
scalar d() const
Definition: cubicEqnI.H:73
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
Definition: Roots.H:68
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar cbrt(const dimensionedScalar &ds)
scalar b() const
Definition: cubicEqnI.H:61
dimensionedScalar sin(const dimensionedScalar &ds)
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
Quadratic equation of the form a*x^2 + b*x + c = 0.
Linear equation of the form a*x + b = 0.
Definition: linearEqn.H:48
Roots< 3 > roots() const
Get the roots.
Definition: cubicEqn.C:32
void type(const direction i, const roots::type t)
Set the type of the i-th root.
Definition: RootsI.H:120
Roots< 1 > roots() const
Get the roots.
Definition: linearEqnI.H:89
scalar a() const
Definition: cubicEqnI.H:55
dimensioned< scalar > mag(const dimensioned< Type > &)
volScalarField & p
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)