数值分析实验(2.3.5章)

试验2.1 多项式插值的振荡现象实验目的：观察多项式插值的振荡现象，了解多项式的次数与逼近效果的关系。

实验内容：问题提出：考虑在一个固定的区间上用插值逼近一个函数。

显然Lagrange 插值中使用的节点越多，插值多项式的次数就越高。

我们自然关心插值多项式的次数增加时，Ln(x)是否也更加靠近被逼近的函数。

Runge 给出的一个例子是极著名并富有启发性的。

设区间[-1,1]上的函数225x11)x (+=f ，考虑区间[-1，1]上的一个等距划分，分点为n2i 1x i +-=，i=0，1，2，…，n则拉格朗日插值多项式为：)x (l 25x11)x (Ln i ni 2i∑=+=，其中的)x (l i ，i=0，1，2，…，n 是n 次拉格朗日插值基函数。

实验要求：1、选择不断增大的分点数目n=2，3，………，画出原函数)x (f 及插值多项式函数)x (Ln 在[-1，1]上的图像，比较并分析试验结果。

2、选择其他的函数，例如定义在区间[-5，5]上的函数4()1x h x x=+，()arctan g x x =，重复上述的实验看其结果如何。

实验步骤及结果分析：1、选择不断增大的分点数目n=2，3，4，5，6，7，8，9，10做)x (f 的拉格朗日插值多项式)x (Ln ，并与原函数值做比较，如下图所示。

观察图像可知：n=2,3时插值函数和原函数差别很大，n=4，5，6时插值函数与原函数的逼近程度相对较好，继续增加插值次数n ，插值函数在插值区域的中间部分收敛，而在这区间外是发散的，此外，n=7,9时在插值中间区域逼近效果不好。

因此，适当提高插值多项式次数，可以提高逼近的精度，但是次数太高反而产生相反的效果。

2、选择其他的函数进行插值。

原函数4()1x h x x=+，区间[-5，5]，插值结果如下图：观察图像可知：低次插值时，插值效果不好。

n=7,8,9,10时，在区间[-2，2]，插值函数与原函数逼近程度好，但在区间外插值函数发散。

其中，n=8,10插值效果比n=7,9好。

原函数()arctan，区间[-5，5]，插值结果如下图：g x x观察图像可知：n=5,6,7,8,9,10时，在区间[-3，3]，插值函数与原函数逼近程度好，但在区间外插值函数发散。

其中，n=7,9在插值区间两端发散的更剧烈。

分析在插值区间两端发散的原因：次数越高，计算量就越大，积累误差也大，在整个区间上做高次多项式，当局部插值节点处的值有微小偏差时，可能引起整个区间上函数值的很大变化，使计算不稳定。

Matlab程序如下：function t_charpt2%数值实验二%输入：实验选择，函数式选择，插值节点数%输出：拟合函数及原函数的图形result=inputdlg({'请选择实验,若选2.1,请输入1,否则输入2:'},'charpt_2',1,{'1'}); Nb=str2num(char(result));if(Nb~=1)&(Nb~=2) errordlg('实验选择错误!');return;endpromps={'请选择实验函数，若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:'};result=inputdlg(promps,'charpt 2',1,{'f'});Nb_f=char(result);if(Nb_f~='f'&Nb_f~='h'&Nb_f~='g')errordlg('实验选择错误！');return;endresult=inputdlg({'请输入插值点数N:'},'charpt_2',1,{'10'});Nd=str2num(char(result));if(Nd<1)errordlg('结点输入错误！');return;endswitch Nb_fcase 'f'f=inline('1./(1+25*x.^2)');a=-1;b=1;case 'h'f=inline('x./(1+x.^4)');a=-5;b=5;case 'g'f=inline('atan(x)');a=-5;b=5;endif(Nb==1)x0=linspace(a,b,Nd+1);y0=feval(f,x0);x=a:0.1:b;y=lagrange(x0,y0,x);% clf;把曲线显示在一张图上fplot(f,[a,b],'co');hold on;plot(x,y,'b--');xlabel('x') ;ylabel('y=f(x)o and y=ln(x)--');elseif(Nb==2)x0=linspace(a,b,Nd+1);y0=feval(f,x0);x=a:0.1:b;cs=spline(x0,y0);y=ppval(cs,x);% clf;把曲线显示在一张图上fplot(f,[a,b],'co');hold on;plot(x,y,'k-');xlabel('x');ylabel('y=f(x) o and y=spline(x)-'); end%function y=lagrange(x0,y0,x)%Lagrange插值n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif (j~=k)p=p*(z-x0(j))/(x0(k)-x0(j));endends=s+p*y0(k);endy(i)=s;end试验3.1 多项式最小二乘拟合实验目的：掌握多项式最小二乘拟合程序，会对数据作多项式最小二乘拟合。

实验内容及要求：编制以函数{}nk k x 0=为基的多项式最小二乘拟合程序，并用于对表中的数据作3次多项式最小拟合。

取权数1≡i w ，求拟合曲线knk kxa∑=**=ϕ中的参数{}k a 、平方误差2δ，并作离散数据{}i i y x ,的拟合函数()x y *=ϕ的图形。

实验步骤及结果分析：拟合函数()x y *=ϕ的图形：拟合曲线knk kxa∑=**=ϕ中的平方误差2δ为2.17619e-005；参数{}k a 为alph 中的四个数值，按降幂排列1.99911，-2.99767，-3.96825e-005，0.549119。

观察图像可知：表中数据均在拟合曲线上，拟合结果很好，平方误差为10-5数量级。

Matlab程序如下：function charpt3%数值实验三：含“实验3.1”和“实验3.2”%子函数调用：dlsa%输入：实验选择%输出：原函数及求得的相应插值多项式的函数的图像以及参数alph和误差rresult=inputdlg({'请选择实验，若选3.1，请输入1，否则输入2：'},'charpt_3',1,{'1'}); Nb=str2num(char(result));if((Nb~=1)&(Nb~=2)) errordlg('实验选择错误！'); return; endx0=-1:0.5:2;y0=[-4.447 -0.452 0.551 0.048 -0.447 0.549 4.552];n=3;%n为拟合阶次if Nb==1alph=polyfit(x0,y0,n);y=polyval(alph,x0);r=(y0-y) * (y0-y)';%平方误差x=-1:0.01:2;y=polyval(alph,x);plot(x,y,'k--');xlabel('x'); ylabel('y0 * and polyfit.y--');hold onplot(x0,y0,'*');title('离散数据的多项式拟合');grid on;elseresult=inputdlg({'请输入权向量w:'},'charpt_3',1,{'[1 1 1 1 1 1 1]'});w=str2num(char(result));[a,b,c,alph,r]=dlsa(x0,y0,w,n);enddisp(['平方误差：',sprintf('%g',r)]);disp(['参数alph：', sprintf('%g\t',alph)])function [a,b,c,alph,r]=dlsa(x,y,w,n)%功能：用正交化方法对离散数据作多项式最小二乘拟合。

%输入：m+1个离散点（x,y,w），x，y，w分别用行向量给出。

% 拟合多项式的次数n，0〈n<m。

%输出：三项递推公式的参数a,b,拟合多项式s(x)的系数c和alph，% 平方误差r=（y-s,y-s），并作离散点列和拟合曲线的图形m=length(x)-1;if (n<1|n>=m) errordlg('错误：n<1或者n>=m!'); return;end%求三项递推公式的参数a，b，拟合多项式s(x)的系数c，其中d(k)=(y,sk);s1=0;s2=ones(1,m+1);v2=sum(w);d(1)=y*w';c(1)=d(1)/v2;for k=1:nxs=x.*s2.^2*w';a(k)=xs/v2;if k==1 b(k)=0;else b(k)=v2/v1;ends3=(x-a(k)).*s2-b(k)*s1;v3=s3.^2*w';d(k+1)=y.*s3*w';c(k+1)=d(k+1)/v3;s1=s2;s2=s3;v1=v2;v2=v3;end%求平方误差rr=y.^2*w'-c*d';%求拟合多项式s(x)的降幂系数alphalph=zeros(1,n+1); T=zeros(n+1,n+2);T(:,2)=ones(n+1,1); T(2,3)=-a(1);if n>=2for k=3:n+1for i=3:k+1T(k,i)=T(k-1,i)-a(k-1)*T(k-1,i-1)-b(k-1)*T(k-2,i-2);endendendfor i=1:n+1for k=i:n+1alph(n+2-i)=alph(n+2-i)+c(k)*T(k,k+2-i);endend%用秦九韶方法计算s(x)的输出序列(t,s)xmin=min(x); xmax=max(x); dx=(xmax-xmin)/(25*m);t=(xmin-dx):dx:(xmax+dx);s=alph(1);for k=2:n+1s=s.*t+alph(k);end%输出点列x-y和拟合曲线t-s的图形plot(t,s,'-',x,y,'*');title('离散数据的正交化多项式拟合');xlabel('x');ylabel('y');grid on;试验5.1 常微分方程性态和R-K 法稳定性试验实验目的：考察下面微分方程右端项中函数y 前面的参数对方程性态的影响（它可使方程为好条件的或坏条件的）和研究计算步长对R-K 法计算稳定性的影响。