我正在尝试使用 Java 计算二阶贝塞尔函数 (Y_v(x))。其中,贝塞尔函数的阶数 v 是非整数。
我知道,Apache 有一个计算一阶贝塞尔函数的函数(J_v(x))
import org.apache.commons.math3.special.BesselJ;
BesselJ bj = new BesselJ(10.5);
double bjr = bj.value(r);
它可以计算非整数一阶贝塞尔函数(但不能计算负阶函数)。我不相信它有二阶的等价物。
我找到了一个计算二阶贝塞尔函数的公式(其他地方的数值计算公式 6.4.2),它从一阶函数计算二阶贝塞尔函数。
然而,这个方程需要一阶贝塞尔函数来接受 apache BesselJ 函数不接受的负 v 值。
我编写了一个程序来计算 BesselJ 函数,该函数使用公式(数字食谱 6.4.1)接受负值
import org.apache.commons.math3.special.Gamma;
import org.apache.commons.math3.util.ArithmeticUtils;
private static int lmax;
private static double[][] Gammas;
private static int besselcount;
public static Double Bessel_J(double m, double x){
double BJ = 0.0;
for (int k = 0; k < besselcount;k++){
double kk = (double) k;
double llmm = (double) lmax;
double A = Math.pow(x/2.,(2.*kk)+m);
double B = Gammas[(int) (m+llmm+0.5)][k];
if (B==0.0||B==-0.0){break;}
BJ = BJ + A*B;
}
return BJ;
}
public static void main(String[] args) throws Exception {
besselcount = 80;
lmax = 20;
Gammas = new double[2*lmax+2][besselcount+1];
for (int n = 0;n <= ((2*lmax)+1);n++){
for (int m = 0;m <=besselcount;m++) {
double nn = (double) n;
double mm = (double) m;
double llmm = (double) lmax;
Gammas[n][m] = Math.pow(-1,m)/(ArithmeticUtils.factorialDouble(m)*Gamma.gamma(1.+nn-llmm-0.5+mm));
}
}
for (double x = 0.02; x < 50.0 ; x = x + 0.02){
double bj = Bessel_J(-10.5, x);
System.out.println(x+" "+bj);
}
}
(将伽马预先计算并存储在数组中只是因为在其余代码中如何使用贝塞尔函数,以避免不必要地重新计算)
我编写的函数适用于低 x 值,但随后会失去稳定性。对于我在x = 30左右之后使用的v = -10.5值。如图所示。
在 wolfram alpha 网站上,mathematica 使用相同的方程计算贝塞尔函数,mathematica 的 BesselJ 函数可以计算
Show[Plot[BesselJ[-10.5, x], {x, 0, 50}], ImageSize -> Full]
不失稳定性。
任何关于如何纠正我的函数或使用 Java 计算二阶贝塞尔函数的替代方法的建议都将非常感激。
请说是否有任何澄清会有用。
谢谢