求高斯函数拟合的算法

[qujc]

我的程序中需要对采集到的数据进行高斯拟合，但找不到相关的资料

[sharkhuang]

高斯？1+到100？

[freshairfly]

《数值计算》里就有吧
/*利用高斯主元消元法求解法方程*/
int Gauess(double a[][MAX], double b[])
{
int i, j, k;

double s[MAX];
int p[MAX];

double fTmp;
int pi, pr, pk, pn, iTmp;

int nYuan = n + 1;

for(i = 0; i < nYuan; i++)
{
for(j = 0; j < nYuan; j++)
{
fTmp = (float)fabs(a[i][j]);
if(j == 0)
s[i] = fTmp;
else if(s[i] < fTmp)
s[i] = fTmp;
}

p[i] = i;
}


for(k = 0; k < nYuan - 1; k++)
{
/*---------------求最小的r----------------*/
int r = k;
for(i = k; i < nYuan; i++)
{
int pr = p[r];
int pi = p[i];

if(fabs(a[pr][k])/s[pr] < fabs(a[pi][k])/s[pi])
{
r = i;
}
}
/*---------------求最小的r----------------*/

/*****************如果是奇异方程，则要停止*****************/
pr = p[r];
if(a[pr][k] == 0)
{
printf("奇异方程组!\n");
return -1;
}

/*---------------方程交换---------------*/
if(k != r)
{
iTmp = p[r];
p[r] = p[k];
p[k] = iTmp;

fTmp = s[r];
s[r] = s[k];
s[k] = fTmp;
}

/*^^^^^^^^^^^^^^^^^^执行Gauss消元^^^^^^^^^^^^^^^^^^^*/
for(i = k + 1; i < nYuan; i++)
{
pi = p[i];
pk = p[k];
a[pi][k] = -a[pi][k]/a[pk][k];
for(j = k + 1; j < nYuan; j++)
{
a[pi][j] = a[pi][j] + a[pi][k]*a[pk][j];
}
b[pi] = b[pi] + a[pi][k]*b[pk];
}
}

/**********************计算结果*********************/
pn = p[nYuan-1];
c[nYuan-1] = b[pn]/a[pn][nYuan-1];

for(i = nYuan - 2; i >= 0; i--)
{
pi = p[i];
fTmp = 0.0;
for(j = i + 1; j < nYuan; j++)
{
fTmp += a[pi][j]*c[j];
}
c[i] = (b[pi] - fTmp)/a[pi][i];
}

return 0;
}

[freshairfly]

不好意思，没看清楼主的要的东西
我的是高斯主元消元法求解线性方程

呵呵，楼主看看Matlab里有没有，如果有可以拿过来用的

[qujc]

1+到100?
什么意思？

[aa1298]

1加到100吗

[libbyliugang]

天那，难道高斯一辈子就只发明了1+...100的简单算发吗？？

[langtuteng]

能不能说清楚点啊

[wangpingpaul]

????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????

[ayingshmily]

一定用高斯吗？我这里倒是有拉格朗日和牛顿的算法。
牛顿的算法：给出一些离散的点（知道这些点的坐标），求出拟合的函数。我这里用sin(x)做的例子，画出了拟合函数的图像和准确函数的图像来做比较，比较直观。
学数值分析的时候代码写的不够规范，撮合着读吧。
#include "stdio.h"
#include "graphics.h"
#include"math.h"
#define maxsize 50
#define e 13
float newton(int x[],float y[],int n,int k)
{ float a[maxsize][maxsize];
float N=0,l=1;
int i,j;
for(j=0;j<n+1;j++)
a[0][j]=y[j];
for(i=1;i<n;i++)
for(j=i;j<n;j++)
a[i][j]=(a[i-1][j-1]-a[i-1][j])/(x[j-i]-x[j]);
for(i=0;i<n;i++)
{N=N+l*a[i][i];
l=l*(k-x[i]);
}
return(N);
}

main()
{ int a[maxsize];float b[maxsize];
int gdriver=DETECT, gmode;
float i;
registerbgidriver(EGAVGA_driver);
initgraph(&gdriver, &gmode, "");
/* setbkcolor(8);*/
setcolor(LIGHTRED);
line(50,220,560,220);
line(100,50,100,400);
outtextxy(90,222,"o");
outtextxy(557,217,">");
outtextxy(97,50,"^");
outtextxy(560,224,"x");
outtextxy(107,50,"y");
setcolor(LIGHTGREEN);
outtextxy(170,30,"y=sin(x)");
i=100;
while(i!=540)
{ putpixel(i,(sin((i-100)/20))*120+220,LIGHTGREEN);
i=i+0.5;
line(i-1,(sin((i-1-100)/20))*120+220,i,(sin((i-100)/20))*120+220);
}
setcolor(YELLOW);
outtextxy(170,50,"Nn=y");
a[0]=100;b[0]=sin((a[0]-100)/20)*120+220;
for(i=0;i<e;i++)
{a[i+1]=a[i]+40;
b[i+1]=sin((a[i+1]-100)/20)*120+220;
}
i=100;
while(i!=540)
{putpixel(i,newton(a,b,e,i),YELLOW);
i=i+0.5;
line(i-1,newton(a,b,e,i-1),i,newton(a,b,e,i));
}
setcolor(LIGHTRED);
outtextxy(129,217,"|");outtextxy(154,225,"PI");
outtextxy(158,217,"|");
outtextxy(190,217,"|");
outtextxy(222,217,"|");outtextxy(220,225,"2PI");
outtextxy(250,217,"|");
outtextxy(285,217,"|");outtextxy(283,225,"3PI");
outtextxy(315,217,"|");
outtextxy(348,217,"|");outtextxy(347,225,"4PI");
outtextxy(378,217,"|");
outtextxy(410,217,"|");outtextxy(406,225,"5PI");
outtextxy(440,217,"|");
outtextxy(473,217,"|");outtextxy(470,225,"6PI");
outtextxy(503,217,"|");
outtextxy(535,217,"|");outtextxy(530,225,"7PI");
outtextxy(100,156,"-");outtextxy(72,156,"0.5");
outtextxy(100,97,"-"); outtextxy(85,97,"1");
getch();
closegraph();
}