je code en C++ et je voudrais réaliser une régression polynomiale d'un certain ordre avec x et y mes tableaux de coordonnées, a mon tableau de coefficients du polynome obtenu, et order_poly l'ordre du polynome souhaité.
voici mon code, si qqn peut me dire ce qui cloche, je n'y arrive pas!!
void polyReg(double *x, double *y, int sizexy, double *a, int order_poly)
{
//realize polynomial interpolation.
//x : array of all points x coordinates.
//y : array of all points y coordinates.
//sizexy : size of x and y arrays.
//a : polynomial coeff returned.
//order_poly : order of polynom wanted.
//-----"Moindres Carrés" Method-----
int ord = order_poly + 1;
int n = sizexy;
double **A = (double **)malloc( ord *sizeof(double *));
for(int i=0; i<ord; i++) A[i] = (double *)calloc(ord, sizeof(double));
//à voir la dimension de S
double *S = (double *)malloc( (2*order_poly +1 ) * sizeof(double));
double *W = (double *)calloc(ord, sizeof(double));
double *sol = (double *)calloc(ord, sizeof(double));
//tableaux x et y OK
//-----computing the Sk-----
for(int k=0; k< 2*order_poly+1 ; k++)
{
S[k]=0;
for(int i=1; i< n+1; i++)
{
S[k] = S[k] + pow(x[i],k);
}//for i
}//for k
//-----computing the Wk-----
for(int k=0; k< order_poly+1; k++)
{
W[k-1]=0;
for(int i=1; i<n+1; i++)
{
W[k-1] = W[k-1] + y[i-1] * pow(x[i-1],k-1);
}//for i
}//for k
//---computing matrix coefficients---
for(int i =1; i<=order_poly+1; i++)
{
for(int j =1; j<=order_poly+1; j++)
{
if(i==1 && j==1)
A[i-1][j-1] = n;
else
A[i-1][j-1] = S[i+j-2]; //à voir pour le -2????
}//for j
}//for i
gauss_solv(A,ord, W,sol);
for(k=0; k<ord; k++)
a[k] = sol[k];
}//end polyReg
void gauss_solv(double **A,int sizeofA, double *W, double* sol)
{
//to solve a system A*sol = W
long n = sizeofA;
long i, j, k, maxi;
// double *X = (double *)malloc( n * sizeof(double));
//M is a matrix with n lines and n+1 columns.
double **M = (double **)malloc( n *sizeof(double *));
for(int i=0; i<n; i++) M[i] = (double *)malloc( (n+1) * sizeof(double));
for( i=0; i<n; i++)
{
for( j=0; j<n;j++)
{
M[i][j]= A[i][j];
}//for j
M[i][n]=W[i];
}//for i
for( i = n-1; i>0 ; i--)
{
maxi = i;
for( j =0 ; j<i-1; j++)
{
if( abs(M[i][j])> abs(M[maxi][i]) )
maxi=j;
}//for j
if (maxi != i)
{
for( k=0; k<n; k++)
{
double temp=M[i][k];
M[i][k]=M[maxi][k];
M[maxi][k]=temp;
}//for k
for( j=0; j<(i-1); j++)
{
M[j][n+1] = W[j] - W[i] * M[j][i] / M[i][i];
for( k=0; k<i; i++)
{
M[j][k] = M[j][k] - M[i][k] * M[j][i] / M[i][i];
}//for k
}//for j
}//end if
}//for i
sol[0]= W[0]/M[0][0];
for( i = 1; i< n; i++)
{
sol[i]=W[i];
for( k = 0; k< i-1; k++)
{
sol[i] = sol[i] - M[i][k] * sol[k];
}//for k
sol[i] = sol[i] - M[i][i];
}//for i
}//end gauss_solv