From 1a5c1fba740dd72835f41912c4bf33147376107e Mon Sep 17 00:00:00 2001 From: diluculo Date: Tue, 6 Nov 2018 19:53:42 +0900 Subject: [PATCH 1/3] Special Functions: cleanup Bessel-type functions - Scale enum is added and replaced expScaled boolean flag. - Scaled Bessel functions are added. - HankelH1 and HankelH2 - allow only complex argument because they return usually complex values. - SphericalBesselJ and SphericalBesselY: - scaling options are removed. --- .../SpecialFunctionsTests/BesselTests.cs | 4 +- src/Numerics/Constants.cs | 3 + src/Numerics/SpecialFunctions/Airy.cs | 214 ++++++++++------ src/Numerics/SpecialFunctions/Bessel.cs | 231 ++++++++++++------ src/Numerics/SpecialFunctions/Hankel.cs | 68 +++--- src/Numerics/SpecialFunctions/Options.cs | 23 ++ .../SpecialFunctions/SphericalBessel.cs | 126 +++++++--- 7 files changed, 454 insertions(+), 215 deletions(-) create mode 100644 src/Numerics/SpecialFunctions/Options.cs diff --git a/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs b/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs index 60ed1416..f8dd3451 100644 --- a/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs +++ b/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs @@ -251,7 +251,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests public void BesselIRatioExact(int n, double zr, double zi, double cyr, double cyi, int decimalPlaces) { var z = new Complex(zr, zi); - var actual = SpecialFunctions.BesselI(n + 1, z, true) / SpecialFunctions.BesselI(n, z, true); + var actual = SpecialFunctions.BesselI(n + 1, z, SpecialFunctions.Scale.Exponential) / SpecialFunctions.BesselI(n, z, SpecialFunctions.Scale.Exponential); AssertHelpers.AlmostEqualRelative(new Complex(cyr, cyi), actual, decimalPlaces); } @@ -267,7 +267,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests public void BesselKRatioExact(int n, double zr, double zi, double cyr, double cyi, int decimalPlaces) { var z = new Complex(zr, zi); - var actual = SpecialFunctions.BesselK(n + 1, z, true) / SpecialFunctions.BesselK(n, z, true); + var actual = SpecialFunctions.BesselK(n + 1, z, SpecialFunctions.Scale.Exponential) / SpecialFunctions.BesselK(n, z, SpecialFunctions.Scale.Exponential); AssertHelpers.AlmostEqualRelative(new Complex(cyr, cyi), actual, decimalPlaces); } diff --git a/src/Numerics/Constants.cs b/src/Numerics/Constants.cs index 5b6fb130..59931864 100644 --- a/src/Numerics/Constants.cs +++ b/src/Numerics/Constants.cs @@ -96,6 +96,9 @@ namespace MathNet.Numerics /// The number sqrt(2pi) public const double Sqrt2Pi = 2.5066282746310005024157652848110452530069867406099d; + /// The number sqrt(pi/2) + public const double SqrtPiOver2 = 1.2533141373155002512078826424055226265034933703050d; + /// The number sqrt(2*pi*e) public const double Sqrt2PiE = 4.1327313541224929384693918842998526494455219169913d; diff --git a/src/Numerics/SpecialFunctions/Airy.cs b/src/Numerics/SpecialFunctions/Airy.cs index 1e69093c..c2f7936a 100644 --- a/src/Numerics/SpecialFunctions/Airy.cs +++ b/src/Numerics/SpecialFunctions/Airy.cs @@ -8,121 +8,195 @@ namespace MathNet.Numerics public static partial class SpecialFunctions { /// - /// Airy function Ai(z). - ///

- /// If expScaled is true, returns Exp(zta) * Ai(z), where zta = (2/3) * z * Sqrt(z). + /// Returns the Airy function Ai. + /// AiryAi(z) is a solution to the Airy equation, y'' - y * z = 0. + /// AiryAi(z, Scale.Exponential) returns Exp(zta) * AiryAi(z), where zta = (2/3) * z * Sqrt(z). ///

/// The value to compute the Airy function of. - /// If true, returns exponentially-scaled Airy function - /// - public static Complex AiryAi(Complex z, bool expScaled = false) + /// The option to set the scaling factor. + /// The Airy function Ai. + public static Complex AiryAi(Complex z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCairy(z) : Amos.Cairy(z); + return (scale == Scale.Exponential) ? Amos.ScaledCairy(z) : Amos.Cairy(z); } /// - /// Airy function Ai(z). - ///

- /// If expScaled is true, returns Exp(zta) * Ai(z), where zta = (2/3) * z * Sqrt(z). + /// Returns the exponentially scaled Airy function Ai. + /// ScaledAiryAi(z) is given by Exp(zta) * AiryAi(z), where zta = (2/3) * z * Sqrt(z). ///

/// The value to compute the Airy function of. - /// If true, returns exponentially-scaled Airy function - /// - public static double AiryAi(double z, bool expScaled = false) + /// The exponentially scaled Airy function Ai. + public static Complex ScaledAiryAi(Complex z) { - if (expScaled) - { - return Amos.ScaledCairy(z); - } - else - { - return AiryAi(new Complex(z, 0), expScaled).Real; - } + return Amos.ScaledCairy(z); } /// - /// Derivative of the Airy function Ai. - ///

- /// If expScaled is true, returns Exp(zta) * d/dz Ai(z), where zta = (2/3) * z * Sqrt(z). + /// Returns the Airy function Ai. + /// AiryAi(z) is a solution to the Airy equation, y'' - y * z = 0. + /// AiryAi(z, Scale.Exponential) returns Exp(zta) * AiryAi(z), where zta = (2/3) * z * Sqrt(z). + ///

+ /// The value to compute the Airy function of. + /// The option to set the scaling factor. + /// The Airy function Ai. + public static double AiryAi(double z, Scale scale = Scale.Unity) + { + return (scale == Scale.Exponential) ? Amos.ScaledCairy(z) : AiryAi(new Complex(z, 0), scale).Real; + } + + /// + /// Returns the exponentially scaled Airy function Ai. + /// ScaledAiryAi(z) is given by Exp(zta) * AiryAi(z), where zta = (2/3) * z * Sqrt(z). + /// + /// The value to compute the Airy function of. + /// The exponentially scaled Airy function Ai. + public static double ScaledAiryAi(double z) + { + return Amos.ScaledCairy(z); + } + + /// + /// Returns the derivative of the Airy function Ai. + /// AiryAiPrime(z) is defined as d/dz AiryAi(z). + /// AiryAiPrime(z, Scale.Exponential) returns Exp(zta) * AiryAiPrime(z), where zta = (2/3) * z * Sqrt(z). + /// + /// The value to compute the derivative of the Airy function of. + /// The option to set the scaling factor. + /// The derivative of the Airy function Ai. + public static Complex AiryAiPrime(Complex z, Scale scale = Scale.Unity) + { + return (scale == Scale.Exponential) ? Amos.ScaledCairyPrime(z) : Amos.CairyPrime(z); + } + + /// + /// Returns the exponentially scaled derivative of Airy function Ai + /// ScaledAiryAiPrime(z) is given by Exp(zta) * AiryAiPrime(z), where zta = (2/3) * z * Sqrt(z). + /// + /// The value to compute the derivative of the Airy function of. + /// The exponentially scaled derivative of Airy function Ai. + public static Complex ScaledAiryAiPrime(Complex z) + { + return Amos.ScaledCairyPrime(z); + } + + /// + /// Returns the derivative of the Airy function Ai. + /// AiryAiPrime(z) is defined as d/dz AiryAi(z). + /// AiryAiPrime(z, Scale.Exponential) returns Exp(zta) * AiryAiPrime(z), where zta = (2/3) * z * Sqrt(z). /// /// The value to compute the derivative of the Airy function of. - /// If true, returns exponentially-scaled Airy function - /// - public static Complex AiryAiPrime(Complex z, bool expScaled = false) + /// The option to set the scaling factor. + /// The derivative of the Airy function Ai. + public static double AiryAiPrime(double z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCairyPrime(z) : Amos.CairyPrime(z); + return (scale == Scale.Exponential) ? Amos.ScaledCairyPrime(z) : AiryAiPrime(new Complex(z, 0), scale).Real; } /// - /// Derivative of the Airy function Ai. - ///

- /// If expScaled is true, returns Exp(zta) * d/dz Ai(z), where zta = (2/3) * z * Sqrt(z). + /// Returns the expoenntially scaled derivative of the Airy function Ai. + /// ScaledAiryAiPrime(z) is given by Exp(zta) * AiryAiPrime(z), where zta = (2/3) * z * Sqrt(z). ///

/// The value to compute the derivative of the Airy function of. - /// If true, returns exponentially-scaled Airy function - /// - public static double AiryAiPrime(double z, bool expScaled = false) + /// The expoenntially scaled derivative of the Airy function Ai. + public static double ScaledAiryAiPrime(double z) { - if (expScaled) - { - return Amos.ScaledCairyPrime(z); - } - else - { - return AiryAiPrime(new Complex(z, 0), expScaled).Real; - } + return Amos.ScaledCairyPrime(z); } /// - /// Airy function Bi(z). - ///

- /// If expScaled is true, returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real). + /// Returns the Airy function Bi. + /// AiryBi(z) is a solution to the Airy equation, y'' - y * z = 0. + /// AiryBi(z, Scale.Exponential) returns Exp(-Abs(zta.Real)) * AiryBi(z) where zta = (2 / 3) * z * Sqrt(z). ///

/// The value to compute the Airy function of. - /// If true, returns exponentially-scaled Airy function - /// - public static Complex AiryBi(Complex z, bool expScaled = false) + /// The option to set the scaling factor. + /// The Airy function Bi. + public static Complex AiryBi(Complex z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCbiry(z) : Amos.Cbiry(z); + return (scale == Scale.Exponential) ? Amos.ScaledCbiry(z) : Amos.Cbiry(z); } /// - /// Airy function Bi(x). - ///

- /// If expScaled is true, returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real). + /// Returns the exponentially scaled Airy function Bi. + /// ScaledAiryBi(z) is given by Exp(-Abs(zta.Real)) * AiryBi(z) where zta = (2 / 3) * z * Sqrt(z). ///

/// The value to compute the Airy function of. - /// If true, returns exponentially-scaled Airy function - /// - public static double AiryBi(double z, bool expScaled = false) + /// The exponentially scaled Airy function Bi(z). + public static Complex ScaledAiryBi(Complex z) + { + return Amos.ScaledCbiry(z); + } + + /// + /// Returns the Airy function Bi. + /// AiryBi(z) is a solution to the Airy equation, y'' - y * z = 0. + /// AiryBi(z, Scale.Exponential) returns Exp(-Abs(zta.Real)) * AiryBi(z) where zta = (2 / 3) * z * Sqrt(z). + /// + /// The value to compute the Airy function of. + /// The option to set the scaling factor. + /// The Airy function Bi. + public static double AiryBi(double z, Scale scale = Scale.Unity) + { + return AiryBi(new Complex(z, 0), scale).Real; + } + + /// + /// Returns the exponentially scaled Airy function Bi. + /// ScaledAiryBi(z) is given by Exp(-Abs(zta.Real)) * AiryBi(z) where zta = (2 / 3) * z * Sqrt(z). + /// + /// The value to compute the Airy function of. + /// The exponentially scaled Airy function Bi. + public static double ScaledAiryBi(double z) + { + return AiryBi(new Complex(z, 0), Scale.Exponential).Real; + } + + /// + /// Returns the derivative of the Airy function Bi. + /// AiryBiPrime(z) is defined as d/dz AiryBi(z). + /// AiryBiPrime(z, Scale.Exponential) returns Exp(-Abs(zta.Real)) * AiryBiPrime(z) where zta = (2 / 3) * z * Sqrt(z). + /// + /// The value to compute the derivative of the Airy function of. + /// The option to set the scaling factor. + /// The derivative of the Airy function Bi. + public static Complex AiryBiPrime(Complex z, Scale scale = Scale.Unity) + { + return (scale == Scale.Exponential) ? Amos.ScaledCbiryPrime(z) : Amos.CbiryPrime(z); + } + + /// + /// Returns the exponentially scaled derivative of the Airy function Bi. + /// ScaledAiryBiPrime(z) is given by Exp(-Abs(zta.Real)) * AiryBiPrime(z) where zta = (2 / 3) * z * Sqrt(z). + /// + /// The value to compute the derivative of the Airy function of. + /// The exponentially scaled derivative of the Airy function Bi. + public static Complex ScaledAiryBiPrime(Complex z) { - return AiryBi(new Complex(z, 0), expScaled).Real; + return Amos.ScaledCbiryPrime(z); } /// - /// Derivative of the Airy function Bi(z). - ///

- /// If expScaled is true, returns Exp(-axzta) * d/dz Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real). + /// Returns the derivative of the Airy function Bi. + /// AiryBiPrime(z) is defined as d/dz AiryBi(z). + /// AiryBiPrime(z, Scale.Exponential) returns Exp(-Abs(zta.Real)) * AiryBiPrime(z) where zta = (2 / 3) * z * Sqrt(z). ///

/// The value to compute the derivative of the Airy function of. - /// If true, returns exponentially-scaled Airy function - /// - public static Complex AiryBiPrime(Complex z, bool expScaled = false) + /// The option to set the scaling factor. + /// The derivative of the Airy function Bi. + public static double AiryBiPrime(double z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCbiryPrime(z) : Amos.CbiryPrime(z); + return AiryBiPrime(new Complex(z, 0), scale).Real; } /// - /// Derivative of the Airy function Bi(z). - ///

- /// If expScaled is true, returns Exp(-axzta) * d/dz Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real). + /// Returns the exponentially scaled derivative of the Airy function Bi. + /// ScaledAiryBiPrime(z) is given by Exp(-Abs(zta.Real)) * AiryBiPrime(z) where zta = (2 / 3) * z * Sqrt(z). ///

/// The value to compute the derivative of the Airy function of. - /// If true, returns exponentially-scaled Airy function - /// - public static double AiryBiPrime(double z, bool expScaled = false) + /// The exponentially scaled derivative of the Airy function Bi. + public static double ScaledAiryBiPrime(double z) { - return AiryBiPrime(new Complex(z, 0), expScaled).Real; + return AiryBiPrime(new Complex(z, 0), Scale.Exponential).Real; } } } diff --git a/src/Numerics/SpecialFunctions/Bessel.cs b/src/Numerics/SpecialFunctions/Bessel.cs index d9c15f7f..0bbb8945 100644 --- a/src/Numerics/SpecialFunctions/Bessel.cs +++ b/src/Numerics/SpecialFunctions/Bessel.cs @@ -8,122 +8,211 @@ namespace MathNet.Numerics public static partial class SpecialFunctions { /// - /// Bessel function of the first kind, J(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(y)) * J(v, z) where y = z.Imaginary. + /// Returns the Bessel function of the first kind. + /// BesselJ(n, z) is a solution to the Bessel differential equation. + /// BesselJ(n, z, Scale.Exponential) returns Exp(-Abs(z.Imaginary)) * BesselJ(n, z). ///

- /// The order of the Bessel function + /// The order of the Bessel function. /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Bessel function - /// - public static Complex BesselJ(double v, Complex z, bool expScaled = false) + /// The option to set the scaling factor. + /// The Bessel function of the first kind. + public static Complex BesselJ(double n, Complex z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCbesj(v, z) : Amos.Cbesj(v, z); + return (scale == Scale.Exponential) ? Amos.ScaledCbesj(n, z) : Amos.Cbesj(n, z); } /// - /// Bessel function of the first kind, J(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(y)) * J(v, z) where y = z.Imaginary. + /// Returns the exponentially scaled Bessel function of the first kind. + /// ScaledBesselJ(n, z) is given by Exp(-Abs(z.Imaginary)) * BesselJ(n, z). ///

- /// The order of the Bessel function + /// The order of the Bessel function. /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Bessel function - /// - public static double BesselJ(double v, double z, bool expScaled = false) + /// The exponentially scaled Bessel function of the first kind. + public static Complex ScaledBesselJ(double n, Complex z) { - return (expScaled) ? Amos.ScaledCbesj(v, z) : Amos.Cbesj(v, z); + return Amos.ScaledCbesj(n, z); } /// - /// Bessel function of the second kind, Y(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary. + /// Returns the Bessel function of the first kind. + /// BesselJ(n, z) is a solution to the Bessel differential equation. + /// BesselJ(n, z, Scale.Exponential) returns Exp(-Abs(z.Imaginary)) * J(n, z). ///

- /// The order of the Bessel function + /// The order of the Bessel function. /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Bessel function - /// - public static Complex BesselY(double v, Complex z, bool expScaled = false) + /// The option to set the scaling factor. + /// The Bessel function of the first kind. + public static double BesselJ(double n, double z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCbesy(v, z) : Amos.Cbesy(v, z); + return (scale == Scale.Exponential) ? Amos.ScaledCbesj(n, z) : Amos.Cbesj(n, z); } /// - /// Bessel function of the second kind, Y(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary. + /// Returns the exponentially scaled Bessel function of the first kind. + /// ScaledBesselJ(n, z) is given by Exp(-Abs(z.Imaginary)) * BesselJ(n, z). ///

- /// The order of the Bessel function + /// The order of the Bessel function. /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Bessel function - /// - public static double BesselY(double v, double z, bool expScaled = false) + /// The exponentially scaled Bessel function of the first kind. + public static double ScaledBesselJ(double n, double z) { - return (expScaled) ? Amos.ScaledCbesy(v, z) : Amos.Cbesy(v, z); + return Amos.ScaledCbesj(n, z); } /// - /// Modified Bessel function of the first kind, I(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(x)) * I(v, z) where x = z.Real. + /// Returns the Bessel function of the second kind. + /// BesselY(n, z) is a solution to the Bessel differential equation. + /// BesselY(n, z, Scale.Exponential) returns Exp(-Abs(z.Imaginary)) * BesselY(n, z). ///

- /// The order of the Bessel function + /// The order of the Bessel function. /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Bessel function - /// - public static Complex BesselI(double v, Complex z, bool expScaled = false) + /// The option to set the scaling factor. + /// The Bessel function of the second kind. + public static Complex BesselY(double n, Complex z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCbesi(v, z) : Amos.Cbesi(v, z); + return (scale == Scale.Exponential) ? Amos.ScaledCbesy(n, z) : Amos.Cbesy(n, z); } /// - /// Modified Bessel function of the first kind, I(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(x)) * I(v, z) where x = z.Real. + /// Returns the exponentially scaled Bessel function of the second kind. + /// ScaledBesselY(n, z) is given by Exp(-Abs(z.Imaginary)) * Y(n, z). ///

- /// The order of the Bessel function + /// The order of the Bessel function. /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Bessel function - /// - public static double BesselI(double v, double z, bool expScaled = false) + /// The exponentially scaled Bessel function of the second kind. + public static Complex ScaledBesselY(double n, Complex z) { - if (expScaled) - { - return Amos.ScaledCbesi(v, z); - } - else - { - return BesselI(v, new Complex(z, 0), expScaled).Real; - } + return Amos.ScaledCbesy(n, z); } /// - /// Modified Bessel function of the second kind, K(v, z). - ///

- /// If expScaled is true, returns Exp(z) * K(v, z). + /// Returns the Bessel function of the second kind. + /// BesselY(n, z) is a solution to the Bessel differential equation. + /// BesselY(n, z, Scale.Exponential) returns Exp(-Abs(z.Imaginary)) * BesselY(n, z). ///

- /// The order of the Bessel function + /// The order of the Bessel function. /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Bessel function - /// - public static Complex BesselK(double v, Complex z, bool expScaled = false) + /// The option to set the scaling factor. + /// The Bessel function of the second kind. + public static double BesselY(double n, double z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCbesk(v, z) : Amos.Cbesk(v, z); + return (scale == Scale.Exponential) ? Amos.ScaledCbesy(n, z) : Amos.Cbesy(n, z); } /// - /// Modified Bessel function of the second kind, K(v, z). - ///

- /// If expScaled is true, returns Exp(z) * K(v, z). + /// Returns the exponentially scaled Bessel function of the second kind. + /// ScaledBesselY(n, z) is given by Exp(-Abs(z.Imaginary)) * BesselY(n, z). ///

- /// The order of the Bessel function + /// The order of the Bessel function. /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Bessel function - /// - public static double BesselK(double v, double z, bool expScaled = false) + /// The exponentially scaled Bessel function of the second kind. + public static double ScaledBesselY(double n, double z) { - return (expScaled) ? Amos.ScaledCbesk(v, z) : Amos.Cbesk(v, z); + return Amos.ScaledCbesy(n, z); + } + + /// + /// Returns the modified Bessel function of the first kind. + /// BesselI(n, z) is a solution to the modified Bessel differential equation. + /// BesselI(n, z, Scale.Exponential) returns Exp(-Abs(z.Real)) * BesselI(n, z). + /// + /// The order of the modified Bessel function. + /// The value to compute the modified Bessel function of. + /// The option to set the scaling factor. + /// The modified Bessel function of the first kind. + public static Complex BesselI(double n, Complex z, Scale scale = Scale.Unity) + { + return (scale == Scale.Exponential) ? Amos.ScaledCbesi(n, z) : Amos.Cbesi(n, z); + } + + /// + /// Returns the exponentially scaled modified Bessel function of the first kind. + /// ScaledBesselI(n, z) is given by Exp(-Abs(z.Real)) * BesselI(n, z). + /// + /// The order of the modified Bessel function. + /// The value to compute the modified Bessel function of. + /// The exponentially scaled modified Bessel function of the first kind. + public static Complex ScaledBesselI(double n, Complex z) + { + return Amos.ScaledCbesi(n, z); + } + + /// + /// Returns the modified Bessel function of the first kind. + /// BesselI(n, z) is a solution to the modified Bessel differential equation. + /// BesselI(n, z, Scale.Exponential) returns Exp(-Abs(z.Real)) * BesselI(n, z). + /// + /// The order of the modified Bessel function. + /// The value to compute the modified Bessel function of. + /// The option to set the scaling factor. + /// The modified Bessel function of the first kind. + public static double BesselI(double n, double z, Scale scale = Scale.Unity) + { + return (scale == Scale.Exponential) ? Amos.ScaledCbesi(n, z) : BesselI(n, new Complex(z, 0), scale).Real; + } + + /// + /// Returns the exponentially scaled modified Bessel function of the first kind. + /// ScaledBesselI(n, z) is given by Exp(-Abs(z.Real)) * BesselI(n, z). + /// + /// The order of the modified Bessel function. + /// The value to compute the modified Bessel function of. + /// The exponentially scaled modified Bessel function of the first kind. + public static double ScaledBesselI(double n, double z) + { + return Amos.ScaledCbesi(n, z); + } + + /// + /// Returns the modified Bessel function of the second kind. + /// BesselK(n, z) is a solution to the modified Bessel differential equation. + /// BesselK(n, z, Scale.Exponential) returns Exp(z) * BesselK(n, z). + /// + /// The order of the modified Bessel function. + /// The value to compute the modified Bessel function of. + /// The option to set the scaling factor. + /// The modified Bessel function of the second kind. + public static Complex BesselK(double n, Complex z, Scale scale = Scale.Unity) + { + return (scale == Scale.Exponential) ? Amos.ScaledCbesk(n, z) : Amos.Cbesk(n, z); + } + + /// + /// Returns the exponentially scaled modified Bessel function of the second kind. + /// ScaledBesselK(n, z) is given by Exp(z) * BesselK(n, z). + /// + /// The order of the modified Bessel function. + /// The value to compute the modified Bessel function of. + /// The exponentially scaled modified Bessel function of the second kind. + public static Complex ScaledBesselK(double n, Complex z) + { + return Amos.ScaledCbesk(n, z); + } + + /// + /// Returns the modified Bessel function of the second kind. + /// BesselK(n, z) is a solution to the modified Bessel differential equation. + /// BesselK(n, z, Scale.Exponential) returns Exp(z) * BesselK(n, z). + /// + /// The order of the modified Bessel function. + /// The value to compute the modified Bessel function of. + /// The option to set the scaling factor. + /// The modified Bessel function of the second kind. + public static double BesselK(double n, double z, Scale scale = Scale.Unity) + { + return (scale == Scale.Exponential) ? Amos.ScaledCbesk(n, z) : Amos.Cbesk(n, z); + } + + /// + /// Returns the exponentially scaled modified Bessel function of the second kind. + /// ScaledBesselK(n, z) is given by Exp(z) * BesselK(n, z). + /// + /// The order of the modified Bessel function. + /// The value to compute the modified Bessel function of. + /// The exponentially scaled modified Bessel function of the second kind. + public static double ScaledBesselK(double n, double z) + { + return Amos.ScaledCbesk(n, z); } } } diff --git a/src/Numerics/SpecialFunctions/Hankel.cs b/src/Numerics/SpecialFunctions/Hankel.cs index 3ad42386..6eaf4156 100644 --- a/src/Numerics/SpecialFunctions/Hankel.cs +++ b/src/Numerics/SpecialFunctions/Hankel.cs @@ -8,59 +8,55 @@ namespace MathNet.Numerics public static partial class SpecialFunctions { /// - /// Hankel function of the first kind, H1(n, z). - ///

- /// If expScaled is true, returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1). + /// Returns the Hankel function of the first kind. + /// HankelH1(n, z) is defined as BesselJ(n, z) + j * BesselY(n, z). + /// HankelH1(n, z, Scale.Exponential) returns Exp(-z * j) * HankelH1(n, z) where j = Sqrt(-1). ///

- /// The order of the Bessel function - /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Hankel function - /// - public static Complex HankelH1(double n, Complex z, bool expScaled = false) + /// The order of the Hankel function. + /// The value to compute the Hankel function of. + /// The option to set the scaling factor. + /// The Hankel function of the first kind. + public static Complex HankelH1(double n, Complex z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCbesh1(n, z) : Amos.Cbesh1(n, z); + return (scale == Scale.Exponential) ? Amos.ScaledCbesh1(n, z) : Amos.Cbesh1(n, z); } /// - /// Hankel function of the first kind, H1(n, z). - ///

- /// If expScaled is true, returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1). + /// Returns the exponentially scaled Hankel function of the first kind. + /// ScaledHankelH1(n, z) is given by Exp(-z * j) * HankelH1(n, z) where j = Sqrt(-1). ///

- /// The order of the Bessel function - /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Hankel function - /// - public static double HankelH1(double n, double z, bool expScaled = false) + /// The order of the Hankel function. + /// The value to compute the Hankel function of. + /// The exponentially scaled Hankel function of the first kind. + public static Complex ScaledHankelH1(double n, Complex z) { - return HankelH1(n, new Complex(z, 0), expScaled).Real; + return Amos.ScaledCbesh1(n, z); } /// - /// Hankel function of the second kind, H2(n, z). - ///

- /// If expScaled is true, returns Exp(z * j) * H2(n, z) where j = Sqrt(-1). + /// Returns the Hankel function of the second kind. + /// HankelH2(n, z) is defined as BesselJ(n, z) - j * BesselY(n, z). + /// HankelH2(n, z, Scale.Exponential) returns Exp(z * j) * HankelH2(n, z) where j = Sqrt(-1). ///

- /// The order of the Hankel function - /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Hankel function - /// - public static Complex HankelH2(double n, Complex z, bool expScaled = false) + /// The order of the Hankel function. + /// The value to compute the Hankel function of. + /// The option to set the scaling factor. + /// The Hankel function of the second kind. + public static Complex HankelH2(double n, Complex z, Scale scale = Scale.Unity) { - return (expScaled) ? Amos.ScaledCbesh2(n, z) : Amos.Cbesh2(n, z); + return (scale == Scale.Exponential) ? Amos.ScaledCbesh2(n, z) : Amos.Cbesh2(n, z); } /// - /// Hankel function of the second kind, H2(n, z). - ///

- /// If expScaled is true, returns Exp(z * j) * H2(n, z) where j = Sqrt(-1). + /// Returns the exponentially scaled Hankel function of the second kind. + /// ScaledHankelH2(n, z) is given by Exp(z * j) * HankelH2(n, z) where j = Sqrt(-1). ///

- /// The order of the Bessel function - /// The value to compute the Bessel function of. - /// If true, returns exponentially-scaled Hankel function - /// - public static double HankelH2(double n, double z, bool expScaled = false) + /// The order of the Hankel function. + /// The value to compute the Hankel function of. + /// The exponentially scaled Hankel function of the second kind. + public static Complex ScaledHankelH2(double n, Complex z) { - return HankelH2(n, new Complex(z, 0), expScaled).Real; + return Amos.ScaledCbesh2(n, z); } } } diff --git a/src/Numerics/SpecialFunctions/Options.cs b/src/Numerics/SpecialFunctions/Options.cs new file mode 100644 index 00000000..6342f0c7 --- /dev/null +++ b/src/Numerics/SpecialFunctions/Options.cs @@ -0,0 +1,23 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics +{ + public static partial class SpecialFunctions + { + public enum Scale + { + /// + /// For Bessel-related functions, no scaling factor is applied. + /// + Unity = 0, + + /// + /// For Bessel-related functions, exponential scaling is applied. + /// + Exponential = 1 + } + } +} diff --git a/src/Numerics/SpecialFunctions/SphericalBessel.cs b/src/Numerics/SpecialFunctions/SphericalBessel.cs index afa1a1b8..2858cc5a 100644 --- a/src/Numerics/SpecialFunctions/SphericalBessel.cs +++ b/src/Numerics/SpecialFunctions/SphericalBessel.cs @@ -9,67 +9,121 @@ namespace MathNet.Numerics public static partial class SpecialFunctions { /// - /// Spherical Bessel function of the first kind, j(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(y)) * j(v, z) where y = z.Imaginary. + /// Returns the spherical Bessel function of the first kind. + /// SphericalBesselJ(n, z) is given by Sqrt(pi/2) / Sqrt(z) * BesselJ(n + 1/2, z). ///

- /// The order of the spherical Bessel function + /// The order of the spherical Bessel function. /// The value to compute the spherical Bessel function of. - /// If true, returns exponentially-scaled spherical Bessel function - /// - public static Complex SphericalBesselJ(double v, Complex z, bool expScaled = false) + /// The spherical Bessel function of the first kind. + public static Complex SphericalBesselJ(double n, Complex z) { - const double rthpi = 1.2533141373155002512; //sqrt(pi/2) + if (double.IsNaN(n) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) + { + return new Complex(double.NaN, double.NaN); + } + + if (double.IsInfinity(z.Real)) + { + return (z.Imaginary == 0) ? Complex.Zero : new Complex(double.PositiveInfinity, double.PositiveInfinity); + } - return rthpi * BesselJ(v + 0.5, z, expScaled) / Complex.Sqrt(z); + if (z.Real == 0 && z.Imaginary == 0) + { + return (n == 0) ? 1 : 0; + } + + return Constants.SqrtPiOver2 * BesselJ(n + 0.5, z, Scale.Unity) / Complex.Sqrt(z); } /// - /// Spherical Bessel function of the first kind, j(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(y)) * j(v, z) where y = z.Imaginary. + /// Returns the spherical Bessel function of the first kind. + /// SphericalBesselJ(n, z) is given by Sqrt(pi/2) / Sqrt(z) * BesselJ(n + 1/2, z). ///

- /// The order of the spherical Bessel function + /// The order of the spherical Bessel function. /// The value to compute the spherical Bessel function of. - /// If true, returns exponentially-scaled spherical Bessel function - /// - public static double SphericalBesselJ(double v, double z, bool expScaled = false) + /// The spherical Bessel function of the first kind. + public static double SphericalBesselJ(double n, double z) { - const double rthpi = 1.2533141373155002512; //sqrt(pi/2) + if (double.IsNaN(n) || double.IsNaN(z)) + { + return double.NaN; + } + + if (n < 0) + { + return double.NaN; + } + + if (double.IsInfinity(z)) + { + return 0; + } - return rthpi * BesselJ(v + 0.5, z, expScaled) / Math.Sqrt(z); + if (z == 0) + { + return (n == 0) ? 1 : 0; + } + + return Constants.SqrtPiOver2 * BesselJ(n + 0.5, z, Scale.Unity) / Math.Sqrt(z); } /// - /// Spherical Bessel function of the second kind, y(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(y)) * y(v, z) where y = z.Imaginary. + /// Returns the spherical Bessel function of the second kind. + /// SphericalBesselY(n, z) is given by Sqrt(pi/2) / Sqrt(z) * BesselY(n + 1/2, z). ///

- /// The order of the spherical Bessel function + /// The order of the spherical Bessel function. /// The value to compute the spherical Bessel function of. - /// If true, returns exponentially-scaled spherical Bessel function - /// - public static Complex SphericalBesselY(double v, Complex z, bool expScaled = false) + /// The spherical Bessel function of the second kind. + public static Complex SphericalBesselY(double n, Complex z) { - const double rthpi = 1.2533141373155002512; //sqrt(pi/2) + if (double.IsNaN(n) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary)) + { + return new Complex(double.NaN, double.NaN); + } + + if (double.IsInfinity(z.Real)) + { + return (z.Imaginary == 0) ? Complex.Zero : new Complex(double.PositiveInfinity, double.PositiveInfinity); + } - return rthpi * BesselY(v + 0.5, z, expScaled) / Complex.Sqrt(z); + if (z.Real == 0 && z.Imaginary == 0) + { + return new Complex(double.NaN, double.NaN); + } + + return Constants.SqrtPiOver2 * BesselY(n + 0.5, z, Scale.Unity) / Complex.Sqrt(z); } /// - /// Spherical Bessel function of the second kind, y(v, z). - ///

- /// If expScaled is true, returns Exp(-Abs(y)) * y(v, z) where y = z.Imaginary. + /// Returns the spherical Bessel function of the second kind. + /// SphericalBesselY(n, z) is given by Sqrt(pi/2) / Sqrt(z) * BesselY(n + 1/2, z). ///

- /// The order of the spherical Bessel function + /// The order of the spherical Bessel function. /// The value to compute the spherical Bessel function of. - /// If true, returns exponentially-scaled spherical Bessel function - /// - public static double SphericalBesselY(double v, double z, bool expScaled = false) + /// The spherical Bessel function of the second kind. + public static double SphericalBesselY(double n, double z) { - const double rthpi = 1.2533141373155002512; //sqrt(pi/2) + if (double.IsNaN(n) || double.IsNaN(z)) + { + return double.NaN; + } + + if (n < 0) + { + return double.NaN; + } + + if (double.IsInfinity(z)) + { + return 0; + } + + if (z == 0) + { + return double.NegativeInfinity; + } - return rthpi * BesselY(v + 0.5, z, expScaled) / Math.Sqrt(z); + return Constants.SqrtPiOver2 * BesselY(n + 0.5, z, Scale.Unity) / Math.Sqrt(z); } } } From 4cd8dbea6b3673f47ccab946ad8f3f51903a8454 Mon Sep 17 00:00:00 2001 From: diluculo Date: Thu, 8 Nov 2018 19:45:59 +0900 Subject: [PATCH 2/3] UnitTests: replaced Evaluate.Polynomial with Polynomial.Evaluate. --- .../SpecialFunctionsTests/BesselTests.cs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs b/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs index f8dd3451..378b9481 100644 --- a/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs +++ b/src/Numerics.Tests/SpecialFunctionsTests/BesselTests.cs @@ -17,7 +17,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests public void BesselJ0Approx([Range(-3, 3, 0.25)] double x) { // Approx by Abramowitz/Stegun 9.4.1 - Assert.AreEqual(Evaluate.Polynomial(x / 3.0, 1.0, 0.0, -2.2499997, 0.0, 1.2656208, 0.0, -0.3163866, 0.0, 0.0444479, 0.0, -0.0039444, 0.0, 0.0002100), SpecialFunctions.BesselJ(0, x), 1e-7); + Assert.AreEqual(Polynomial.Evaluate(x / 3.0, 1.0, 0.0, -2.2499997, 0.0, 1.2656208, 0.0, -0.3163866, 0.0, 0.0444479, 0.0, -0.0039444, 0.0, 0.0002100), SpecialFunctions.BesselJ(0, x), 1e-7); } [TestCase(0, 0.0, 0.0, 1.0000000000000000000, 0.0000000000000000000, 14)] @@ -57,7 +57,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests { // Approx by Abramowitz/Stegun 9.4.2 Assert.AreEqual( - Evaluate.Polynomial(x / 3.0, 2.0 / Math.PI * Math.Log(x / 2.0) * SpecialFunctions.BesselJ(0, x) + 0.36746691, 0.0, 0.60559366, 0.0, -0.74350384, 0.0, 0.25300117, 0.0, -0.04261214, 0.0, 0.00427916, 0.0, -0.00024846), + Polynomial.Evaluate(x / 3.0, 2.0 / Math.PI * Math.Log(x / 2.0) * SpecialFunctions.BesselJ(0, x) + 0.36746691, 0.0, 0.60559366, 0.0, -0.74350384, 0.0, 0.25300117, 0.0, -0.04261214, 0.0, 0.00427916, 0.0, -0.00024846), SpecialFunctions.BesselY(0, x), 1e-7); } @@ -91,7 +91,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests public void BesselI0Approx([Range(-3.75, 3.75, 0.25)] double x) { // Approx by Abramowitz/Stegun 9.8.1 - Assert.AreEqual(Evaluate.Polynomial(x / 3.75, 1.0, 0.0, 3.5156229, 0.0, 3.0899424, 0.0, 1.2067492, 0.0, 0.2659732, 0.0, 0.0360768, 0.0, 0.0045813), SpecialFunctions.BesselI(0, x), 1e-7); + Assert.AreEqual(Polynomial.Evaluate(x / 3.75, 1.0, 0.0, 3.5156229, 0.0, 3.0899424, 0.0, 1.2067492, 0.0, 0.2659732, 0.0, 0.0360768, 0.0, 0.0045813), SpecialFunctions.BesselI(0, x), 1e-7); } [TestCase(0.0, 1.0)] @@ -111,7 +111,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests public void BesselI1Approx([Range(-3.75, 3.75, 0.25)] double x) { // Approx by Abramowitz/Stegun 9.8.3 - Assert.AreEqual(Evaluate.Polynomial(x / 3.75, 0.5, 0.0, 0.87890594, 0.0, 0.51498869, 0.0, 0.15084934, 0.0, 0.02658733, 0.0, 0.00301532, 0.0, 0.00032411) * x, SpecialFunctions.BesselI(1, x), 1e-8); + Assert.AreEqual(Polynomial.Evaluate(x / 3.75, 0.5, 0.0, 0.87890594, 0.0, 0.51498869, 0.0, 0.15084934, 0.0, 0.02658733, 0.0, 0.00301532, 0.0, 0.00032411) * x, SpecialFunctions.BesselI(1, x), 1e-8); } [TestCase(0.0, 0.0)] @@ -177,7 +177,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests public void BesselK0Approx([Range(0.20, 2.0, 0.20)] double x) { // Approx by Abramowitz/Stegun 9.8.5 - Assert.AreEqual(Evaluate.Polynomial(x / 2.0, -Math.Log(x / 2.0) * SpecialFunctions.BesselI(0, x) - 0.57721566, 0.0, 0.42278420, 0.0, 0.23069756, 0.0, 0.03488590, 0.0, 0.00262698, 0.0, 0.00010750, 0.0, 0.00000740), SpecialFunctions.BesselK(0, x), 1e-8); + Assert.AreEqual(Polynomial.Evaluate(x / 2.0, -Math.Log(x / 2.0) * SpecialFunctions.BesselI(0, x) - 0.57721566, 0.0, 0.42278420, 0.0, 0.23069756, 0.0, 0.03488590, 0.0, 0.00262698, 0.0, 0.00010750, 0.0, 0.00000740), SpecialFunctions.BesselK(0, x), 1e-8); } [TestCase(1e-10, 23.14178244559887)] @@ -196,7 +196,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests public void BesselK1Approx([Range(0.20, 2.0, 0.20)] double x) { // Approx by Abramowitz/Stegun 9.8.7 - Assert.AreEqual(Evaluate.Polynomial(x / 2.0, x * Math.Log(x / 2.0) * SpecialFunctions.BesselI(1, x) + 1.0, 0.0, 0.15443144, 0.0, -0.67278579, 0.0, -0.18156897, 0.0, -0.01919402, 0.0, -0.00110404, 0.0, -0.00004686), SpecialFunctions.BesselK(1, x) * x, 1e-8); + Assert.AreEqual(Polynomial.Evaluate(x / 2.0, x * Math.Log(x / 2.0) * SpecialFunctions.BesselI(1, x) + 1.0, 0.0, 0.15443144, 0.0, -0.67278579, 0.0, -0.18156897, 0.0, -0.01919402, 0.0, -0.00110404, 0.0, -0.00004686), SpecialFunctions.BesselK(1, x) * x, 1e-8); } [TestCase(1e-10, 1.0e+10)] From 284ca1ee5b387023e3f6576a4c492972ba1011fe Mon Sep 17 00:00:00 2001 From: diluculo Date: Fri, 9 Nov 2018 14:37:20 +0900 Subject: [PATCH 3/3] Special Functions: added Kelvin functions --- .../SpecialFunctionsTests/KelvinTests.cs | 135 +++++++++ src/Numerics/SpecialFunctions/Kelvin.cs | 264 ++++++++++++++++++ 2 files changed, 399 insertions(+) create mode 100644 src/Numerics.Tests/SpecialFunctionsTests/KelvinTests.cs create mode 100644 src/Numerics/SpecialFunctions/Kelvin.cs diff --git a/src/Numerics.Tests/SpecialFunctionsTests/KelvinTests.cs b/src/Numerics.Tests/SpecialFunctionsTests/KelvinTests.cs new file mode 100644 index 00000000..96857265 --- /dev/null +++ b/src/Numerics.Tests/SpecialFunctionsTests/KelvinTests.cs @@ -0,0 +1,135 @@ +using NUnit.Framework; +using System; + +namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests +{ + /// + /// Kelvin functions tests. + /// + [TestFixture, Category("Functions")] + public class KelvinTests + { + [Test] + public void KelvinBerApprox([Range(-8, 8, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.11.1 + Assert.AreEqual(Polynomial.Evaluate(x / 8.0, + 1.0, + 0.0, 0.0, 0.0, -64.0, 0.0, + 0.0, 0.0, 113.77777774, 0.0, 0.0, + 0.0, -32.36345652, 0.0, 0.0, 0.0, + 2.64191397, 0.0, 0.0, 0.0, -0.08349609, + 0.0, 0.0, 0.0, 0.00122552, 0.0, + 0.0, 0.0, -0.00000901), SpecialFunctions.KelvinBer(x), 1e-9); + } + + [Test] + public void KelvinBeiApprox([Range(-8, 8, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.11.2 + Assert.AreEqual(Polynomial.Evaluate(x / 8.0, + 0.0, + 0.0, 16.0, 0.0, 0.0, 0.0, + -113.77777774, 0.0, 0.0, 0.0, 72.81777742, + 0.0, 0.0, 0.0, -10.56765779, 0.0, + 0.0, 0.0, 0.52185615, 0.0, 0.0, + 0.0, -0.01103667, 0.0, 0.0, 0.0, + 0.00011346), + SpecialFunctions.KelvinBei(x), 6e-9); + } + + [Test] + public void KelvinKerApprox([Range(0.25, 8, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.11.3 + Assert.AreEqual( + Polynomial.Evaluate(x / 8.0, + -Math.Log(x / 2.0) * SpecialFunctions.KelvinBer(x) + SpecialFunctions.KelvinBei(x) * Constants.PiOver4 - 0.57721566, + 0.0, 0.0, 0.0, -59.05819744, 0.0, + 0.0, 0.0, 171.36272133, 0.0, 0.0, + 0.0, -60.60977451, 0.0, 0.0, 0.0, + 5.65539121, 0.0, 0.0, 0.0, -0.19636347, + 0.0, 0.0, 0.0, 0.00309699, 0.0, + 0.0, 0.0, -0.00002458), + SpecialFunctions.KelvinKer(x), 1e-8); + } + + [Test] + public void KelvinKeiApprox([Range(0.25, 8, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.11.4 + Assert.AreEqual( + -Math.Log(x / 2.0) * SpecialFunctions.KelvinBei(x) - Constants.PiOver4 * SpecialFunctions.KelvinBer(x) + + Polynomial.Evaluate(x / 8.0, + 0.0, + 0.0, 6.76454936, 0.0, 0.0, 0.0, + -142.91827687, 0.0, 0.0, 0.0, 124.23569650, + 0.0, 0.0, 0.0, -21.30060904, 0.0, + 0.0, 0.0, 1.17509064, 0.0, 0.0, + 0.0, -0.02695875, 0.0, 0.0, 0.0, + 0.00029532), + SpecialFunctions.KelvinKei(x), 3e-9); + } + + [Test] + public void KelvinBerPrimeApprox([Range(-8, 8, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.11.5 + Assert.AreEqual(x * Polynomial.Evaluate(x / 8.0, + 0.0, + 0.0, -4.0, 0.0, 0.0, 0.0, + 14.22222222, 0.0, 0.0, 0.0, -6.06814810, + 0.0, 0.0, 0.0, 0.66047849, 0.0, + 0.0, 0.0, -0.02609253, 0.0, 0.0, + 0.0, 0.00045957, 0.0, 0.0, 0.0, + -0.00000394), SpecialFunctions.KelvinBerPrime(x), 2.1e-8); + } + + [Test] + public void KelvinBeiPrimeApprox([Range(-8, 8, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.11.6 + Assert.AreEqual(x * Polynomial.Evaluate(x / 8.0, + 0.5, + 0.0, 0.0, 0.0, -10.66666666, 0.0, + 0.0, 0.0, 11.37777772, 0.0, 0.0, + 0.0, -2.31167514, 0.0, 0.0, 0.0, + 0.14677204, 0.0, 0.0, 0.0, -0.00379386, + 0.0, 0.0, 0.0, 0.00004609), + SpecialFunctions.KelvinBeiPrime(x), 7e-8); + } + + [Test] + public void KelvinKerPrimeApprox([Range(0.25, 8, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.11.7 + Assert.AreEqual( + -Math.Log(x / 2.0) * SpecialFunctions.KelvinBerPrime(x) - SpecialFunctions.KelvinBer(x) / x + Constants.PiOver4 * SpecialFunctions.KelvinBeiPrime(x) + + x * Polynomial.Evaluate(x / 8.0, + 0.0, + 0.0, -3.69113734, 0.0, 0.0, 0.0, + 21.42034017, 0.0, 0.0, 0.0, -11.36433272, + 0.0, 0.0, 0.0, 1.41384780, 0.0, + 0.0, 0.0, -0.06136358, 0.0, 0.0, + 0.0, 0.00116137, 0.0, 0.0, 0.0, + -0.00001075), + SpecialFunctions.KelvinKerPrime(x), 8e-8); + } + + [Test] + public void KelvinKeiPrimeApprox([Range(0.25, 8, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.11.8 + Assert.AreEqual( + -Math.Log(x / 2.0) * SpecialFunctions.KelvinBeiPrime(x) - SpecialFunctions.KelvinBei(x) / x - Constants.PiOver4 * SpecialFunctions.KelvinBerPrime(x) + + x * Polynomial.Evaluate(x / 8.0, + 0.21139217, + 0.0, 0.0, 0.0, -13.39858846, 0.0, + 0.0, 0.0, 19.41182758, 0.0, 0.0, + 0.0, -4.65950823, 0.0, 0.0, 0.0, + 0.33049424, 0.0, 0.0, 0.0, -0.00926707, + 0.0, 0.0, 0.0, 0.00011997), + SpecialFunctions.KelvinKeiPrime(x), 7e-8); + } + } +} diff --git a/src/Numerics/SpecialFunctions/Kelvin.cs b/src/Numerics/SpecialFunctions/Kelvin.cs new file mode 100644 index 00000000..9c6e467a --- /dev/null +++ b/src/Numerics/SpecialFunctions/Kelvin.cs @@ -0,0 +1,264 @@ +using System; +using System.Numerics; + +namespace MathNet.Numerics +{ + /// + /// This partial implementation of the SpecialFunctions class contains all methods related to the modified Bessel function. + /// + public static partial class SpecialFunctions + { + /// + /// Returns the Kelvin function of the first kind. + /// KelvinBe(nu, x) is given by BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1). + /// KelvinBer(nu, x) and KelvinBei(nu, x) are the real and imaginary parts of the KelvinBe(nu, x) + /// + /// the order of the the Kelvin function. + /// The value to compute the Kelvin function of. + /// The Kelvin function of the first kind. + public static Complex KelvinBe(double nu, double x) + { + Complex ISqrtI = new Complex(-Constants.Sqrt1Over2, Constants.Sqrt1Over2); // j * sqrt(j) = (-1)^(3/4) = (-1 + j)/sqrt(2) + return BesselJ(nu, ISqrtI * x); + } + + /// + /// Returns the Kelvin function ber. + /// KelvinBer(nu, x) is given by the real part of BesselJ(nu, j * sqrt(j) * x) where j = sqrt(-1). + /// + /// the order of the the Kelvin function. + /// The value to compute the Kelvin function of. + /// The Kelvin function ber. + public static double KelvinBer(double nu, double x) + { + return KelvinBe(nu, x).Real; + } + + /// + /// Returns the Kelvin function ber. + /// KelvinBer(x) is given by the real part of BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1). + /// KelvinBer(x) is equivalent to KelvinBer(0, x). + /// + /// The value to compute the Kelvin function of. + /// The Kelvin function ber. + public static double KelvinBer(double x) + { + return KelvinBe(0, x).Real; + } + + /// + /// Returns the Kelvin function bei. + /// KelvinBei(nu, x) is given by the imaginary part of BesselJ(nu, j * sqrt(j) * x) where j = sqrt(-1). + /// + /// the order of the the Kelvin function. + /// The value to compute the Kelvin function of. + /// The Kelvin function bei. + public static double KelvinBei(double nu, double x) + { + return KelvinBe(nu, x).Imaginary; + } + + /// + /// Returns the Kelvin function bei. + /// KelvinBei(x) is given by the imaginary part of BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1). + /// KelvinBei(x) is equivalent to KelvinBei(0, x). + /// + /// The value to compute the Kelvin function of. + /// The Kelvin function bei. + public static double KelvinBei(double x) + { + return KelvinBe(0, x).Imaginary; + } + + /// + /// Returns the derivative of the Kelvin function ber. + /// + /// The order of the Kelvin function. + /// The value to compute the derivative of the Kelvin function of. + /// the derivative of the Kelvin function ber + public static double KelvinBerPrime(double nu, double x) + { + const double inv2Sqrt2 = 0.35355339059327376220042218105242451964241796884424; // 1/(2 * sqrt(2)) + return inv2Sqrt2 * (-KelvinBer(nu - 1, x) + KelvinBer(nu + 1, x) - KelvinBei(nu - 1, x) + KelvinBei(nu + 1, x)); + } + + /// + /// Returns the derivative of the Kelvin function ber. + /// + /// The value to compute the derivative of the Kelvin function of. + /// The derivative of the Kelvin function ber. + public static double KelvinBerPrime(double x) + { + return KelvinBerPrime(0, x); + } + + /// + /// Returns the derivative of the Kelvin function bei. + /// + /// The order of the Kelvin function. + /// The value to compute the derivative of the Kelvin function of. + /// the derivative of the Kelvin function bei. + public static double KelvinBeiPrime(double nu, double x) + { + const double inv2Sqrt2 = 0.35355339059327376220042218105242451964241796884424; // 1/(2 * sqrt(2)) + return inv2Sqrt2 * (KelvinBer(nu - 1, x) - KelvinBer(nu + 1, x) - KelvinBei(nu - 1, x) + KelvinBei(nu + 1, x)); + } + + /// + /// Returns the derivative of the Kelvin function bei. + /// + /// The value to compute the derivative of the Kelvin function of. + /// The derivative of the Kelvin function bei. + public static double KelvinBeiPrime(double x) + { + return KelvinBeiPrime(0, x); + } + + /// + /// Returns the Kelvin function of the second kind + /// KelvinKe(nu, x) is given by Exp(-nu * pi * j / 2) * BesselK(nu, x * sqrt(j)) where j = sqrt(-1). + /// KelvinKer(nu, x) and KelvinKei(nu, x) are the real and imaginary parts of the KelvinBe(nu, x) + /// + /// The order of the Kelvin function. + /// The value to calculate the kelvin function of, + /// + public static Complex KelvinKe(double nu, double x) + { + Complex PiIOver2 = new Complex(0.0, Constants.PiOver2); // pi * I / 2 + Complex SqrtI = new Complex(Constants.Sqrt1Over2, Constants.Sqrt1Over2); // sqrt(j) = (-1)^(1/4) = (1 + j)/sqrt(2) + return Complex.Exp(-nu * PiIOver2) * BesselK(nu, SqrtI * x); + } + + /// + /// Returns the Kelvin function ker. + /// KelvinKer(nu, x) is given by the real part of Exp(-nu * pi * j / 2) * BesselK(nu, sqrt(j) * x) where j = sqrt(-1). + /// + /// the order of the the Kelvin function. + /// The non-negative real value to compute the Kelvin function of. + /// The Kelvin function ker. + public static double KelvinKer(double nu, double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + + return KelvinKe(nu, x).Real; + } + + /// + /// Returns the Kelvin function ker. + /// KelvinKer(x) is given by the real part of Exp(-nu * pi * j / 2) * BesselK(0, sqrt(j) * x) where j = sqrt(-1). + /// KelvinKer(x) is equivalent to KelvinKer(0, x). + /// + /// The non-negative real value to compute the Kelvin function of. + /// The Kelvin function ker. + public static double KelvinKer(double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + + return KelvinKe(0, x).Real; + } + + /// + /// Returns the Kelvin function kei. + /// KelvinKei(nu, x) is given by the imaginary part of Exp(-nu * pi * j / 2) * BesselK(nu, sqrt(j) * x) where j = sqrt(-1). + /// + /// the order of the the Kelvin function. + /// The non-negative real value to compute the Kelvin function of. + /// The Kelvin function kei. + public static double KelvinKei(double nu, double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + + return KelvinKe(nu, x).Imaginary; + } + + /// + /// Returns the Kelvin function kei. + /// KelvinKei(x) is given by the imaginary part of Exp(-nu * pi * j / 2) * BesselK(0, sqrt(j) * x) where j = sqrt(-1). + /// KelvinKei(x) is equivalent to KelvinKei(0, x). + /// + /// The non-negative real value to compute the Kelvin function of. + /// The Kelvin function kei. + public static double KelvinKei(double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + + return KelvinKe(0, x).Imaginary; + } + + /// + /// Returns the derivative of the Kelvin function ker. + /// + /// The order of the Kelvin function. + /// The non-negative real value to compute the derivative of the Kelvin function of. + /// The derivative of the Kelvin function ker. + public static double KelvinKerPrime(double nu, double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + + const double inv2Sqrt2 = 0.35355339059327376220042218105242451964241796884424; // 1/(2 * sqrt(2)) + return inv2Sqrt2 * (-KelvinKer(nu - 1, x) + KelvinKer(nu + 1, x) - KelvinKei(nu - 1, x) + KelvinKei(nu + 1, x)); + } + + /// + /// Returns the derivative of the Kelvin function ker. + /// + /// The value to compute the derivative of the Kelvin function of. + /// The derivative of the Kelvin function ker. + public static double KelvinKerPrime(double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + + return KelvinKerPrime(0, x); + } + + /// + /// Returns the derivative of the Kelvin function kei. + /// + /// The order of the Kelvin function. + /// The value to compute the derivative of the Kelvin function of. + /// The derivative of the Kelvin function kei. + public static double KelvinKeiPrime(double nu, double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + + const double inv2Sqrt2 = 0.35355339059327376220042218105242451964241796884424; // 1/(2 * sqrt(2)) + return inv2Sqrt2 * (KelvinKer(nu - 1, x) - KelvinKer(nu + 1, x) - KelvinKei(nu - 1, x) + KelvinKei(nu + 1, x)); + } + + /// + /// Returns the derivative of the Kelvin function kei. + /// + /// The value to compute the derivative of the Kelvin function of. + /// The derivative of the Kelvin function kei. + public static double KelvinKeiPrime(double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + + return KelvinKeiPrime(0, x); + } + } +}