diff --git a/MathNet.Numerics.sln.DotSettings b/MathNet.Numerics.sln.DotSettings index 104d6356..ba81b902 100644 --- a/MathNet.Numerics.sln.DotSettings +++ b/MathNet.Numerics.sln.DotSettings @@ -1,9 +1,20 @@  + True True + SOLUTION + DoShow + DoShow + DoShow + DoHide + DoShow + DoShow + DoShow + DoShow + DoShow HINT HINT <?xml version="1.0" encoding="utf-16"?><Profile name="Full Cleanup (Math.NET)"><CSArrangeThisQualifier>True</CSArrangeThisQualifier><CSRemoveCodeRedundancies>True</CSRemoveCodeRedundancies><CSUseAutoProperty>True</CSUseAutoProperty><CSMakeFieldReadonly>True</CSMakeFieldReadonly><CSUpdateFileHeader>True</CSUpdateFileHeader><CSOptimizeUsings><OptimizeUsings>True</OptimizeUsings><EmbraceInRegion>False</EmbraceInRegion><RegionName></RegionName></CSOptimizeUsings><CSShortenReferences>True</CSShortenReferences><CSReformatCode>True</CSReformatCode><XMLReformatCode>True</XMLReformatCode><CssAlphabetizeProperties>True</CssAlphabetizeProperties><CssReformatCode>True</CssReformatCode><JsReformatCode>True</JsReformatCode><JsInsertSemicolon>True</JsInsertSemicolon><VBFormatDocComments>True</VBFormatDocComments><VBReformatCode>True</VBReformatCode><VBShortenReferences>True</VBShortenReferences><VBOptimizeImports>True</VBOptimizeImports><HtmlReformatCode>True</HtmlReformatCode><AspOptimizeRegisterDirectives>True</AspOptimizeRegisterDirectives><CSReorderTypeMembers>True</CSReorderTypeMembers><CSUseVar><BehavourStyle>CAN_CHANGE_BOTH</BehavourStyle><LocalVariableStyle>IMPLICIT_WHEN_INITIALIZER_HAS_TYPE</LocalVariableStyle><ForeachVariableStyle>ALWAYS_EXPLICIT</ForeachVariableStyle></CSUseVar></Profile> - Full Cleanup (Math.NET) + Default: Reformat Code False False False @@ -18,7 +29,20 @@ True False False + True + False + False False + True + False + True + True + True + True + False + True + True + True True True False @@ -51,6 +75,7 @@ WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. </copyright> + True AVX CDF DFT @@ -79,4 +104,5 @@ OTHER DEALINGS IN THE SOFTWARE. WH True <data /> - <data><IncludeFilters /><ExcludeFilters /></data> \ No newline at end of file + <data><IncludeFilters /><ExcludeFilters /></data> + \ No newline at end of file diff --git a/src/Numerics/FindRoots.cs b/src/Numerics/FindRoots.cs index 5dfeec0f..f02cf07a 100644 --- a/src/Numerics/FindRoots.cs +++ b/src/Numerics/FindRoots.cs @@ -50,6 +50,11 @@ namespace MathNet.Numerics { double root; + if (!ZeroCrossingBracketing.Expand(f, ref lowerBound, ref upperBound, 1.6, 100)) + { + throw new NonConvergenceException(Resources.RootFindingFailed); + } + if (Brent.TryFindRoot(f, lowerBound, upperBound, accuracy, maxIterations, out root)) { return root; @@ -79,12 +84,7 @@ namespace MathNet.Numerics return root; } - if (Bisection.TryFindRoot(f, lowerBound, upperBound, accuracy, maxIterations, out root)) - { - return root; - } - - throw new NonConvergenceException(Resources.RootFindingFailed); + return OfFunction(f, lowerBound, upperBound, accuracy, maxIterations); } /// diff --git a/src/Numerics/RootFinding/Bisection.cs b/src/Numerics/RootFinding/Bisection.cs index 0539bde2..384e5786 100644 --- a/src/Numerics/RootFinding/Bisection.cs +++ b/src/Numerics/RootFinding/Bisection.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2013 Math.NET -// +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -34,7 +34,7 @@ using MathNet.Numerics.Properties; namespace MathNet.Numerics.RootFinding { /// - /// Bisection root-finding algorithm without any recovery measures in case of lacking bracketing. + /// Bisection root-finding algorithm. /// public static class Bisection { diff --git a/src/Numerics/RootFinding/Brent.cs b/src/Numerics/RootFinding/Brent.cs index a2bb24e8..f94b748f 100644 --- a/src/Numerics/RootFinding/Brent.cs +++ b/src/Numerics/RootFinding/Brent.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2013 Math.NET -// +// +// Copyright (c) 2009-2014 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -39,6 +39,22 @@ namespace MathNet.Numerics.RootFinding /// public static class Brent { + /// Find a solution of the equation f(x)=0. + /// The function to find roots from. + /// Guess for the low value of the range where the root is supposed to be. Will be expanded if needed. + /// Guess for the high value of the range where the root is supposed to be. Will be expanded if needed. + /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Default 1e-8. + /// Maximum number of iterations. Default 100. + /// Factor at which to expand the bounds, if needed. Default 1.6. + /// Maximum number of expand iterations. Default 100. + /// Returns the root with the specified accuracy. + /// + public static double FindRootExpand(Func f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, int maxIterations = 100, double expandFactor = 1.6, int maxExpandIteratons = 100) + { + ZeroCrossingBracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons); + return FindRoot(f, guessLowerBound, guessUpperBound, accuracy, maxIterations); + } + /// Find a solution of the equation f(x)=0. /// The function to find roots from. /// The low value of the range where the root is supposed to be. @@ -103,9 +119,9 @@ namespace MathNet.Numerics.RootFinding } // convergence check - double xAcc1 = 2.0 * Precision.DoublePrecision * Math.Abs(root) + 0.5 * accuracy; + double xAcc1 = 2.0*Precision.DoublePrecision*Math.Abs(root) + 0.5*accuracy; double xMidOld = xMid; - xMid = (upperBound - root) / 2.0; + xMid = (upperBound - root)/2.0; if (Math.Abs(xMid) <= xAcc1 || froot.AlmostEqualNormRelative(0, froot, accuracy)) { @@ -121,20 +137,20 @@ namespace MathNet.Numerics.RootFinding if (Math.Abs(e) >= xAcc1 && Math.Abs(fmin) > Math.Abs(froot)) { // Attempt inverse quadratic interpolation - double s = froot / fmin; + double s = froot/fmin; double p; double q; if (lowerBound.AlmostEqualRelative(upperBound)) { - p = 2.0 * xMid * s; + p = 2.0*xMid*s; q = 1.0 - s; } else { - q = fmin / fmax; - double r = froot / fmax; - p = s * (2.0 * xMid * q * (q - r) - (root - lowerBound) * (r - 1.0)); - q = (q - 1.0) * (r - 1.0) * (s - 1.0); + q = fmin/fmax; + double r = froot/fmax; + p = s*(2.0*xMid*q*(q - r) - (root - lowerBound)*(r - 1.0)); + q = (q - 1.0)*(r - 1.0)*(s - 1.0); } if (p > 0.0) @@ -144,11 +160,11 @@ namespace MathNet.Numerics.RootFinding } p = Math.Abs(p); - if (2.0 * p < Math.Min(3.0 * xMid * q - Math.Abs(xAcc1 * q), Math.Abs(e * q))) + if (2.0*p < Math.Min(3.0*xMid*q - Math.Abs(xAcc1*q), Math.Abs(e*q))) { // Accept interpolation e = d; - d = p / q; + d = p/q; } else { diff --git a/src/Numerics/RootFinding/ZeroCrossingBracketing.cs b/src/Numerics/RootFinding/ZeroCrossingBracketing.cs index aa23b477..d07af4d2 100644 --- a/src/Numerics/RootFinding/ZeroCrossingBracketing.cs +++ b/src/Numerics/RootFinding/ZeroCrossingBracketing.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2013 Math.NET -// +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -100,12 +100,12 @@ namespace MathNet.Numerics.RootFinding if (Math.Abs(fmin) < Math.Abs(fmax)) { - lowerBound += factor * (lowerBound - upperBound); + lowerBound += factor*(lowerBound - upperBound); fmin = f(lowerBound); } else { - upperBound += factor * (upperBound - lowerBound); + upperBound += factor*(upperBound - lowerBound); fmax = f(upperBound); } } diff --git a/src/UnitTests/RootFindingTests/BisectionTest.cs b/src/UnitTests/RootFindingTests/BisectionTest.cs index 7a7edc97..ee58b44a 100644 --- a/src/UnitTests/RootFindingTests/BisectionTest.cs +++ b/src/UnitTests/RootFindingTests/BisectionTest.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2013 Math.NET -// +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -41,7 +41,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void MultipleRoots() { // Roots at -2, 2 - Func f1 = x => x * x - 4; + Func f1 = x => x*x - 4; Assert.AreEqual(0, f1(Bisection.FindRoot(f1, 0, 5, 1e-14))); Assert.AreEqual(0, f1(Bisection.FindRootExpand(f1, 3, 4, 1e-14))); Assert.AreEqual(-2, Bisection.FindRoot(f1, -5, -1, 1e-14)); @@ -51,7 +51,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests Assert.AreEqual(2, Bisection.FindRoot(x => -f1(x), 1, 4, 1e-14)); // Roots at 3, 4 - Func f2 = x => (x - 3) * (x - 4); + Func f2 = x => (x - 3)*(x - 4); Assert.AreEqual(0, f2(Bisection.FindRoot(f2, 3.5, 5, 1e-14)), 1e-14); Assert.AreEqual(3, Bisection.FindRoot(f2, -5, 3.5, 1e-14)); Assert.AreEqual(4, Bisection.FindRoot(f2, 3.2, 5, 1e-14)); @@ -62,7 +62,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void LocalMinima() { - Func f1 = x => x * x * x - 2 * x + 2; + Func f1 = x => x*x*x - 2*x + 2; Assert.AreEqual(0, f1(Bisection.FindRoot(f1, -5, 5, 1e-14)), 1e-14); Assert.AreEqual(0, f1(Bisection.FindRoot(f1, -2, 4, 1e-14)), 1e-14); } @@ -70,7 +70,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void NoRoot() { - Func f1 = x => x * x + 4; + Func f1 = x => x*x + 4; Assert.That(() => Bisection.FindRoot(f1, -5, 5, 1e-14), Throws.TypeOf()); } @@ -78,7 +78,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq1() { // Test case from http://www.polymath-software.com/library/nle/Oneeq1.htm - Func f1 = z => 8 * Math.Pow((4 - z) * z, 2) / (Math.Pow(6 - 3 * z, 2) * (2 - z)) - 0.186; + Func f1 = z => 8*Math.Pow((4 - z)*z, 2)/(Math.Pow(6 - 3*z, 2)*(2 - z)) - 0.186; double x = Bisection.FindRoot(f1, 0.1, 0.9, accuracy: 1e-9, maxIterations: 80); Assert.AreEqual(0.277759543089215, x, 1e-9); Assert.AreEqual(0, f1(x), 1e-16); @@ -94,15 +94,15 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double x2 = 1 - x1; const double G12 = 0.07858889; const double G21 = 0.30175355; - double P2 = Math.Pow(10, 6.87776 - 1171.53 / (224.366 + T)); - double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + T)); - const double t1 = x1 + x2 * G12; - const double t2 = x2 + x1 * G21; - double gamma2 = Math.Exp(-Math.Log(t2) - (x1 * (G12 * t2 - G21 * t1)) / (t1 * t2)); - double gamma1 = Math.Exp(-Math.Log(t1) + (x2 * (G12 * t2 - G21 * t1)) / (t1 * t2)); - double k1 = gamma1 * P1 / 760; - double k2 = gamma2 * P2 / 760; - return 1 - k1 * x1 - k2 * x2; + double P2 = Math.Pow(10, 6.87776 - 1171.53/(224.366 + T)); + double P1 = Math.Pow(10, 8.04494 - 1554.3/(222.65 + T)); + const double t1 = x1 + x2*G12; + const double t2 = x2 + x1*G21; + double gamma2 = Math.Exp(-Math.Log(t2) - (x1*(G12*t2 - G21*t1))/(t1*t2)); + double gamma1 = Math.Exp(-Math.Log(t1) + (x2*(G12*t2 - G21*t1))/(t1*t2)); + double k1 = gamma1*P1/760; + double k2 = gamma2*P2/760; + return 1 - k1*x1 - k2*x2; }; double x = Bisection.FindRoot(f1, 0, 100); @@ -120,13 +120,13 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double x2 = 1 - x1; const double A = 0.75807689; const double B = 1.1249035; - double P2 = Math.Pow(10, 6.9024 - 1268.115 / (216.9 + T)); - double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + T)); - double gamma2 = Math.Pow(10, B * Math.Pow(x1, 2) / Math.Pow(x1 + B * x2 / A, 2)); - double gamma1 = Math.Pow(10, A * Math.Pow(x2, 2) / Math.Pow(A * x1 / B + x2, 2)); - double k1 = gamma1 * P1 / 760; - double k2 = gamma2 * P2 / 760; - return 1 - k1 * x1 - k2 * x2; + double P2 = Math.Pow(10, 6.9024 - 1268.115/(216.9 + T)); + double P1 = Math.Pow(10, 8.04494 - 1554.3/(222.65 + T)); + double gamma2 = Math.Pow(10, B*Math.Pow(x1, 2)/Math.Pow(x1 + B*x2/A, 2)); + double gamma1 = Math.Pow(10, A*Math.Pow(x2, 2)/Math.Pow(A*x1/B + x2, 2)); + double k1 = gamma1*P1/760; + double k2 = gamma2*P2/760; + return 1 - k1*x1 - k2*x2; }; double x = Bisection.FindRoot(f1, 0, 100); @@ -138,7 +138,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq3() { // Test case from http://www.polymath-software.com/library/nle/Oneeq3.htm - Func f1 = T => Math.Exp(21000 / T) / (T * T) - 1.11e11; + Func f1 = T => Math.Exp(21000/T)/(T*T) - 1.11e11; double x = Bisection.FindRoot(f1, 550, 560, 1e-2); Assert.AreEqual(551.773822885233, x, 1e-5); Assert.AreEqual(0, f1(x), 1e-2); @@ -152,7 +152,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests { const double A = 0.38969; const double B = 0.55954; - return A * B * (B * Math.Pow(1 - x, 2) - A * x * x) / Math.Pow(x * (A - B) + B, 2) + 0.14845; + return A*B*(B*Math.Pow(1 - x, 2) - A*x*x)/Math.Pow(x*(A - B) + B, 2) + 0.14845; }; double r = Bisection.FindRoot(f1, 0, 1); @@ -170,7 +170,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double BETA = 2.298E-3; const double GAMMA = 0.283E-6; const double DH = -57798.0; - return DH + TR * (ALPHA + TR * (BETA / 2 + TR * GAMMA / 3)) - 298.0 * (ALPHA + 298.0 * (BETA / 2 + 298.0 * GAMMA / 3)); + return DH + TR*(ALPHA + TR*(BETA/2 + TR*GAMMA/3)) - 298.0*(ALPHA + 298.0*(BETA/2 + 298.0*GAMMA/3)); }; double x = Bisection.FindRoot(f1, 3000, 5000); @@ -192,17 +192,17 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double A = 0.01855; const double B = -0.01587; const double P = 100.0; - const double Beta = R * T * B0 - A0 - R * C / (T * T); - const double Gama = -R * T * B0 * B + A0 * A - R * C * B0 / (T * T); - const double Delta = R * B0 * B * C / (T * T); + const double Beta = R*T*B0 - A0 - R*C/(T*T); + const double Gama = -R*T*B0*B + A0*A - R*C*B0/(T*T); + const double Delta = R*B0*B*C/(T*T); - return R * T / V + Beta / (V * V) + Gama / (V * V * V) + Delta / (V * V * V * V) - P; + return R*T/V + Beta/(V*V) + Gama/(V*V*V) + Delta/(V*V*V*V) - P; }; double x = Bisection.FindRoot(f1, 0.1, 1); Assert.AreEqual(0.174749531708621, x, 1e-5); Assert.AreEqual(0, f1(x), 1e-12); - } + } [Test] public void Oneeq6b() @@ -218,23 +218,23 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double A = 0.01855; const double B = -0.01587; const double P = 100.0; - const double Beta = R * T * B0 - A0 - R * C / (T * T); - const double Gama = -R * T * B0 * B + A0 * A - R * C * B0 / (T * T); - const double Delta = R * B0 * B * C / (T * T); + const double Beta = R*T*B0 - A0 - R*C/(T*T); + const double Gama = -R*T*B0*B + A0*A - R*C*B0/(T*T); + const double Delta = R*B0*B*C/(T*T); - return R * T * V * V * V + Beta * V * V + Gama * V + Delta - P * V * V * V * V; + return R*T*V*V*V + Beta*V*V + Gama*V + Delta - P*V*V*V*V; }; double x = Bisection.FindRoot(f1, 0.1, 1); Assert.AreEqual(0.174749600398905, x, 1e-5); Assert.AreEqual(0, f1(x), 1e-13); - } - + } + [Test] public void Oneeq7() { // Test case from http://www.polymath-software.com/library/nle/Oneeq7.htm - Func f1 = x => x / (1 - x) - 5 * Math.Log(0.4 * (1 - x) / (0.4 - 0.5 * x)) + 4.45977; + Func f1 = x => x/(1 - x) - 5*Math.Log(0.4*(1 - x)/(0.4 - 0.5*x)) + 4.45977; double r = Bisection.FindRoot(f1, 0, 0.79, 1e-2); Assert.AreEqual(0.757396293891, r, 1e-5); Assert.AreEqual(0, f1(r), 1e-13); @@ -250,7 +250,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double b = 40; const double c = 200; - return a * v * v + b * Math.Pow(v, 7 / 4) - c; + return a*v*v + b*Math.Pow(v, 7/4) - c; }; double x = Bisection.FindRoot(f1, 0.01, 1); @@ -269,7 +269,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double gama = 1.41; const double Pd = 100; - return Math.Pow((gama + 1) / 2, (gama + 1) / (gama - 1)) * (2 / (gama - 1)) * (Math.Pow(P / Pd, 2 / gama) - Math.Pow(P / Pd, (gama + 1) / gama)) - Math.Pow(At / A, 2); + return Math.Pow((gama + 1)/2, (gama + 1)/(gama - 1))*(2/(gama - 1))*(Math.Pow(P/Pd, 2/gama) - Math.Pow(P/Pd, (gama + 1)/gama)) - Math.Pow(At/A, 2); }; double x = Bisection.FindRoot(f1, 1, 50); @@ -284,7 +284,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq10() { // Test case from http://www.polymath-software.com/library/nle/Oneeq10.htm - Func f1 = x => (1.0 / 63.0) * Math.Log(x) + (64.0 / 63.0) * Math.Log(1 / (1 - x)) + Math.Log(0.95 - x) - Math.Log(0.9); + Func f1 = x => (1.0/63.0)*Math.Log(x) + (64.0/63.0)*Math.Log(1/(1 - x)) + Math.Log(0.95 - x) - Math.Log(0.9); double r = Bisection.FindRoot(f1, .01, 0.35); Assert.AreEqual(0.036210083704, r, 1e-5); @@ -306,7 +306,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double kp = 604500; const double P = 0.00243; - return a - Math.Pow((c + 2 * x) * (a + b + c - 2 * x), 2) / (kp * P * P * Math.Pow(b - 3 * x, 3)) - x; + return a - Math.Pow((c + 2*x)*(a + b + c - 2*x), 2)/(kp*P*P*Math.Pow(b - 3*x, 3)) - x; }; double r = Bisection.FindRoot(f1, 0, 0.25); @@ -329,13 +329,13 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double kp = 604500; const double P = 0.00243; - return (b - Math.Pow(Math.Pow((c + 2 * x) * (a + b + c - 2 * x), 2) / (kp * P * P * (a - x)) , 1.0 / 3.0)) / 3.0 - x; + return (b - Math.Pow(Math.Pow((c + 2*x)*(a + b + c - 2*x), 2)/(kp*P*P*(a - x)), 1.0/3.0))/3.0 - x; }; double r = Bisection.FindRoot(f1, 0.01, 0.45); Assert.AreEqual(0.058654571042804, r, 1e-5); Assert.AreEqual(0, f1(r), 1e-14); - // Could not test the following root, since Math.Pow(y,1/3) does not work for y<0 + // Could not test the following root, since Math.Pow(y,1/3) does not work for y<0 // r = Bisection.FindRoot(f1, 0.55, 0.7); // Assert.AreEqual(0.600323117488527, r, 1e-5); // Assert.AreEqual(0, f1(r), 1e-14); @@ -353,7 +353,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double T = 430.85; const double R = 82.05; - return (R * T / P) * (1 + b / v + c / (v * v)) - v; + return (R*T/P)*(1 + b/v + c/(v*v)) - v; }; double x = Bisection.FindRoot(f1, 100, 300); @@ -372,10 +372,10 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double P = 100; const double T = 210.0; const double R = 0.08205; - double A = 0.42748 * (R * R * Math.Pow(Tc, 2.5)) / Pc; - const double B = 0.08664 * R * Tc / Pc; + double A = 0.42748*(R*R*Math.Pow(Tc, 2.5))/Pc; + const double B = 0.08664*R*Tc/Pc; - return R * Math.Pow(T, 1.5) * V * (V + B) / (P * Math.Sqrt(T) * V * (V + B) + A) + B - V; + return R*Math.Pow(T, 1.5)*V*(V + B)/(P*Math.Sqrt(T)*V*(V + B) + A) + B - V; }; double x = Bisection.FindRoot(f1, .01, .3); @@ -397,10 +397,10 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double T1 = 390.0; const double t1 = 100.0; const double cp = 0.49; - double DT1 = T1 - t1 - Q / (M * Cp); - double DT2 = T1 - t1 - Q / (m * cp); + double DT1 = T1 - t1 - Q/(M*Cp); + double DT2 = T1 - t1 - Q/(m*cp); - return U * A * (DT2 - DT1) / Math.Log(DT2 / DT1) - Q; + return U*A*(DT2 - DT1)/Math.Log(DT2/DT1) - Q; }; double x = Bisection.FindRoot(f1, 1000000, 7000000); @@ -412,7 +412,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq15a() { // Test case from http://www.polymath-software.com/library/nle/Oneeq15a.htm - Func f1 = h => h * (1 + h * (3 - h)) - 2.4 - h; + Func f1 = h => h*(1 + h*(3 - h)) - 2.4 - h; double x = Bisection.FindRoot(f1, -2, 0); Assert.AreEqual(-0.795219754733581, x, 1e-5); @@ -429,7 +429,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq15b() { // Test case from http://www.polymath-software.com/library/nle/Oneeq15b.htm - Func f1 = h => 2.4 / (h * (3 - h)) - h; + Func f1 = h => 2.4/(h*(3 - h)) - h; double x = Bisection.FindRoot(f1, -2, -0.5); Assert.AreEqual(-0.795219745444464, x, 1e-5); @@ -453,18 +453,18 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double z = 1 - y - 0.02; const double CH4 = 0; const double C2H6 = 0; - const double COtwo = y + 2 * z; - const double H2O = 2 * y + 3 * z; - const double N2 = 0.02 + 3.76 * (2 * y + 7 * z / 2) * x; - const double O2 = (2 * y + 7 * z / 2) * (x - 1); - const double alp = 3.381 * CH4 + 2.247 * C2H6 + 6.214 * COtwo + 7.256 * H2O + 6.524 * N2 + 6.148 * O2; - const double bet = 18.044 * CH4 + 38.201 * C2H6 + 10.396 * COtwo + 2.298 * H2O + 1.25 * N2 + 3.102 * O2; - const double gam = -4.3 * CH4 - 11.049 * C2H6 - 3.545 * COtwo + 0.283 * H2O - 0.001 * N2 - 0.923 * O2; - const double H0 = alp * 298 + bet * 0.001 * 298 * 298 / 2 + gam * 1e-6 * 298 * 298 * 298 / 3; - double Hf = alp * T + bet * 0.001 * T * T / 2 + gam * 1e-6 * T * T * T / 3; + const double COtwo = y + 2*z; + const double H2O = 2*y + 3*z; + const double N2 = 0.02 + 3.76*(2*y + 7*z/2)*x; + const double O2 = (2*y + 7*z/2)*(x - 1); + const double alp = 3.381*CH4 + 2.247*C2H6 + 6.214*COtwo + 7.256*H2O + 6.524*N2 + 6.148*O2; + const double bet = 18.044*CH4 + 38.201*C2H6 + 10.396*COtwo + 2.298*H2O + 1.25*N2 + 3.102*O2; + const double gam = -4.3*CH4 - 11.049*C2H6 - 3.545*COtwo + 0.283*H2O - 0.001*N2 - 0.923*O2; + const double H0 = alp*298 + bet*0.001*298*298/2 + gam*1e-6*298*298*298/3; + double Hf = alp*T + bet*0.001*T*T/2 + gam*1e-6*T*T*T/3; const double xx = 1; - return 212798 * y * xx + 372820 * z * xx + H0 - Hf; + return 212798*y*xx + 372820*z*xx + H0 - Hf; }; double r = Bisection.FindRoot(f1, 1000, 3000); @@ -476,7 +476,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq17() { // Test case from http://www.polymath-software.com/library/nle/Oneeq17.htm - Func f1 = rp => rp - 0.327 * Math.Pow(0.06 - 161 * rp, 0.804) * Math.Exp(-5230 / (1.987 * (373 + 1.84e6 * rp))); + Func f1 = rp => rp - 0.327*Math.Pow(0.06 - 161*rp, 0.804)*Math.Exp(-5230/(1.987*(373 + 1.84e6*rp))); double x = Bisection.FindRoot(f1, 0, 0.00035); Assert.AreEqual(0.000340568862275, x, 1e-5); @@ -489,12 +489,12 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Test case from http://www.polymath-software.com/library/nle/Oneeq18a.htm Func f1 = T => { - double Tp = (269.267 * T - 1752) / 266.667; - double k = 0.0006 * Math.Exp(20.7 - 15000 / Tp); - double Pp = -(0.1 / (321 * (320.0 / 321.0 - (1 + k)))); - double P = (320 * Pp + 0.1) / 321; + double Tp = (269.267*T - 1752)/266.667; + double k = 0.0006*Math.Exp(20.7 - 15000/Tp); + double Pp = -(0.1/(321*(320.0/321.0 - (1 + k)))); + double P = (320*Pp + 0.1)/321; - return 1.296 * (T - Tp) + 10369 * k * Pp; + return 1.296*(T - Tp) + 10369*k*Pp; }; double x = Bisection.FindRoot(f1, 500, 725); @@ -514,12 +514,12 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Test case from http://www.polymath-software.com/library/nle/Oneeq18b.htm Func f1 = T => { - double Tp = (269 * T - 1752) / 267; - double k = 0.0006 * Math.Exp(20.7 - 15000 / Tp); - double Pp = -(0.1 / (321 * (320.0 / 321.0 - (1 + k)))); - double P = (320 * Pp + 0.1) / 321; + double Tp = (269*T - 1752)/267; + double k = 0.0006*Math.Exp(20.7 - 15000/Tp); + double Pp = -(0.1/(321*(320.0/321.0 - (1 + k)))); + double P = (320*Pp + 0.1)/321; - return 1.3 * (T - Tp) + 1.04e4 * k * Pp; + return 1.3*(T - Tp) + 1.04e4*k*Pp; }; double x = Bisection.FindRoot(f1, 500, 1500); @@ -537,25 +537,25 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double Pc = 33.5; const double T = 631*2; const double Tc = 126.2; - const double Pr = P / Pc; - const double Tr = T / Tc; - double Asqr = 0.4278 * Pr / Math.Pow(Tr, 2.5); - const double B = 0.0867 * Pr / Tr; - double Q = B * B + B - Asqr; - double r = Asqr * B; - double p = (-3 * Q - 1) / 3; - double q = (-27 * r - 9 * Q - 2) / 27; - double R = Math.Pow(p / 3, 3) + Math.Pow(q / 2, 2); - double V = (R > 0) ? Math.Pow(-q / 2 + Math.Sqrt(R), 1 / 3) : 0; - double WW = (R > 0) ? (-q / 2 - Math.Sqrt(R)) : (0); - double psi1 = (R < 0) ? (Math.Acos(Math.Sqrt((q * q / 4) / (-p * p * p / 27)))) : (0); - double W = (R > 0) ? (Math.Sign(WW) * Math.Pow(Math.Abs(WW), 1 / 3)) : (0); - double z1 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 0 / 3) + 1 / 3) : (0); - double z2 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 1 / 3) + 1 / 3) : (0); - double z3 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 2 / 3) + 1 / 3) : (0); - double z0 = (R > 0) ? (V + W + 1 / 3) : (0); - - return z * z * z - z * z - Q * z - r; + const double Pr = P/Pc; + const double Tr = T/Tc; + double Asqr = 0.4278*Pr/Math.Pow(Tr, 2.5); + const double B = 0.0867*Pr/Tr; + double Q = B*B + B - Asqr; + double r = Asqr*B; + double p = (-3*Q - 1)/3; + double q = (-27*r - 9*Q - 2)/27; + double R = Math.Pow(p/3, 3) + Math.Pow(q/2, 2); + double V = (R > 0) ? Math.Pow(-q/2 + Math.Sqrt(R), 1/3) : 0; + double WW = (R > 0) ? (-q/2 - Math.Sqrt(R)) : (0); + double psi1 = (R < 0) ? (Math.Acos(Math.Sqrt((q*q/4)/(-p*p*p/27)))) : (0); + double W = (R > 0) ? (Math.Sign(WW)*Math.Pow(Math.Abs(WW), 1/3)) : (0); + double z1 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*0/3) + 1/3) : (0); + double z2 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*1/3) + 1/3) : (0); + double z3 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*2/3) + 1/3) : (0); + double z0 = (R > 0) ? (V + W + 1/3) : (0); + + return z*z*z - z*z - Q*z - r; }; double x = Bisection.FindRoot(f1, -0.5, -0.02); @@ -579,25 +579,25 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double Pc = 33.5; const double T = 631; const double Tc = 126.2; - const double Pr = P / Pc; - const double Tr = T / Tc; - double Asqr = 0.4278 * Pr / Math.Pow(Tr, 2.5); - const double B = 0.0867 * Pr / Tr; - double Q = B * B + B - Asqr; - double r = Asqr * B; - double p = (-3 * Q - 1) / 3; - double q = (-27 * r - 9 * Q - 2) / 27; - double R = Math.Pow(p / 3, 3) + Math.Pow(q / 2, 2); - double V = (R > 0) ? Math.Pow(-q / 2 + Math.Sqrt(R), 1 / 3) : 0; - double WW = (R > 0) ? (-q / 2 - Math.Sqrt(R)) : (0); - double psi1 = (R < 0) ? (Math.Acos(Math.Sqrt((q * q / 4) / (-p * p * p / 27)))) : (0); - double W = (R > 0) ? (Math.Sign(WW) * Math.Pow(Math.Abs(WW), 1 / 3)) : (0); - double z1 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 0 / 3) + 1 / 3) : (0); - double z2 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 1 / 3) + 1 / 3) : (0); - double z3 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 2 / 3) + 1 / 3) : (0); - double z0 = (R > 0) ? (V + W + 1 / 3) : (0); - - return z * z * z - z * z - Q * z - r; + const double Pr = P/Pc; + const double Tr = T/Tc; + double Asqr = 0.4278*Pr/Math.Pow(Tr, 2.5); + const double B = 0.0867*Pr/Tr; + double Q = B*B + B - Asqr; + double r = Asqr*B; + double p = (-3*Q - 1)/3; + double q = (-27*r - 9*Q - 2)/27; + double R = Math.Pow(p/3, 3) + Math.Pow(q/2, 2); + double V = (R > 0) ? Math.Pow(-q/2 + Math.Sqrt(R), 1/3) : 0; + double WW = (R > 0) ? (-q/2 - Math.Sqrt(R)) : (0); + double psi1 = (R < 0) ? (Math.Acos(Math.Sqrt((q*q/4)/(-p*p*p/27)))) : (0); + double W = (R > 0) ? (Math.Sign(WW)*Math.Pow(Math.Abs(WW), 1/3)) : (0); + double z1 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*0/3) + 1/3) : (0); + double z2 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*1/3) + 1/3) : (0); + double z3 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*2/3) + 1/3) : (0); + double z0 = (R > 0) ? (V + W + 1/3) : (0); + + return z*z*z - z*z - Q*z - r; }; double x = Bisection.FindRoot(f1, -0.5, 1.2); @@ -613,14 +613,14 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests { const double T = 313; const double P0 = 10; - const double FA0 = 20 * P0 / (0.082 * 450); - double k1 = 1.277 * 1.0e9 * Math.Exp(-90000 / (8.31 * T)); - double k2 = 1.29 * 1.0e11 * Math.Exp(-135000 / (8.31 * T)); - double den = 1 + 7 * P0 * (1 - xa) / (1 + xa); + const double FA0 = 20*P0/(0.082*450); + double k1 = 1.277*1.0e9*Math.Exp(-90000/(8.31*T)); + double k2 = 1.29*1.0e11*Math.Exp(-135000/(8.31*T)); + double den = 1 + 7*P0*(1 - xa)/(1 + xa); double xa1 = 1 + xa; - double ra = (-k1 * P0 * (1 - xa) / xa1 + k2 * P0 * P0 * xa * xa / (xa1 * xa1)) / den; + double ra = (-k1*P0*(1 - xa)/xa1 + k2*P0*P0*xa*xa/(xa1*xa1))/den; - return -ra / FA0; + return -ra/FA0; }; double x = Bisection.FindRoot(f1, 0.75, 1.02); @@ -636,13 +636,13 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests { const double T = 313; const double P0 = 10; - const double FA0 = 20 * P0 / (0.082 * 450); - double k1 = 1.277 * 1.0e9 * Math.Exp(-90000 / (8.31 * T)); - double k2 = 1.29 * 1.0e11 * Math.Exp(-135000 / (8.31 * T)); + const double FA0 = 20*P0/(0.082*450); + double k1 = 1.277*1.0e9*Math.Exp(-90000/(8.31*T)); + double k2 = 1.29*1.0e11*Math.Exp(-135000/(8.31*T)); double xa1 = 1 + xa; - double ra = (-k1 * P0 * (1 - xa) / xa1 + k2 * P0 * P0 * xa * xa / (xa1 * xa1)); + double ra = (-k1*P0*(1 - xa)/xa1 + k2*P0*P0*xa*xa/(xa1*xa1)); - return -ra / FA0; + return -ra/FA0; }; double x = Bisection.FindRoot(f1, 0.75, 1.02); @@ -656,37 +656,37 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Test case from http://www.polymath-software.com/library/nle/Oneeq21a.htm Func f1 = h => { - double sg = Math.Acos(2 * h - 1); - double si = Math.Pow(1 - Math.Pow(2 * h - 1, 2), 0.5); - double ag = (sg - (2 * h - 1) * si) / 4; - double al = 3.1415927 / 4 - ag; + double sg = Math.Acos(2*h - 1); + double si = Math.Pow(1 - Math.Pow(2*h - 1, 2), 0.5); + double ag = (sg - (2*h - 1)*si)/4; + double al = 3.1415927/4 - ag; const double d = 1.44; - double dg = 4 * ag / (sg + si) * d; + double dg = 4*ag/(sg + si)*d; const double ugs = -0.0295; - double ug = ugs * 3.1415927 / 4 / ag; + double ug = ugs*3.1415927/4/ag; const double uls = -0.00001; - double ul = uls * 3.1415927 / 4 / al; + double ul = uls*3.1415927/4/al; double sl = 3.1415927 - sg; - double dl = 4 * al / (sl) * d; + double dl = 4*al/(sl)*d; const double rol = 1; const double rog = 0.835; const double bt = .6; - double g = 980 * Math.Sin(bt * 3.1415927 / 180); - double dpg = (rol * al + rog * ag) * g / 3.1415927 * 4; + double g = 980*Math.Sin(bt*3.1415927/180); + double dpg = (rol*al + rog*ag)*g/3.1415927*4; const double vig = 1.0; - double reg = Math.Abs(ug) * rog * dg / vig; - double tg = -16 / reg * (.5 * rog * ug * Math.Abs(ug)); + double reg = Math.Abs(ug)*rog*dg/vig; + double tg = -16/reg*(.5*rog*ug*Math.Abs(ug)); const double vil = .01; - double rel = Math.Abs(ul) * rol * dl / vil; - double tl = -16 / rel * (.5 * rol * ul * Math.Abs(ul)); + double rel = Math.Abs(ul)*rol*dl/vil; + double tl = -16/rel*(.5*rol*ul*Math.Abs(ul)); const double it = 1; - double ti = -16 / reg * rog * (.5 * (ug - ul) * Math.Abs(ug - ul)) * it; - double dpf = (tl * sl + tg * sg) * d / (al + ag) / (d * d) * 4; + double ti = -16/reg*rog*(.5*(ug - ul)*Math.Abs(ug - ul))*it; + double dpf = (tl*sl + tg*sg)*d/(al + ag)/(d*d)*4; double dp = dpg + dpf; - double tic = (Math.Abs(ul) < Math.Abs(ug)) ? (rog / reg) : (rol / rel); - double y = (rol - rog) * g * Math.Pow(d / 2, 2) / (8 * vig * ugs); + double tic = (Math.Abs(ul) < Math.Abs(ug)) ? (rog/reg) : (rol/rel); + double y = (rol - rog)*g*Math.Pow(d/2, 2)/(8*vig*ugs); - return -tg * sg * al + tl * sl * ag - ti * si * (al + ag) + d * (rol - rog) * al * ag * g; + return -tg*sg*al + tl*sl*ag - ti*si*(al + ag) + d*(rol - rog)*al*ag*g; }; double x = Bisection.FindRoot(f1, 0, 0.3); @@ -700,37 +700,37 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Test case from http://www.polymath-software.com/library/nle/Oneeq21b.htm Func f1 = h => { - double sg = Math.Acos(2 * h - 1); - double si = Math.Pow(1 - Math.Pow(2 * h - 1, 2), 0.5); - double ag = (sg - (2 * h - 1) * si) / 4; - double al = 3.1415927 / 4 - ag; + double sg = Math.Acos(2*h - 1); + double si = Math.Pow(1 - Math.Pow(2*h - 1, 2), 0.5); + double ag = (sg - (2*h - 1)*si)/4; + double al = 3.1415927/4 - ag; const double d = 1.44; - double dg = 4 * ag / (sg + si) * d; + double dg = 4*ag/(sg + si)*d; const double ugs = -0.029; - double ug = ugs * 3.1415927 / 4 / ag; + double ug = ugs*3.1415927/4/ag; const double uls = -0.00001; - double ul = uls * 3.1415927 / 4 / al; + double ul = uls*3.1415927/4/al; double sl = 3.1415927 - sg; - double dl = 4 * al / (sl) * d; + double dl = 4*al/(sl)*d; const double rol = 1; const double rog = 0.835; const double bt = .6; - double g = 980 * Math.Sin(bt * 3.1415927 / 180); - double dpg = (rol * al + rog * ag) * g / 3.1415927 * 4; + double g = 980*Math.Sin(bt*3.1415927/180); + double dpg = (rol*al + rog*ag)*g/3.1415927*4; const double vig = 1.0; - double reg = Math.Abs(ug) * rog * dg / vig; - double tg = -16 / reg * (.5 * rog * ug * Math.Abs(ug)); + double reg = Math.Abs(ug)*rog*dg/vig; + double tg = -16/reg*(.5*rog*ug*Math.Abs(ug)); const double vil = .01; - double rel = Math.Abs(ul) * rol * dl / vil; - double tl = -16 / rel * (.5 * rol * ul * Math.Abs(ul)); + double rel = Math.Abs(ul)*rol*dl/vil; + double tl = -16/rel*(.5*rol*ul*Math.Abs(ul)); const double it = 1; - double ti = -16 / reg * rog * (.5 * (ug - ul) * Math.Abs(ug - ul)) * it; - double dpf = (tl * sl + tg * sg) * d / (al + ag) / (d * d) * 4; + double ti = -16/reg*rog*(.5*(ug - ul)*Math.Abs(ug - ul))*it; + double dpf = (tl*sl + tg*sg)*d/(al + ag)/(d*d)*4; double dp = dpg + dpf; - double tic = (Math.Abs(ul) < Math.Abs(ug)) ? (rog / reg) : (rol / rel); - double y = (rol - rog) * g * Math.Pow(d / 2, 2) / (8 * vig * ugs); + double tic = (Math.Abs(ul) < Math.Abs(ug)) ? (rog/reg) : (rol/rel); + double y = (rol - rog)*g*Math.Pow(d/2, 2)/(8*vig*ugs); - return -tg * sg * al + tl * sl * ag - ti * si * (al + ag) + d * (rol - rog) * al * ag * g; + return -tg*sg*al + tl*sl*ag - ti*si*(al + ag) + d*(rol - rog)*al*ag*g; }; double x = Bisection.FindRoot(f1, 0, 0.1); diff --git a/src/UnitTests/RootFindingTests/BrentTest.cs b/src/UnitTests/RootFindingTests/BrentTest.cs index dadaba80..5fe59ab3 100644 --- a/src/UnitTests/RootFindingTests/BrentTest.cs +++ b/src/UnitTests/RootFindingTests/BrentTest.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2013 Math.NET -// +// +// Copyright (c) 2009-2014 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -41,8 +41,9 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void MultipleRoots() { // Roots at -2, 2 - Func f1 = x => x * x - 4; + Func f1 = x => x*x - 4; Assert.AreEqual(0, f1(Brent.FindRoot(f1, 1, 5, 1e-14, 100))); + Assert.AreEqual(0, f1(Brent.FindRootExpand(f1, 3, 5, 1e-14, 100))); Assert.AreEqual(-2, Brent.FindRoot(f1, -5, -1, 1e-14, 100)); Assert.AreEqual(2, Brent.FindRoot(f1, 1, 4, 1e-14, 100)); Assert.AreEqual(0, f1(Brent.FindRoot(x => -f1(x), 1, 5, 1e-14, 100))); @@ -50,7 +51,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests Assert.AreEqual(2, Brent.FindRoot(x => -f1(x), 1, 4, 1e-14, 100)); // Roots at 3, 4 - Func f2 = x => (x - 3) * (x - 4); + Func f2 = x => (x - 3)*(x - 4); Assert.AreEqual(0, f2(Brent.FindRoot(f2, -5, 3.5, 1e-14, 100))); Assert.AreEqual(3, Brent.FindRoot(f2, -5, 3.5, 1e-14, 100)); Assert.AreEqual(4, Brent.FindRoot(f2, 3.2, 5, 1e-14, 100)); @@ -61,7 +62,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void LocalMinima() { - Func f1 = x => x * x * x - 2 * x + 2; + Func f1 = x => x*x*x - 2*x + 2; Assert.AreEqual(0, f1(Brent.FindRoot(f1, -5, 5, 1e-14, 100)), 1e-14); Assert.AreEqual(0, f1(Brent.FindRoot(f1, -2, 4, 1e-14, 100)), 1e-14); } @@ -83,7 +84,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void NoRoot() { - Func f1 = x => x * x + 4; + Func f1 = x => x*x + 4; Assert.That(() => Brent.FindRoot(f1, -5, 5, 1e-14, 50), Throws.TypeOf()); } @@ -91,7 +92,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq1() { // Test case from http://www.polymath-software.com/library/nle/Oneeq1.htm - Func f1 = z => 8 * Math.Pow((4 - z) * z, 2) / (Math.Pow(6 - 3 * z, 2) * (2 - z)) - 0.186; + Func f1 = z => 8*Math.Pow((4 - z)*z, 2)/(Math.Pow(6 - 3*z, 2)*(2 - z)) - 0.186; double x = Brent.FindRoot(f1, 0.1, 0.9, 1e-14, 100); Assert.AreEqual(0.277759543089215, x, 1e-9); Assert.AreEqual(0, f1(x), 1e-16); @@ -107,15 +108,15 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double x2 = 1 - x1; const double G12 = 0.07858889; const double G21 = 0.30175355; - double P2 = Math.Pow(10, 6.87776 - 1171.53 / (224.366 + T)); - double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + T)); - const double t1 = x1 + x2 * G12; - const double t2 = x2 + x1 * G21; - double gamma2 = Math.Exp(-Math.Log(t2) - (x1 * (G12 * t2 - G21 * t1)) / (t1 * t2)); - double gamma1 = Math.Exp(-Math.Log(t1) + (x2 * (G12 * t2 - G21 * t1)) / (t1 * t2)); - double k1 = gamma1 * P1 / 760; - double k2 = gamma2 * P2 / 760; - return 1 - k1 * x1 - k2 * x2; + double P2 = Math.Pow(10, 6.87776 - 1171.53/(224.366 + T)); + double P1 = Math.Pow(10, 8.04494 - 1554.3/(222.65 + T)); + const double t1 = x1 + x2*G12; + const double t2 = x2 + x1*G21; + double gamma2 = Math.Exp(-Math.Log(t2) - (x1*(G12*t2 - G21*t1))/(t1*t2)); + double gamma1 = Math.Exp(-Math.Log(t1) + (x2*(G12*t2 - G21*t1))/(t1*t2)); + double k1 = gamma1*P1/760; + double k2 = gamma2*P2/760; + return 1 - k1*x1 - k2*x2; }; double x = Brent.FindRoot(f1, 0, 100, 1e-14, 100); @@ -133,13 +134,13 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double x2 = 1 - x1; const double A = 0.75807689; const double B = 1.1249035; - double P2 = Math.Pow(10, 6.9024 - 1268.115 / (216.9 + T)); - double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + T)); - double gamma2 = Math.Pow(10, B * Math.Pow(x1, 2) / Math.Pow(x1 + B * x2 / A, 2)); - double gamma1 = Math.Pow(10, A * Math.Pow(x2, 2) / Math.Pow(A * x1 / B + x2, 2)); - double k1 = gamma1 * P1 / 760; - double k2 = gamma2 * P2 / 760; - return 1 - k1 * x1 - k2 * x2; + double P2 = Math.Pow(10, 6.9024 - 1268.115/(216.9 + T)); + double P1 = Math.Pow(10, 8.04494 - 1554.3/(222.65 + T)); + double gamma2 = Math.Pow(10, B*Math.Pow(x1, 2)/Math.Pow(x1 + B*x2/A, 2)); + double gamma1 = Math.Pow(10, A*Math.Pow(x2, 2)/Math.Pow(A*x1/B + x2, 2)); + double k1 = gamma1*P1/760; + double k2 = gamma2*P2/760; + return 1 - k1*x1 - k2*x2; }; double x = Brent.FindRoot(f1, 0, 100, 1e-14); @@ -151,7 +152,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq3() { // Test case from http://www.polymath-software.com/library/nle/Oneeq3.htm - Func f1 = T => Math.Exp(21000 / T) / (T * T) - 1.11e11; + Func f1 = T => Math.Exp(21000/T)/(T*T) - 1.11e11; Assert.AreEqual(551.773822885233, Brent.FindRoot(f1, 550, 560, 1e-2), 1e-2); } @@ -163,7 +164,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests { const double A = 0.38969; const double B = 0.55954; - return A * B * (B * Math.Pow(1 - x, 2) - A * x * x) / Math.Pow(x * (A - B) + B, 2) + 0.14845; + return A*B*(B*Math.Pow(1 - x, 2) - A*x*x)/Math.Pow(x*(A - B) + B, 2) + 0.14845; }; double r = Brent.FindRoot(f1, 0, 1, 1.0e-14); @@ -182,7 +183,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double BETA = 2.298E-3; const double GAMMA = 0.283E-6; const double DH = -57798.0; - return DH + TR * (ALPHA + TR * (BETA / 2 + TR * GAMMA / 3)) - 298.0 * (ALPHA + 298.0 * (BETA / 2 + 298.0 * GAMMA / 3)); + return DH + TR*(ALPHA + TR*(BETA/2 + TR*GAMMA/3)) - 298.0*(ALPHA + 298.0*(BETA/2 + 298.0*GAMMA/3)); }; double x = Brent.FindRoot(f1, 3000, 5000, 1.0e-3); @@ -204,11 +205,11 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double A = 0.01855; const double B = -0.01587; const double P = 100.0; - const double Beta = R * T * B0 - A0 - R * C / (T * T); - const double Gama = -R * T * B0 * B + A0 * A - R * C * B0 / (T * T); - const double Delta = R * B0 * B * C / (T * T); + const double Beta = R*T*B0 - A0 - R*C/(T*T); + const double Gama = -R*T*B0*B + A0*A - R*C*B0/(T*T); + const double Delta = R*B0*B*C/(T*T); - return R * T / V + Beta / (V * V) + Gama / (V * V * V) + Delta / (V * V * V * V) - P; + return R*T/V + Beta/(V*V) + Gama/(V*V*V) + Delta/(V*V*V*V) - P; }; Assert.AreEqual(0.174749531708621, Brent.FindRoot(f1, 0.1, 1), 1e-8); @@ -228,11 +229,11 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double A = 0.01855; const double B = -0.01587; const double P = 100.0; - const double Beta = R * T * B0 - A0 - R * C / (T * T); - const double Gama = -R * T * B0 * B + A0 * A - R * C * B0 / (T * T); - const double Delta = R * B0 * B * C / (T * T); + const double Beta = R*T*B0 - A0 - R*C/(T*T); + const double Gama = -R*T*B0*B + A0*A - R*C*B0/(T*T); + const double Delta = R*B0*B*C/(T*T); - return R * T * V * V * V + Beta * V * V + Gama * V + Delta - P * V * V * V * V; + return R*T*V*V*V + Beta*V*V + Gama*V + Delta - P*V*V*V*V; }; double x = Brent.FindRoot(f1, 0.1, 1, 1e-14); @@ -244,7 +245,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq7() { // Test case from http://www.polymath-software.com/library/nle/Oneeq7.htm - Func f1 = x => x / (1 - x) - 5 * Math.Log(0.4 * (1 - x) / (0.4 - 0.5 * x)) + 4.45977; + Func f1 = x => x/(1 - x) - 5*Math.Log(0.4*(1 - x)/(0.4 - 0.5*x)) + 4.45977; Assert.AreEqual(0.757396293891, Brent.FindRoot(f1, 0, 0.79, 1e-2), 1e-2); } @@ -258,7 +259,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double b = 40; const double c = 200; - return a * v * v + b * Math.Pow(v, 7 / 4) - c; + return a*v*v + b*Math.Pow(v, 7/4) - c; }; Assert.AreEqual(0.842524411168525, Brent.FindRoot(f1, 0.01, 1, 1e-2), 1e-2); @@ -275,7 +276,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double gama = 1.41; const double Pd = 100; - return Math.Pow((gama + 1) / 2, (gama + 1) / (gama - 1)) * (2 / (gama - 1)) * (Math.Pow(P / Pd, 2 / gama) - Math.Pow(P / Pd, (gama + 1) / gama)) - Math.Pow(At / A, 2); + return Math.Pow((gama + 1)/2, (gama + 1)/(gama - 1))*(2/(gama - 1))*(Math.Pow(P/Pd, 2/gama) - Math.Pow(P/Pd, (gama + 1)/gama)) - Math.Pow(At/A, 2); }; double x = Brent.FindRoot(f1, 1, 50, 1e-14); @@ -290,7 +291,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq10() { // Test case from http://www.polymath-software.com/library/nle/Oneeq10.htm - Func f1 = x => (1.0 / 63.0) * Math.Log(x) + (64.0 / 63.0) * Math.Log(1 / (1 - x)) + Math.Log(0.95 - x) - Math.Log(0.9); + Func f1 = x => (1.0/63.0)*Math.Log(x) + (64.0/63.0)*Math.Log(1/(1 - x)) + Math.Log(0.95 - x) - Math.Log(0.9); double r = Brent.FindRoot(f1, .01, 0.35, 1e-14); Assert.AreEqual(0.036210083704, r, 1e-5); @@ -312,7 +313,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double kp = 604500; const double P = 0.00243; - return a - Math.Pow((c + 2 * x) * (a + b + c - 2 * x), 2) / (kp * P * P * Math.Pow(b - 3 * x, 3)) - x; + return a - Math.Pow((c + 2*x)*(a + b + c - 2*x), 2)/(kp*P*P*Math.Pow(b - 3*x, 3)) - x; }; double r = Brent.FindRoot(f1, 0, 0.25, 1e-14); @@ -335,7 +336,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double kp = 604500; const double P = 0.00243; - return (b - Math.Pow(Math.Pow((c + 2 * x) * (a + b + c - 2 * x), 2) / (kp * P * P * (a - x)), 1.0 / 3.0)) / 3.0 - x; + return (b - Math.Pow(Math.Pow((c + 2*x)*(a + b + c - 2*x), 2)/(kp*P*P*(a - x)), 1.0/3.0))/3.0 - x; }; double r = Brent.FindRoot(f1, 0.01, 0.45, 1e-14); @@ -359,7 +360,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double T = 430.85; const double R = 82.05; - return (R * T / P) * (1 + b / v + c / (v * v)) - v; + return (R*T/P)*(1 + b/v + c/(v*v)) - v; }; double x = Brent.FindRoot(f1, 100, 300, 1e-14); @@ -378,10 +379,10 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double P = 100; const double T = 210.0; const double R = 0.08205; - double A = 0.42748 * (R * R * Math.Pow(Tc, 2.5)) / Pc; - const double B = 0.08664 * R * Tc / Pc; + double A = 0.42748*(R*R*Math.Pow(Tc, 2.5))/Pc; + const double B = 0.08664*R*Tc/Pc; - return R * Math.Pow(T, 1.5) * V * (V + B) / (P * Math.Sqrt(T) * V * (V + B) + A) + B - V; + return R*Math.Pow(T, 1.5)*V*(V + B)/(P*Math.Sqrt(T)*V*(V + B) + A) + B - V; }; double x = Brent.FindRoot(f1, .01, .3, 1e-14); @@ -403,10 +404,10 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double T1 = 390.0; const double t1 = 100.0; const double cp = 0.49; - double DT1 = T1 - t1 - Q / (M * Cp); - double DT2 = T1 - t1 - Q / (m * cp); + double DT1 = T1 - t1 - Q/(M*Cp); + double DT2 = T1 - t1 - Q/(m*cp); - return U * A * (DT2 - DT1) / Math.Log(DT2 / DT1) - Q; + return U*A*(DT2 - DT1)/Math.Log(DT2/DT1) - Q; }; double x = Brent.FindRoot(f1, 1000000, 7000000, 1e-9); @@ -418,7 +419,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq15a() { // Test case from http://www.polymath-software.com/library/nle/Oneeq15a.htm - Func f1 = h => h * (1 + h * (3 - h)) - 2.4 - h; + Func f1 = h => h*(1 + h*(3 - h)) - 2.4 - h; double x = Brent.FindRoot(f1, -2, 0, 1e-14); Assert.AreEqual(-0.795219754733581, x, 1e-5); @@ -435,7 +436,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq15b() { // Test case from http://www.polymath-software.com/library/nle/Oneeq15b.htm - Func f1 = h => 2.4 / (h * (3 - h)) - h; + Func f1 = h => 2.4/(h*(3 - h)) - h; double x = Brent.FindRoot(f1, -2, -0.5, 1e-14); Assert.AreEqual(-0.795219745444464, x, 1e-5); @@ -459,18 +460,18 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double z = 1 - y - 0.02; const double CH4 = 0; const double C2H6 = 0; - const double COtwo = y + 2 * z; - const double H2O = 2 * y + 3 * z; - const double N2 = 0.02 + 3.76 * (2 * y + 7 * z / 2) * x; - const double O2 = (2 * y + 7 * z / 2) * (x - 1); - const double alp = 3.381 * CH4 + 2.247 * C2H6 + 6.214 * COtwo + 7.256 * H2O + 6.524 * N2 + 6.148 * O2; - const double bet = 18.044 * CH4 + 38.201 * C2H6 + 10.396 * COtwo + 2.298 * H2O + 1.25 * N2 + 3.102 * O2; - const double gam = -4.3 * CH4 - 11.049 * C2H6 - 3.545 * COtwo + 0.283 * H2O - 0.001 * N2 - 0.923 * O2; - const double H0 = alp * 298 + bet * 0.001 * 298 * 298 / 2 + gam * 1e-6 * 298 * 298 * 298 / 3; - double Hf = alp * T + bet * 0.001 * T * T / 2 + gam * 1e-6 * T * T * T / 3; + const double COtwo = y + 2*z; + const double H2O = 2*y + 3*z; + const double N2 = 0.02 + 3.76*(2*y + 7*z/2)*x; + const double O2 = (2*y + 7*z/2)*(x - 1); + const double alp = 3.381*CH4 + 2.247*C2H6 + 6.214*COtwo + 7.256*H2O + 6.524*N2 + 6.148*O2; + const double bet = 18.044*CH4 + 38.201*C2H6 + 10.396*COtwo + 2.298*H2O + 1.25*N2 + 3.102*O2; + const double gam = -4.3*CH4 - 11.049*C2H6 - 3.545*COtwo + 0.283*H2O - 0.001*N2 - 0.923*O2; + const double H0 = alp*298 + bet*0.001*298*298/2 + gam*1e-6*298*298*298/3; + double Hf = alp*T + bet*0.001*T*T/2 + gam*1e-6*T*T*T/3; const double xx = 1; - return 212798 * y * xx + 372820 * z * xx + H0 - Hf; + return 212798*y*xx + 372820*z*xx + H0 - Hf; }; double r = Brent.FindRoot(f1, 1000, 3000); @@ -482,7 +483,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq17() { // Test case from http://www.polymath-software.com/library/nle/Oneeq17.htm - Func f1 = rp => rp - 0.327 * Math.Pow(0.06 - 161 * rp, 0.804) * Math.Exp(-5230 / (1.987 * (373 + 1.84e6 * rp))); + Func f1 = rp => rp - 0.327*Math.Pow(0.06 - 161*rp, 0.804)*Math.Exp(-5230/(1.987*(373 + 1.84e6*rp))); double x = Brent.FindRoot(f1, 0, 0.00035, 1e-13); Assert.AreEqual(0.000340568862275, x, 1e-5); @@ -495,12 +496,12 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Test case from http://www.polymath-software.com/library/nle/Oneeq18a.htm Func f1 = T => { - double Tp = (269.267 * T - 1752) / 266.667; - double k = 0.0006 * Math.Exp(20.7 - 15000 / Tp); - double Pp = -(0.1 / (321 * (320.0 / 321.0 - (1 + k)))); - double P = (320 * Pp + 0.1) / 321; + double Tp = (269.267*T - 1752)/266.667; + double k = 0.0006*Math.Exp(20.7 - 15000/Tp); + double Pp = -(0.1/(321*(320.0/321.0 - (1 + k)))); + double P = (320*Pp + 0.1)/321; - return 1.296 * (T - Tp) + 10369 * k * Pp; + return 1.296*(T - Tp) + 10369*k*Pp; }; double x = Brent.FindRoot(f1, 500, 725); @@ -520,12 +521,12 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Test case from http://www.polymath-software.com/library/nle/Oneeq18b.htm Func f1 = T => { - double Tp = (269 * T - 1752) / 267; - double k = 0.0006 * Math.Exp(20.7 - 15000 / Tp); - double Pp = -(0.1 / (321 * (320.0 / 321.0 - (1 + k)))); - double P = (320 * Pp + 0.1) / 321; + double Tp = (269*T - 1752)/267; + double k = 0.0006*Math.Exp(20.7 - 15000/Tp); + double Pp = -(0.1/(321*(320.0/321.0 - (1 + k)))); + double P = (320*Pp + 0.1)/321; - return 1.3 * (T - Tp) + 1.04e4 * k * Pp; + return 1.3*(T - Tp) + 1.04e4*k*Pp; }; double x = Brent.FindRoot(f1, 500, 1500); @@ -541,27 +542,27 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests { const double P = 200; const double Pc = 33.5; - const double T = 631 * 2; + const double T = 631*2; const double Tc = 126.2; - const double Pr = P / Pc; - const double Tr = T / Tc; - double Asqr = 0.4278 * Pr / Math.Pow(Tr, 2.5); - const double B = 0.0867 * Pr / Tr; - double Q = B * B + B - Asqr; - double r = Asqr * B; - double p = (-3 * Q - 1) / 3; - double q = (-27 * r - 9 * Q - 2) / 27; - double R = Math.Pow(p / 3, 3) + Math.Pow(q / 2, 2); - double V = (R > 0) ? Math.Pow(-q / 2 + Math.Sqrt(R), 1 / 3) : 0; - double WW = (R > 0) ? (-q / 2 - Math.Sqrt(R)) : (0); - double psi1 = (R < 0) ? (Math.Acos(Math.Sqrt((q * q / 4) / (-p * p * p / 27)))) : (0); - double W = (R > 0) ? (Math.Sign(WW) * Math.Pow(Math.Abs(WW), 1 / 3)) : (0); - double z1 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 0 / 3) + 1 / 3) : (0); - double z2 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 1 / 3) + 1 / 3) : (0); - double z3 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 2 / 3) + 1 / 3) : (0); - double z0 = (R > 0) ? (V + W + 1 / 3) : (0); - - return z * z * z - z * z - Q * z - r; + const double Pr = P/Pc; + const double Tr = T/Tc; + double Asqr = 0.4278*Pr/Math.Pow(Tr, 2.5); + const double B = 0.0867*Pr/Tr; + double Q = B*B + B - Asqr; + double r = Asqr*B; + double p = (-3*Q - 1)/3; + double q = (-27*r - 9*Q - 2)/27; + double R = Math.Pow(p/3, 3) + Math.Pow(q/2, 2); + double V = (R > 0) ? Math.Pow(-q/2 + Math.Sqrt(R), 1/3) : 0; + double WW = (R > 0) ? (-q/2 - Math.Sqrt(R)) : (0); + double psi1 = (R < 0) ? (Math.Acos(Math.Sqrt((q*q/4)/(-p*p*p/27)))) : (0); + double W = (R > 0) ? (Math.Sign(WW)*Math.Pow(Math.Abs(WW), 1/3)) : (0); + double z1 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*0/3) + 1/3) : (0); + double z2 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*1/3) + 1/3) : (0); + double z3 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*2/3) + 1/3) : (0); + double z0 = (R > 0) ? (V + W + 1/3) : (0); + + return z*z*z - z*z - Q*z - r; }; double x = Brent.FindRoot(f1, -0.5, -0.02, 1e-14); @@ -585,25 +586,25 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double Pc = 33.5; const double T = 631; const double Tc = 126.2; - const double Pr = P / Pc; - const double Tr = T / Tc; - double Asqr = 0.4278 * Pr / Math.Pow(Tr, 2.5); - const double B = 0.0867 * Pr / Tr; - double Q = B * B + B - Asqr; - double r = Asqr * B; - double p = (-3 * Q - 1) / 3; - double q = (-27 * r - 9 * Q - 2) / 27; - double R = Math.Pow(p / 3, 3) + Math.Pow(q / 2, 2); - double V = (R > 0) ? Math.Pow(-q / 2 + Math.Sqrt(R), 1 / 3) : 0; - double WW = (R > 0) ? (-q / 2 - Math.Sqrt(R)) : (0); - double psi1 = (R < 0) ? (Math.Acos(Math.Sqrt((q * q / 4) / (-p * p * p / 27)))) : (0); - double W = (R > 0) ? (Math.Sign(WW) * Math.Pow(Math.Abs(WW), 1 / 3)) : (0); - double z1 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 0 / 3) + 1 / 3) : (0); - double z2 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 1 / 3) + 1 / 3) : (0); - double z3 = (R < 0) ? (2 * Math.Sqrt(-p / 3) * Math.Cos((psi1 / 3) + 2 * 3.1416 * 2 / 3) + 1 / 3) : (0); - double z0 = (R > 0) ? (V + W + 1 / 3) : (0); - - return z * z * z - z * z - Q * z - r; + const double Pr = P/Pc; + const double Tr = T/Tc; + double Asqr = 0.4278*Pr/Math.Pow(Tr, 2.5); + const double B = 0.0867*Pr/Tr; + double Q = B*B + B - Asqr; + double r = Asqr*B; + double p = (-3*Q - 1)/3; + double q = (-27*r - 9*Q - 2)/27; + double R = Math.Pow(p/3, 3) + Math.Pow(q/2, 2); + double V = (R > 0) ? Math.Pow(-q/2 + Math.Sqrt(R), 1/3) : 0; + double WW = (R > 0) ? (-q/2 - Math.Sqrt(R)) : (0); + double psi1 = (R < 0) ? (Math.Acos(Math.Sqrt((q*q/4)/(-p*p*p/27)))) : (0); + double W = (R > 0) ? (Math.Sign(WW)*Math.Pow(Math.Abs(WW), 1/3)) : (0); + double z1 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*0/3) + 1/3) : (0); + double z2 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*1/3) + 1/3) : (0); + double z3 = (R < 0) ? (2*Math.Sqrt(-p/3)*Math.Cos((psi1/3) + 2*3.1416*2/3) + 1/3) : (0); + double z0 = (R > 0) ? (V + W + 1/3) : (0); + + return z*z*z - z*z - Q*z - r; }; double x = Brent.FindRoot(f1, -0.5, 1.2, 1e-14); @@ -619,14 +620,14 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests { const double T = 313; const double P0 = 10; - const double FA0 = 20 * P0 / (0.082 * 450); - double k1 = 1.277 * 1.0e9 * Math.Exp(-90000 / (8.31 * T)); - double k2 = 1.29 * 1.0e11 * Math.Exp(-135000 / (8.31 * T)); - double den = 1 + 7 * P0 * (1 - xa) / (1 + xa); + const double FA0 = 20*P0/(0.082*450); + double k1 = 1.277*1.0e9*Math.Exp(-90000/(8.31*T)); + double k2 = 1.29*1.0e11*Math.Exp(-135000/(8.31*T)); + double den = 1 + 7*P0*(1 - xa)/(1 + xa); double xa1 = 1 + xa; - double ra = (-k1 * P0 * (1 - xa) / xa1 + k2 * P0 * P0 * xa * xa / (xa1 * xa1)) / den; + double ra = (-k1*P0*(1 - xa)/xa1 + k2*P0*P0*xa*xa/(xa1*xa1))/den; - return -ra / FA0; + return -ra/FA0; }; double x = Brent.FindRoot(f1, 0.75, 1.02); @@ -642,13 +643,13 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests { const double T = 313; const double P0 = 10; - const double FA0 = 20 * P0 / (0.082 * 450); - double k1 = 1.277 * 1.0e9 * Math.Exp(-90000 / (8.31 * T)); - double k2 = 1.29 * 1.0e11 * Math.Exp(-135000 / (8.31 * T)); + const double FA0 = 20*P0/(0.082*450); + double k1 = 1.277*1.0e9*Math.Exp(-90000/(8.31*T)); + double k2 = 1.29*1.0e11*Math.Exp(-135000/(8.31*T)); double xa1 = 1 + xa; - double ra = (-k1 * P0 * (1 - xa) / xa1 + k2 * P0 * P0 * xa * xa / (xa1 * xa1)); + double ra = (-k1*P0*(1 - xa)/xa1 + k2*P0*P0*xa*xa/(xa1*xa1)); - return -ra / FA0; + return -ra/FA0; }; double x = Brent.FindRoot(f1, 0.75, 1.02, 1e-10); @@ -662,37 +663,37 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Test case from http://www.polymath-software.com/library/nle/Oneeq21a.htm Func f1 = h => { - double sg = Math.Acos(2 * h - 1); - double si = Math.Pow(1 - Math.Pow(2 * h - 1, 2), 0.5); - double ag = (sg - (2 * h - 1) * si) / 4; - double al = 3.1415927 / 4 - ag; + double sg = Math.Acos(2*h - 1); + double si = Math.Pow(1 - Math.Pow(2*h - 1, 2), 0.5); + double ag = (sg - (2*h - 1)*si)/4; + double al = 3.1415927/4 - ag; const double d = 1.44; - double dg = 4 * ag / (sg + si) * d; + double dg = 4*ag/(sg + si)*d; const double ugs = -0.0295; - double ug = ugs * 3.1415927 / 4 / ag; + double ug = ugs*3.1415927/4/ag; const double uls = -0.00001; - double ul = uls * 3.1415927 / 4 / al; + double ul = uls*3.1415927/4/al; double sl = 3.1415927 - sg; - double dl = 4 * al / (sl) * d; + double dl = 4*al/(sl)*d; const double rol = 1; const double rog = 0.835; const double bt = .6; - double g = 980 * Math.Sin(bt * 3.1415927 / 180); - double dpg = (rol * al + rog * ag) * g / 3.1415927 * 4; + double g = 980*Math.Sin(bt*3.1415927/180); + double dpg = (rol*al + rog*ag)*g/3.1415927*4; const double vig = 1.0; - double reg = Math.Abs(ug) * rog * dg / vig; - double tg = -16 / reg * (.5 * rog * ug * Math.Abs(ug)); + double reg = Math.Abs(ug)*rog*dg/vig; + double tg = -16/reg*(.5*rog*ug*Math.Abs(ug)); const double vil = .01; - double rel = Math.Abs(ul) * rol * dl / vil; - double tl = -16 / rel * (.5 * rol * ul * Math.Abs(ul)); + double rel = Math.Abs(ul)*rol*dl/vil; + double tl = -16/rel*(.5*rol*ul*Math.Abs(ul)); const double it = 1; - double ti = -16 / reg * rog * (.5 * (ug - ul) * Math.Abs(ug - ul)) * it; - double dpf = (tl * sl + tg * sg) * d / (al + ag) / (d * d) * 4; + double ti = -16/reg*rog*(.5*(ug - ul)*Math.Abs(ug - ul))*it; + double dpf = (tl*sl + tg*sg)*d/(al + ag)/(d*d)*4; double dp = dpg + dpf; - double tic = (Math.Abs(ul) < Math.Abs(ug)) ? (rog / reg) : (rol / rel); - double y = (rol - rog) * g * Math.Pow(d / 2, 2) / (8 * vig * ugs); + double tic = (Math.Abs(ul) < Math.Abs(ug)) ? (rog/reg) : (rol/rel); + double y = (rol - rog)*g*Math.Pow(d/2, 2)/(8*vig*ugs); - return -tg * sg * al + tl * sl * ag - ti * si * (al + ag) + d * (rol - rog) * al * ag * g; + return -tg*sg*al + tl*sl*ag - ti*si*(al + ag) + d*(rol - rog)*al*ag*g; }; double x = Brent.FindRoot(f1, 0, 0.3, 1e-14); @@ -706,37 +707,37 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Test case from http://www.polymath-software.com/library/nle/Oneeq21b.htm Func f1 = h => { - double sg = Math.Acos(2 * h - 1); - double si = Math.Pow(1 - Math.Pow(2 * h - 1, 2), 0.5); - double ag = (sg - (2 * h - 1) * si) / 4; - double al = 3.1415927 / 4 - ag; + double sg = Math.Acos(2*h - 1); + double si = Math.Pow(1 - Math.Pow(2*h - 1, 2), 0.5); + double ag = (sg - (2*h - 1)*si)/4; + double al = 3.1415927/4 - ag; const double d = 1.44; - double dg = 4 * ag / (sg + si) * d; + double dg = 4*ag/(sg + si)*d; const double ugs = -0.029; - double ug = ugs * 3.1415927 / 4 / ag; + double ug = ugs*3.1415927/4/ag; const double uls = -0.00001; - double ul = uls * 3.1415927 / 4 / al; + double ul = uls*3.1415927/4/al; double sl = 3.1415927 - sg; - double dl = 4 * al / (sl) * d; + double dl = 4*al/(sl)*d; const double rol = 1; const double rog = 0.835; const double bt = .6; - double g = 980 * Math.Sin(bt * 3.1415927 / 180); - double dpg = (rol * al + rog * ag) * g / 3.1415927 * 4; + double g = 980*Math.Sin(bt*3.1415927/180); + double dpg = (rol*al + rog*ag)*g/3.1415927*4; const double vig = 1.0; - double reg = Math.Abs(ug) * rog * dg / vig; - double tg = -16 / reg * (.5 * rog * ug * Math.Abs(ug)); + double reg = Math.Abs(ug)*rog*dg/vig; + double tg = -16/reg*(.5*rog*ug*Math.Abs(ug)); const double vil = .01; - double rel = Math.Abs(ul) * rol * dl / vil; - double tl = -16 / rel * (.5 * rol * ul * Math.Abs(ul)); + double rel = Math.Abs(ul)*rol*dl/vil; + double tl = -16/rel*(.5*rol*ul*Math.Abs(ul)); const double it = 1; - double ti = -16 / reg * rog * (.5 * (ug - ul) * Math.Abs(ug - ul)) * it; - double dpf = (tl * sl + tg * sg) * d / (al + ag) / (d * d) * 4; + double ti = -16/reg*rog*(.5*(ug - ul)*Math.Abs(ug - ul))*it; + double dpf = (tl*sl + tg*sg)*d/(al + ag)/(d*d)*4; double dp = dpg + dpf; - double tic = (Math.Abs(ul) < Math.Abs(ug)) ? (rog / reg) : (rol / rel); - double y = (rol - rog) * g * Math.Pow(d / 2, 2) / (8 * vig * ugs); + double tic = (Math.Abs(ul) < Math.Abs(ug)) ? (rog/reg) : (rol/rel); + double y = (rol - rog)*g*Math.Pow(d/2, 2)/(8*vig*ugs); - return -tg * sg * al + tl * sl * ag - ti * si * (al + ag) + d * (rol - rog) * al * ag * g; + return -tg*sg*al + tl*sl*ag - ti*si*(al + ag) + d*(rol - rog)*al*ag*g; }; double x = Brent.FindRoot(f1, 0, 0.1, 1e-14); diff --git a/src/UnitTests/RootFindingTests/FindRootsTest.cs b/src/UnitTests/RootFindingTests/FindRootsTest.cs index 8c54d708..3e66d011 100644 --- a/src/UnitTests/RootFindingTests/FindRootsTest.cs +++ b/src/UnitTests/RootFindingTests/FindRootsTest.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2013 Math.NET -// +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -37,6 +37,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif [TestFixture, Category("RootFinding")] @@ -46,7 +47,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void MultipleRoots() { // Roots at -2, 2 - Func f1 = x => x * x - 4; + Func f1 = x => x*x - 4; Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, 0, 5, 1e-14))); Assert.AreEqual(-2, FindRoots.OfFunction(f1, -5, -1, 1e-14)); Assert.AreEqual(2, FindRoots.OfFunction(f1, 1, 4, 1e-14)); @@ -55,7 +56,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests Assert.AreEqual(2, FindRoots.OfFunction(x => -f1(x), 1, 4, 1e-14)); // Roots at 3, 4 - Func f2 = x => (x - 3) * (x - 4); + Func f2 = x => (x - 3)*(x - 4); Assert.AreEqual(0, f2(FindRoots.OfFunction(f2, 3.5, 5, 1e-14)), 1e-14); Assert.AreEqual(3, FindRoots.OfFunction(f2, -5, 3.5, 1e-14)); Assert.AreEqual(4, FindRoots.OfFunction(f2, 3.2, 5, 1e-14)); @@ -66,7 +67,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void LocalMinima() { - Func f1 = x => x * x * x - 2 * x + 2; + Func f1 = x => x*x*x - 2*x + 2; Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, -5, 5, 1e-14)), 1e-14); Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, -2, 4, 1e-14)), 1e-14); } @@ -74,7 +75,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void NoRoot() { - Func f1 = x => x * x + 4; + Func f1 = x => x*x + 4; Assert.That(() => FindRoots.OfFunction(f1, -5, 5, 1e-14), Throws.TypeOf()); } @@ -82,8 +83,8 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq3() { // Test case from http://www.polymath-software.com/library/nle/Oneeq3.htm - Func f1 = T => Math.Exp(21000 / T) / (T * T) - 1.11e11; - double x = FindRoots.OfFunction(f1, 550, 560, 1e-2); + Func f1 = T => Math.Exp(21000/T)/(T*T) - 1.11e11; + double x = FindRoots.OfFunction(f1, 550, 560, 1e-12); Assert.AreEqual(551.773822885233, x, 1e-5); Assert.AreEqual(0, f1(x), 1e-2); } @@ -98,12 +99,12 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double BETA = 2.298E-3; const double GAMMA = 0.283E-6; const double DH = -57798.0; - return DH + TR * (ALPHA + TR * (BETA / 2 + TR * GAMMA / 3)) - 298.0 * (ALPHA + 298.0 * (BETA / 2 + 298.0 * GAMMA / 3)); + return DH + TR*(ALPHA + TR*(BETA/2 + TR*GAMMA/3)) - 298.0*(ALPHA + 298.0*(BETA/2 + 298.0*GAMMA/3)); }; double x = FindRoots.OfFunction(f1, 3000, 5000, 1e-10); Assert.AreEqual(4305.30991366774, x, 1e-5); - Assert.AreEqual(0, f1(x), 1e-10); + Assert.AreEqual(0, f1(x), 1e-8); } [Test] @@ -120,26 +121,26 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double A = 0.01855; const double B = -0.01587; const double P = 100.0; - const double Beta = R * T * B0 - A0 - R * C / (T * T); - const double Gama = -R * T * B0 * B + A0 * A - R * C * B0 / (T * T); - const double Delta = R * B0 * B * C / (T * T); + const double Beta = R*T*B0 - A0 - R*C/(T*T); + const double Gama = -R*T*B0*B + A0*A - R*C*B0/(T*T); + const double Delta = R*B0*B*C/(T*T); - return R * T / V + Beta / (V * V) + Gama / (V * V * V) + Delta / (V * V * V * V) - P; + return R*T/V + Beta/(V*V) + Gama/(V*V*V) + Delta/(V*V*V*V) - P; }; double x = FindRoots.OfFunction(f1, 0.1, 1); Assert.AreEqual(0.174749531708621, x, 1e-5); - Assert.AreEqual(0, f1(x), 1e-12); + Assert.AreEqual(0, f1(x), 1e-5); } - + [Test] public void Oneeq7() { // Test case from http://www.polymath-software.com/library/nle/Oneeq7.htm - Func f1 = x => x / (1 - x) - 5 * Math.Log(0.4 * (1 - x) / (0.4 - 0.5 * x)) + 4.45977; - double r = FindRoots.OfFunction(f1, 0, 0.79, 1e-2); - Assert.AreEqual(0.757396293891, r, 1e-5); - Assert.AreEqual(0, f1(r), 1e-13); + Func f1 = x => x/(1 - x) - 5*Math.Log(0.4*(1 - x)/(0.4 - 0.5*x)) + 4.45977; + double r = FindRoots.OfFunction(f1, 0, 0.79, 1e-10); + Assert.AreEqual(0.757396293891, r, 1e-6); + Assert.AreEqual(0, f1(r), 1e-6); } [Test] @@ -152,15 +153,15 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests const double b = 40; const double c = 200; - return a * v * v + b * Math.Pow(v, 7 / 4) - c; + return a*v*v + b*Math.Pow(v, 7/4) - c; }; double x = FindRoots.OfFunction(f1, 0.01, 1); Assert.AreEqual(0.842524411168525, x, 1e-2); - Assert.AreEqual(0, f1(x), 1e-12); + Assert.AreEqual(0, f1(x), 1e-5); } - private void AssertComplexEqual(Complex expected, Complex actual, double delta) + void AssertComplexEqual(Complex expected, Complex actual, double delta) { Assert.AreEqual(expected.Real, actual.Real, delta); Assert.AreEqual(expected.Imaginary, actual.Imaginary, delta); @@ -180,7 +181,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Expected Roots: -u, -v Complex r1 = new Complex(-u, 0d), r2 = new Complex(-v, 0d); Assert.IsTrue(x1.AlmostEqualRelative(r1, 1e-14) && x2.AlmostEqualRelative(r2, 1e-14) - || x1.AlmostEqualRelative(r2, 1e-14) && x2.AlmostEqualRelative(r1, 1e-14)); + || x1.AlmostEqualRelative(r2, 1e-14) && x2.AlmostEqualRelative(r1, 1e-14)); // Verify they really are roots AssertComplexEqual(Complex.Zero, c + b*x1 + a*x1*x1, 1e-14); @@ -200,7 +201,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Expected roots? Complex r1 = new Complex(x1R, x1I), r2 = new Complex(x2R, x2I); Assert.IsTrue(x1.AlmostEqualRelative(r1, 1e-14) && x2.AlmostEqualRelative(r2, 1e-14) - || x1.AlmostEqualRelative(r2, 1e-14) && x2.AlmostEqualRelative(r1, 1e-14)); + || x1.AlmostEqualRelative(r2, 1e-14) && x2.AlmostEqualRelative(r1, 1e-14)); // Verify they really are roots AssertComplexEqual(Complex.Zero, c + b*x1 + a*x1*x1, 1e-14);