三次样条拟合范例.docx
- 文档编号:28143448
- 上传时间:2023-07-08
- 格式:DOCX
- 页数:12
- 大小:102.63KB
三次样条拟合范例.docx
《三次样条拟合范例.docx》由会员分享,可在线阅读,更多相关《三次样条拟合范例.docx(12页珍藏版)》请在冰豆网上搜索。
三次样条拟合范例
三次样条拟合范例
1设计目的、要求
对龙格函数
在区间[-1,1]上取
的等距节点,分别作多项式插值、三次样条插值和三次曲线拟合,画出
及各逼近函数的图形,比较各结果。
2设计原理
(1)多项式插值:
利用拉格朗日多项式插值的方法,其主要原理是拉格朗日多项式,即:
表示待插值函数的
个节点,
,其中
;
(2)三次样条插值:
三次样条插值有三种方法,在本例中,我们选择第一边界条件下的样条插值,即两端一阶导数已知的插值方法:
(3)三次曲线拟合:
本题中采用最小二乘法的三次多项式拟合。
最小二乘拟合是利用已知的数据得出一条直线或者曲线,使之在坐标系上与已知数据之间的距离的平方和最小。
在本题中,n=10,故有11个点,以这11个点的
和
值为已知数据,进行三次多项式拟合,设该多项式为
,该拟合曲线只需
的值最小即可。
3采用软件、设备
计算机、matlab软件
4设计内容
1、多项式插值:
在区间
上取
的等距节点,带入拉格朗日插值多项式中,求出各个节点的插值,并利用matlab软件建立m函数,画出其图形。
在matlab中建立一个lagrange.m文件,里面代码如下:
%lagrange函数
functiony=lagrange(x0,y0,x)
n=length(x0);m=length(x);
fori=1:
m
z=x(i);
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
建立一个polynomial.m文件,用于多项式插值的实现,代码如下:
%lagrange插值
x=[-1:
0.2:
1];
y=1./(1+25*x.^2);
x0=[-1:
0.02:
1];
y0=lagrange(x,y,x0);
y1=1./(1+25*x0.^2);
plot(x0,y0,'--r')
%插值曲线
holdon
%原曲线
plot(x0,y1,'-b')
运行duoxiangshi.m文件,得到如下图形:
2、三次样条插值:
所谓三次样条插值多项式
是一种分段函数,它在节点
分成的每个小区间
上是3次多项式,其在此区间上的表达式如下:
因此,只要确定了
的值,就确定了整个表达式,
的计算方法如下:
令:
则
满足如下
个方程:
对于第一种边界条件下有
如果令
那么解就可以为
求函数的二阶导数:
>>symsx
>>f=sym(1/(1+25*x^2))
f=
1/(1+25*x^2)
>>diff(f)
ans=
-(50*x)/(25*x^2+1)^2
将函数的两个端点,代入上面的式子中:
f’(-1)=0.0740
f’
(1)=-0.0740
求出从-1到1的n=10的等距节点,对应的x,y值
对应m文件代码如下:
forx=-1:
0.2:
1
y=1/(1+25*x^2)
end
y=
得出
x=-1-0.8-0.6-0.4-0.200.20.40.60.81
y=0.03850.05880.10.20.510.50.20.10.05880.0385
编写m文件Three_Spline.m
x=linspace(-1,1,11);
y=1./(1+25*x.^2);
[m,p]=scyt1(x,y,0.0740,-0.0740);
holdon
x0=-1:
0.01:
1;
y0=1./(1+25*x0.^2);
plot(x0,y0,'--b')
得到如下图像:
.
其中蓝色曲线为原图,红色曲线为拟合后的图像。
3、三次曲线拟合:
这里我们使用最小二乘法的3次拟合
建立一个Three_fitting.m文件,代码如下:
%主要代码
x=[-1-0.8-0.6-0.4-0.200.20.40.60.81];
y=[0.03850.05880.10.20.510.50.20.10.05880.0385];
a=polyfit(x,y,3);
x1=[-1:
0.01:
1];
y1=a(4)+a(3)*x1+a
(2)*x1.^2+a
(1)*x1.^3;
x0=[-1:
0.01:
1];
y0=1./(1+25*x0.^2)
%原曲线
plot(x0,y0,'-r')
holdon
%三次拟合曲线
plot(x1,y1,'-b')
上图中,蓝色部分为三次拟合曲线,红色部分为原曲线
6结果分析
拉格朗日插值的优点是对于某一区域,不限于被估计点周围,公式简单易实施。
一般认为n的次数越高,逼近
的精度就越好,但在本题中,对龙格函数
,中间部分插值效果比较好,而对于两端,插值结果是非常不理想的,即龙格现象。
样条函数可以给出光滑的插值曲线,从本题中就能体现出来。
从以上图形可以看出,三次样条插值的图形是比较逼近于原图的,收敛性相对而言是非常好的,但在本题中,仅将原区间分成10个等距区间,因此,逼近效果还不是特别理想,当我们将n增大时,插值后的曲线越逼近于原曲线。
总的来说,三次样条插值的稳定性比较好,收敛性比较强。
在这三种方法中,三次曲线拟合的效果是最差的,所得的图形与原曲线差距甚远。
最小二乘法中,并不要求拟合后的曲线经过所有已知的点,只需要拟合多项式上的点在某种标准上与定点之间的差距最小即可,因此与原曲线的逼近程度是最差的。
最小二乘法的多项式拟合只适用于多项式,而本题中的函数并不是一个多项式,因此,不建议使用最小二乘法拟合。
参考文献:
[1]李庆扬王能超等.数值分析[M].清华大学出版社
[2]吴振远.科学计算实验指导书基于MATLAB数值分析[M].中国地质大学出版社
[3]宋叶志.MATLAB数值分析与应用[M].机械工业出版社,2009.07
附录
三次样条插值主要代码:
function[m,p]=scyt1(x,y,df0,dfn)
n=length(x);
r=ones(n-1,1);
u=ones(n-1,1);
d=ones(n,1);
r
(1)=1;
d
(1)=6*((y
(2)-y
(1))/(x
(2)-x
(1))-df0)/(x
(2)-x
(1));
u(n-1)=1;
d(n)=6*(dfn-(y(n)-y(n-1))/(x(n)-x(n-1)))/(x(n)-x(n-1));
fork=2:
n-1
u(k-1)=(x(k)-x(k-1))/(x(k+1)-x(k-1));r(k)=(x(k+1)-x(k))/(x(k+1)-x(k-1));
d(k)=6*((y(k+1)-y(k))/(x(k+1)-x(k))-(y(k)-y(k-1))/(x(k)-x(k-1)))/(x(k+1)-x(k-1));
end
A=eye(n,n)*2;
fork=1:
n-1
A(k,k+1)=r(k);
A(k+1,k)=u(k);
end
m=A\d;
ft=d
(1);
symst
fork=1:
n-1%求s(x)即插值多项式
p(k,1)=m(k)/(6*(x(k+1)-x(k)));
p(k,2)=m(k+1)/(6*(x(k+1)-x(k)));
p(k,3)=(y(k)-m(k)*(x(k+1)-x(k))^2/6)/(x(k+1)-x(k));
p(k,4)=(y(k+1)-m(k+1)*(x(k+1)-x(k))^2/6)/(x(k+1)-x(k));
sx(k)=p(k,1)*(x(k+1)-t)^3+p(k,2)*(t-x(k))^3+p(k,3)*(x(k+1)-t)+p(k,4)*(t-x(k));
end
kmax=1000;
xt=linspace(x
(1),x(n),kmax);
fori=1:
n-1%出点xt对应的y值
fork=1:
kmax
ifx(i)<=xt(k)&xt(k)<=x(i+1)
fx(k)=subs(sx(i),xt(k));
end
end
end
plot(xt,fx,'r');xlabel('x');
ylabel('y');title('f');
text(x(fix(n/2)),y(fix(n/2)),'f')
holdon
plot(x,y,'*')
holdoff
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 三次 拟合 范例