数值分析实验报告

1、 实验2.1 多项式插值的振荡现象实验目的：在一个固定的区间上用插值逼近一个函数，显然Lagrange插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，Ln(x)是否也更加靠近被逼近的函数。Runge给出的一个例子是极著名并富有启发性的。实验容：设区间-1，1上函数 f(x)=1/(1+25x2)。考虑区间-1,1的一个等距划分，分点为 xi= -1 + 2i/n，i=0，1，2，n，则拉格朗日插值多项式为.其中，li(x)，i=0，1，2，n是n次Lagrange插值基函数。实验步骤与结果分析：实验源程序function Chap2Interpolation%

2、 数值实验二：“实验2.1：多项式插值的震荡现象”% 输入：函数式选择，插值结点数% 输出：拟合函数及原函数的图形promps = '请选择实验函数，若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:'titles = 'charpt_2'result = inputdlg(promps,'charpt 2',1,'f');Nb_f = char(result);if(Nb_f = 'f' & Nb_f = 'h' & Nb_f = 'g')er

3、rordlg('实验函数选择错误！');return;endresult = inputdlg('请输入插值结点数N:','charpt_2',1,'10');Nd = str2num(char(result);if(Nd <1)errordlg('结点输入错误！');return;endswitch Nb_f case 'f' f=inline('1./(1+25*x.2)'); a = -1;b = 1; case 'h' f=inline('x./

4、(1+x.4)'); a = -5; b = 5; case 'g' f=inline('atan(x)'); a = -5; b= 5;end x0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b; y = Lagrange(x0, y0, x); fplot(f, a b, 'co'); hold on; plot(x, y, 'b-'); xlabel('x'); ylabel('y = f(x) o and y = Ln(x)-&

5、#39;);%-function y=Lagrange(x0, y0, x);n= length(x0); m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if(j = k) p = p*(z - x0(j)/(x0(k) - x0(j); end end s = s + p*y0(k); end y(i) = s;end实验结果分析(1)增大分点n=2，3，时，拉格朗日插值函数曲线如图所示。 n=6 n=7 n=8 n=9 n=10从图中可以看出，随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f

6、(x)，但是却在两端x= -1和x=1处出现了很大的振荡现象。并且，仔细分析图形，可以看出，当n为奇数时，虽然有振荡，但振荡的幅度不算太大，n为偶数时，其振荡幅度变得很大。通过思考分析，我认为，可能的原因是f(x)本身是偶函数，如果n为奇数，那么Lagrange插值函数Ln(x)的最高次项xn-1是偶次幂，比较符合f(x)本身是偶函数的性质；如果n为偶数，那么Lagrange插值函数Ln(x)的最高次项xn-1是奇次幂，与f(x)本身是偶函数的性质相反，因此振荡可能更剧烈。(2)将原来的f(x)换为其他函数如h(x)、g(x)，结果如图所示。其中h(x), g(x)均定义在-5，5区间上，h(

7、x)=x/(1+x4)，g(x)=arctan x。h(x), n=7 h(x), n=8h(x), n=9 h(x), n=10g(x), n=7g(x), n=8g(x), n=9g(x), n=10分析两个函数的插值图形，可以看出：随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x)，但是却在两端x= -5和x=5处出现了很大的振荡现象。并且，仔细分析图形，可以看出，当n为偶数时，虽然有振荡，但振荡的幅度不算太大，n为奇数时，其振荡幅度变得很大。原因和上面f(x)的插值类似，h(x)、g(x)本身是奇函数，如果n为偶数，那么Lagrange插值函数Ln(x)的最高次项

8、xn-1是奇次幂，比较符合h(x)、g(x)本身是奇函数的性质；如果n为奇数，那么Lagrange插值函数Ln(x)的最高次项xn-1是偶次幂，与h(x)、g(x)本身是奇函数的性质相反，因此振荡可能更剧烈。实验3.1 多项式最小二乘拟合实验目的： 编制以函数xkk=0,n；为基的多项式最小二乘拟合程序。实验容： 对表中的数据作三次多项式最小二乘拟合。xi-1.0-0.50.00.51.01.52.0yi-4.447-0.4520.5510.048-0.4470.5494.552取权函数wi1，求拟合曲线中的参数k、平方误差2，并作离散据xi, yi的拟合函数的图形。实验源程序function

9、 Chap3CurveFitting% 数值实验三:“实验3.1”% 输出:原函数及求得的相应插值多项式的函数的图像以及参数alph和误差rx0 = -1:0.5:2;y0 = -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552;n = 3; % n为拟合阶次alph = 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('

10、x'); ylabel('y0 * and polyfit.y-');hold onplot(x0, y0, '*')grid on;disp('平方误差:', num2str(r)disp('参数alph:', num2str(alph)实验结果平方误差:2.1762e-005参数alph:1.9991 -2.9977 -3.9683e-005 0.54912实验4.1 实验目的：复化求积公式计算定积分. 实验题目：数值计算下列各式右端定积分的近似值. 实验要求： （1）若用复化梯形公式、复化Simpson公式和复化Ga

11、uss-Legendre I 型公式做计算，要求绝对误差限为，分别利用它们的余项对每种算法做出步长的事前估计. （2）分别用复化梯形公式，复化Simpson 公式和复化Gauss-Legendre I 型公式作计算. （3）将计算结果与精确解做比较，并比较各种算法的计算量. 实验程序：1.事前估计的Matlab程序如下： （1）用复化梯形公式进行事前估计的Matlab程序 format long g x=2:0.01:3; f=-4*(3*x.2+1)./(x.2-1).3; %二阶导函数 %plot(x,f) %画出二阶导函数图像 x=2.0; %计算导函数最大值 f=-4*(3*x2+1)

12、/(x2-1)3; h2=0.5*10(-7)*12/f; h=sqrt(abs(h2) %步长 n=1/h; n=ceil(1/h)+1 %选取的点数 format long g x=0:0.01:1; f=8.*(3*x.2-1)./(x.2+1).3; %二阶导函数 %plot(x,f) %画出二阶导函数图像 x=1; %计算导函数最大值 f=8.*(3*x.2-1)./(x.2+1).3; h2=0.5*10(-7)*12/f; h=sqrt(abs(h2) %步长 n=1/h n=ceil(1/h)+1 %选取的点数 format long g x=0:0.01:1; f=log(3

13、).*log(3).*3.x; %二阶导函数 %plot(x,f); %画出二阶导函数图像 x=1; %计算导函数最大值 f=log(3)*log(3)*3x; h2=0.5*10(-7)*12/f; h=sqrt(abs(h2) %步长 n=1/h n=ceil(1/h)+1 %选取的点数 format long g x=1:0.01:2; f=2.*exp(x)+x.*exp(x);%二阶导函数 %plot(x,f) %画出二阶导函数图像 x=2; %计算导函数最大值 f=2.*exp(x)+x.*exp(x); h2=0.5*10(-7)*12/f; h=sqrt(abs(h2) %步长

14、 n=1/h n=ceil(1/h)+1 %选取的点数估计结果步长h及结点数n分别为 h = 0.8 n =1793 h =0.6 n =1827 h =0.9 n =2458 h =0.9 n =7020 （2）用复化simpson公式进行事前估计的Matlab程序 format long g x=2:0.01:3; f=-2*(-72*x.2-24).*(x.2-1)-192*x.2.*(x.2+1)./(x.2-1).5;%四阶导函数 x=2.0; f=-2*(-72*x2-24)*(x2-1)-192*x2*(x2+1)/(x2-1)5; %计算导函数最大值 h4=0.5*10(-7)

15、*180*16/f; h=sqrt(sqrt(abs(h4) %步长 n=1/h; %求分段区间个数 n=2*ceil(1/h)+1 %选取的点数 format long g x=0:0.01:1; f=4*(-72*x.2+24).*(x.2+1)-192*x.2.*(-x.2+1)./(x.2+1).5;%四阶导函数 x=1; f=4*(-72*x2+24)*(x2+1)-192*x2*(-x2+1)/(x2+1)5; %计算导函数最大值 h4=0.5*10(-7)*180*16/f; h=sqrt(sqrt(abs(h4)%步长 n=1/h; %求分段区间个数 n=2*ceil(1/h)

16、+1 %选取的点数 format long g x=0:0.01:1; f=log(3)4*3.x;%四阶导函数 x=1; f=log(3)4*3.x;%计算导函数最大值 h4=0.5*10(-7)*180*16/f; h=sqrt(sqrt(abs(h4)%步长 n=1/h; %求分段区间个数 n=2*ceil(1/h)+1 %选取的点数 format long g x=1:0.01:2; f=4*exp(x)+x.*exp(x);%四阶导函数 plot(x,f) %画出原函数 x=2; f=4*exp(x)+x.*exp(x); %计算导函数最大值 h4=0.5*10(-7)*180*16

17、/f; h=sqrt(sqrt(abs(h4) n=1/h; %求分段区间个数 n=2*ceil(1/h)+1 %选取的点数 估计结果步长h及结点数n分别为 h =0.13411 n =47 h =0.76542 n =35 h =0.18433 n =29 h =0.18546 n =49 2. 积分计算的Matlab程序： format long g promps='请选择积分公式，若用复化梯形，请输入T，用复化simpson，输入S，用复化Gauss_Legendre，输入GL：' result=inputdlg(promps,'charpt 4',1,&

18、#39;T'); Nb=char(result); if(Nb='T'&Nb='S'&Nb='GL') errordlg('积分公式选择错误'); return; end result=inputdlg('请输入积分式题号1-4：','实验4.1',1,'1'); Nb_f=str2num(char(result); if(Nb_f<1|Nb_f>4) errordlg('没有该积分式'); return; end switch N

19、b_f case 1 fun=inline('-2./(x.2-1)');a=2;b=3; case 2 fun=inline('4./(x.2+1)');a=0;b=1; case 3 fun=inline('3.x');a=0;b=1; case 4 fun=inline('x.*exp(x)');a=1;b=2; end if(Nb='T')%用复化梯形公式 promps='请输入用复化梯形公式应取的步长：' result=inputdlg(promps,'实验4.2',1,&

20、#39;0.01'); h=str2num(char(result); if(h<=0) errordlg('请输入正确的步长！'); return; end tic; N=floor(b-a)/h); detsum=0; for i=1:N-1 xk=a+i*h; detsum=detsum+fun(xk); end t=h*(fun(a)+fun(b)+2*detsum)/2; time=toc; t end if(Nb='S')%用复化Simpson公式 promps='请输入用复化Simpson公式应取的步长：' resul

21、t=inputdlg(promps,'实验4.2',1,'0.01'); h=str2num(char(result); if(h<=0) errordlg('请输入正确的步长！'); return; end tic; N=floor(b-a)/h); detsum_1=0; detsum_2=0; for i=1:N-1 xk_1=a+i*h; detsum_1=detsum_1+fun(xk_1); end for i=1:N xk_2=a+h*(2*i-1)/2; detsum_2=detsum_2+fun(xk_2); end t=

22、h*(fun(a)+fun(b)+2*detsum_1+4*detsum_2)/6; time=toc; t end if(Nb='GL')%用复化Gauss_Legendre I %先根据复化Gauss_Legendre I公式的余项估计步长 promps='请输入用复化Gauss_Legendre I 公式应取的步长：' result=inputdlg(promps,'实验4.2',1,'0.01'); h=str2num(char(result); if(h<=0) errordlg('请输入正确的步长!

23、9;); return; end tic; N=floor(b-a)/h);t=0; for k=0:N-1 xk=a+k*h+h/2; t=t+fun(xk-h/(2*sqrt(3)+fun(xk+h/(2*sqrt(3); end t=t*h/2; time=toc; t end switch Nb_f case 1 disp('精确解：ln2-ln3=-0.4054651081') disp('绝对误差：',num2str(abs(t+0.4054651081); disp('运行时间：',num2str(time); case 2 dis

24、p('精确解：pi=3.979') disp('绝对误差：',num2str(abs(t-pi); disp('运行时间：',num2str(time); case 3 disp('精确解：2/ln3=1.368') disp('绝对误差：',num2str(abs(t-1.368); disp('运行时间：',num2str(time); case 4 disp('精确解：e2=7.065') disp('绝对误差：',num2str(abs(t-7.065); d

25、isp('运行时间：',num2str(time); end1. 当选用复化梯形公式时： （1）式运行结果为： t =-0.351 精确解：ln2-ln3=-0.4054651081 绝对误差：1.3944e-008 运行时间：0.003 （2）式运行结果为： t =3.336 精确解：pi=3.979 绝对误差：3.9736e-008 运行时间：0.005 （3）式运行结果为： t = 1.861 精确解：2/ln3=1.368 绝对误差：4.3655e-008 运行时间：0.016 （4）式运行结果为： t =7.610 精确解：e2=7.065 绝对误差：2.0775e-

26、008 运行时间：0.007 2. 当选用复化Simpson公式进行计算时： （1）式运行结果为： t =-0.7519 精确解：ln2-ln3=-0.4054651081 绝对误差：2.7519e-011 运行时间：0.022 （2）式运行结果为： t =3.979 精确解：pi=3.979 绝对误差：0 运行时间：0.021 （3）式运行结果为： t =1.288 精确解：2/ln3=1.368 绝对误差：9.2018e-012 运行时间：0.019 （4）式运行结果为： t =7.118 精确解：e2=7.065 绝对误差：9.0528e-011 运行时间：0.021 3. 当选用复化G

27、auss-Legendre I型公式进行计算时： （1）式运行结果为： t =-0.5262 精确解：ln2-ln3=-0.4054651081 绝对误差：4.7385e-012 运行时间：0.023 （2）式运行结果为： t =3.979 精确解：pi=3.979 绝对误差：1.3323e-015 运行时间：0.021 （3）式运行结果为： t =1.754 精确解：2/ln3=1.368 绝对误差：6.1431e-012 运行时间：0.019 （4）式运行结果为： t =7.046 精确解：e2=7.065 绝对误差：1.441e-012 运行时间：0.021 结果分析： 当选用复化梯形公

28、式时，对步长的事前估计所要求的步长很小，选取的节点很多，误差绝对限要达到时，对不同的函数n 的取值需达到1000-10000之间，计算量是很大。 用复化simpson公式对步长的事前估计所要求的步长相对大些，选取的节点较少，误差绝对限要达到时，对不同的函数n的取值只需在10-100之间，计算量相对小了很多，可满足用较少的节点达到较高的精度，比复化梯形公式的计算量小了很多。用复化simpson公式计算所得的结果比用复化梯形公式计算所得的结果精度高很多，而且计算量小。 当选用Gauss-Lagrange I型公式进行计算时，选用较少的节点就可以达到很高的精度。实验5.1 常微分方程性态和R-K法稳

29、定性试验实验目的：考察下面微分方程右端项中函数y前面的参数对方程性态的影响（它可使方程为好条件的或坏条件的）和研究计算步长对R-K法计算稳定性的影响。实验容及要求： 实验题目：常微分方程初值问题其中，。其精确解为。实验要求：本实验题都用4阶经典R-K法计算。（1）对参数a分别取4个不同的数值：一个大的正值，一个小的正值，一个绝对值小的负值和一个绝对值大的负值。取步长h=0.01，分别用经典的R-K法计算，将四组计算结果画在同一图上，进行比较并说明相应初值问题的性态。（2）取参数a为一个绝对值不大的负值和两个计算步长，一个步长使参数ah在经典R-K法的稳定域，另一个步长在经典R-K法的稳定域外。

30、分别用经典R-K法计算并比较计算结果。取全域等距的10个点上的计算值，列表说明。实验程序：Matlab程序如下：function charp5RK%数值试验5.1：常微分方程性态和R-K法稳定性试验%输入：参数a，步长h%输出：精确解和数值解图形对比%clf;result=inputdlg('请输入-50，50间的参数a:','实验5.1',1,'-40');a=str2num(char(result); if (a<-50|a>50) errordlg('请输入正确的参数a!'); return;endresult=inputdlg('请输入（0 1）之间的步长:','实验5.1'