integrating a function 2

[aaadetos]

The following code below is from a text on numerical recipes in C++...Its purpose is to evaluate the integral of a function, "func", using the Gaussian Quadrature formula, viz:

double qgaus(double func, double a, double b)
{
	//elements of the integral
	double xr,xm,dx,s;	//s unrelated to Laplace
	//at abscissas x and with weights w,
	double x[] = {0.1488743389816312,0.4333953941292472,
		0.6794095682990244,0.8650633666889845,0.9739065285171717};
	double w[] = {0.2955242247147529,0.2692667193099963,
		0.2190863625159821,0.1494513491505806,0.0666713443086881};

	xm = 0.5*(b+a);
	xr = 0.5*(b-a);
	s = 0;
	
	//performing the integration on func
	for(int j=0;j<5;j++)
	{
		dx = xr*x[j];
		s += w[j]*(func(xm+dx)+func(xm-dx));
	}
	return s *= xr;
}

Im trying to adapt the code to integrate my function, "func", which returns an array of values, viz:

void func(double *func,CKuchukDlg *pDialog)
{
	enum{nSize = 100};
	double Ip[nSize],p[nSize],gh[1][nSize];
	double alpha[nSize],beta[nSize];	//this 's' is unrelated to Laplace's

	for(int t=0;t<nSize;t++)
	{
		if(t==0)
		{
				alpha[0] = beta[0] = 1;
				p[0] = sqrt(pow(alpha[0],2.0)+pow(beta[0],2.0));

				homo(gh,nSize,pDialog);
				Integrals(Ip,20);

				func[0] = (gh[0][0] * Ip[0])/p[0];
		}
		else
		{
			alpha[t] = beta[t] = 2*t*PI;

			homo(gh,nSize,pDialog);
			Integrals(Ip,20);
		
			p[t] = sqrt(pow(alpha[t],2.0)+pow(beta[t],2.0));
		
			func[t] = (gh[0][t] * Ip[t])/p[t];
		}
	}
}

Can anyone please show me how I may do this? I have tried a number of ways, but I keep running into errors. Thank you!

 

[aaadetos]

I modified my function to:

void func(double (*func)[nSize],CKuchukDlg *pDialog)
{
    enum{nSize = 100};
    double Ip[nSize],p[nSize],gh[1][nSize];
    double alpha[nSize],beta[nSize];    //this 's' is unrelated to Laplace's

    for(int t=0;t<nSize;t++)
    {
        if(t==0)
        {
                alpha[0] = beta[0] = 1;
                p[0] = sqrt(pow(alpha[0],2.0)+pow(beta[0],2.0));

                homo(gh,nSize,pDialog);
                Integrals(Ip,nSize);

                func[0][0] = (gh[0][0] * Ip[0])/p[0];
        }
        else
        {
            alpha[t] = beta[t] = 2*t*PI;

            homo(gh,nSize,pDialog);
            Integrals(Ip,nSize);
        
            p[t] = sqrt(pow(alpha[t],2.0)+pow(beta[t],2.0));
        
            func[0][t] = (gh[0][t] * Ip[t])/p[t];
        }
    }
}

 

[Salem]

> void func(double (*func)[nSize],CKuchukDlg *pDialog)
Well having a parameter the same name as the function isn't a good idea.

Also, you've posted several different versions now, which has the problem?

> but I keep running into errors.
Well posting how you call the function, and your error messages would be helpful.

You should be able to show how you are calling your function from a very minimal main() (like this example), which just calls the function you're having problems with.

void testcode ( char *p ) {
}
int main ( ) {
  char arr[10];
  testcode(arr);
  return 0;
}

That way, it's easy for us to look at. All we need to do is copy/paste the code to our own compilers, and see what happens.
It's all too easy to create a plausable call to your function which would work, but that doesn't help you.

--

 

[aaadetos]

Thank you Salem. Here's my function:
I tried implementing it in the function below, but I have error C2664 at line (*): Function cannot convert parameter 2 from 'void(_cdecl*)(double*,double,class CKuchukDlg*)to double* and at line (**):Function cannot convert parameter 1 from double to double*.

void Integrals(double *Ip,int nSize);
void homo(double (*gh)[nSize], int nSize,CKuchukDlg *pDialog);
void function(double *p,double *func,int nSize,CKuchukDlg *pDialog);


void function(double *p,double *func,int nSize,CKuchukDlg *pDialog)
{
	enum{nSize = 10};
	double Ip[nSize],gh[1][nSize];
	double alpha[nSize],beta[nSize];	//this 's' is unrelated to Laplace's

	for(int t=0;t<nSize;t++)
	{
		if(t==0)
		{
				alpha[0] = beta[0] = 1;
				p[0] = sqrt(pow(alpha[0],2.0)+pow(beta[0],2.0));

				homo(gh,nSize,pDialog);
				Integrals(Ip,nSize);

				func[0] = (gh[0][0] * Ip[0])/p[0];
		}
		else
		{
			alpha[t] = beta[t] = 2*t*PI;

			homo(gh,nSize,pDialog);
			Integrals(Ip,nSize);
		
			p[t] = sqrt(pow(alpha[t],2.0)+pow(beta[t],2.0));
		
			func[t] = (gh[0][t] * Ip[t])/p[t];
		}
	}
}
//-------------------------------------------------
//Implementing the function in the integral
void qgaus(double *ans, int nSize, CKuchukDlg *pDialog)
{
	enum{nSize = 10};
	//elements of the integral
	double xr,xm,dx;	//s unrelated to Laplace
	double *p = new double[nSize];
	double a = 0;
	double b = 1e+20;
	//at abscissas x and with weights w,
	double x[] = {0.1488743389816312,0.4333953941292472,
		0.6794095682990244,0.8650633666889845,0.9739065285171717};
	double w[] = {0.2955242247147529,0.2692667193099963,
		0.2190863625159821,0.1494513491505806,0.0666713443086881};

	xm = 0.5*(b+a);
	xr = 0.5*(b-a);

	//performing the integration on func
	for(int t=0;t<nSize;t++)
	{
		ans[t] = 0;
		for(int j=0;j<5;j++)
		{
			function(p,func,nSize,pDialog);//(*)
			dx = xr*x[j];
			ans[t] += w[j]*(function(xm+dx,func,nSize,pDialog)[t]+function(xm-dx,func,nSize,pDialog)[t]);//(**)
		}
		ans[t] *= xr;
	}
}

 

[timmay3141]

At (*), you pass func as the second parameter, which is apparently the name of a void function when you are expecting a pointer to a double.

At (**), you call function where the first parameter xm + dx is a double instead of a pointer to a double (it looks like it should be an array).

You seem to be confused about how you want to integrate the function. Using arrays is pretty cumbersome here. Unless there's a very good reason for doing this, I'd forget the arrays and recommend having your integrate function look like this:

double Integrate(double (*func)(double x), double stepsize, double start, double stop);

It shouldn't be difficult to write the implimentation of this function...just use a for loop going from start to stop in intervals of stepsize and evaluating func at the intervals. I've never used the Gaussian Quadrature formula, but it shouldn't be too hard to adapt this to that.

 

[globos]

From your first question, here is a sample source code from which you can restart :

#include <functional>// Give some facilities on pointers on functions.
#include <iostream>
using namespace std;

// Type alias for pointer on unary functions returning reals.
typedef pointer_to_unary_function<double, double> pointer_to_real_function;

// A dummy identity function (implement your own one).
double id_function (double x)
{
  return x;
}

// Pass a pointer to the function to integrate (see the main).
double qgaus (pointer_to_real_function func_ptr,
	      double a, double b)
{
    //elements of the integral
    double xr,xm,dx,s;    //s unrelated to Laplace
    //at abscissas x and with weights w,
    double x[] = {0.1488743389816312,0.4333953941292472,
        0.6794095682990244,0.8650633666889845,0.9739065285171717};
    double w[] = {0.2955242247147529,0.2692667193099963,
        0.2190863625159821,0.1494513491505806,0.0666713443086881};

    xm = 0.5*(b+a);
    xr = 0.5*(b-a);
    s = 0;
    
    //performing the integration on func
    for(int j=0;j<5;j++)
    {
        dx = xr * x[j];
        s += w[j] * (func_ptr (xm + dx) + func_ptr (xm - dx));
    }
    return s *= xr;
}


int main ()
{
  cout << qgaus (ptr_fun (id_function), 0, 15) << endl;
  return 0;
}

You should follow the above code, and try to not mix GUI entities (like CKuchukDlg *) with the computation. Or make an dedicated abstraction for that (if you want to see the progress of the computation for example).

--
Globos

 

[aaadetos]

Thanks, timmay3141 and globos; I eventually wrote the integral thus, but have a problem calling it:

double* function(double *p,double *func,int nSize,CKuchukDlg *pDialog);

void qgaus(double *ans, int nSize, CKuchukDlg *pDialog)
{
	enum{nSize = 6};
	//elements of the integral
	double xr,xm,dx;	//s unrelated to Laplace
	double *p = new double[nSize];
	double a = 0;
	double b = 1e+20;
	double *func = NULL;

	func = new double[nSize];
	double deltapos = 0;
	double deltaneg = 0;
	//at abscissas x and with weights w,
	double x[] = {0.1488743389816312,0.4333953941292472,
		0.6794095682990244,0.8650633666889845,0.9739065285171717};
	double w[] = {0.2955242247147529,0.2692667193099963,
		0.2190863625159821,0.1494513491505806,0.0666713443086881};

	xm = 0.5*(b+a);
	xr = 0.5*(b-a);

	//performing the integration on func
	for(int t=0;t<nSize;t++)
	{
		ans[t] = 0;
		for(int j=0;j<5;j++)
		{
			function(p,func,nSize,pDialog);
			dx = xr*x[j];
			deltapos = xm + dx;
			deltaneg = xm - dx;
			ans[t] += w[j]*(function(&deltapos,func,nSize,pDialog)[t]+function(&deltaneg,func,nSize,pDialog)[t]);
		}
		ans[t] *= xr;
	}
}

But when i attempt to call it here, the debugger doesn't respond:

extern double PI;
double *kh_md = NULL,*kv_md = NULL,*por = NULL,*visc_cp = NULL,*ct_psi = NULL;
void unbound(double *unbd,int nSize,CKuchukDlg *pDialog);
void integral(double *ans,int nSize,CKuchukDlg *pDialog);//from cplxintegral.cpp
void qgaus(double *ans,int nSize,CKuchukDlg *pDialog);//from cplxintegral.cpp
void Time(double tStart, double tStop, double *time_hr, const int ntime);

void gwkfinal(double *gwk, int nSize, CKuchukDlg *pDialog)	
{
	enum{nSize = 6};
	double unbd[nSize];
	double ans[nSize];//from cplxintegral.cpp

	for (int i=0; i < nSize; i++)
	{
		unbd[i] = 0;
		ans[i] = 0;
	}

	unbound(unbd,nSize,pDialog);
	//integral(ans,nSize,pDialog);
	qgaus(ans,nSize,pDialog);

	for (int t=0;t<nSize;t++)
	{
		gwk[t] = unbd[t] + ans[t];	//equation (6)
	}
}

 

[globos]

Why don't you follow the advices that were given?
Again, separate concerns in your code.
There is really several things that are wrong with your code. Why have you written this :

void qgaus(double *ans, int nSize, CKuchukDlg *pDialog)
{
    enum{nSize = 6};
    //...

The signature you made of qgaus is pretty weird. Why do you use a pointer to return the result? What is nSize for ? (and what is the enum "enum {nSize = 6}" for).
There is also allocated memory that seems to be never freed ("double *p = new double[nSize];" and the same for "func").

Again should simplify your code, this is all too cryptic.



--
Globos