гаусс-легенда в с ++

Question

Я пытаюсь создать код Гаусс-Легендра по следующему алгоритму:

для n баллов

То есть создается система уравнений 2n (если мы требуем точности для многочленов порядка 2n-1,

ti - корни легендарных полиномов порядка n. Легендарные полиномы даны:

и 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" для двоичного оператора "*" -> at result = w * f (* t);

2) Я не уверен, правильно ли я все сделал. Я имею в виду, собрал ли я все вместе и правильно ли реализовал алгоритм.

Answer 1

Я не знаю алгоритм, но ваш код неверен.
Первый:

        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 (увеличение, уменьшение и т. Д.), Иначе это бесконечный цикл.

Второй:

//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);

Поскольку вы используете указатели w и t в качестве массивов, вы должны предоставить какой-то индекс, например result=w[ind] * f(t[ind]);. Здесь вы просто умножаете значение указателя w, а не значение, на которое указывает w (значение w равно 0, образуйте ваш код, кстати) на первое значение двойного массива, на которое указывает t.

Также я не смог получить какую-либо связь между вашим кодом и формулами в вопросе. Поэтому мой скромный совет: не используйте C или C ++ для этого. Если вы должны, то не используйте указатели, потому что кажется, что вы не знакомы с ними. Вы можете легко использовать std :: vector вместо указателей.

Answer 2

w - указатель, и вы пытаетесь умножить его на что-то ... вы должны использовать индекс

w[index] * f(*t)

также *t - первый элемент массива t. Это то, что вы имеете в виду?

Answer 3

Оба alglib и gsl имеют реализации гауссовской квадратуры. Оба лицензированы по лицензии GPL.

Answer 4

Что касается вашего алгоритма, предполагается, что x (значения абсциссы) являются нулями полинома Лежандра. Я не видел, чтобы вы их нигде не определяли. Их немного сложно определить. Я делал что-то похожее и нашел это (это файл Matlab, а не файл C ++), который определяет значения N xi и wi. Алгоритм работает отлично: http://www.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes