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

Question

我使用gnu科学库（gsl）来定义抛物线柱U函数。Here是这个U函数的定义。但是我的定义和数学世界有一点不同，因为我想用我的定义来计算无穷区间积分。所以我忽略了指数部分。这是我的密码：

#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