计算方法实验三
一、实验名称：
实验三、插值逼近
二、实验目的：
1． 掌握Lagrange插值、Newton插值的概念；
2． 编写程序实现；
3． 观察Runge现象。
三、实验内容及要求：
（a）取，以等分节点为插值节点，构造Lagrange插值公式，用程序计算各等分区间中点处的值，列表显示（给出插值多项式及函数在等分区间中点处值），并分别作出和插值函数的图形（将5个图形画在同一张图里，并用不同颜色表示，四个插值多项式及一
个图形）；
（b）取，以等分节点为插值节点，构造牛顿插值函数，用程序计算各等分区间中点处的值，列表显示（给出插值多项式及函数在等分区间中点处值），并分别作出和插值函数的图形（将5个图形画在同一张图里，并用不同颜色表示，四个插值多项式及一个图形）；
（c）取，在剖分的基础上简历分片线性Lagrange插值函数，用程序计算各等分区间中点处的值，列表显示（给出插值多项式及函数在等分区间中点处值），并分别作出和插值函数的图形（将5个图形画在同一张图里，并用不同颜色表示，四个插值多项式及一个图形）；
2.实现课本39页例2.12，取，分别作出和插值函数的图形（将5个图形画在同一张图里，并用不同颜色表示，四个插值多项式及一个图形）；观察Runge现象。
三、实验步骤（或记录）
1、（a）
插值多项式及函数在等分区间中点处值
Y 0.038461538 0.042440318 0.047058824 0.052459016 0.058823529 0.066390041 0.075471698 0.086486486 0.1 0.116788321 0.137931034 0.164948454 0.2 0.246153846 0.307692308 0.390243902 0.5 0.64 0.8 0.941176471 1 0.941176471 0.8 0.64 0.5 0.390243902 0.307692308 0.246153846 0.2 0.164948454 0.137931034 0.116788321 0.1 0.086486486 0.075471698 0.066390041 0.058823529 0.052459016 0.047058824 0.042440318 0.038461538
y1 -0.048076923 0.321153846 0.567307692 0.321153846 -0.048076923
y2 1.57872099 -0.226196289 0.253755457 0.235346591 0.84340743 0.84340743 0.235346591 0.253755457 -0.226196289 1.57872099
y3 -39.95244903 3.4549578 -0.447051961 0.202422616 0.080659993 0.17976263 0.238445934 0.395093054 0.636755336 0.94249038 0.94249038 0.636755336 0.395093054 0.238445934 0.17976263 0.080659993 0.202422616 -0.447051961 3.4549578 -39.95244903
y4 -57409.17974 2287.728499 -156.169717 15.42434982 -1.939791158 0.39938813 0.014910807 0.10858616 0.103537161 0.128151723 0.150063036 0.181523667 0.221349033 0.274733128 0.345913865 0.441399666 0.566357999 0.719110364 0.876706771 0.984617272 0.984617272 0.876706771 0.719110364 0.566357999 0.441399666 0.345913865 0.274733128 0.221349033 0.181523667 0.150063036 0.128151723 0.103537161 0.10858616 0.014910807 0.39938813 -1.939791158 15.42434982 -156.169717 2287.728499 -57409.17974
%构造lagrange插值函数，参考：https://blog.csdn.net/qq_44692057/article/details/108314856
function y=lagr(X,Y,x)
n=length(X);
m=length(x);
for i=1:m
z=x(i);
s=0;
for k=1:n
p=1;
for j=1:n
if j~=k
p=p*(z-X(j))/(X(k)-X(j));
end
end
s=p*Y(k)+s;
end
y(i)=s;
end
%输入已知量(第2次a问）
clear;
n=[5 10 20 40];
for a=1:4
b=10/n(a);
i=1:n(a)+1;
j=-5:b:5;
x=[i;j];
for k=1:n(a)+1
X(k)=x(2,k)
Y(k)=1/(1+X(k)*X(k));
end
for l=1:n(a)
x0(l)=(x(2,l)+x(2,l+1))/2;
end
if a==1
y1=lagr(X,Y,x0);
x1=x0;
elseif a==2
y2=lagr(X,Y,x0);
x2=x0;
elseif a==3
y3=lagr(X,Y,x0);
x3=x0;
elseif a==4
y4=lagr(X,Y,x0);
x4=x0;
end
end
%画图比较精确值
figure(1)
xc=[-5:0.05:5]
plot(xc,1./(1+xc.*xc),'k')
hold on
plot(x1,y1,'r')
hold on
plot(x2,y2,'y')
hold on
plot(x3,y3,'b')
hold on
plot(x4,y4,'g')
hold off
figure(2)
xc=[-5:0.05:5]
plot(xc,1./(1+xc.*xc),'k')
hold on
plot(x1,y1,'r')
hold on
plot(x2,y2,'y')
hold on
plot(x3,y3,'b')
hold on
plot(x4,y4,'g')
hold off
axis([-5,5,-1,1])
得出图像
（b）
%构造newton插值参考https://zhuanlan.zhihu.com/p/34883522
function N = Newton( x,y,t )
syms p ; %定义符号变量
N = y(1);
dd = 0;
dxs = 1;
n = length(x);
for(i = 1:n-1)
for(j = i+1:n)
dd(j) = (y(j)-y(i))/(x(j)-x(i));
end
temp1(i) = dd(i+1);
dxs = dxs*(p-x(i));
N = N + temp1(i)*dxs;
y = dd;
end
simplify(N);
%以上为计算部分，下面是输出规则；
if(nargin == 2)
N = subs(N,'p','x');
N = collect(N);
N = vpa(N,4);
else
m = length(t);
for i = 1:m
temp(i) = subs(N,'p',t(i));
end
N = temp;
end
%输入已知量(第2次a问）
clear;
n=[5 10 20 40];
for a=1:4
b=10/n(a);
i=1:n(a)+1;
j=-5:b:5;
x=[i;j];
for k=1:n(a)+1
X(k)=x(2,k)
Y(k)=1/(1+X(k)*X(k));
end
for l=1:n(a)
x0(l)=(x(2,l)+x(2,l+1))/2;
end
if a==1
y1=Newton(X,Y,x0);
x1=x0;
elseif a==2
y2=Newton(X,Y,x0);
x2=x0;
elseif a==3
y3=Newton(X,Y,x0);
x3=x0;
elseif a==4
y4=Newton(X,Y,x0);
x4=x0;
end
end
%画图比较精确值
figure(1)
xc=[-5:0.05:5]
plot(xc,1./(1+xc.*xc),'k')
hold on
plot(x1,y1,'r')
hold on
plot(x2,y2,'y')
hold on
plot(x3,y3,'b')
hold on
plot(x4,y4,'g')
hold off
figure(2)
xc=[-5:0.05:5]
plot(xc,1./(1+xc.*xc),'k')
hold on
plot(x1,y1,'r')
hold on
plot(x2,y2,'y')
hold on
plot(x3,y3,'b')
hold on
plot(x4,y4,'g')
hold off
axis([-5,5,-1,1])
得出图像
（c）
%构造分段线性插值函数
function y=fenduan(X,Y,x)
n=length(X);
m=length(x);
for j=1:m
for i=1:n-1
if x(j)>X(i)&&x(j)<=X(i+1)
y(j)=((x(j)-X(i+1))/(X(i)-X(i+1)))*Y(i)+(((x(j)-X(i))/(X(i+1)-X(i)))*Y(i+1));
end
end
end
%第二次c问
clear;
n=[5 10 20 40];
for a=1:4
b=10/n(a);
i=1:n(a)+1;
j=-5:b:5;
x=[i;j];
for k=1:n(a)+1
X(k)=x(2,k)
Y(k)=1/(1+X(k)*X(k));
end
for l=1:n(a)
x0(l)=(x(2,l)+x(2,l+1))/2;
end
if a==1
y1=fenduan(X,Y,x0);
x1=x0;
elseif a==2
y2=fenduan(X,Y,x0);
x2=x0;
elseif a==3
y3=fenduan(X,Y,x0);
x3=x0;
elseif a==4
y4=fenduan(X,Y,x0);
x4=x0;
end
end
%画图比较精确值
xc=[-5:0.05:5]
plot(xc,1./(1+xc.*xc),'k')
hold on
plot(x1,y1,'r')
hold on
plot(x2,y2,'y')
hold on
plot(x3,y3,'b')
hold on
plot(x4,y4,'g')
hold off
得出图像
2、
%d2
clear;
n=[4 6 8 10];
x0=[-1:0.005:1];
for a=1:4
b=2/n(a);
X=-1:b:1;
for i=1:n(a)+1
Y(i)=1/(1+25*X(i)*X(i))
end
if a==1
y1=lagr(X,Y,x0);
elseif a==2
y2=lagr(X,Y,x0);
elseif a==3
y3=lagr(X,Y,x0);
elseif a==4
y4=lagr(X,Y,x0);
end
end
%画图比较精确值
plot(x0,1./(1+25*x0.*x0),'k')
hold on
plot(x0,y1,'r')
hold on
plot(x0,y2,'y')
hold on
plot(x0,y3,'b')
hold on
plot(x0,y4,'g')
hold off
得出图像
