The open source CFD toolbox
cubicEqn.H
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2017 OpenFOAM Foundation
9 Copyright (C) 2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Class
28 Foam::cubicEqn
29
30Description
31 Container to encapsulate various operations for
32 cubic equation of the forms with real coefficients:
33
34 \f[
35 a*x^3 + b*x^2 + c*x + d = 0
36 x^3 + B*x^2 + C*x + D = 0
37 \f]
38
39 The following two substitutions into the above forms are used:
40
41 \f[
42 x = t - B/3
43 t = w - P/3/w
44 \f]
45
46 This reduces the problem to a quadratic in \c w^3:
47
48 \f[
49 w^6 + Q*w^3 - P = 0
50 \f]
51
52 where \c Q and \c P are given in the code.
53
54 The properties of the cubic can be related to the properties of this
56
57 If the quadratic eqn has two identical real roots at zero, three identical
58 real roots exist in the cubic eqn.
59
60 If the quadratic eqn has two identical real roots at non-zero, two identical
61 and one distinct real roots exist in the cubic eqn.
62
63 If the quadratic eqn has two complex roots (a complex conjugate pair),
64 three distinct real roots exist in the cubic eqn.
65
66 If the quadratic eqn has two distinct real roots, one real root and two
67 complex roots (a complex conjugate pair) exist in the cubic eqn.
68
69 The quadratic eqn is solved for the most numerically accurate value
71 pick a value. This single value of \c w^3 can yield up to three cube
72 roots for \c w, which relate to the three solutions for \c x.
73
74 Only a single root, or pair of conjugate roots, is directly evaluated; the
75 one, or ones with the lowest relative numerical error. Root identities are
76 then used to recover the remaining roots, possibly utilising a quadratic
77 and/or linear solution. This seems to be a good way of maintaining the
78 accuracy of roots at very different magnitudes.
79
80 Reference:
81 \verbatim
82 Kahan's algo. to compute 'p' using fused multiply-adds (tag:JML):
83 Jeannerod, C. P., Louvet, N., & Muller, J. M. (2013).
84 Further analysis of Kahan's algorithm for the accurate
85 computation of 2× 2 determinants.
86 Mathematics of Computation, 82(284), 2245-2264.
87 DOI:10.1090/S0025-5718-2013-02679-8
88 \endverbatim
89
91 Test-cubicEqn.C
92
93SourceFiles
94 cubicEqnI.H
95 cubicEqn.C
96
97\*---------------------------------------------------------------------------*/
98
99#ifndef cubicEqn_H
100#define cubicEqn_H
101
102#include "Roots.H"
103
104// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
105
106namespace Foam
107{
108
109/*---------------------------------------------------------------------------*\
110 Class cubicEqn Declaration
111\*---------------------------------------------------------------------------*/
113class cubicEqn
114:
115 public VectorSpace<cubicEqn, scalar, 4>
116{
117public:
118
119 //- Component labeling enumeration
120 enum components { A, B, C, D };
121
122
123 // Constructors
124
125 //- Construct null
126 inline cubicEqn();
127
128 //- Construct initialized to zero
129 inline cubicEqn(const Foam::zero);
130
131 //- Construct from components
132 inline cubicEqn
133 (
134 const scalar a,
135 const scalar b,
136 const scalar c,
137 const scalar d
138 );
139
140
141 // Member Functions
142
143 // Access
144
145 inline scalar a() const;
146 inline scalar b() const;
147 inline scalar c() const;
148 inline scalar d() const;
149
150 inline scalar& a();
151 inline scalar& b();
152 inline scalar& c();
153 inline scalar& d();
154
155 // Evaluate
156
157 //- Evaluate the cubic equation at x
158 inline scalar value(const scalar x) const;
159
160 //- Evaluate the derivative of the cubic equation at x
161 inline scalar derivative(const scalar x) const;
162
163 //- Estimate the error of evaluation of the cubic equation at x
164 inline scalar error(const scalar x) const;
165
166 //- Return the roots of the cubic equation with no particular order
167 // if discriminant > 0: return three distinct real roots
168 // if discriminant < 0: return one real root and one complex root
169 // (one member of the complex conjugate pair)
170 // if discriminant = 0: return two identical roots,
171 // and one distinct real root
172 // if identical zero Hessian: return three identical real roots
173 // where the discriminant = - 4p^3 - 27q^2.
174 Roots<3> roots() const;
175};
176
177
178// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
179
180} // End namespace Foam
181
182// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
183
184#include "cubicEqnI.H"
185
186// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
187
188#endif
189
190// ************************************************************************* //
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
Definition: Roots.H:73
Templated vector space.
Definition: VectorSpace.H:79
Container to encapsulate various operations for cubic equation of the forms with real coefficients:
Definition: cubicEqn.H:115
cubicEqn()
Construct null.
Definition: cubicEqnI.H:30
components
Component labeling enumeration.
Definition: cubicEqn.H:119
scalar d() const
Definition: cubicEqnI.H:75
scalar value(const scalar x) const
Evaluate the cubic equation at x.
Definition: cubicEqnI.H:105
scalar c() const
Definition: cubicEqnI.H:69
scalar derivative(const scalar x) const
Evaluate the derivative of the cubic equation at x.
Definition: cubicEqnI.H:111
Roots< 3 > roots() const
Return the roots of the cubic equation with no particular order.
Definition: cubicEqn.C:35
scalar a() const
Definition: cubicEqnI.H:57
scalar b() const
Definition: cubicEqnI.H:63
Class to handle errors and exceptions in a simple, consistent stream-based manner.
Definition: error.H:77
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
Definition: zero.H:63
Namespace for OpenFOAM.