[求助]哪位大侠知道曲线拟合的源程序？？
我想知道怎么样拟合一个曲线，就是给一组数据（xi，yi）求出直线、多项式、乘幂等曲线的方程式。还望各位大侠多多指教。 找本数值方法（机械工业出版社的浙大出的就可以） 这是几年前写的一个东东，应该能满足你的需要，稍做修改还能实现更强大的函数曲线的绘制。; DRAW.LSP by XBBCAD
；运行环境：AutoCAD R12或以上版本
；程序一：绘制一元函数曲线
(defun C:DW(/ x x1 x2 n y p1 p2 xdis b osm)
(setq s (getstring T"输入一元函数表达式(以x为X坐标自变量):"))
(if (/= s "") (setq dwexp s))
(princ (strcat "现在的表达式是-> " dwexp "\n"))
(setq x1 (getreal"X起始坐标="))
(setq x2 (getreal"X终点坐标="))
(setq n (getint"区间数量n="))
(setq x x1)
(setq xdis (/ (- x2 x1) n)) ;计算步长
(setq p1 (list x (cal dwexp))) ;起点坐标
(setq b (ssadd)) ;初始化选择集
(setq osm (getvar "OSMODE"));保存OSNAP模式
(setvar "OSMODE" 0); 关闭OSNAP模式
(repeat n
(progn
(setq x (+ x xdis)) ;x坐标增量
(setq y (cal dwexp)) ;计算y坐标
(setq p2 (list x y)) ;下一点
(command "PLINE" p1 p2 "") ;绘制PLINE线
(setq p1 p2) ;重设起点
(ssadd (entlast) b) ;加入选择集
);end progn
);end repeat
(command "PEDIT" (entlast) "j" b "" "") ;连接成一条PLINE线
(setvar "OSMODE" osm);恢复OSNAP模式
(princ "\nDraw end.\n")
(princ)
);end DW.
；程序二：绘制极坐标函数曲线
(defun C:PW(/ x x1 x2 n p1 p2 xdis len b osm)
(setq s (getstring T "输入极坐标函数式(以x为极角自变量): "))
(if (/= s "") (setq polarexp s))
(princ (strcat "现在的表达式是-> " polarexp "\n"))
(setq x1 (getreal"起始角="))
(setq x2 (getreal"终止角="))
(setq n (getint"区间数量n="))
(setq basepoint (getpoint"输入基点<0,0> :"))
(if (null basepoint) (setq basepoint '(0 0)))
(setq x x1 len 0.0)
(setq xdis (/ (- x2 x1) n)) ;计算极角的增量
(setq len (cal polarexp))
(setq p1 (polar basepoint x len));起点
(setq b (ssadd)) ;初始化选择集
(setq osm (getvar "OSMODE")) ;保存OSNAP模式
(setvar "OSMODE" 0) ;关闭OSNAP模式
(repeat n
(progn
(setq x (+ x xdis)) ;极角增量
(setq len (cal polarexp)) ;计算极径
(setq p2 (polar basepoint x len)) ;下一点坐标
(command "PLINE" p1 p2 "") ;绘制PLINE线
(setq p1 p2) ;重设起点
(ssadd (entlast) b) ;加入选择集
);end progn
);end repeat
(command "PEDIT" (entlast) "j" b "" "") ; 连接成一条PLINE线
(setvar "OSMODE" osm) ;恢复目标捕捉模式
(princ "\nDraw end.\n")
(princ)
);end PW. 我这有个用VC写的读ASCII文件，绘二次样条曲线的程序，看你需要吗？
void CTttView::OnDraw(CDC* pDC)
{
if(flag==1)
{
pDC->SetMapMode (MM_LOMETRIC);
pDC->SetWindowOrg (-150,1700);
pDC->MoveTo (0,0);
pDC->LineTo (0,1500);
pDC->LineTo (-35,1455);
pDC->MoveTo (0,1500);
pDC->LineTo (35,1455);
pDC->MoveTo (0,0);
pDC->LineTo (2000,0);//DRAW X AXIS
pDC->LineTo (1955,35);
pDC->MoveTo (2000,0);
pDC->LineTo (1955,-35);//DRAW X AXIS
pDC->TextOut (-50,1450,"Y");
pDC->TextOut (-50,-50,"(0,0)");
pDC->TextOut (2000,-50,"X");
pDC->Ellipse (x-10,y-10,x+10,y+10);
for(int i=1;i<10;i++)
{
pDC->Ellipse (x-10,y-10,x+10,y+10);
}
CSe pse;
double step=10;
double sx=x-step;
double sy=pse.calstart (x,y,x,y,step);
double ex=x+step;
double ey=pse.calend (x,y,x,y,step);
CCal pcal;
double tstep;
double tx,ty;
i=0;
for(tstep=0.0;tstep<=1.0;tstep=tstep+0.1)
{
pcal.input (sx,sy,x,y,x,y,x,y,tstep);
tx=pcal.calpx ();
ty=pcal.calpy ();
i++;
}
pDC->MoveTo (x,y);
for(i=1;i<=10;i++)
{
if((tx-floor(tx))>=0.5)
{
tx=ceil(tx);
}
else
{
tx=floor(tx);
}
if((ty-floor(ty))>=0.5)
{
ty=ceil(ty);
}
else
{
ty=floor(ty);
}
pDC->LineTo ((int)(tx),(int)(ty));
}
for(int j=1;j<=7;j++)
{
i=0;
for(tstep=0.0;tstep<1.0;tstep=tstep+0.1)
{
pcal.input (x,y,x,y,x,y,x,y,tstep);
tx=pcal.calpx ();
ty=pcal.calpy ();
i++;
}
pDC->MoveTo (x,y);
for(i=1;i<=10;i++)
{
if((tx-floor(tx))>=0.5)
{
tx=ceil(tx);
}
else
{
tx=floor(tx);
}
if((ty-floor(ty))>=0.5)
{
ty=ceil(ty);
}
else
{
ty=floor(ty);
}
pDC->LineTo ((int)(tx),(int)(ty));
}
}
i=0;
for(tstep=0.0;tstep<1.0;tstep=tstep+0.1)
{
pcal.input (x,y,x,y,x,y,ex,ey,tstep);
tx=pcal.calpx ();
ty=pcal.calpy ();
i++;
}
pDC->MoveTo (x,y);
for(i=1;i<=10;i++)
{
if((tx-floor(tx))>=0.5)
{
tx=ceil(tx);
}
else
{
tx=floor(tx);
}
if((ty-floor(ty))>=0.5)
{
ty=ceil(ty);
}
else
{
ty=floor(ty);
}
pDC->LineTo ((int)(tx),(int)(ty));
}
}
}
///////////////////////////
void CTttView::OnFileOpen()
{
CFileDialog input(TRUE);
CString filename;
FILE *fp;
float dx,dy;
if(input.DoModal ()!= IDOK)
return;
filename=input.GetPathName ();
if((fp=fopen(filename,"r++"))==NULL)
{
MessageBox("Failed");
return;
}
for(int i=0;i<=9;i++)
{
fscanf(fp,"%f,%f;",&dx,&dy);
}
//////////////////////////////////////////////////////////////////////////////
for(i=0;i<10;i++)
{
if((dx-floor(dx))>=0.5)
{
x=(int)(150*dx+1);
}
else
{
x=(int)(150*dx-1);
}
if((dy-floor(dy))>=0.5)
{
y=(int)(10*dy+1);
}
else
{
y=(int)(10*dy-1);
}
}
flag=1;
Invalidate(TRUE);
}
页:
[1]