三种插值多项式程序

程序如下：
% y=1./(1+25*x.^2)原始函数图形
clear
format long
x0=-1:0.2:1; y0=1./(1+25*x0.^2);
subplot(2,2,1), plot(x0,y0,'-ko')
grid; xlabel('x轴'),ylabel('y轴'),title('y=1/(1+25*x.^2)原始函数图形');
legend('y=1/(1+25*x.^2)',6)
% (1) P(n)等距节点插值多项式
syms f x df s s1 s2 s3 s4;
f=1./(1+25*x^2);
df=diff(f);
N=10;
h=2/N;
xx=-1:2/N:1;
p=1; L=0;
ff=zeros(1,length(xx));
for i=1:(N+1)
x=xx(i); ff(i)=eval(f); dff(i)=eval(df);
end
syms x
for i=1:N
s1=(x-xx(i+1))^2*(h+2*(x-xx(i)))*ff(i)/h^3; s2=(x-xx(i))^3*(h+2*(xx(i+1)-x))*ff(i+1)/h^3;
s3=(x-xx(i+1))^2*(x-xx(i))*dff(i)/h^2; s4=(x-xx(i))^2*(x-xx(i+1))*dff(i+1)/h^2;
s(i)=s1+s2+s3+s4;
end
hold on
subplot(2,2,2), plot(xx,ff,'--gd')
grid; xlabel('x轴'),ylabel('y轴'),title('P(n)等距节点插值多项式');
legend('y=1/(1+25*x.^2)',6)
% (2) n+1次chebyshev多项式的零点作插值多项式
clear
format long
n=10; %结点个数
lb=-1; %下界
ub=1; %上界
step=0.2; %作图点步长
% 插值函数
for i=1:n+1
xi(i)=(lb+ub)/2+((ub-lb)/2)*cos((2*i-1)*pi/(2*(n+1)));
yi(i)=1/(1+25*xi(i)^2);
end
count=1;
for x=lb:step:ub
fl=0;
%求出pn(xk)
for k=1:n+1
up=1; dn=1;
%求出f(xk)
for i=1:n+1
if k~=i
up=up*(x-xi(i)); dn=dn*(xi(k)-xi(i));
end
end
fl=fl+yi(k)*up/dn;
end
pn(count)=fl;
fi(count)=1./(1+25*x^2);%求原函数的值
count=count+1;
end
% L插值函数图
x=lb:step:ub;
subplot(2,2,3), plot(x,pn,'-.rh')
grid; xlabel('x轴'),ylabel('y轴'),title('n+1次chebyshev多项式的零点作插值多项式'); legend('y=1/(1+25*x.^2)',6)
% (3) 三次样条函数Sn(x):
clear;
format long
x = -1:0.2:1; y =1./(1+25*x.^2);
xx = -1:0.2:1;
yy =spline(x,y,xx);
subplot(2,2,4), plot(x,y,':bp',xx,yy)
grid; xlabel('x轴'),ylabel('y轴'),title('三次样条函数Sn(x)');
legend('y=1/(1+25*x.^2)',6)
程序如下：
(1) 程序：确定n,并用复合梯形公式Tn
syms x
L= inline('exp(x)'); ;
[QS,FCNTS] =quad8(L,0, 1,1.e-8)
结果：n=18
(2) 程序：确定n,并用复合辛普森公式Sn
syms x
L= inline('exp(x)');
[QS,FCNTS] =quad(L,0, 1,1.e-8)
结果：n=33
(3) 程序：龙贝格积分：
syms f x;
f=input('请输入积分函数f=');
A=input('请输入积分区间[a,b]=');
e=input('请输入误差限e=');
x=A(1);
fa=eval(f);
x=A(2);
fb=eval(f);
T(1)=(A(2)-A(1))*(fa+fb)/2;
m=1;
x=(A(1)+A(2));
t=T(1)/2+(A(2)-A(1))*eval(f)/2;
while abs(t-T(1))>e
t=T(1);
new=0;
for i=1:(2^m)
x=A(1)+(2*i-1)*(A(2)-A(1))/(2^(m+1));
new=new+eval(f);
end
T(m+1)=T(m)/2+(A(2)-A(1))*new/(2^(m+1));
for i=m:(-1):1
L(m,i)=T(i);
end
for p=m:(-1):1
T(p)=((4^(m+1-p))*T(p+1)-T(p))/(4^(m+1-p)-1);
end
m=m+1;
end
fprintf('数值积分结果为T(0)=%.8f\n',T(1));
fprintf('经过了%d次迭代\n',m);
结果：
请输入积分函数f=exp(x)
请输入积分区间[a,b]=[0,1]
请输入误差限e=1.0000e-008
数值积分结果为T(0)=1.71828184
经过了24次迭代。