Christoph Ruegg
8 years ago
9 changed files with 859 additions and 221 deletions
@ -0,0 +1,135 @@ |
|||
using NUnit.Framework; |
|||
using System; |
|||
|
|||
namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests |
|||
{ |
|||
/// <summary>
|
|||
/// Kelvin functions tests.
|
|||
/// </summary>
|
|||
[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); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,264 @@ |
|||
using System; |
|||
using System.Numerics; |
|||
|
|||
namespace MathNet.Numerics |
|||
{ |
|||
/// <summary>
|
|||
/// This partial implementation of the SpecialFunctions class contains all methods related to the modified Bessel function.
|
|||
/// </summary>
|
|||
public static partial class SpecialFunctions |
|||
{ |
|||
/// <summary>
|
|||
/// Returns the Kelvin function of the first kind.
|
|||
/// <para>KelvinBe(nu, x) is given by BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// <para>KelvinBer(nu, x) and KelvinBei(nu, x) are the real and imaginary parts of the KelvinBe(nu, x)</para>
|
|||
/// </summary>
|
|||
/// <param name="nu">the order of the the Kelvin function.</param>
|
|||
/// <param name="x">The value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function of the first kind.</returns>
|
|||
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); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function ber.
|
|||
/// <para>KelvinBer(nu, x) is given by the real part of BesselJ(nu, j * sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// </summary>
|
|||
/// <param name="nu">the order of the the Kelvin function.</param>
|
|||
/// <param name="x">The value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function ber.</returns>
|
|||
public static double KelvinBer(double nu, double x) |
|||
{ |
|||
return KelvinBe(nu, x).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function ber.
|
|||
/// <para>KelvinBer(x) is given by the real part of BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// <para>KelvinBer(x) is equivalent to KelvinBer(0, x).</para>
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function ber.</returns>
|
|||
public static double KelvinBer(double x) |
|||
{ |
|||
return KelvinBe(0, x).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function bei.
|
|||
/// <para>KelvinBei(nu, x) is given by the imaginary part of BesselJ(nu, j * sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// </summary>
|
|||
/// <param name="nu">the order of the the Kelvin function.</param>
|
|||
/// <param name="x">The value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function bei.</returns>
|
|||
public static double KelvinBei(double nu, double x) |
|||
{ |
|||
return KelvinBe(nu, x).Imaginary; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function bei.
|
|||
/// <para>KelvinBei(x) is given by the imaginary part of BesselJ(0, j * sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// <para>KelvinBei(x) is equivalent to KelvinBei(0, x).</para>
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function bei.</returns>
|
|||
public static double KelvinBei(double x) |
|||
{ |
|||
return KelvinBe(0, x).Imaginary; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the derivative of the Kelvin function ber.
|
|||
/// </summary>
|
|||
/// <param name="nu">The order of the Kelvin function.</param>
|
|||
/// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
|
|||
/// <returns>the derivative of the Kelvin function ber</returns>
|
|||
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)); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the derivative of the Kelvin function ber.
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
|
|||
/// <returns>The derivative of the Kelvin function ber.</returns>
|
|||
public static double KelvinBerPrime(double x) |
|||
{ |
|||
return KelvinBerPrime(0, x); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the derivative of the Kelvin function bei.
|
|||
/// </summary>
|
|||
/// <param name="nu">The order of the Kelvin function.</param>
|
|||
/// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
|
|||
/// <returns>the derivative of the Kelvin function bei.</returns>
|
|||
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)); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the derivative of the Kelvin function bei.
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
|
|||
/// <returns>The derivative of the Kelvin function bei.</returns>
|
|||
public static double KelvinBeiPrime(double x) |
|||
{ |
|||
return KelvinBeiPrime(0, x); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function of the second kind
|
|||
/// <para>KelvinKe(nu, x) is given by Exp(-nu * pi * j / 2) * BesselK(nu, x * sqrt(j)) where j = sqrt(-1).</para>
|
|||
/// <para>KelvinKer(nu, x) and KelvinKei(nu, x) are the real and imaginary parts of the KelvinBe(nu, x)</para>
|
|||
/// </summary>
|
|||
/// <param name="nu">The order of the Kelvin function.</param>
|
|||
/// <param name="x">The value to calculate the kelvin function of,</param>
|
|||
/// <returns></returns>
|
|||
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); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function ker.
|
|||
/// <para>KelvinKer(nu, x) is given by the real part of Exp(-nu * pi * j / 2) * BesselK(nu, sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// </summary>
|
|||
/// <param name="nu">the order of the the Kelvin function.</param>
|
|||
/// <param name="x">The non-negative real value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function ker.</returns>
|
|||
public static double KelvinKer(double nu, double x) |
|||
{ |
|||
if (x <= 0.0) |
|||
{ |
|||
throw new ArithmeticException(); |
|||
} |
|||
|
|||
return KelvinKe(nu, x).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function ker.
|
|||
/// <para>KelvinKer(x) is given by the real part of Exp(-nu * pi * j / 2) * BesselK(0, sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// <para>KelvinKer(x) is equivalent to KelvinKer(0, x).</para>
|
|||
/// </summary>
|
|||
/// <param name="x">The non-negative real value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function ker.</returns>
|
|||
public static double KelvinKer(double x) |
|||
{ |
|||
if (x <= 0.0) |
|||
{ |
|||
throw new ArithmeticException(); |
|||
} |
|||
|
|||
return KelvinKe(0, x).Real; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function kei.
|
|||
/// <para>KelvinKei(nu, x) is given by the imaginary part of Exp(-nu * pi * j / 2) * BesselK(nu, sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// </summary>
|
|||
/// <param name="nu">the order of the the Kelvin function.</param>
|
|||
/// <param name="x">The non-negative real value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function kei.</returns>
|
|||
public static double KelvinKei(double nu, double x) |
|||
{ |
|||
if (x <= 0.0) |
|||
{ |
|||
throw new ArithmeticException(); |
|||
} |
|||
|
|||
return KelvinKe(nu, x).Imaginary; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the Kelvin function kei.
|
|||
/// <para>KelvinKei(x) is given by the imaginary part of Exp(-nu * pi * j / 2) * BesselK(0, sqrt(j) * x) where j = sqrt(-1).</para>
|
|||
/// <para>KelvinKei(x) is equivalent to KelvinKei(0, x).</para>
|
|||
/// </summary>
|
|||
/// <param name="x">The non-negative real value to compute the Kelvin function of.</param>
|
|||
/// <returns>The Kelvin function kei.</returns>
|
|||
public static double KelvinKei(double x) |
|||
{ |
|||
if (x <= 0.0) |
|||
{ |
|||
throw new ArithmeticException(); |
|||
} |
|||
|
|||
return KelvinKe(0, x).Imaginary; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the derivative of the Kelvin function ker.
|
|||
/// </summary>
|
|||
/// <param name="nu">The order of the Kelvin function.</param>
|
|||
/// <param name="x">The non-negative real value to compute the derivative of the Kelvin function of.</param>
|
|||
/// <returns>The derivative of the Kelvin function ker.</returns>
|
|||
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)); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the derivative of the Kelvin function ker.
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
|
|||
/// <returns>The derivative of the Kelvin function ker.</returns>
|
|||
public static double KelvinKerPrime(double x) |
|||
{ |
|||
if (x <= 0.0) |
|||
{ |
|||
throw new ArithmeticException(); |
|||
} |
|||
|
|||
return KelvinKerPrime(0, x); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the derivative of the Kelvin function kei.
|
|||
/// </summary>
|
|||
/// <param name="nu">The order of the Kelvin function.</param>
|
|||
/// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
|
|||
/// <returns>The derivative of the Kelvin function kei.</returns>
|
|||
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)); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the derivative of the Kelvin function kei.
|
|||
/// </summary>
|
|||
/// <param name="x">The value to compute the derivative of the Kelvin function of.</param>
|
|||
/// <returns>The derivative of the Kelvin function kei.</returns>
|
|||
public static double KelvinKeiPrime(double x) |
|||
{ |
|||
if (x <= 0.0) |
|||
{ |
|||
throw new ArithmeticException(); |
|||
} |
|||
|
|||
return KelvinKeiPrime(0, x); |
|||
} |
|||
} |
|||
} |
|||
@ -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 |
|||
{ |
|||
/// <summary>
|
|||
/// For Bessel-related functions, no scaling factor is applied.
|
|||
/// </summary>
|
|||
Unity = 0, |
|||
|
|||
/// <summary>
|
|||
/// For Bessel-related functions, exponential scaling is applied.
|
|||
/// </summary>
|
|||
Exponential = 1 |
|||
} |
|||
} |
|||
} |
|||
Loading…
Reference in new issue