Re: about QAWF integration -- Fourier integrals

[xie qiong <xieq at yahoo dot com>]

Hello,

Thank you for your reply. 

The following is my code. Even though I changed the
tolerance to 0.1, core dump still happens.


double myfunc(double x, void *params)
{
        double a = *(double *)params;
        
	return exp(- log(x) * log(x)/(2.0 * a * a)) /
(sqrt(2.0 * PI) * a * x);
}

int main ()
{
        int i, m;
	double w0, error1 = 0.00, error2 = 0.00;
	double imag[ARRAYSIZE] = {0.0}; 
	double real[ARRAYSIZE] = {0.0}; 
	double a = 0.25;
	FILE *fp;

	gsl_integration_workspace * wp; 
	gsl_integration_workspace * c_wp;
	gsl_integration_qawo_table * wf = 
	      gsl_integration_qawo_table_alloc(0.0,          
    
					       10000.0,           
					       GSL_INTEG_COSINE,  
					       10000);            
	
	gsl_function F;
	F.function = &myfunc;
	F.params = &a;
	
	for ( w0 = 0.0, m = 0; w0 < OMEGAMAX; w0 +=
OMEGASPACE, m++ )  {
	       wp = gsl_integration_workspace_alloc(20000);
	       c_wp = gsl_integration_workspace_alloc(20000);

	       gsl_integration_qawo_table_set(wf,            
    
					      w0,                 
					      10000.0,            
					      GSL_INTEG_COSINE);  

	       gsl_integration_qawf(&F,       
				    1e-30,    
				    1e-4,      
				    20000,    
				    wp,       
				    c_wp,     
				    wf,       
				    real+m,   
				    &error1); 

	       printf("\t\t\t\t%e\n", error1);
	       printf("%e", real[m]);
	       gsl_integration_workspace_free(wp);
	       gsl_integration_workspace_free(c_wp);

	       wp = gsl_integration_workspace_alloc(20000);
	       c_wp = gsl_integration_workspace_alloc(20000);

	       gsl_integration_qawo_table_set(wf, w0,
10000.0, GSL_INTEG_SINE);
	       gsl_integration_qawf(&F, 1e-30, 1e-4, 20000,
wp, c_wp, wf, imag+m, &error2);
	       printf("      %e\n", imag[m]);
	       gsl_integration_workspace_free(wp);
	       gsl_integration_workspace_free(c_wp);
	}

	gsl_integration_qawo_table_free(wf);



Thanks,

Qiong



> rectangular or normal pdf
>  > using this routine. But when I changed the
> integrand to be
>  > log-normal pdf, core dump happens due to the
> ERROR: cannot reach
>  > tolerance because of roundoff error Aborted.
> 
> Only a larger tolerance will prevent that.  Send a
> small example program
> to the list if you still have problems.
> 
>  > I have ever tried to make changes to some
> parameters, such as table
>  > level n, workspace size and desired tolerance
> limit, but still
>  > cannot provent core dump from happening.
>  >  Any solution, help or suggestion would be
> greatly appreciated.
>  >  Regards, Qiong


__________________________________________________
Do You Yahoo!?
Try FREE Yahoo! Mail - the world's greatest free email!
http://mail.yahoo.com/