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*' に与えます --> at result=w*f(*t);

2）すべてを正しく行ったかどうかはわかりません。つまり、すべてを組み合わせて、アルゴリズムを正しく実装したかどうかです。

4

4 に答える 4

2

w はポインターであり、それを何かで乗算しようとしています... index を使用する必要があります

w[index] * f(*t)

また*t、t 配列の最初の要素です。そうですか？

于 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(インクリメント、デクリメントなど)の値を変更する必要があります。変更しないと、無限ループになります。

2番 ： 

//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);Do not make it equal to 0 を使用してください。

最後：

result=w*f(*t);

wとポインターを配列として使用してtいるため、 のような何らかのインデックスを提供する必要がありますresult=w[ind] * f(t[ind]);。ここでは、pointerwの値ではなく、pointedwの値 (ちなみにコードからの値wは 0 です) に double 配列の最初の値を掛けていますt。

また、あなたのコードと質問の数式との関係を得ることができませんでした。したがって、私の謙虚なアドバイスは、これには C または C++ を使用しないことです。必要に応じて、ポインターを使用しないでください。ポインターに慣れていないようです。ポインターの代わりに std::vector を簡単に使用できます。

于 2011-03-17T12:51:20.983 に答える
1

アルゴリズムに関して、x (横座標値) は、ルジャンドル多項式のゼロであると想定されています。あなたがそれらをどこでも定義しているのを見たことがありません。それらを定義するのは少し面倒です。私は似たようなことをしていて、N xi と wi の値を定義する this (C++ ファイルではなく、Matlab ファイルです) を見つけました。アルゴリズムは正常に動作します: 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 に答える