# Polynomial Fitting – C PROGRAM

Okay, so here I am sharing a code for fitting a polynomial to a given set of data-points using the Least Squares Approximation Method(Wikipedia).

Let’s say we have $N$ data-point pairs and we are trying to fit them using a polynomial of degree $n$. If N=n+1 then the polynomial will pass exactly through each point and it will correspond to the interpolating polynomial that I wrote about earlier.

Let’s say the polynomial we are using is given as: $a_0+a_1x+a_2x^2+a_3x^3+.....+a_nx^n$

with errors given by $e_i=Y_i-y_i=Y_i-a_0-a_1x_i-a_2x_i^2-....-a_nx_i^n$

Here, we are using $Y_i$ to represent the observed data-points corresponding to $x_i$. We now minimize the following quantity $S=\Sigma_{i=1}^Ne_i^2=\Sigma_{i=1}^N(Y_i-y_i=Y_i-a_0-a_1x_i-a_2x_i^2-....-a_nx_i^n)^2$

At the minimum all the partial derivatives with respect to the coefficients will vanish. This will give us the following $n+1$ equations: $\frac{\partial S}{\partial a_0}=0=\Sigma_{i=1}^N2(Y_i-a_0-a_1x_i-....-a_ix_i^n)(-1),$ $\frac{\partial S}{\partial a_1}=0=\Sigma_{i=1}^N2(Y_i-a_0-a_1x_i-....-a_ix_i^n)(-x_i),$
.
.
. $\frac{\partial S}{\partial a_n}=0=\Sigma_{i=1}^N2(Y_i-a_0-a_1x_i-....-a_nx_i^n)(-x_i^n),$

Dividing each by -2 and rearranging gives the $n+1$ normal equations to be solved simultaneously: $\begin{bmatrix} N& \Sigma x_i& \Sigma x_i^2& \dots& \Sigma x_i^n& \\ \Sigma x_i& \Sigma x_i^2& \Sigma x_i^3& \dots& \Sigma x_i^{n+1}& \\ \Sigma x_i^2& \Sigma x_i^3& \Sigma x_i^4& \dots& \Sigma x_i^{n+2}& \\ & \vdots& & & \vdots& \\ \Sigma x_i^n& \Sigma x_i^{n+1}& \Sigma x_i^{n+2}& \dots& \Sigma x_i^{2n}& \end{bmatrix}\begin{bmatrix} a \end{bmatrix}=\begin{bmatrix} \Sigma Y_i\\ \Sigma x_iY_i\\ \Sigma x_i^2Y_i\\ \vdots\\ \Sigma x_i^nY_i \end{bmatrix}$
where $x_i$ and $Y_i$ are the data-points entered by the user and $[a]=[a_0,a_1,..a_n]^T$ which are the required coefficients.

So we just need to build up the above system of equations and then solve it using Gaussian elimination to get the coefficients.

The following program illustrates the process.

### CODE:

/******************************************************
*************Chi-square fitting**************
Polynomial Fitting
******************************************************/
#include<stdio.h>
#include<math.h>
/*******
Function that performs Gauss-Elimination and returns the Upper triangular matrix and solution of equations:
There are two options to do this in C.
1. Pass the augmented matrix (a) as the parameter, and calculate and store the upperTriangular(Gauss-Eliminated Matrix) in it.
2. Use malloc and make the function of pointer type and return the pointer.
This program uses the first option.
********/
void gaussEliminationLS(int m, int n, double a[m][n], double x[n-1]){
int i,j,k;
for(i=0;i<m-1;i++){
//Partial Pivoting
for(k=i+1;k<m;k++){
//If diagonal element(absolute vallue) is smaller than any of the terms below it
if(fabs(a[i][i])<fabs(a[k][i])){
//Swap the rows
for(j=0;j<n;j++){
double temp;
temp=a[i][j];
a[i][j]=a[k][j];
a[k][j]=temp;
}
}
}
//Begin Gauss Elimination
for(k=i+1;k<m;k++){
double  term=a[k][i]/ a[i][i];
for(j=0;j<n;j++){
a[k][j]=a[k][j]-term*a[i][j];
}
}

}
//Begin Back-substitution
for(i=m-1;i>=0;i--){
x[i]=a[i][n-1];
for(j=i+1;j<n-1;j++){
x[i]=x[i]-a[i][j]*x[j];
}
x[i]=x[i]/a[i][i];
}

}
/*******
Function that prints the elements of a matrix row-wise
Parameters: rows(m),columns(n),matrix[m][n]
*******/
void printMatrix(int m, int n, double matrix[m][n]){
int i,j;
for(i=0;i<m;i++){
for(j=0;j<n;j++){
printf("%lf\t",matrix[i][j]);
}
printf("\n");
}
}
main(){
//no. of data-points
int N;
//degree of polynomial
int n;
printf("Enter the no. of data-points:\n");
scanf("%d",&N);
//arrays to store the c and y-axis data-points
double x[N], y[N];
printf("Enter the x-axis values:\n");
int i,j;
for(i=0;i<N;i++){
scanf("%lf",&x[i]);
}
printf("Enter the y-axis values:\n");
for(i=0;i<N;i++){
scanf("%lf",&y[i]);
}
printf("Enter the degree of polynomial to be used:\n");
scanf("%d",&n);
// an array of size 2*n+1 for storing N, Sig xi, Sig xi^2, ...., etc. which are the independent components of the normal matrix
double X[2*n+1];
for(i=0;i<=2*n;i++){
X[i]=0;
for(j=0;j<N;j++){
X[i]=X[i]+pow(x[j],i);
}
}
//the normal augmented matrix
double B[n+1][n+2];
// rhs
double Y[n+1];
for(i=0;i<=n;i++){
Y[i]=0;
for(j=0;j<N;j++){
Y[i]=Y[i]+pow(x[j],i)*y[j];
}
}
for(i=0;i<=n;i++){
for(j=0;j<=n;j++){
B[i][j]=X[i+j];
}
}
for(i=0;i<=n;i++){
B[i][n+1]=Y[i];
}
double A[n+1];
printf("The polynomial fit is given by the equation:\n");
printMatrix(n+1,n+2,B);
gaussEliminationLS(n+1,n+2,B,A);
for(i=0;i<=n;i++){
printf("%lfx^%d+",A[i],i);
}

}



### OUTPUT: So that’s it! That’s how you perform a polynomial fit to a given set of data.

I had written a C++ code for this a long time ago, and coincidentally it got very popular for some reason. But then I felt the need to make an Android app that does the same.

So I ported my code to JAVA so that it works in my Android App.

So if you want you can check out those posts too.

Hope you guys find it useful!
If you have any questions/doubts, hit me up in the comments section below.

Linear Fitting – Lab Write-Up
Linear Fitting – C++ Program
Linear Fitting – Scilab Code
Curve Fit Tools – Android App (using the above code)
Curve Fit Tools – Documentation
Curve Fit Tools – Play Store
Curve Fit Tools – GitHub Repository
Curve Fitters – Scilab Toolbox [wpedon id="7041" align="center"]

## 13 thoughts on “Polynomial Fitting – C PROGRAM”

1. Well done, this is a bit cleaner implementation. Also, this balances out to a database implementation that I have developed.

1. I’m glad you found it useful!

2. I already have a couple of .dat files, and do not want to type the X and Y values, because it would take too long. Basically, let a single program give the polynomial fit equations for different sets of data.

1. That’s extremely easy.
If you’re looking for some code go through some codes on this blog, for ex: plotting exercises

3. Thanks for leading me here, this is quite handy! So I am still trying to incorporating weights to perform a weighted least squares fit. Considering the theoretical equation I would incorporate the weight in the following position as ‘w’:

Y[i]=Y[i]+ w * pow(x[j],i)*y[j];

Since their it should be reside after the derivations.

Thank you in advance for you input!

Greetings
Eric

1. Hi Eric,
I’m looking at doing the same thing, weighted LS fit. Did you ever try this? Do you have hints?
Thanks,
Randy

4. Wonderful program! It just worked without any problems.
And the explanation too is really easy to understand.
Greetings,
Avinash

1. Thanks!

5. Hello Manas,

can you share a surface fit polynomial c/c++ code or the logic to write one?

Thanks,
Rajeev

6. Hello Manas,
I am trying to write a linear surface fit algorithm f(x,y)=A+Bx+Cy+Dxy+Ex^2……..; having the order of 3, 4; 4, 4 and/or 4, 3; Can you suggest me how does the augmented matrix looks like? what is the best approach to get PseudoInverse? what changes to the above code, I should be doing to address a surface fit on z; Any pointers would be much appreciated.

7. I do not think you even compiled it 🙂
a) gaussEliminationLS(int m, int n, double a[m][n], double x[n-1])
<== cannot have variable array sizez

b) Y[i]=Y[i]+pow(x[j],i)*y[j];<== 'y' is not even declared

c) indecis i,j are not declared anywhere …

1. These do work on certain compilers. For example if you use Dev C++ it will run guaranteed.

8. It does not even compile

a) cannot declare the variable array size like function (m,n, a[m][n] ..)
b) indecis i,j are never declared
c) “y” is not declared ( should be capital Y