一、第一题：
线性插值：
-
两点式线性插值
function liner() x=input("请输入插值点：");
y=x*((sin(0.6)-sin(0.5))/(0.6-0.5));
Re=y-sin(x);
fprintf("插值点近似值y=%2f，截断误差Re=%2f\n",y,Re)
运行结果： >> liner 请输入插值点：0.57891 插值点近似值y=0.493329，截断误差Re=-0.053783-
调用Matlab库函数
function liner()
x0=0.5:0.1:0.6;
y0=sin(x0);
x=input("请输入插值点：");
y=interp1(x0,y0,x,"linear"); %没有指定的情况下默认选择内插
Re=y-sin(x);
fprintf("插值点近似值y=%2f，截断误差Re=%2f\n",y,Re)
运行结果 >> liner 请输入插值点：0.57891 插值点近似值y=0.546670，截断误差Re=-0.000442二次插值：列举了两种方法（拉格朗日二次插值、牛顿二次插值）
-
拉格朗日二次插值：
function lagrangen_Quadratic()
x0=0.4:0.1:0.6;
y0=sin(x0);
n=length(x0); %记录插值节点个数 x=input("请输入插值点：");
y=0.0;
for k=1:n p=1.0;
for j=1:n
if(j~=k)
p=p*(x-x0(j))/(x0(k)-x0(j)); %连乘
end
end
y=y+p*y0(k); %求和
end
fprintf('插值点近似值y=%f\n',y)
end
运行结果 >> lagrangen_Quadratic 请输入插值点：0.57891 插值点近似值y=0.547069，截断误差Re=-0.000043-
牛顿二次插值
function newton_Quadratic()
x0=0.4:0.1:0.6;
y0=sin(x0);
n=length(x0);
x=input("请输入插值点：");
A=zeros(n,n);A(:,1)=y0'; % A是差商表
y=0.0;
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)- A(i-1,j-1))/(x0(i)-x0(i-j+1));
end
end
for k=1:n
p=1.0;
for l=1:k-1
p=p*(x-x0(l));
end
y=y+A(k,k)*p;
end
fprintf("插值点近似值y=%f\n",y)
end
运行结果
>> newton_Quadratic
请输入插值点：0.57891
插值点近似值y=0.547069，截断误差Re=-0.000043三次插值
function Cubic()
%三次插值
%% 三次拉格朗日插值
x0=0.3:0.1:0.6;
y0=sin(x0);
n=length(x0); %记录插值节点个数
x=input("请输入插值点：");
y1=0.0;
for k=1:n
p=1.0;
for j=1:n
if(j~=k)
p=p*(x-x0(j))/(x0(k)-x0(j)); %连乘
end
end
y1=y1+p*y0(k); %求和
end
%误差计算
Re=y1-sin(x);
fprintf("三次拉格朗日插值法，插值点近似值y1=%2f，截断误差Re=%2f\n",y1,Re)
%% 三次牛顿插值
y2=0.0;
A=zeros(n,n);A(:,1)=y0'; % A是差商表
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)- A(i-1,j-1))/(x0(i)-x0(i-j+1));
end
end
for k=1:n
p=1.0;
for l=1:k-1
p=p*(x-x0(l));
end
y2=y2+A(k,k)*p;
end
%误差计算
Re=y2-sin(x);
fprintf("三次牛顿插值法，插值点近似值y2=%2f，截断误差Re=%2f\n",y2,Re)
%% 分段三次插值
y3=interp1(x0,y0,x,"pchip");
%误差计算
Re=y3-sin(x);
fprintf("三次分段插值，插值点近似值y3=%2f，截断误差Re=%2f\n",y3,Re)
%% 三次样条插值
y4=interp1(x0,y0,x,"spline");
%误差计算
Re=y4-sin(x);
fprintf("三次样条插值，插值点近似值y4=%2f，截断误差Re=%2f\n",y4,Re)
运行结果
>> Cubic
请输入插值点：0.57891
三次拉格朗日插值法，插值点近似值y1=0.547113，截断误差Re=0.000002
三次牛顿插值法，插值点近似值y2=0.547113，截断误差Re=0.000002
三次分段插值，插值点近似值y3=0.547067，截断误差Re=-0.000045
三次样条插值，插值点近似值y4=0.547113，截断误差Re=0.000002结果分析：通过对比不同插值方法，可以看到在一定范围内（高次会出现龙格现象），插值次数越高，截断误差越小（插值结果越接近于真实函数值）；同时，对于相同次数的插值，由于不同的插值方法它们算法不尽相同，故得到的结果也有所差异。
二、第二题
拉格朗日插值法得到插值多项式
function lagrange_Interplotion()
%% lagrange插值多项式
x0=1:0.5:5;
y0=[2.71828,4.48169,7.38906,12.0855,20.0855,33.1155,54.5982,54.5982,90.0171];
n=length(x0); %记录插值节点个数
for k=1:n
c_lx=1.0; %基函数多项式系数
for j=1:n
if(j~=k)
c_x0_j=poly(x0(j)); %得到以节点x0(j)为根的多项式的系数
c_lx=conv(c_lx,c_x0_j)/(x0(k)-x0(j)); %卷积，得到下一个节点的基函数多项式系数
end
end
c(k,:)=c_lx; %将得到的系数写成列向量，方便后面图形绘制
l(k,:)=poly2sym(c_lx); %生成多项式
end
C=y0*c; %拉格朗日插值函数系数向量
fprintf("基函数:\n")
for k=1:n
fprintf("q%d(x)=%s\n",k,l(k))
end
L=vpa(y0*l);
fprintf('拉格朗日插值多项式为：\nP(x)=%s\n',L)
plot(x0,y0,"d",LineWidth=2)
hold on;
plot(x0,y0,"LineWidth",3,"LineStyle","-.")
hold on;
x=input("请输入估计点值x=");
y=polyval(C,x);
fprintf("由插值多项式计算得到估计点函数值y=%2f\n",y)
x1=1:0.1:5;
plot(x1,polyval(C,x1),"-g","LineWidth",2)
title("拉格朗日法插值多项式")
legend("原始数据点","线性插值曲线","拉格朗日插值曲线")
end
L=vpa
运行结果：
>> Interplotion
基函数:
q1(x)=(43903*x^2)/280 - (3601*x)/28 - (19201*x^3)/180 + (15967*x^4)/360 - (104*x^5)/9 + (83*x^6)/45 - (52*x^7)/315 + (2*x^8)/315 + 45
q2(x)=(16084*x)/21 - (8897680915588293*x^2)/8796093022208 + (3654*x^3)/5 - (4774*x^4)/15 + 86*x^5 - (212*x^6)/15 + (136*x^7)/105 - (16*x^8)/315 - 240
q3(x)=(29379*x^2)/10 - (4231*x)/2 - (39953*x^3)/18 + (90409*x^4)/90 - (2524*x^5)/9 + (2132*x^6)/45 - (40*x^7)/9 + (8*x^8)/45 + 630
q4(x)=(17428*x)/5 - (24994*x^2)/5 + (4290489840315779*x^3)/1099511627776 - (81994*x^4)/45 + (23558*x^5)/45 - (4084*x^6)/45 + (392*x^7)/45 - (16*x^8)/45 - 1008
q5(x)=(5965736298384953*x^2)/1099511627776 - (22205*x)/6 - 4339*x^3 + (24931*x^4)/12 - 612*x^5 + (326*x^6)/3 - (32*x^7)/3 + (4*x^8)/9 + 1050
q6(x)=2572*x - (19158*x^2)/5 + (857826755194971*x^3)/274877906944 - (68614*x^4)/45 + (4130*x^5)/9 - (3748*x^6)/45 + (376*x^7)/45 - (16*x^8)/45 - 720
q7(x)=(17137*x^2)/10 - (2273*x)/2 - (6226998586338327*x^3)/4398046511104 + (63289*x^4)/90 - (1940*x^5)/9 + (1796*x^6)/45 - (184*x^7)/45 + (8*x^8)/45 + 315
q8(x)=(2036*x)/7 - (3894463204556383*x^2)/8796093022208 + 370*x^3 - (2794*x^4)/15 + 58*x^5 - (164*x^6)/15 + (8*x^7)/7 - (16*x^8)/315 - 80
q9(x)=(14139*x^2)/280 - (4609*x)/140 - (7667*x^3)/180 + (7807*x^4)/360 - (308*x^5)/45 + (59*x^6)/45 - (44*x^7)/315 + (2*x^8)/315 + 9
拉格朗日插值多项式为：
P(x)=13219.923975626989297054402568969*x^2 - 8717.4714292857121464782161573177*x - 10997.479117888890432702431275579*x^3 + 5508.5890782777766034013937870542*x^4 - 1705.6709622222218751682652913991*x^5 + 319.55153377777771566815722002907*x^6 - 33.174790603174597008322933337666*x^7 + 1.4632923174603171998627196502177*x^8 + 2406.9866999999993839111311899615
请输入估计点值x=4.25
由插值多项式计算得到估计点函数值y=59.771361
牛顿法插值多项式
function nweton_Interplotion()
%% 牛顿插值多项式
x0=1:0.5:5;
y0=[2.71828,4.48169,7.38906,12.0855,20.0855,33.1155,54.5982,54.5982,90.0171];
n=length(x0); %记录插值节点个数
A=zeros(n,n);A(:,1)=y0'; % A是差商表
%先求差商表
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)- A(i-1,j-1))/(x0(i)-x0(i-(j-1)));
if i==j
C(1,1)=y0(1); %补上初始值y0
C(:,j)=A(i,j); %n阶均差
end
end
end
%求基函数
for k=1:n
c_lx=1.0;
for j=1:k-1
c_lx=conv(c_lx,poly(x0(j)));
end
l(k,:)=poly2sym(c_lx);
end
fprintf("基函数为：\n")
for k=1:n
fprintf("q(x)=%s\n",l(k))
end
L=vpa(C*l);
fprintf('牛顿插值多项式为：\nP(x)=%s\n',L)
plot(x0,y0,"d",LineWidth=2)
hold on;
plot(x0,y0,"LineWidth",3,"LineStyle","-.")
hold on;
x=input("请输入估计点值x=");
Coff=sym2poly(C*l); %得到插值多项式系数
y=polyval(Coff,x);
fprintf("由插值多项式计算得到估计点函数值y=%2f\n",y)
%绘制插值曲线
x1=1:0.01:5;
plot(x1,polyval(Coff,x1),"-g","LineWidth",2)
title("牛顿法插值多项式")
legend("原始数据","线性插值曲线","牛顿插值曲线")
end
运行结果：
>> nweton_Interplotion
基函数为：
q(x)=1
q(x)=x - 1
q(x)=x^2 - (5*x)/2 + 3/2
q(x)=(13*x)/2 - (9*x^2)/2 + x^3 - 3
q(x)=(71*x^2)/4 - (77*x)/4 - 7*x^3 + x^4 + 15/2
q(x)=(261*x)/4 - (145*x^2)/2 + (155*x^3)/4 - 10*x^4 + x^5 - 45/2
q(x)=319*x^2 - (2007*x)/8 - (1665*x^3)/8 + (295*x^4)/4 - (27*x^5)/2 + x^6 + 315/4
q(x)=(4329*x)/4 - (12215*x^2)/8 + (2303*x^3)/2 - (4025*x^4)/8 + (511*x^5)/4 - (35*x^6)/2 + x^7 - 315
q(x)=(127251*x^2)/16 - (41481*x)/8 - (53669*x^3)/8 + (54649*x^4)/16 - 1078*x^5 + (413*x^6)/2 - 22*x^7 + x^8 + 2835/2
牛顿插值多项式为：
P(x)=13219.923975626990916816144346531*x^2 - 8717.4714292857188404819035465465*x - 10997.479117888894418232977878702*x^3 + 5508.5890782777804758410492939428*x^4 - 1705.6709622222230317151181466784*x^5 + 319.55153377777792386704724947322*x^6 - 33.174790603174617698023496359383*x^7 + 1.4632923174603180704878013784764*x^8 + 2406.9867000000012735332943769606
请输入估计点值x=4.25
由插值多项式计算得到估计点函数值y=59.771361
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
