gauss-legendre in c++
i am trying to create gauss-legendre code according to the following algorithm:
for n points
That is,it is created a 2n equation system (if we demand to be accurate for polynominals of order 2n-1 ,
ti are roots of the legendre polynominals of order n.The legendre poynominals are given :
and wi :
My code is :
#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;
}
My problems are:
1) The compiler gives me :error: invalid operands of types ‘double*’ and ‘double’ to binary ‘operator*’ --> at result=w*f(*t);
2) I am not su开发者_StackOverflow社区re if i have done the whole thing right.I mean ,if i combined all the things together and if i implemented right the algorithm.
I do not know the algorithm but your code is wrong.
First :
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);
}
You must change the value of n
(increment, decrement, etc.) otherwise this is an endless loop.
Second :
//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
is a pointer. If you want to reset it, use memset(w,0,sizeof(double)*n);
Do not make it equal to 0.
Last:
result=w*f(*t);
Since you are using the w
and t
pointers as arrays, you have to provide some sort of index like result=w[ind] * f(t[ind]);
. Here you are simply multiplying the value of pointer w
, not the value that is pointed by w
(the value of w
is 0 form your code by the way) with the first value of the double array pointed by t
.
Also I could not get any relation between your code and the formulas in the question. So my humble advise is do not use C or C++ for this. If you must, then do not use pointers, because it seems you are not familiar with them. You can easily have std::vector instead of pointers.
w is a pointer, and you are trying to multiply it with something... you must use index
w[index] * f(*t)
also *t
is the first element of t array. Is that what you mean?
Regarding your algorithm, the x (the abscissa values) are supposed to be the zeros of the Legendre polynomial. I didn't see you define them anywhere. They're a bit of a pain to define. I was doing something similar and found this (it's a Matlab file, not a C++ file) that defines N xi and wi values. The algorithm works fine: http://www.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes
Both alglib and gsl have implementations of gaussian quadrature. Both are GPL licensed, though.
精彩评论