数值分析实验报告.docx

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 Chap2Inter

2、polation% 数值实验二：“实验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)errordlg(实验函数选择错误！);return;endresult = inputdlg(请输入插值结点数N:,charpt_2,1

3、,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./(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

4、, a b, co); hold on; plot(x, y, b-); xlabel(x); ylabel(y = f(x) o and y = Ln(x)-);%-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，时，拉

5、格朗日插值函数曲线如图所示。 n=6 n=7 n=8 n=9 n=10 从图中可以看出，随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(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)本身是偶函数

6、的性质相反，因此振荡可能更剧烈。(2)将原来的f(x)换为其他函数如h(x)、g(x)，结果如图所示。其中h(x), g(x)均定义在-5，5区间上，h(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

7、为奇数时，其振荡幅度变得很大。原因和上面f(x)的插值类似，h(x)、g(x)本身是奇函数，如果n为偶数，那么Lagrange插值函数Ln(x)的最高次项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

8、.048-0.4470.5494.552取权函数wi1，求拟合曲线中的参数k、平方误差2，并作离散据xi, yi的拟合函数的图形。实验源程序function 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); % 平方

9、误差x = -1:0.01:2;y = polyval(alph, x);plot(x, y, k-);xlabel(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）若用复化梯形公式、复化

10、Simpson公式和复化Gauss-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; %计算导函数最大值

11、 f=-4*(3*x2+1)/(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

12、.01:1; f=log(3).*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=sq

13、rt(abs(h2) %步长 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; %计算导函数最大值

14、 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=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; %求分段区间个数

15、 n=2*ceil(1/h)+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

16、*10(-7)*180*16/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,T); Nb

17、=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_f4) errordlg(没有该积分式); return; end switch Nb_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

18、 4 fun=inline(x.*exp(x);a=1;b=2; end if(Nb=T)%用复化梯形公式 promps=请输入用复化梯形公式应取的步长：; result=inputdlg(promps,实验4.2,1,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=

19、toc; t end if(Nb=S)%用复化Simpson公式 promps=请输入用复化Simpson公式应取的步长：; result=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; d

20、etsum_2=detsum_2+fun(xk_2); end t=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(请输入正确的步长!); return; end ti

21、c; 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 disp(精确解：pi=3.979) disp(绝对误差：,num2str(abs(t-pi); dis

22、p(运行时间：,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); disp(运行时间：,num2str(time); end1. 当选用复化梯形公式时： （1）式运行结果为： t =-0.351 精确解：ln2-ln3=-0.4054651081 绝对误差：1.3944e-008 运行时间：0.003 （2）式运行结果为：

23、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-008 运行时间：0.007 2. 当选用复化Simpson公式进行计算时： （1）式运行结果为： t =-0.7519 精确解：ln2-ln3=-0.4054651081 绝对误差：2.7519e-011 运行时间：0.022 （2）式运行结果为： t =3.979 精确解：

24、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. 当选用复化Gauss-Legendre I型公式进行计算时： （1）式运行结果为： t =-0.5262 精确解：ln2-ln3=-0.4054651081 绝对误差：4.7385e-012 运行时间：0.023 （2）式运行结果为： t =3.979 精确解：pi=3.979 绝对误差：

25、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 结果分析： 当选用复化梯形公式时，对步长的事前估计所要求的步长很小，选取的节点很多，误差绝对限要达到时，对不同的函数n 的取值需达到1000-10000之间，计算量是很大。 用复化simpson公式对步长的事前估计所要求的步长相对大些，选取的节点较少，误差绝对限要达到时，对不同的函数n的取值只需在10-1

26、00之间，计算量相对小了很多，可满足用较少的节点达到较高的精度，比复化梯形公式的计算量小了很多。用复化simpson公式计算所得的结果比用复化梯形公式计算所得的结果精度高很多，而且计算量小。 当选用Gauss-Lagrange I型公式进行计算时，选用较少的节点就可以达到很高的精度。实验5.1 常微分方程性态和R-K法稳定性试验实验目的：考察下面微分方程右端项中函数y前面的参数对方程性态的影响（它可使方程为好条件的或坏条件的）和研究计算步长对R-K法计算稳定性的影响。实验内容及要求： 实验题目：常微分方程初值问题其中，。其精确解为。实验要求：本实验题都用4阶经典R-K法计算。（1）对参数a分别

27、取4个不同的数值：一个大的正值，一个小的正值，一个绝对值小的负值和一个绝对值大的负值。取步长h=0.01，分别用经典的R-K法计算，将四组计算结果画在同一X图上，进行比较并说明相应初值问题的性态。（2）取参数a为一个绝对值不大的负值和两个计算步长，一个步长使参数ah在经典R-K法的稳定域内，另一个步长在经典R-K法的稳定域外。分别用经典R-K法计算并比较计算结果。取全域等距的10个点上的计算值，列表说明。实验程序：Matlab程序如下：function charp5RK%数值试验5.1：常微分方程性态和R-K法稳定性试验%输入：参数a，步长h%输出：精确解和数值解图形对比%clf;result

28、=inputdlg(请输入-50，50间的参数a:,实验5.1,1,-40);a=str2num(char(result); if (a50) errordlg(请输入正确的参数a!); return;endresult=inputdlg(请输入（0 1）之间的步长:,实验5.1,1,0.01);h=str2num(char(result);if (h=1) errordlg(请输入正确的（0 1）间的步长!); return;endx=0:h:1;y=x;N=length(x);y(1)=1;func=inline(1+(y-x).*a);for n=1:N-1 k1=func(a,x(n),y(n); k2=func(a,x(n)+h/2,y(n)+k1*h/2); k3=func(a,x(n)+h/2,y(n)+k2*h/2); k4=func(a,x(n)+h,y(n)+k3*h); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endy0=exp(