home | career | drupal | java | mac | mysql | perl | scala | uml | unix  

Java example source code file (BesselJ.java)

This example Java source code file (BesselJ.java) is included in the alvinalexander.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" TM.

Learn more about this Java project at its project page.

Java - Java tags/keywords

besseljresult, convergenceexception, enmten, ensig, enten, fact, mathillegalargumentexception, rtnsig, towpi1, twopi2, univariatefunction, x_max, x_min

The BesselJ.java Java example source code

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math3.special;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;

/**
 * This class provides computation methods related to Bessel
 * functions of the first kind. Detailed descriptions of these functions are
 * available in <a
 * href="http://en.wikipedia.org/wiki/Bessel_function">Wikipedia</a>, DLMF (Ch. 10).
 * <p>
 * This implementation is based on the rjbesl Fortran routine at
 * <a href="http://www.netlib.org/specfun/rjbesl">Netlib.

* <p> * From the Fortran code: </p> * <p> * This program is based on a program written by David J. Sookne (2) that * computes values of the Bessel functions J or I of real argument and integer * order. Modifications include the restriction of the computation to the J * Bessel function of non-negative real argument, the extension of the * computation to arbitrary positive order, and the elimination of most * underflow.</p> * <p> * References:</p> * <ul> * <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne, * D. J., Math. Comp. 26, 1972, pp 941-947.</li> * <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS * Jour. of Res. B. 77B, 1973, pp 125-132.</li> * </ul>

* @since 3.4 */ public class BesselJ implements UnivariateFunction { // --------------------------------------------------------------------- // Mathematical constants // --------------------------------------------------------------------- /** -2 / pi */ private static final double PI2 = 0.636619772367581343075535; /** first few significant digits of 2pi */ private static final double TOWPI1 = 6.28125; /** 2pi - TWOPI1 to working precision */ private static final double TWOPI2 = 1.935307179586476925286767e-3; /** TOWPI1 + TWOPI2 */ private static final double TWOPI = TOWPI1 + TWOPI2; // --------------------------------------------------------------------- // Machine-dependent parameters // --------------------------------------------------------------------- /** * 10.0^K, where K is the largest integer such that ENTEN is * machine-representable in working precision */ private static final double ENTEN = 1.0e308; /** * Decimal significance desired. Should be set to (INT(log_{10}(2) * (it)+1)). * Setting NSIG lower will result in decreased accuracy while setting * NSIG higher will increase CPU time without increasing accuracy. * The truncation error is limited to a relative error of * T=.5(10^(-NSIG)). */ private static final double ENSIG = 1.0e16; /** 10.0 ** (-K) for the smallest integer K such that K >= NSIG/4 */ private static final double RTNSIG = 1.0e-4; /** Smallest ABS(X) such that X/4 does not underflow */ private static final double ENMTEN = 8.90e-308; /** Minimum acceptable value for x */ private static final double X_MIN = 0.0; /** * Upper limit on the magnitude of x. If abs(x) = n, then at least * n iterations of the backward recursion will be executed. The value of * 10.0 ** 4 is used on every machine. */ private static final double X_MAX = 1.0e4; /** First 25 factorials as doubles */ private static final double[] FACT = { 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, 1.307674368e12, 2.0922789888e13, 3.55687428096e14, 6.402373705728e15, 1.21645100408832e17, 2.43290200817664e18, 5.109094217170944e19, 1.12400072777760768e21, 2.585201673888497664e22, 6.2044840173323943936e23 }; /** Order of the function computed when {@link #value(double)} is used */ private final double order; /** * Create a new BesselJ with the given order. * * @param order order of the function computed when using {@link #value(double)}. */ public BesselJ(double order) { this.order = order; } /** * Returns the value of the constructed Bessel function of the first kind, * for the passed argument. * * @param x Argument * @return Value of the Bessel function at x * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order} * @throws ConvergenceException if the algorithm fails to converge */ public double value(double x) throws MathIllegalArgumentException, ConvergenceException { return BesselJ.value(order, x); } /** * Returns the first Bessel function, \(J_{order}(x)\). * * @param order Order of the Bessel function * @param x Argument * @return Value of the Bessel function of the first kind, \(J_{order}(x)\) * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order} * @throws ConvergenceException if the algorithm fails to converge */ public static double value(double order, double x) throws MathIllegalArgumentException, ConvergenceException { final int n = (int) order; final double alpha = order - n; final int nb = n + 1; final BesselJResult res = rjBesl(x, alpha, nb); if (res.nVals >= nb) { return res.vals[n]; } else if (res.nVals < 0) { throw new MathIllegalArgumentException(LocalizedFormats.BESSEL_FUNCTION_BAD_ARGUMENT,order, x); } else if (FastMath.abs(res.vals[res.nVals - 1]) < 1e-100) { return res.vals[n]; // underflow; return value (will be zero) } throw new ConvergenceException(LocalizedFormats.BESSEL_FUNCTION_FAILED_CONVERGENCE, order, x); } /** * Encapsulates the results returned by {@link BesselJ#rjBesl(double, double, int)}. * <p> * {@link #getVals()} returns the computed function values. * {@link #getnVals()} is the number of values among those returned by {@link #getnVals()} * that can be considered accurate. * </p>

* <ul> * <li>nVals < 0: An argument is out of range. For example, nb <= 0, alpha * < 0 or > 1, or x is too large. In this case, b(0) is set to zero, the * remainder of the b-vector is not calculated, and nVals is set to * MIN(nb,0) - 1 so that nVals != nb.</li> * <li>nb > nVals > 0: Not all requested function values could be calculated * accurately. This usually occurs because nb is much larger than abs(x). In * this case, b(n) is calculated to the desired accuracy for n < nVals, but * precision is lost for nVals < n <= nb. If b(n) does not vanish for n > * nVals (because it is too small to be represented), and b(n)/b(nVals) = * \(10^{-k}\), then only the first NSIG-k significant figures of b(n) can be * trusted.</li>

*/ public static class BesselJResult { /** Bessel function values */ private final double[] vals; /** Valid value count */ private final int nVals; /** * Create a new BesselJResult with the given values and valid value count. * * @param b values * @param n count of valid values */ public BesselJResult(double[] b, int n) { vals = MathArrays.copyOf(b, b.length); nVals = n; } /** * @return the computed function values */ public double[] getVals() { return MathArrays.copyOf(vals, vals.length); } /** * @return the number of valid function values (normally the same as the * length of the array returned by {@link #getnVals()}) */ public int getnVals() { return nVals; } } /** * Calculates Bessel functions \(J_{n+alpha}(x)\) for * non-negative argument x, and non-negative order n + alpha. * <p> * Before using the output vector, the user should check that * nVals = nb, i.e., all orders have been calculated to the desired accuracy. * See BesselResult class javadoc for details on return values. * </p> * @param x non-negative real argument for which J's are to be calculated * @param alpha fractional part of order for which J's or exponentially * scaled J's (\(J\cdot e^{x}\)) are to be calculated. 0 <= alpha < 1.0. * @param nb integer number of functions to be calculated, nb > 0. The first * function calculated is of order alpha, and the last is of order * nb - 1 + alpha. * @return BesselJResult a vector of the functions * \(J_{alpha}(x)\) through \(J_{nb-1+alpha}(x)\), or the corresponding exponentially * scaled functions and an integer output variable indicating possible errors */ public static BesselJResult rjBesl(double x, double alpha, int nb) { final double[] b = new double[nb]; int ncalc = 0; double alpem = 0; double alp2em = 0; // --------------------------------------------------------------------- // Check for out of range arguments. // --------------------------------------------------------------------- final int magx = (int) x; if ((nb > 0) && (x >= X_MIN) && (x <= X_MAX) && (alpha >= 0) && (alpha < 1)) { // --------------------------------------------------------------------- // Initialize result array to zero. // --------------------------------------------------------------------- ncalc = nb; for (int i = 0; i < nb; ++i) { b[i] = 0; } // --------------------------------------------------------------------- // Branch to use 2-term ascending series for small X and asymptotic // form for large X when NB is not too large. // --------------------------------------------------------------------- double tempa; double tempb; double tempc; double tover; if (x < RTNSIG) { // --------------------------------------------------------------------- // Two-term ascending series for small X. // --------------------------------------------------------------------- tempa = 1; alpem = 1 + alpha; double halfx = 0; if (x > ENMTEN) { halfx = 0.5 * x; } if (alpha != 0) { tempa = FastMath.pow(halfx, alpha) / (alpha * Gamma.gamma(alpha)); } tempb = 0; if (x + 1 > 1) { tempb = -halfx * halfx; } b[0] = tempa + (tempa * tempb / alpem); if ((x != 0) && (b[0] == 0)) { ncalc = 0; } if (nb != 1) { if (x <= 0) { for (int n = 1; n < nb; ++n) { b[n] = 0; } } else { // --------------------------------------------------------------------- // Calculate higher order functions. // --------------------------------------------------------------------- tempc = halfx; tover = tempb != 0 ? ENMTEN / tempb : 2 * ENMTEN / x; for (int n = 1; n < nb; ++n) { tempa /= alpem; alpem += 1; tempa *= tempc; if (tempa <= tover * alpem) { tempa = 0; } b[n] = tempa + (tempa * tempb / alpem); if ((b[n] == 0) && (ncalc > n)) { ncalc = n; } } } } } else if ((x > 25.0) && (nb <= magx + 1)) { // --------------------------------------------------------------------- // Asymptotic series for X > 25 // --------------------------------------------------------------------- final double xc = FastMath.sqrt(PI2 / x); final double mul = 0.125 / x; final double xin = mul * mul; int m = 0; if (x >= 130.0) { m = 4; } else if (x >= 35.0) { m = 8; } else { m = 11; } final double xm = 4.0 * m; // --------------------------------------------------------------------- // Argument reduction for SIN and COS routines. // --------------------------------------------------------------------- double t = (double) ((int) ((x / TWOPI) + 0.5)); final double z = x - t * TOWPI1 - t * TWOPI2 - (alpha + 0.5) / PI2; double vsin = FastMath.sin(z); double vcos = FastMath.cos(z); double gnu = 2 * alpha; double capq; double capp; double s; double t1; double xk; for (int i = 1; i <= 2; i++) { s = (xm - 1 - gnu) * (xm - 1 + gnu) * xin * 0.5; t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0)); capp = (s * t) / FACT[2 * m]; t1 = (gnu - (xm + 1)) * (gnu + (xm + 1)); capq = (s * t1) / FACT[2 * m + 1]; xk = xm; int k = 2 * m; t1 = t; for (int j = 2; j <= m; j++) { xk -= 4.0; s = (xk - 1 - gnu) * (xk - 1 + gnu); t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0)); capp = (capp + 1 / FACT[k - 2]) * s * t * xin; capq = (capq + 1 / FACT[k - 1]) * s * t1 * xin; k -= 2; t1 = t; } capp += 1; capq = (capq + 1) * ((gnu * gnu) - 1) * (0.125 / x); b[i - 1] = xc * (capp * vcos - capq * vsin); if (nb == 1) { return new BesselJResult(MathArrays.copyOf(b, b.length), ncalc); } t = vsin; vsin = -vcos; vcos = t; gnu += 2.0; } // --------------------------------------------------------------------- // If NB > 2, compute J(X,ORDER+I) I = 2, NB-1 // --------------------------------------------------------------------- if (nb > 2) { gnu = 2 * alpha + 2.0; for (int j = 2; j < nb; ++j) { b[j] = gnu * b[j - 1] / x - b[j - 2]; gnu += 2.0; } } } else { // --------------------------------------------------------------------- // Use recurrence to generate results. First initialize the // calculation of P*S. // --------------------------------------------------------------------- final int nbmx = nb - magx; int n = magx + 1; int nstart = 0; int nend = 0; double en = 2 * (n + alpha); double plast = 1; double p = en / x; double pold; // --------------------------------------------------------------------- // Calculate general significance test. // --------------------------------------------------------------------- double test = 2 * ENSIG; boolean readyToInitialize = false; if (nbmx >= 3) { // --------------------------------------------------------------------- // Calculate P*S until N = NB-1. Check for possible // overflow. // --------------------------------------------------------------------- tover = ENTEN / ENSIG; nstart = magx + 2; nend = nb - 1; en = 2 * (nstart - 1 + alpha); double psave; double psavel; for (int k = nstart; k <= nend; k++) { n = k; en += 2.0; pold = plast; plast = p; p = (en * plast / x) - pold; if (p > tover) { // --------------------------------------------------------------------- // To avoid overflow, divide P*S by TOVER. Calculate // P*S until // ABS(P) > 1. // --------------------------------------------------------------------- tover = ENTEN; p /= tover; plast /= tover; psave = p; psavel = plast; nstart = n + 1; do { n += 1; en += 2.0; pold = plast; plast = p; p = (en * plast / x) - pold; } while (p <= 1); tempb = en / x; // --------------------------------------------------------------------- // Calculate backward test and find NCALC, the // highest N such that // the test is passed. // --------------------------------------------------------------------- test = pold * plast * (0.5 - 0.5 / (tempb * tempb)); test /= ENSIG; p = plast * tover; n -= 1; en -= 2.0; nend = FastMath.min(nb, n); for (int l = nstart; l <= nend; l++) { pold = psavel; psavel = psave; psave = (en * psavel / x) - pold; if (psave * psavel > test) { ncalc = l - 1; readyToInitialize = true; break; } } ncalc = nend; readyToInitialize = true; break; } } if (!readyToInitialize) { n = nend; en = 2 * (n + alpha); // --------------------------------------------------------------------- // Calculate special significance test for NBMX > 2. // --------------------------------------------------------------------- test = FastMath.max(test, FastMath.sqrt(plast * ENSIG) * FastMath.sqrt(2 * p)); } } // --------------------------------------------------------------------- // Calculate P*S until significance test passes. // --------------------------------------------------------------------- if (!readyToInitialize) { do { n += 1; en += 2.0; pold = plast; plast = p; p = (en * plast / x) - pold; } while (p < test); } // --------------------------------------------------------------------- // Initialize the backward recursion and the normalization sum. // --------------------------------------------------------------------- n += 1; en += 2.0; tempb = 0; tempa = 1 / p; int m = (2 * n) - 4 * (n / 2); double sum = 0; double em = (double) (n / 2); alpem = em - 1 + alpha; alp2em = 2 * em + alpha; if (m != 0) { sum = tempa * alpem * alp2em / em; } nend = n - nb; boolean readyToNormalize = false; boolean calculatedB0 = false; // --------------------------------------------------------------------- // Recur backward via difference equation, calculating (but not // storing) B(N), until N = NB. // --------------------------------------------------------------------- for (int l = 1; l <= nend; l++) { n -= 1; en -= 2.0; tempc = tempb; tempb = tempa; tempa = (en * tempb / x) - tempc; m = 2 - m; if (m != 0) { em -= 1; alp2em = 2 * em + alpha; if (n == 1) { break; } alpem = em - 1 + alpha; if (alpem == 0) { alpem = 1; } sum = (sum + tempa * alp2em) * alpem / em; } } // --------------------------------------------------------------------- // Store B(NB). // --------------------------------------------------------------------- b[n - 1] = tempa; if (nend >= 0) { if (nb <= 1) { alp2em = alpha; if (alpha + 1 == 1) { alp2em = 1; } sum += b[0] * alp2em; readyToNormalize = true; } else { // --------------------------------------------------------------------- // Calculate and store B(NB-1). // --------------------------------------------------------------------- n -= 1; en -= 2.0; b[n - 1] = (en * tempa / x) - tempb; if (n == 1) { calculatedB0 = true; } else { m = 2 - m; if (m != 0) { em -= 1; alp2em = 2 * em + alpha; alpem = em - 1 + alpha; if (alpem == 0) { alpem = 1; } sum = (sum + (b[n - 1] * alp2em)) * alpem / em; } } } } if (!readyToNormalize && !calculatedB0) { nend = n - 2; if (nend != 0) { // --------------------------------------------------------------------- // Calculate via difference equation and store B(N), // until N = 2. // --------------------------------------------------------------------- for (int l = 1; l <= nend; l++) { n -= 1; en -= 2.0; b[n - 1] = (en * b[n] / x) - b[n + 1]; m = 2 - m; if (m != 0) { em -= 1; alp2em = 2 * em + alpha; alpem = em - 1 + alpha; if (alpem == 0) { alpem = 1; } sum = (sum + b[n - 1] * alp2em) * alpem / em; } } } } // --------------------------------------------------------------------- // Calculate b[0] // --------------------------------------------------------------------- if (!readyToNormalize) { if (!calculatedB0) { b[0] = 2.0 * (alpha + 1) * b[1] / x - b[2]; } em -= 1; alp2em = 2 * em + alpha; if (alp2em == 0) { alp2em = 1; } sum += b[0] * alp2em; } // --------------------------------------------------------------------- // Normalize. Divide all B(N) by sum. // --------------------------------------------------------------------- if (FastMath.abs(alpha) > 1e-16) { sum *= Gamma.gamma(alpha) * FastMath.pow(x * 0.5, -alpha); } tempa = ENMTEN; if (sum > 1) { tempa *= sum; } for (n = 0; n < nb; n++) { if (FastMath.abs(b[n]) < tempa) { b[n] = 0; } b[n] /= sum; } } // --------------------------------------------------------------------- // Error return -- X, NB, or ALPHA is out of range. // --------------------------------------------------------------------- } else { if (b.length > 0) { b[0] = 0; } ncalc = FastMath.min(nb, 0) - 1; } return new BesselJResult(MathArrays.copyOf(b, b.length), ncalc); } }

Other Java examples (source code examples)

Here is a short list of links related to this Java BesselJ.java source code file:



my book on functional programming

 

new blog posts

 

Copyright 1998-2019 Alvin Alexander, alvinalexander.com
All Rights Reserved.

A percentage of advertising revenue from
pages under the /java/jwarehouse URI on this website is
paid back to open source projects.