2016-09-01 5 views
1

J'utilise gnu bibliothèque scientifique (gsl) pour définir cylindre parabolique U fonction. Here est la définition de cette fonction U. Mais il y a une petite différence entre ma définition et celle de Mathworld, parce que je veux utiliser la mienne pour utiliser la méthode d'intégration de Gauss-Hermite pour calculer l'intégration de l'intervalle infini. Donc je rejette la partie exponentielle. Voici mon code:pourquoi cette fonction se trompe résultat à une plus grande valeur variable

#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; 
} 

J'utilise scipy.special.pbdv d'examiner si ma définition est correcte. Pour la plus petite valeur x, le résultat est d'accord avec le pbdv, mais lorsque la valeur augmente, le résultat est plus étrange. Comment puis-je résoudre ce problème. Il y a une partie de la production:

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 

la première colonne est la valeur « x », la seconde colonne est la valeur de fonction parabolique, il arrive de se tromper à environ 5.83999999, parce que si nous saisissons cette ParablicCylinderU(-0.999999-1./2.,sqrt(2)*x), cela est presque un ligne droite. Ce qui suit est mon code python d'examen.

#!/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") 
+1

[Ce] (https://msdn.microsoft.com /en-us/library/c151dt3s.aspx) peut aider. –

Répondre

1

Vous pouvez obtenir pbdv code source à partir specfun scipy code Rechercher SUBROUTINE PBDV, si vous comprenez pas Fortran, vous pouvez le convertir en C via f2c