种插值多项式程序.doc
文本预览下载声明
程序如下：
% 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=);
显示全部