请问谁有曲线拟合的原程序。(100分)

X

xjlqs

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-01

• #1

请问谁有曲线拟合的原程序，比如最小二乘法的。

 

J

joysun

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-02

• #2

#include <stdio.h>
#include <iostream.h>
#include <math.h>

//列主元消元法求解线性方程组
int gs1(int n,float a[],float b[],float x[])
{
int i,j,k,r;
int temp; //交换数据时的临时变量
double *s;
double ar; //a[P[r]][k]
int *P;
double ax;
s=new double[n];
P=new int[n];
for(i=0;i<n;i++)
P=i;
for(k=0;k<n-1;k++)
{
ar=fabs(a[P[k]*n+k]);
r=k;
for(i=k+1;i<n;i++)
if(fabs(a[P*n+k])>ar)
{
ar=fabs(a[i*n+k]);
r=i;
}
if(a[P[r]*n+k]==0)
return(-1); //消元错误 奇异
if(k!=r)
{
temp=P[k];
P[k]=P[r];
P[r]=temp;
}
for(i=k+1;i<n;i++)
{
a[P*n+k]=-a[P*n+k]/a[P[k]*n+k];
for(j=k+1;j<n;j++)
a[P*n+j]+=a[P*n+k]*a[P[k]*n+j];
b[P]+=a[P*n+k]*b[P[k]];
}
}
x[n-1]=b[P[n-1]]/a[P[n-1]*n+n-1];
for(i=n-2;i>=0;i--)
{
ax=0;
for(j=i+1;j<n;j++)
ax+=a[P*n+j]*x[j];
x=(b[P]-ax)/a[P*n+i];
}
delete[]P;
delete[]s;
return 0;
}

//求矩阵乘积
int mul(int m,int n,float A[],float B[],float C[],float D[])
{
int i,j,k;
float temp;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
temp=0;
for(k=0;k<m;k++)
temp+=A[k*n+i]*A[k*n+j];
C[i*n+j]=temp;
}
for(i=0;i<n;i++)
{
temp=0;
for(k=0;k<m;k++)
temp+=A[k*n+i]*B[k];
D=temp;
}
return 0;
}

void main()
{
float *amn,*bm;
FILE *fp;
int m,n;
int i,j,rtn;
float *a,*b;
float *x;
if((fp=fopen("data.txt","r"))==NULL)
{
cout<<"无法打开数据文件！"<<endl;
exit(0);
}
fscanf(fp,"%d%d",&amp;m,&amp;n);
a=new float[m*m]; //a[i*n+j] 表示a[j] a+i*n+j 表示地址
b=new float[n];
x=new float[n];
amn=new float[m*n];
bm=new float[m];
for(i=0;i<m;i++)
for(j=0;j<n;j++)
fscanf(fp,"%f",amn+i*n+j);
for(i=0;i<m;i++)
fscanf(fp,"%f",bm+i);
fclose(fp);
rtn=mul(m,n,amn,bm,a,b);
rtn=gs1(n,a,b,x);
if((fp=fopen("result.txt","w"))==NULL)
{
cout<<"无法打开输出数据文件"<<endl;
exit(1);
}
switch(rtn)
{
case 0:
fprintf(fp,"方程组的解为/n",n);
for(i=0;i<n;i++)
fprintf(fp,"x%d=%g/t",i+1,x);
break;
case -1:
fprintf(fp,"列主元消元法失败");
}
fclose(fp);
delete[]a;
delete[]b;
delete[]x;
delete[]amn;
delete[]bm;
}

 

X

xjlqs

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-02

• #3

to joysun:请问你的程序是delphi吗？

 

X

xjlqs

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-02

• #4

to joysun:如果你能用delphi改写上面的程序我再加100分。

 

J

joysun

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-03

• #5

//上面是c++ builder的代码，已经经过测试！
//以下给出了delphi的代码，两个函数的调用仿照上面即可！

unit Unit2;

interface
uses math;
implementation
//列主元消元法求解线性方程组
function gs1(n : integer;var a :array of real ;var b:array of real;var x:array of real):integer;
var i,j,k,r : integer;
temp : integer;
s : array of real;
ar : real;
P : array of integer;
ax : real;
begin
SetLength(s,n);
SetLength(P,n);

for i := 0 to n do
P := i;
for k := 0 to n do
begin
ar := abs(a[P[k]*n+k]);
r := k;
for i := k+1 to n do
begin
if (abs(a[P*n+k])>ar) then
begin
ar := abs(a[P*n+k]);
r := i;
end;
end;
if (a[P[r]*n+k]=0) then result := -1;
if (k <> r) then
begin
temp := P[k];
P[k] := P[r];
P[r] := temp;
end;
for i := k+1 to n do
begin
a[P*n+k] := - a[P*n+k]/a[P[k]*n+k];
for j := k+1 to n do
a[P*n+j] := a[P*n+k] * a[P[k]*n+j] + a[P*n+j];
b[P] := a[P*n+k] * b[P[k]] + b[P];
end;
end;
x[n-1] :=b[P[n-1]]/a[P[n-1]*n+n-1];

for i := n-2 to 0 do
begin
ax := 0;
for j := i+1 to n do
ax := ax + a[P*n+j]*x[j];
x := (b[P]-ax)/a[P*n+i];
end;
end;
//求矩阵乘积
function mul(m : integer;n : integer;var A: array of real;var B: array of real;var C:array of real;var D:array of real):integer;
var i,j,k : integer;
temp : real;
begin
for i := 0 to n do
for j := 0 to n do
begin
temp := 0;
for k := 0 to m do
temp := temp + A[k*n+i]*A[k*n+j];
C[i*n+j] := temp;
end;
for i := 0 to n do
begin
temp := 0;
for k := 0 to n do
temp := temp + A[k*n+i]*B[k];
D := temp;
end;
result := 0;
end;
end.

 

X

xjlqs

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-03

• #6

假设曲线为Y=1+x^2 我的主程序为：
procedure TForm1.FormCreate(Sender: TObject);
var
m,n,i,j,rtn:integer;
amn: array of real;
bm:array of real;
b,x:array of real;
a: array of real;
begin
m:=6;
n:=2;
setlength(a,(m+1)*(m+1)); //a[i*n+j] 表示a[j] a+i*n+j 表示地址
setlength(b,n+1);
setlength(x,n+1);
setlength(bm,m+1);
setlength(amn,(m+1)*(n+1));
for i:=0 to m do
for j:=0 to n do
amn[i*n+j]:=i;
for i:=0 to m do
bm:=i*i+1;
rtn:=mul(m,n,amn,bm,a,b);
rtn:=gs1(n,a,b,x);
// showmessage(inttostr(rtn));
memo1.clear;
case rtn of
0: begin
memo1.lines.add('方程组的解为/n:'+inttostr);
for i:=0 to n do
memo1.lines.add(inttostr(i)+': '+formatfloat('0.00',x));
end;
-1:memo1.lines.add('列主元消元法失败');
end;




end;

结果为：
方程组的解为/n:2
0: -0.06
1: 0.19
2: 0.00

好象不正确，我的数学有点抱歉，还请多多指教

 

J

joysun

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-03

• #7

对于你的问题需要重新写mul（）函数！稍后我给你贴出来！

 

X

xjlqs

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-03

• #8

to joysun:我主要是不知道你的输入数据的内容和排列方式

 

X

xjlqs

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-01-04

• #9

to joysun:问题我自己解决了，不过分我还会给你，请到
http://www.delphibbs.com/delphibbs/dispq.asp?lid=1503171
处发言，我再给你100分

 

M

miniking

Unregistered / Unconfirmed

GUEST, unregistred user!

• 2003-06-27

• #10

joysum先生你好，我现在遇到了问题，你可不可以把用最小二乘法的那段代码（正确的那段）给我一份呢？我现在也遇到了这个问题