1

我正在尝试根据以下算法创建 gauss-legendre 代码：

对于 n 个点 对于 n 个点

也就是说，它创建了一个 2n 方程系统（如果我们要求对 2n-1 阶的多项式准确，

ti是n阶legendre多项式的根。legendre多项式给出：

和 wi ：

我的代码是：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <cmath>


using namespace std;

const double pi=3.14;


//my function with limits (-1,1)
double f(double x){
double y;
y=(pi/4.0)*(log((pi*(x+1.0))/4.0 +1.0));
return y;

}

double legendre (int n){

    double *L,*w,*t;
    double x,sum1,sum2,result;
    L=new double [n];
    w=new double [n];
    t=new double [n];


        while(n<10){

         L[0]=1;
         L[1]=x;


        //legendre coef
        for (int i=1;i<=10;i++){
        L[i+1]=((2.0*i+1.0)*x*L[i] - i*L[i-1])/(i+1.0);


        }

        //weights w
        w=0;
        for (int i=1;i<=10;i++){
        w[i]+=(2.0*(1.0-x*x))/(i*i*(L[i-1]*L[i-1]));
        }


        //sums  w*t
        for (int i=1;i<=10;i++){
            sum1=0.0; //for k=1,3,5,2n-1
            for (int k=1;k<=2*n-1;k+=2){
                sum1+=w[i]*(pow(t[i],k));
            }
                sum1=0;
                sum2=0.0;//for k=0,2,4,2n-2
                for(int k=0;k<=2*n-2;k+=2){
                    sum2+=w[i]*(pow(t[i],k));
                }
                sum2=2.0/n;
        }


    }

    result=w*f(*t);

    return result;

}


int main()
{
    double eps=1e-8;//accuracy
    double exact=0.8565899396;//exact solution for the integral
    double error=1.0;
    double result;

    int n=1;//initial point


    while (fabs(result-exact)>eps) {
        result=legendre(n);
        cout <<"\nFor n = "<<n<<",error = "<<fabs(error-exact)<<",value = "<<result;

    n++;
    }

    return 0;
}

我的问题是：

1) 编译器给了我 :error: 'double*' 和 'double' 类型的无效操作数到二进制 'operator*' --> 在 result=w*f(*t);

2）我不确定我是否做对了整个事情。我的意思是，如果我将所有事情组合在一起并且我是否正确实施了算法。

4

4 回答 4

2

w 是一个指针，你正试图将它与某事相乘......你必须使用 index

w[index] * f(*t)

也是*tt 数组的第一个元素。你是这个意思吗？

于 2011-03-16T18:23:10.280 回答
2

我不知道算法，但你的代码是错误的。
第一的 ：

        while(n<10)
        {
         L[0]=1;
         L[1]=x;
        //legendre coef
        for (int i=1;i<=10;i++){
        L[i+1]=((2.0*i+1.0)*x*L[i] - i*L[i-1])/(i+1.0);
        }

您必须更改n(increment, decrement, etc.) 的值，否则这是一个无限循环。

第二 ： 

//weights w
    w=0;
    for (int i=1;i<=10;i++){
    w[i]+=(2.0*(1.0-x*x))/(i*i*(L[i-1]*L[i-1]));
    }

w是一个指针。如果要重置它，请使用memset(w,0,sizeof(double)*n);不要使其等于 0。

最后的：

result=w*f(*t);

由于您将wandt指针用作数组，因此您必须提供某种索引，例如result=w[ind] * f(t[ind]);. 在这里，您只是将指针w的值，而不是所指向的w值（顺便说一下，您的代码中的值w是 0）与所指向的双精度数组的第一个值相乘t。

此外，我无法得到您的代码与问题中的公式之间的任何关系。所以我的谦虚建议是不要为此使用 C 或 C++。如果必须，请不要使用指针，因为您似乎不熟悉它们。您可以轻松地使用 std::vector 而不是指针。

于 2011-03-17T12:51:20.983 回答
1

关于您的算法，x（横坐标值）应该是勒让德多项式的零点。我没有看到你在任何地方定义它们。定义它们有点痛苦。我正在做类似的事情，发现这个（它是一个 Matlab 文件，而不是一个 C++ 文件）定义了 N xi 和 wi 值。该算法工作正常： http: //www.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes

于 2011-04-11T16:31:45.143 回答
1

alglib和gsl都有高斯正交的实现。不过，两者都是 GPL 许可的。

于 2011-06-09T07:27:42.540 回答