龙格函数.doc

对龙格函数在区间[-1,1]上取n=10,20等距节点，分别作多项式插值、三次样条插值比较各结果。
1、 多项式插值：
在区间[-1,1]上取n=10,20等距节点，带入拉格朗日插值多项式中，求出各个节点的插值，并利用matlab软件建立m函数，画出其图形。
在matlab中建立一个lagrange.m文件，里面代码如下：
%lagrange 函数
function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=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插值
n=10
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')
%插值曲线
hold on
%原曲线
plot(x0,y1,'-b')

取n=20
x=[-1:0.1: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')
%插值曲线
hold on
%原曲线
plot(x0,y1,'-b')

2三次样条插值：
所谓三次样条插值多项式是一种分段函数，它在节点分成的每个小区间上是3次多项式，其在此区间上的表达式如下：

因此，只要确定了的值，就确定了整个表达式，的计算方法如下：
令：

则满足如下个方程：

对于第一种边界条件下有

如果令那么解就可以为

求函数的二阶导数：

>> syms x
>> 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文件代码如下：

for x=-1:0.2:1
y=1/(1+25*x^2)
end
y =

得出
x=-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
y=0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385

编写m文件Three_Spline.m
n=10
x=linspace(-1,1,11);
y=1./(1+25*x.^2);
[m,p]=scyt1(x,y,0.0740,-0.0740);
hold on
x0=-1:0.01:1;
y0=1./(1+25*x0.^2);
plot(x0,y0,'--b')
得到如下图像
：

n=20
x=linspace(-1,1,21);
y=1./(1+25*x.^2);
[m,p]=scyt1(x,y,0.0740,-0.0740);
hold on
x0=-1:0.01:1;
y0=1./(1+25*x0.^2);
plot(x0,y0,'--b')
得到如下图像

附录
三次样条插值主要代码：
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));
for k=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;
for k=1:n-1
A(k,k+1)=r(k);
A(k+1,k)=u(k);
end
m=A\d;
ft=d(1);
syms t
for k=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);
for i=1:n-1 %出点xt对应的y值
for k=1:kmax
if x(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')
hold on
plot(x,y,'*')
hold off