曲线拟合的数值计算方法实验Word格式.docx
- 文档编号:4275863
- 上传时间:2023-05-03
- 格式:DOCX
- 页数:24
- 大小:172.07KB
曲线拟合的数值计算方法实验Word格式.docx
《曲线拟合的数值计算方法实验Word格式.docx》由会员分享,可在线阅读,更多相关《曲线拟合的数值计算方法实验Word格式.docx(24页珍藏版)》请在冰点文库上搜索。
3.幂函数拟合:
函数曲线为:
设
有N个点,其中横坐标是确定的。
最小二乘幂函数拟合曲线的系数A为:
、
4.对数函数拟合:
对数函数(lograrithmicfunction)的标准式形式为
b>
0时,Y随X增大而增大,先快后慢;
b<
0时,Y随X增大而减少,先快后慢,见图(c)、(d)。
当以Y和lnX绘制的散点图呈直线趋势时,可考虑采用对数函数描述Y与X之间的非线性关系,式中的b和a分别为斜率和截距。
更一般的对数函数
Y=a+bln(X+k)
式中k为一常量,往往未知。
(a)lnY=lna+bX
(b)lnY=lna-bX
(c)Y=a+blnX
(d)Y=a-blnX
5.线性插值:
在代数插值中,为了提高插值多项式对函数的逼近程度一般是增加节点的个数,即提高多项式的次数,但这样做往往不能达到预想的效果。
如下图所示:
f(x)=1/(1+x2)如果在区间[-5,5]上取7个等距节点:
xk=5*(k/3-1)(k=0,1,2,...,6),由lagrange插值公式可得到f(x)的次L7(x)。
如图所示:
L7(x)仅在区间的中部能较好的逼近函数f(x),在其它部位差异较大,而且越接近端点,逼近效果越差。
可以证明,当节点无限加密时,Ln(x)也只能在很小的范围内收敛,这一现象称为Runge现象。
它表明通过增加节点来提高逼近程度是不适宜的,因而不采用高次多项式插值。
如果我们把以上的点用直线连接起来,显然比L7(x)要更逼近f(x)。
这就是分段线性插值。
而在实际应用中通常采用分段低次插值来提高近似程度,比如可用分段线性插值或分段三次埃尔米特插值来逼近已知函数,但它们的总体光滑性较差。
为了克服这一缺点,一种全局化的分段插值方法——三次样条插值成为比较理想的工具。
6.三次样条插值:
设
有N+1个点,其中
。
如果存在N个三次多项式
,系数为
满足如下性质:
则成函数S(x)为三次样条函数。
7.端点约束:
紧压样条:
存在唯一的三次样条曲线,其一阶导数的边界条件是:
natural样条:
存在唯一的三次样条曲线,它的自由边界条件是:
外推样条:
存在唯一的三次样条曲线,其中通过对点x1和x2进行外推得到
,同时通过对点X(n-1)和X(N-2)进行外推得到
端点曲率调整:
存在唯一的三次样条曲线,其中二阶导数的边界条件
和
是确定的。
抛物线终结样条:
存在唯一的三次样条曲线,其中二阶在区间[X0,X1]内
,而在[Xn-1,Xn]内
三、实验内容
1.P2021
胡克定律指出F=kx,其中F是拉伸弹簧的拉力(单位为盎司),x为拉伸的长度(单位为英寸)。
根据下列试验数据,求解拉伸常量k的近似值。
xk
Fk
Xk
(a)(b)
2.P2151
洛杉矶(美国城市)郊区11月8日的温度记录入下表所示,其中共有24个数据点。
(a)根据例中的处理过程(使用fmins命令),对给定的数据求解最小二乘曲线
(b)求
(c)在同一坐标系下画出这些点集和(a)得出的最小二乘曲线。
时间,.
温度
1
66
58
2
3
65
4
64
5
63
57
6
7
62
8
61
9
60
10
11
59
67
午夜
正午
68
3.P2291
一个轿车在时间Tk时经过的距离dk,如下表所示。
使用程序,并根据一阶导数边界条件
求这些数据的三次紧压样条插值。
时间,tk
距离,dk
40
160
300
480
4.P2385
美国洛杉矶郊区11月9日的温度(华氏温度)如表所示。
采用24小时制。
(a)求三角多项式
(b)在同一坐标系下,画出图
和24个数据点。
(c)使用本地的温度情况重新求解问题(a)和问题(b)。
时间,
5.P2461
编写Matlab程序,生成并绘制组合贝塞尔曲线。
利用该程序生成和绘制过3个控制点集{(0,0),(1,2),(1,1),(3,0)},{(3,0),(4,-1),(5,-2),(6,1),(7,0)},{(7,0),(4,-3),(2,-1),(0,0)}的贝塞尔曲线。
四、实验结果及分析
1
实验描述:
由题意可知,此题需要用最小二乘法进行计算,因为已知函数的5个插点并且知道它们的x,y的值。
且函数的表达式为F=kx,所以只需用方程中
便可计算出k的数值。
定义一个动态数组
用来依次存取x和y的插值。
其中x,y的插值通过键盘手动输入并赋予给a中的元素。
然后通过相应的求和和基本运算便可以得到相应插值下的最小二乘法的函数表达式。
(其中k保留小数点后4位)
具体函数运行效果如下:
(a)请在此输入x的各插值
请在此输入y的各插值
此函数的最小二乘法曲线表达式为
Y=*x
请按任意键继续。
(b)请在此输入x的各插值
结果分析:
易得,最小二乘法多项式计算可以很好的做出较为精确的最小二乘法拟合曲线,并且程序通用性高,只要输入相应的插值便可以计算出结果,结果系数的小数点有效位同时也可以自己确定。
给出的最小二乘曲线表达式为:
其中变量有5个,而给出的数据点有24个,在C语言中可以用牛顿-拉夫森算法迭代计算分别得出5个变量的值,但是方法繁琐,且迭代计算量庞大,因此这里采用Matlab进行相关的计算,调用fminsearch函数,求得当5个参量都为1附近时候多项式的最小值,此时便可求出此5个参变量的值.然后继续通过Matlab,将得到的公式以及各点,用plot函数,便可以求得。
实验结果:
运行matlab结果如下:
>
fminsearch('
fiveOne'
[11111])
ans=
此时的所求值便为函数的待定系数。
所以可得最小二乘曲线的表达式为:
然后进行相应的绘制图形便可以求出所要求出的结果。
通过最小二乘法多项式同样可以绘制出三角函数的曲线。
并且程序通用性高,只要输入相应的插值便可以计算出结果,结果系数的小数点有效位同时也可以自己确定。
由题意可知,由于这里涉及到了样条线的运算,计算较为复杂。
且要涉及到画图的部分,所以此处采用matlab计算较为方便快捷。
而书本上给出了一个用来计算三次紧压样条曲线的可调用函数,现在将其引用,并根据已知点得出相应的曲线。
代入后得出的结果如下所示:
X=[02468];
Y=[040160300480];
S=csfit(X,Y,0,98)
S=
00
由结果可知此插值为在区间[0,8]中有三个分别划分了[0,2],[2,4],[4,6],[6,8]四个区间的插点。
且多项式分别为
在matlab中通过polyval作出相应的曲线,再用plot函数便可绘制图线。
绘制后图线如下:
此时便可以直观的看到一个平滑的样条线。
通过给定的一些数据点,便可以绘制出过这些点的相应的样条线,通过观察能发现样条线的平滑度与你选择的样条线类型以及数据点的分布有一定关联。
不仅仅是紧压样条线,相关其它的线也可以用同样的方法一一给出。
5
由题意可知,题目是想求出所绘制数据点的三角多项式的逼近,三角多项式逼近的公式为:
此外,x为离散傅里叶级数,满足条件:
而所给的点x为[1,24]的离散数值点,所以无法直接对其作出逼近公式,需要进行尺度变换,将x点转换为:
(1)通过matlab绘制出相关的三角多项式曲线,然后同样通过matlab的绘制点,将点绘制到这个曲线之中,具体的matlab代码如下:
holdon;
%保持图形不动,绘制新的图形入曲线中
X=[123456789101112131415161718192021222324];
%数据点的x的取值
Y=[585858585757575860646768666665646363626160605958];
%数据点的y的取值
plot(X,Y,’xk’)%绘制出数据点
(2)然后便可以画出如图所示的插值数据点。
三角多项式的曲线拟合度非常高,能很好的绘制出图像的具体形式而且曲线平滑,但是它需要满足x属于-pi到pi的区间内。
由题意可知,首先以3阶贝塞尔函数为带求函数。
其求解格式如下:
其中b为伯恩斯坦多项式,其定义如下:
用matlab先设置一个参数t,然后再根据公式,通过所给的数据点将伯恩斯坦多项式,以及x,y的关于参数t的多项式求出来。
然后再把t设置为精度足够大的单位序列,绘制图线即可得出贝塞尔的效果。
其matlab运行后结果如下:
【以第一组控制点为例】
X=[0113]
X=
0113
Y=[0210]
Y=
0210
fiveTwo(X,Y)
x=
3*t*(t-1)^2-3*t^2*(t-1)+3*t^3
y=
6*t*(t-1)^2-3*t^2*(t-1)
且绘制出第一组控制点所在位置以及三阶贝塞尔曲线如下所示:
同理,第2,3组控制点所作的图形如下所示:
可以通过Matlab函数绘制出三阶的贝塞尔函数。
只要给出控制点,便可自动绘制出所控制的三阶贝塞尔函数以及控制点的位标。
由此可以观察到贝塞尔与控制点的约束作用,以及所求得贝塞尔函数是个相对平滑的曲线。
五、实验结论
曲线拟合对于实际的工程以及理论推导问题都有着重大的作用。
在具体的实验,数据分析上,往往我们只有巨量的已知离散点,想从离散点中得出函数表达式则就需要曲线拟合进行,根据不同的要求,我们可以选择最小二乘法的幂函数拟合或者是一次函数,二次函数拟合,抑或是精度非常高的傅里叶变换的三角函数拟合。
同时在建模方面,贝塞尔函数的引用也从数学层面解决了如何用计算机绘制出光滑圆润的曲线,在一些设计软件中,例如3dmax和maya的三维建模,AutoCAD的工程建模,贝塞尔运算对于曲线曲面的设计有着举足轻重的作用。
附件(代码):
#include<
voidmain()
{
externfloatafunction(floatb[]);
externfloatbfunction(floatz,floatv);
float*a,y,q,c;
inti;
a=newfloat[5];
cos(i*B)+C.*sin(i*D)+E-y(i)).^2+z;
%函数的线性最小二乘法的多项式
end
之后在Matlab的命令窗口输入如下语句:
>
fminsearch(‘fiveOne’,[11111])
(绘制图形方程组)
x1=[0123456789101112131415161718192021222324];
%数值点的x取值
y1=[585858585757575860646768666665646363626160605958];
plot(x1,y1,'
x'
)%绘制出24个数值点的图形
x=linspace(0,25,100);
holdon
plot(x,y,'
k'
)%绘制出此函数的二乘法后的函数曲线
(三次紧压样条线计算函数)
functionS=csfit(x,y,dx0,dxn)
N=length(x)-1;
H=diff(x);
D=diff(y)./H;
%d(k)的值
A=H(2:
N-1);
B=2*(H(1:
N-1)+H(2:
N));
C=H(2:
N);
U=6*diff(D);
%U(k)的值
B
(1)=B
(1)-H
(1)/2;
U
(1)=U
(1)-3*(D
(1)-dx0);
B(N-1)=B(N-1)-H(N)/2;
U(N-1)=U(N-1)-3*(dxn-D(N));
fork=2:
N-1
temp=A(k-1)/B(k-1);
B(k)=B(k)-temp*C(k-1);
U(k)=U(k)-temp*U(k-1);
M(N)=U(N-1)/B(N-1);
fork=N-2:
-1:
M(k+1)=(U(k)-C(k)*M(k+2))/B(k);
M
(1)=3*(D
(1)-dx0)/H
(1)-M
(2)/2;
%三次紧压约束中S"
(0)的值
M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;
(k)的值
fork=0:
S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));
%S(k,3)的值
S(k+1,2)=M(k+1)/2;
%S(k,2)的值
S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;
%S(k,1)的值
S(k+1,4)=y(k+1);
%S(k,0)的值
End
(曲线绘制函数程序)
x1=0:
.01:
2;
y1=polyval(S(1,:
),x1-X
(1));
%第一段样条线
x2=2:
4;
y2=polyval(S(2,:
),x2-X
(2));
%第二段样条线
x3=4:
6;
y3=polyval(S(3,:
),x3-X(3));
%第三段样条线
x4=6:
8;
y4=polyval(S(4,:
),x4-X(4));
%第四段样条线
plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,'
)%绘制连接成完整曲线,并且标注数据点的位置。
externfloatafunction(floatX[],floatY[],inttemp);
externfloatbfunction(floatX[],floatY[],inttemp);
floata[24]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24};
floatb[24]={58,58,58,58,57,57,57,58,60,64,67,68,66,66,65,64,63,63,62,61,60,60,59,58};
float*c;
float*d;
c=newfloat[7],d=newfloat[7];
for(i=0;
i<
=24;
i++)
{
c[i]=afunction(a,b,i);
d[i]=bfunction(a,b,i);
}
printf("
由所给的数据点可以得出\n三角多项式T7(x)中a的系数分别为\n"
);
%f%f%f%f%f%f%f\n"
c[0],c[1],c[2],c[3],c[4],c[5],c[6]);
由所给的数据点可以得出\n三角多项式T7(x)中b的系数分别为\n"
d[0],d[1],d[2],d[3],d[4],d[5],d[6]);
system("
pause"
return;
}
floatafunction(floatX[],floatY[],inttemp)
floatA=0;
A=Y[i]*cos(temp*X[i])/6+A;
returnA;
floatbfunction(floatX[],floatY[],inttemp)
A=Y[i]*sin(temp*X[i])/6+A;
function[]=fiveTwo()
aa=0;
bb=0;
y=0;
symsx
fori=1:
24
X(i)=+2.*i*24;
forj=1:
aa=Y(j).*cos(j.*X(i))+aa;
end
a(i)=12*aa;
%求出多项式aj;
bb=Y(j).*sin(j.*X(i))+aa;
b(i)=12*bb;
%求出多项式bj;
y=a(i)*cos(i.*x)+b(i)*sin(i.*x)+y;
%求出三角多项式逼近
x=:
:
;
plot(x,y,'
g'
%绘制图形
functionz=fiveTwo(X,Y)
symst%定义t为参数变量
a=0;
b=0;
fori=0:
1:
B(i+1)=nchoosek(3,i)*t^i*(1-t)^(3-i);
%计算伯恩斯坦多项式
a=X(i).*B(i)+a;
%计算三阶贝塞尔函数下x关于t的参数方程
b=Y(i).*B(i)+b;
t=0:
1;
x=a;
y=b;
plot(X,Y,'
)%在图形中绘制出已知的4个控制点
X%输出三阶贝塞尔函数下x关于t的参数方程
Y%输出三阶贝塞尔函数下x关于t的参数方程
为了绘制出最终的贝塞尔函数,需要在matlab命令面板中输入如下指令:
x=X(转换成元素的点乘除形式)
y=Y(转换成元素的点乘除形式)
holdon%保持在一张图内绘制
plot(x,y)%在图形中绘制出所求得贝塞尔函数曲线
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 曲线拟合 数值 计算方法 实验
![提示](https://static.bingdoc.com/images/bang_tan.gif)