为什么这个函数在一些更大的变量值上得到错误的结果

Question

我正在使用 gnu 科学库 (gsl) 来定义 抛物柱面 U 函数。 Here就是这个U函数的定义。但是我的定义和 Mathworld 的定义有点不同，因为我想用我的定义来使用 Gauss-Hermite 积分方法来计算无穷大区间积分。所以我忽略了指数部分。这是我的代码：

#include <gsl/gsl_errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf.h>
#include <math.h>
#include <string.h>

double ParabolicCylinderU(double a, double x)
{
   double res;
   double tmp1,tmp2;
   tmp1=sqrt(M_PI)/(gsl_sf_gamma(3./4.+a/2.)*pow(2,a/2.+1./4.))*gsl_sf_hyperg_1F1(a/2.+1./4.,1./2.,pow(x,2)/2.);
   tmp2=sqrt(M_PI)/(gsl_sf_gamma(1./4.+a/2.)*pow(2,a/2.-1./4.))*x*gsl_sf_hyperg_1F1(a/2.+3./4.,3./2.,pow(x,2)/2.);
   res=tmp1-tmp2;
   return res;
}

我用scipy.special.pbdv来检验我的定义是否正确。对于较小的x值，结果与pbdv一致，但当值变大时，结果更奇怪。我怎么解决这个问题。有部分输出：

5.3599999999999994        10.720004342260813
5.4000000000000004        10.801040023779478     
5.4399999999999995        10.879249134730191     
5.4800000000000004        10.961176175487582     
5.5199999999999996        11.036212437780495     
5.5600000000000005        11.115912986845069     
5.5999999999999996        11.189683877125612     
5.6400000000000006        11.289942688341265     
5.6799999999999997        11.385711218186625     
5.7200000000000006        11.379394963160701     
5.7599999999999998        11.532254417763568     
5.8000000000000007        11.575985086141165     
5.8399999999999999        11.881533311150061     
5.8800000000000008        11.896642911301599     
5.9199999999999999        11.639708082791794     
5.9600000000000009        11.550981603033073     
6.0000000000000000        10.455916671990396     
6.0399999999999991        9.3887694230651952 
...
9.6799999999999997        5.1646711481042656E+024
9.7199999999999989       -1.5183962001815768E+025
9.7600000000000016       -3.5383625277571119E+025
9.8000000000000007        4.2306292738245936E+025
9.8399999999999999       -1.3554993237535422E+026
9.8799999999999990       -3.9599339761206319E+026
9.9200000000000017        4.4052629528738374E+026
9.9600000000000009        3.7902904410228328E+027
10.000000000000000       -1.6696086530684606E+021

第一列是'x'值，第二列是抛物线函数值，刚好在5.83999999左右是错误的，因为如果我们输入这个ParablicCylinderU(-0.999999-1./2.,sqrt(2)*x)，这就差不多了一条直线。以下是我的考试python代码。

#!/usr/bin/env python                                                                                                                                                                          

from math import *
import numpy as np
from scipy import special
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt

nz= 0.999998779431

x=np.linspace(0,10,250,endpoint=True)
d,dv=special.pbdv(nz,sqrt(2.)*(x))

plt.figure()
d=d*np.exp(np.power(x,2)/2)*np.power(2,nz/2.)
plt.plot(x,d)
plt.plot(x,np.zeros(x.shape[0]))
plt.savefig("scipy_parabolic.png")

Answer 1

您可以从 specfun scipy code 获得 pbdv 源代码寻找SUBROUTINE PBDV，如果你不懂fortran，你可以通过f2c

将其转换为C

为什么这个函数在一些更大的变量值上得到错误的结果

why this function get wrong result at some larger variable value

c

python

gsl