gauss-legendre в С++

я пытаюсь создать код gauss-legendre в соответствии со следующим алгоритмом:

для n точек для n точек

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

ti являются корнями многочленов Лежандра порядка n. Заданы многочлены Лежандра:

и:

Мой код:

#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) Я не уверен, что все сделал правильно. Я имею в виду, если я объединил все вместе и правильно ли реализовал алгоритм.


person George    schedule 16.03.2011    source источник


Ответы (4)


Я не знаю алгоритма, но ваш код неверен.
Во-первых:

        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 вместо указателей.

person ali_bahoo    schedule 17.03.2011
comment
Спасибо за советы. Я проверяю их. Вы правы. Но я хочу использовать указатели, чтобы учиться. - person George; 17.03.2011

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

w[index] * f(*t)

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

person Armen Tsirunyan    schedule 16.03.2011
comment
@George: понятия не имею: я указываю на синтаксические/языковые проблемы. Я не знаком с алгоритмом - person Armen Tsirunyan; 16.03.2011

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

person ariadne3s    schedule 11.04.2011

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

person Dan    schedule 09.06.2011