想用C语言编写多项式拟合的程序

#include <stdio.h>

#include <conio.h>

#include <stdlib.h>

#include <math.h>

main()

{

int i,j,m,n=7,poly_n=2

double x[7]={1,2,3,4,6,7,8},y[7]={2,3,6,7,5,3,2}

double a[3]

void polyfit(int n,double *x,double *y,int poly_n,double a[])

system("cls")

polyfit(n,x,y,poly_n,a)

for (i=0i<poly_n+1i++)/*这里是升序排列，Matlab是降序排列*/

printf("a[%d]=%g\n",i,a[i])

getch()

}

/*==================polyfit(n,x,y,poly_n,a)===================*/

/*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/

/*=====n是数据个数 xy是数据值 poly_n是多项式的项数======*/

/*===返回a0,a1,a2,……a[poly_n]，系数比项数多一（常数项）=====*/

void polyfit(int n,double x[],double y[],int poly_n,double a[])

{

int i,j

double *tempx,*tempy,*sumxx,*sumxy,*ata

void gauss_solve(int n,double A[],double x[],double b[])

tempx=calloc(n,sizeof(double))

sumxx=calloc(poly_n*2+1,sizeof(double))

tempy=calloc(n,sizeof(double))

sumxy=calloc(poly_n+1,sizeof(double))

ata=calloc((poly_n+1)*(poly_n+1),sizeof(double))

for (i=0i<ni++)

{

tempx[i]=1

tempy[i]=y[i]

}

for (i=0i<2*poly_n+1i++)

for (sumxx[i]=0,j=0j<nj++)

{

sumxx[i]+=tempx[j]

tempx[j]*=x[j]

}

for (i=0i<poly_n+1i++)

for (sumxy[i]=0,j=0j<nj++)

{

sumxy[i]+=tempy[j]

tempy[j]*=x[j]

}

for (i=0i<poly_n+1i++)

for (j=0j<poly_n+1j++)

ata[i*(poly_n+1)+j]=sumxx[i+j]

gauss_solve(poly_n+1,ata,a,sumxy)

free(tempx)

free(sumxx)

free(tempy)

free(sumxy)

free(ata)

}

void gauss_solve(int n,double A[],double x[],double b[])

{

int i,j,k,r

double max

for (k=0k<n-1k++)

{

max=fabs(A[k*n+k])/*find maxmum*/

r=k

for (i=k+1i<n-1i++)

if (max<fabs(A[i*n+i]))

{

max=fabs(A[i*n+i])

r=i

}

if (r!=k)

for (i=0i<ni++) /*change array:A[k]&A[r] */

{

max=A[k*n+i]

A[k*n+i]=A[r*n+i]

A[r*n+i]=max

}

max=b[k] /*change array:b[k]&b[r] */

b[k]=b[r]

b[r]=max

for (i=k+1i<ni++)

{

for (j=k+1j<nj++)

A[i*n+j]-=A[i*n+k]*A[k*n+j]/A[k*n+k]

b[i]-=A[i*n+k]*b[k]/A[k*n+k]

}

}

定义一个myfun.m的m文件，保存：

function f = myfun(par)

data = [20 0.6813

37 0.6027

54 0.4236

71 0.2771

105 0.199

122 0.1665

139 0.1354

156 0.1172

173 0.0973

190 0.0772

207 0.061

224 0.0392

241 0.0364

258 0.0235

275 0.0106]

x = data(:, 1)

y = data(:, 2)

for i = 1:length(x),

wow = par(1)*x(i)

yfit(i,1) = 1 - exp(wow^-1)

end

temp = y - yfit

f = temp'*temp

temp = y - yfit

f = temp'*temp

主窗口里面运行：

options = optimset('MaxFunEvals', 1e6, 'Maxiter', 1e6)

[par, fval] = fminsearch(@myfun, -0.05)

验证：

x = data(:, 1)

y = data(:, 2)

for i = 1:length(x),

wow = par(1)*x(i)

yfit(i,1) = 1 - exp(wow^-1)

end

plot(yfit)hold onplot(y, 'r.')

add:

我之前写的答案不是很好，我建议你固定b = -1, -2, ...，然后比较每次出来的结果，找最优解。因为这个exp((-a*x)^b)函数很容易就变成复数域了，很麻烦.看你的x, y取值，定b为负整数，这是很显然的事情。。。

【代码】

x=[200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000]

y=[0.1 0.25 0.49 0.65 0.7 0.91 1.15 1.26 1.37 1.46 1.52 1.60 1.65 1.67 1.68 1.68 1.69 1.69 1.71]

cftool(x,y)

【拟合方式一：指数拟合】

General model Power2:

     f(x) = a*x^b+c

Coefficients (with 95% confidence bounds):

       a =      -44.95  (-570, 480.1)

       b =    -0.02049  (-0.297, 0.2561)

       c =        40.3  (-490, 570.6)

Goodness of fit:

  SSE: 0.1527

  R-square: 0.9708

  Adjusted R-square: 0.9672

  RMSE: 0.0977

【拟合方式二：最高三次多项式】

Linear model Poly3:

     f(x) = p1*x^3 + p2*x^2 + p3*x + p4

Coefficients (with 95% confidence bounds):

       p1 = -1.208e-011  (-1.778e-010, 1.536e-010)

       p2 = -6.613e-007  (-1.214e-006, -1.088e-007)

       p3 =    0.002397  (0.001855, 0.00294)

       p4 =      -0.388  (-0.5376, -0.2384)

Goodness of fit:

  SSE: 0.02784

  R-square: 0.9947

  Adjusted R-square: 0.9936

  RMSE: 0.04308

【拟合方式三：最高四次多项式】

Linear model Poly4:

     f(x) = p1*x^4 + p2*x^3 + p3*x^2 + p4*x + p5

Coefficients (with 95% confidence bounds):

       p1 =  4.099e-013  (1.256e-013, 6.942e-013)

       p2 = -1.815e-009  (-3.073e-009, -5.575e-010)

       p3 =  2.001e-006  (1.018e-007, 3.9e-006)

       p4 =   0.0009045  (-0.0002188, 0.002028)

       p5 =     -0.1391  (-0.3494, 0.07117)

Goodness of fit:

  SSE: 0.01654

  R-square: 0.9968

  Adjusted R-square: 0.9959

  RMSE: 0.03437

【小结】

实际上多项式拟合经度（SSE,RMSE）四次高于三次，因为四次包含了三次，需要根据物理模型来确定选择多项式最高次数。

指数模型精度低于三次多项式。

本文希望你能认识一个新的有用的函数，曲线拟合工具箱

希望你学习进步。

---