Re: Funzione ipergeometrica di Gauss in C
Inviato: lunedì 23 aprile 2012, 17:47
Vedi allegato.
Il forum della comunità italiana di Ubuntu.
https://forum.ubuntu-it.org/
La funzione ipergeometrica è definita come la soluzione dell'equazione differenziale ipergeometrica di Eulero, e non è stupefacente che un'equazione differenziale del secondo ordine a coefficienti reali abbia una soluzione esprimibile in forma complessa: pensa ad esempio all'oscillatore armonico. La serie ipergeometrica risolve l'equazione ipergeometrica ma ha raggio di convergenza 1 (non 0.5 come erroneamente avevo scritto). Per z reale compreso tra 0 e 1 il valore di 2F1 dovrebbe sempre essere reale.giuseppe morello ha scritto:Di nuovo grazie, ottimo lavoro.
Ho fatto un po' di prove confrontando con i risultati di wolphram alpha: i risultati sono coincidenti, variando al più l'ultima cifra decimale.
Due osservazioni:
1, di carattere matematico: mi aspetterei che se a, b, c, z sono tutti reali, la parte immaginaria della funzione deve essere nulla, ma non è così nemmeno per wolphram (per esempio, plot per a=1, b=2, c=4, z tra 0 e 10)
Vero, ma non so come rimediare. hmin viene passato come argomento alla odeint, ma cambiarlo non risolve.
2, sull'algoritmo: per a=-1, b=5, c=3, z=4 e anche in alcuni altri casi, il codice non converge, ritenendo che hnext risulti troppo piccolo. Per gli stessi valori wolphram non mostra problemi. In effetti, l'ho fatto stampare e risulta nullo almeno fino alla 15^a cifra, mentre hmin=0.0001 (anche se non ho ancora trovato dove viene assegnato questo valore).
Codice: Seleziona tutto
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<gsl/gsl_sf_hyperg.h>
double a, b, c, x, iper;
int main()
{
printf("Inserisci i valori a, b, c, x:\n");
scanf("%lf %lf %lf %lf", &a, &b, &c, &x);
iper=hyperg_2F1_series(a,b,c,x);
printf("%18.15f\n", iper);
return 0;
}
Codice: Seleziona tutto
gcc -c -lm ipergeom.c -o ipergeom.out
Codice: Seleziona tutto
pkg-config --libs gsl_sf_hyperg.h
Codice: Seleziona tutto
gcc -lm -lgsl -lgslcblas -o eseguibile sorgente.c
Non so cosa facesse, ovviamente il con il tuo comando invece funziona. Devo assolutamente studiare un po' o seguire qualche corso, perché a cercare come si compila sul web spesso mi confondo.bite ha scritto:Se si produce un eseguibile significa che gli header sono stati trovati, e anche le librerie.
Mi sembra strano, però, che abbia compilato con il comando che hai riportato. Io ho dovuto dare:
Codice: Seleziona tutto
gcc -lm -lgsl -lgslcblas -o eseguibile sorgente.c
Svista.bite ha scritto: Da dove ti salta fuori hyperg_2F1_series? I nomi di tutte le funzioni della libreria gsl cominciano con gsl_
Gli esempi che facevo erano del tutto casuali, non lavoro solo con interi. Potrebbe anche darsi che il codice che mi hai inviato mi basti (sulla base delle poche prove fatte). Per verificarlo devo prima terminare il programma, per il quale mi mancano ancora una funzione abbastanza facile da computare e in più l'ipergeometrica F1 di Appell, per la quale ancora non ho trovato moltissimo (non mi risulta né su GSL né su Numerical Recipes).bite ha scritto: Io ho usato gsl_sf_hyperg_2F1 (che però converge solo per |x|<1).
... e ahimé, a ben guardare il manuale _tutte_ le funzioni ipergeometriche della gsl hanno questa limitazione.
Comunque casi come quello che hai citato nel tuo messaggio precedente, con a o b interi negativi, li puoi gestire facilmente come casi speciale, perché i simboli di Pochhammer per interi negativi diventano 0 dopo un certo indice e quindi la serie ipergeometrica risulta avere un numero finito di termini.
Temo sia impossibile. odeint sta per Ordinary DIfferential Equation INTegrator, mentre le equazioni differenziali costitutive delle funzioni di Appell sono alle derivate parziali.giuseppe morello ha scritto:Lo sviluppo in serie e' abbastanza facile, l'ho fatto sulla falsariga di quello di 2F1. Ora devo sperare di riuscire ad adattare qualcosa tipo odeint
Il secondo paper che ti ho linkato parla di "higher order radiative corrections to scattering amplitudes"; per questa applicazione sembra che sia sufficiente uno sviluppo della F1 in un parametro piccolo attorno ad un valore noto. Se pensi che possa andar bene anche per te, dai uno sguardo ai paragrafi 1 e 4.2.Quel che sto facendo non e' niente di eccezionale: semplicemente sto cercando di tradurre in C e adattare meglio ai miei scopi un codice che simula curve di luce di transiti, solo che l'originale non ricorre a queste funzioni.
Mi sa che questa e' l'idea giusta . Il dominio di convergenza della serie non era sufficiente, quello dell'integrale pare di si' invece.bite ha scritto: In alternativa, visto che la F1 ha anche una rappresentazione in termini di un integrale semplice (nel senso di non doppio ) può darsi che ti convenga seguire questa strada.
Attenzione però che l'integrando può avere delle singolarità a seconda dei valori dei parametri.giuseppe morello ha scritto: Al momento l'ho calcolato semplicemente con il metodo di Simpson, devo solo perfezionarlo per avere una convergenza piu' rapida. Ti aggiornero'.
Codice: Seleziona tutto
#include<math.h>
//#include<complex.h>
#include<stdio.h>
#include"beta.h"
#include"appellf1_integrand.h"
int main()
{
double integ[30], a, b, c, d, x, y, t, f, h, lim;
int i, n;
//eps=pow(10.,-15.);
//printf("eps=%18.15lf\n", eps);
printf("Inserisci i valori di a, b, c, d, x, y:\n");
scanf("%lf %lf %lf %lf %lf %lf", &a, &b, &c, &d, &x, &y);
if(a<=0.)
{
printf("Can't accept a<=0 in appell_integ\n");
}
else if(d-a<=0.)
{
//nerror("Can't accept d<=a in appell_integ\n");
printf("Can't accept d<=a in appell_integ\n");
}
else if(x>1. || y>1.)
{
//nerror("1 or 2 poles along real axis in appell_integ\n");
printf("Can't accept x>1 or y>1 in appell_integ\n");
}
/*else
{
integ=0.;
//if(a==1.) integ+=1.;
for(i=1; i<1000000; i++)
{
t=0.000001*i;
f=appellf1_integrand(t,a,b,c,d,x,y);
if(i%2==1) integ+=4.*f;
else integ+=2.*f;
}
if(a-1.<0.)
{
f=appellf1_integrand(0.000001,a,b,c,d,x,y);
}
else
{
f=appellf1_integrand(0.,a,b,c,d,x,y);
}
integ+=f;
if(x!=1. && y!=1. && d-a-1.>=0.)
{
f=appellf1_integrand(1.,a,b,c,d,x,y);
}
if(x!=1. && y!=1. && d-a-1.<0.)
{
f=appellf1_integrand(0.999999,a,b,c,d,x,y);
}
if(x==1. && y!=1. && d-a-1.-b>=0.)
{
f=appellf1_integrand(1.,a,b,c,d,x,y);
}
if(x==1. && y!=1. && d-a-1.-b<0.)
{
f=appellf1_integrand(0.999999,a,b,c,d,x,y);
}
if(x!=1. && y==1. && d-a-1.-c>=0.)
{
f=appellf1_integrand(1.,a,b,c,d,x,y);
}
if(x!=1. && y==1. && d-a-1.-c<0.)
{
f=appellf1_integrand(0.999999,a,b,c,d,x,y);
}
if(x==1. && y==1. && d-a-1.-b-c>=0)
{
f=appellf1_integrand(1.,a,b,c,d,x,y);
}
if(x==1. && y==1. && d-a-1.-b-c<0)
{
f=appellf1_integrand(0.999999,a,b,c,d,x,y);
}
integ+=f;
integ*=0.000001/3;
integ/=beta(a,d-a);
printf("appell=%18.15lf\n", integ);
}*/
else
{
integ[0]=0.;
for(n=10; n<=27; n++)
{
integ[n]=0.;
h=pow(10.,-n);
for(i=1; i<=10; i++)
{
t=h*i;
f=appellf1_integrand(t,a,b,c,d,x,y);
//printf("f(%lf)=%lf\n", t, f);
if(i==1 || i==10) integ[n]+=f;
else if(i%2==1) integ[n]+=2.*f;
else if(i%2==0) integ[n]+=4.*f;
t=1.-(10.-i)*h;
f=appellf1_integrand(t,a,b,c,d,x,y);
//printf("f(%lf)=%lf\n", t, f);
if(i==1 || i==10) integ[n]+=f;
else if(i%2==1) integ[n]+=2.*f;
else if(i%2==0) integ[n]+=4.*f;
}
integ[n]*=h/3.;
integ[0]+=integ[n];
}
h=pow(10.,-9);
lim=pow(10.,9);
integ[9]=0.;
for(i=1; i<lim; i++)
{
t=h*i;
f=appellf1_integrand(t,a,b,c,d,x,y);
if(i==1 || i==lim-1) integ[9]+=f;
else if(i%2==1) integ[9]+=2.*f;
else if(i%2==0) integ[9]+=4.*f;
}
integ[9]*=h/3.;
integ[0]+=integ[9];
integ[0]/=beta(a,d-a);
printf("appell=%18.15lf\n", integ[0]);
}
return 0;
}
Codice: Seleziona tutto
#include<math.h>
double appellf1_integrand(double t, double a, double b, double c, double d, double x, double y)
{
double f;
if(t==1. && x==1. && y!=1.) f=pow(t,a-1.)*pow(1.-t,d-a-1.-b)*pow(1.-y*t,-c);
if(t==1. && x!=1. && y==1.) f=pow(t,a-1.)*pow(1.-t,d-a-1.-c)*pow(1.-x*t,-b);
if(t==1. && x==1. && y==1.) f=pow(t,a-1.)*pow(1.-t,d-a-1.-b-c);
else f=pow(t,a-1.)*pow(1.-t,d-a-1.)*pow(1.-x*t,-b)*pow(1.-y*t,-c);
return f;
}
Codice: Seleziona tutto
#include<math.h>
double beta(double x, double y)
{
double B;
B=exp(log(gamma(x))+log(gamma(y))-log(gamma(x+y)));
return B;
}
Sì, qualcuno lo ha recentemente convertito dal Fortran in R (pacchetto statistico), con la collaborazione degli autori. Ecco il linkbite ha scritto:... non si capisce se abbiano reso pubblico il loro codice.
No, il tuo link si riferisce alla traduzione (o meglio wrapping) in R del codice di Colavecchia, che è il primo dei due link che ho postato io.BlueEyes ha scritto:Sì, qualcuno lo ha recentemente convertito dal Fortran in R (pacchetto statistico), con la collaborazione degli autori. Ecco il linkbite ha scritto:... non si capisce se abbiano reso pubblico il loro codice.
Ciao
Sì, è vero. C'è però dentro (il link) un file pdf, da dove ho potuto ricavare il seguente codice Fortran, che ritengo possa essere utile allo studente che ha aperto il topic.bite ha scritto:... il tuo link si riferisce alla traduzione (o meglio wrapping) in R del codice di Colavecchia...
Codice: Seleziona tutto
c***********************************************************************
c
c program : complex hypergeometric function
c
c notation : F(a,b;c;z)
c
c usage: The program has been written to be used as a function
c where the user supplies the hypergeometric parameters
c (a,b, and c), the independent variable (z).
c
c warning: This code is under construction. Warning messages
c will appear in regions where the code has not yet
c been completed.
c
c written by:
c
c Robert C. Forrey
c Institute for Theoretical Atomic and Molecular Physics
c Harvard-Smithsonian Center for Astrophysics
c 60 Garden Street Cambridge, MA 02138
c rforrey@cfa.harvard.edu
c
c***********************************************************************
c
c function name - chyp(a,b,c,z)
c
c computation
c performed - calculates the hypergeometric function.
c
c arguments z - the independent variable of the hypergeometric
c function.
c
c a,b,c - parameters of the hypergeometric function.
c
c precision - complex*16
c
c language - fortran
c
c***********************************************************************
complex*16 function chyp(a,b,c,z)
implicit none
real*8 zero,one,pi,tmp1,tmp2,x,y,r,phi
parameter (zero=0.d0,one=1.d0)
integer flag,n
complex*16 i,a,b,c,z,w,f1,f2,cgamm,csum,
# ctmp,ctmp1,ctmp2,ctmp3,ctmp4,ogamm
pi=dacos(-1.d0)
i=dcmplx(0.d0,1.d0)
c Set error messages
tmp1=dreal(a-b)
tmp2=dimag(a-b)
n=nint(tmp1)
tmp1=tmp1+tmp2-n
if(tmp1.eq.zero)then
write(*,*)'transformation error in chyp'
return
endif
tmp1=dreal(c-a-b)
tmp2=dimag(c-a-b)
n=nint(tmp1)
tmp1=tmp1+tmp2-n
if(tmp1.eq.zero)then
write(*,*)'transformation error in chyp'
return
endif
c Handle the special cases when z=0 or z=1
tmp1=cdabs(z)
if (tmp1.eq.zero) then
chyp=dcmplx(one)
return
endif
ctmp=dcmplx(one)
if (z.eq.ctmp) then
c chyp=cgamm(c)*cgamm(c-a-b)/cgamm(c-a)/cgamm(c-b)
chyp=cgamm(c)*cgamm(c-a-b)*ogamm(c-a)*ogamm(c-b)
return
endif
c Transform to a new variable w whose magnitude is the
c lowest possible value that lies between 0 and 1.
tmp1=cdabs(z)
flag=3
w=1/(1-z)
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=1
endif
w=z/(z-1)
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=2
endif
w=1-z
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=4
endif
w=1-1/z
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=5
endif
w=1/z
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=6
endif
if(tmp1.ge.one)then
write(*,*)'error in transformation'
return
endif
write(20,*)'flag=',flag
c Compute the hypergeometric function in z via the transformation
c theory and the series representation in w.
if (flag .eq. 1)then
w=1/(1-z)
f1=csum(a,c-b,a-b+1,w)
f2=csum(b,c-a,b-a+1,w)
x=dreal(1-z)
y=dimag(1-z)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.lt.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 1'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp1=ctmp**(-a)
ctmp2=ctmp**(-b)
c chyp=f1*ctmp1*cgamm(c)*cgamm(b-a)/cgamm(b)/cgamm(c-a)
c # +f2*ctmp2*cgamm(c)*cgamm(a-b)/cgamm(a)/cgamm(c-b)
chyp=f1*ctmp1*cgamm(c)*cgamm(b-a)*ogamm(b)*ogamm(c-a)
# +f2*ctmp2*cgamm(c)*cgamm(a-b)*ogamm(a)*ogamm(c-b)
c write(*,'(2(2 e13.6,2x))')ctmp1,ctmp1-w**a
c write(*,'(2(2 e13.6,2x))')ctmp2,ctmp2-w**b
elseif (flag .eq. 2)then
w=z/(z-1)
chyp=(1-w)**a*csum(a,c-b,c,w)
elseif (flag .eq. 3)then
w=z
chyp=csum(a,b,c,w)
elseif (flag .eq. 4)then
w=1-z
f1=csum(a,b,a+b-c+1,w)
f2=csum(c-a,c-b,c-a-b+1,w)
x=dreal(w)
y=dimag(w)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.lt.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 4'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp2=ctmp**(c-a-b)
c chyp=f1*cgamm(c)*cgamm(c-a-b)/cgamm(c-a)/cgamm(c-b)
c # +f2*ctmp2*cgamm(c)*cgamm(a+b-c)/cgamm(a)/cgamm(b)
chyp=f1*cgamm(c)*cgamm(c-a-b)*ogamm(c-a)*ogamm(c-b)
# +f2*ctmp2*cgamm(c)*cgamm(a+b-c)*ogamm(a)*ogamm(b)
c write(*,'(2(2 e13.6,2x))')ctmp2,ctmp2-w**(c-a-b)
elseif (flag .eq. 5)then
w=1-1/z
f1=csum(a,a-c+1,a+b-c+1,w)
f2=csum(c-a,1-a,c-a-b+1,w)
x=dreal(z)
y=dimag(z)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.lt.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 5'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp1=ctmp**(-a)
ctmp2=ctmp**(a-c)
x=dreal(1-z)
y=dimag(1-z)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.le.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 5'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp2=ctmp2*ctmp**(c-a-b)
c chyp=f1*ctmp1*cgamm(c)*cgamm(c-a-b)/cgamm(c-a)/cgamm(c-b)
c # +f2*ctmp2*cgamm(c)*cgamm(a+b-c)/cgamm(a)/cgamm(b)
chyp=f1*ctmp1*cgamm(c)*cgamm(c-a-b)*ogamm(c-a)*ogamm(c-b)
# +f2*ctmp2*cgamm(c)*cgamm(a+b-c)*ogamm(a)*ogamm(b)
c write(*,'(2(2 e13.6,2x))')ctmp1,ctmp1-z**(-a)
c write(*,'(2(2 e13.6,2x))')ctmp2,ctmp2-z**(a-c)*(1-z)**(c-a-b)
elseif (flag .eq. 6)then
w=1/z
f1=csum(a,a-c+1,a-b+1,w)
f2=csum(b-c+1,b,b-a+1,w)
x=dreal(-z)
y=dimag(-z)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.le.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 6'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp1=ctmp**(-a)
ctmp2=ctmp**(-b)
c chyp=f1*ctmp1*cgamm(c)*cgamm(b-a)/cgamm(b)/cgamm(c-a)
c # +f2*ctmp2*cgamm(c)*cgamm(a-b)/cgamm(a)/cgamm(c-b)
chyp=f1*ctmp1*cgamm(c)*cgamm(b-a)*ogamm(b)*ogamm(c-a)
# +f2*ctmp2*cgamm(c)*cgamm(a-b)*ogamm(a)*ogamm(c-b)
c write(*,'(2(2 e13.6,2x))')ctmp1,ctmp1-(-z)**(-a)
c write(*,'(2(2 e13.6,2x))')ctmp2,ctmp2-(-z)**(-b)
endif
c write(20,'(i2,2x,2(e13.6,2x))')flag,cdabs(z),cdabs(w)
return
end
c*******************************************************************
c
c234567
complex*16 function csum(a,b,c,z)
implicit none
integer n
real*8 test,eps
complex*16 a,b,c,z,num,den,tmp,sum
eps=1.d-14
tmp=dcmplx(1.d0)
sum=dcmplx(0.d0)
do 10 n=1,300
sum=sum+tmp
test=cdabs(tmp/sum)
if(test.gt.eps)then
num=(a+n-1)*(b+n-1)
den=(c+n-1)*n
tmp=tmp*num*z/den
else
go to 20
endif
10 continue
write(*,*)'warning in csum: sum not completely converged,
#|z|=',cdabs(z)
20 continue
csum=sum
return
end
c*******************************************************************
c
c function name - ogamm
c
c computation
c performed - calculates one over the gamma function
c
c usage - ogamm(z)
c
c argument z - any complex number
c
c precision - double
c
c language - fortran
c
c*******************************************************************
complex*16 function ogamm(z)
implicit none
real*8 one
parameter (one=1.d0)
complex*16 z,ztmp,coeff,cgamm
ztmp=z
coeff=dcmplx(one)
1 if (dreal(ztmp).lt.one) then
coeff=coeff*ztmp
ztmp=ztmp+one
go to 1
endif
ogamm=coeff/cgamm(ztmp)
return
end
c*******************************************************************
c
c function name - cgamm
c
c computation
c performed - calculates the complex gamma function
c
c usage - cgamm(z)
c
c argument z - any complex number
c
c precision - double
c
c language - fortran
c
c*******************************************************************
complex*16 function cgamm(z)
implicit none
integer k
real*8 c(0:15),pi,one
parameter (one=1.d0)
complex*16 z,ztmp,sum,tmp,const
pi=acos(-1.d0)
call coeff(c,15)
ztmp=z
const=dcmplx(one)
1 if (dreal(ztmp).lt.one) then
const=const*ztmp
ztmp=ztmp+one
go to 1
endif
sum=cmplx(0.d0)
tmp=cmplx(1.d0)
do 10 k=0,15
sum=sum+c(k)*tmp
tmp=tmp*(ztmp-k-1)/(ztmp+k)
10 continue
cgamm=sum*dsqrt(2*pi)*(ztmp+9.d0/2)**(ztmp-.5d0)
# *exp(-ztmp-9.d0/2)/const
return
end
c*******************************************************************
c
c subroutine name - coeff
c
c computation
c performed - tabulates the gamma function coefficients
c (see luke vol.2 p.304)
c
c usage - call coeff(c,n)
c
c arguments c - the array (output) which contains the
c coefficients.
c
c n - the dimension (input) of the array 'c'.
c
c precision - double
c
c language - fortran 77
c
c*******************************************************************
subroutine coeff(c,n)
real*8 c(0:n)
c(0) = 41.624436916439068d0
c(1) =-51.224241022374774d0
c(2) = 11.338755813488977d0
c(3) = -0.747732687772388d0
c(4) = 0.008782877493061d0
c(5) = -0.000001899030264d0
c(6) = 0.000000001946335d0
c(7) = -0.000000000199345d0
c(8) = 0.000000000008433d0
c(9) = 0.000000000001486d0
c(10)= -0.000000000000806d0
c(11)= 0.000000000000293d0
c(12)= -0.000000000000102d0
c(13)= 0.000000000000037d0
c(14)= -0.000000000000014d0
c(15)= 0.000000000000006d0
return
end
Quella che hai riportato però mi sembra la 2F1, non la F1 di Appell.BlueEyes ha scritto:C'è però dentro (il link) un file pdf, da dove ho potuto ricavare il seguente codice Fortran, che ritengo possa essere utile allo studente che ha aperto il topic.
giuseppe morello ha scritto: 2, sull'algoritmo: per a=-1, b=5, c=3, z=4 e anche in alcuni altri casi, il codice non converge, ritenendo che hnext risulti troppo piccolo. Per gli stessi valori wolphram non mostra problemi. In effetti, l'ho fatto stampare e risulta nullo almeno fino alla 15^a cifra, mentre hmin=0.0001 (anche se non ho ancora trovato dove viene assegnato questo valore).
Codice: Seleziona tutto
$ python
Python 2.7.1+ (r271:86832, Apr 11 2011, 18:05:24)
[GCC 4.5.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import mpmath
>>> float(mpmath.hyp2f1(float(-1), float(5), float(3), float(4)))
-5.666666666666667
>>>
Codice: Seleziona tutto
>>> import mpmath
>>> mpmath.hyp2f1(-1, 5, 3, 4+3j)
mpc(real='-5.666666666666667', imag='-5.0')