数值逼近实验报告1共17页Word文档下载推荐.docx
《数值逼近实验报告1共17页Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《数值逼近实验报告1共17页Word文档下载推荐.docx(20页珍藏版)》请在冰豆网上搜索。
【实验环境】
(使用的软硬件)
软件:
MATLAB2012a
硬件:
电脑型号:
联想Lenovo昭阳E46A笔记本电脑
操作系统:
Windows8专业版
处理器:
Intel(R)Core(TM)i3CPUM350@2.27GHz2.27GHz
实验内容:
【实验方案设计】
第一步,将书上关于三种插值方法的内容转化成程序语言,用MATLAB实现;
第二步,分别用牛顿多项式插值,三次样条插值,拉格朗日插值求解不同的问题。
【实验过程】
(实验步骤、记录、数据、分析)
实验的主要步骤是:
首先分析问题,根据分析设计MATLAB程序,利用程序算出问题答案,分析所得答案结果,再得出最后结论。
实验一:
已知函数在下列各点的值为
xi
0.2
0.4
0.6
.0.8
1.0
f(xi)
0.98
0.92
0.81
0.64
0.38
试用4次牛顿插值多项式P4(x)及三次样条函数S(x)(自然边界条件)对数据进行插值。
用图给出{(xi,yi),xi=0.2+0.08i,i=0,1,11,10},P4(x)及S(x)。
(1)首先我们先求牛顿插值多项式,此处要用4次牛顿插值多项式处理数据。
已知n次牛顿插值多项式如下:
Pn=f(x0)+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+·
·
+f[x0,x1,·
xn](x-x0)·
(x-xn-1)
我们要知道牛顿插值多项式的系数,即均差表中得部分均差。
在MATLAB的Editor中输入程序代码,计算牛顿插值中多项式系数的程序如下:
functionvarargout=newtonliu(varargin)
clear,clc
x=[0.20.40.60.81.0];
fx=[0.980.920.810.640.38];
newtonchzh(x,fx);
functionnewtonchzh(x,fx)
%由此函数可得差分表
n=length(x);
fprintf('
*****************差分表*****************************\n'
);
FF=ones(n,n);
FF(:
1)=fx'
;
fori=2:
n
forj=i:
FF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1));
end
end
fori=1:
fprintf('
%4.2f'
x(i));
forj=1:
i
%10.5f'
FF(i,j));
\n'
由MATLAB计算得:
xi
f(xi)
一阶差商
二阶差商
三阶差商
四阶差商
0.20
0.980
0.40
0.920
-0.30000
0.60
0.810
-0.55000
-0.62500
0.80
0.640
-0.85000
-0.75000
-0.20833
1.00
0.380
-1.30000
-1.12500
-0.52083
所以有四次插值牛顿多项式为:
P4(x)=0.98-0.3(x-0.2)-0.62500(x-0.2)(x-0.4)-0.20833(x-0.2)(x-0.4)(x-0.6)-0.52083(x-0.2)(x-0.4)(x-0.6)(x-0.8)
(2)接下来我们求三次样条插值函数。
用三次样条插值函数由上题分析知,要求各点的M值:
三次样条插值函数计算的程序如下:
functiontgsanci(n,s,t)%n代表元素数,s,t代表端点的一阶导。
y=[0.980.920.810.640.38];
n=5
forj=1:
1:
n-1
h(j)=x(j+1)-x(j);
forj=2:
r(j)=h(j)/(h(j)+h(j-1));
u(j)=1-r(j);
f(j)=(y(j+1)-y(j))/h(j);
d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));
d
(1)=0
d(n)=0
a=zeros(n,n);
a(j,j)=2;
r
(1)=0;
u(n)=0;
a(j+1,j)=u(j+1);
a(j,j+1)=r(j);
b=inv(a)
m=b*d'
p=zeros(n-1,4);
%p矩阵为S(x)函数的系数矩阵
p(j,1)=m(j)/(6*h(j));
p(j,2)=m(j+1)/(6*h(j));
p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);
p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);
p
得到m=(0-1.6071-1.0714-3.10710)T
即M0=0;
M1=-1.6071;
M2=-1.0714;
M3=-3.1071;
M4=0
则根据三次样条函数定义,可得:
S(x)=
接着,在CommandWindow里输入画图的程序代码,
下面是画牛顿插值以及三次样条插值图形的程序:
plot(x,y)
holdon
5
y(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)
k=[011011]
x0=0.2+0.08*k
4
y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)
plot(x0,y0,'
o'
x0,y0)
holdon
y1=spline(x,y,x0)
plot(x0,y1,'
)
s=csape(x,y,'
variational'
fnplt(s,'
r'
gtext('
三次样条自然边界'
原图像'
4次牛顿插值'
运行上述程序可知:
给出的{(xi,yi),xi=0.2+0.08i,i=0,1,11,10}点,S(x)及P4(x)图形如下所示:
实验二:
在区间[-1,1]上分别取用两组等距节点对龙格函数作多项式插值及三次样条插值,对每个值,分别画出插值函数即的图形。
我们先求多项式插值:
在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令(如牛顿插值公式):
function[C,D]=newpoly(X,Y)
n=length(X);
D=zeros(n,n)
D(:
1)=Y'
fork=j:
D(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1));
C=D(n,n);
fork=(n-1):
-1:
1
C=conv(C,poly(X(k)))
m=length(C);
C(m)=C(m)+D(k,k);
当n=10时,我们在CommandWindow中输入以下的命令:
clear,clf,holdon;
X=-1:
0.2:
1;
Y=1./(1+25*X.^2);
[C,D]=newpoly(X,Y);
x=-1:
0.01:
y=polyval(C,x);
plot(x,y,X,Y,'
.'
gridon;
xp=-1:
z=1./(1+25*xp.^2);
plot(xp,z,'
得到插值函数和f(x)图形:
当n=20时,我们在CommandWindow中输入以下的命令:
0.1:
下面再求三次样条插值函数,在MATLAB的Editor中建立一个多项式的M-file,
输入下列程序代码:
functionS=csfit(X,Y,dx0,dxn)
N=length(X)-1;
H=diff(X);
D=diff(Y)./H;
A=H(2:
N-1);
B=2*(H(1:
N-1)+H(2:
N));
C=H(2:
N);
U=6*diff(D);
B
(1)=B
(1)-H
(1)/2;
U
(1)=U
(1)-3*(D
(1));
B(N-1)=B(N-1)-H(N)/2;
U(N-1)=U(N-1)-3*(-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);
end
M(N)=U(N-1)/B(N-1);
fork=N-2:
M(k+1)=(U(k)-C(k)*M(k+2))/B(k);
M
(1)=3*(D
(1)-dx0)/H
(1)-M
(2)/2;
M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;
fork=0: