2022年数值分析实验报告07158.pdf

《2022年数值分析实验报告07158.pdf》由会员分享，可在线阅读，更多相关《2022年数值分析实验报告07158.pdf（19页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、与结果分析：实验源程序function Chap2Interpolation% 数值实验二： “实验：多项式插值的震荡现象”% 输入：函数式选择，插值结点数% 输出：拟合函数及原函数的图形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;endresul

3、t = inputdlg(请输入插值结点数N:,charpt_2,1,10);Nd = str2num(char(result);精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 1 页，共 19 页 - - - - - - - - - - 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;

4、 case g f=inline(atan(x); a = -5; b= 5;end x0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a: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)-);%-function y=Lagrange(x0, y0, x);n= length(x0); m=length(x);for i=1:m z=x(i); s=; for k=

5、1:n p=; for j=1:n if(j = k) p = p*(z - x0(j)/(x0(k) - x0(j); end精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 2 页，共 19 页 - - - - - - - - - - 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 附近较好地逼近了原来

6、的函数 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) 本身是偶函数的性质相反，因此振荡可能更剧烈。(2) 将原来的f(x)换为其他函数如h(x) 、g(x) ，结果如

7、图所示。其中 h(x), g(x)均定义在 -5 ，5 区间上， h(x)=x/(1+x4) ，g(x)=arctan x。 h(x), n=7 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 3 页，共 19 页 - - - - - - - - - - h(x), n=8 h(x), n=9 h(x), n=10 g(x), n=7 g(x), n=8 g(x), n=9 g(x), n=10分析两个函数的插值图形，可以看出：随着 n 的增大，拉格朗日插值函数在x=0 附近较好地逼近了原来的函数f(x)

8、，但是却在两端 x= -5和 x=5 处出现了很大的振荡现象。并且，仔细分析图形， 可以看出， 当 n 为偶数时， 虽然有振荡， 但振荡的幅度不算太大，n 为奇数时， 其振荡幅度变得很大。原因和上面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) 本身是奇函数的性质相反，因此振荡可能更剧烈。实验多项式最小二乘拟合实验目的：编制以函数

9、 xkk=0,n ；为基的多项式最小二乘拟合程序。实验内容：对表中的数据作三次多项式最小二乘拟合。xiyi精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 4 页，共 19 页 - - - - - - - - - - 取权函数 wi 1，求拟合曲线*0nkkkx中的参数 k、平方误差2，并作离散据xi, yi的拟合函数的图形。实验源程序function Chap3CurveFitting% 数值实验三 : “实验 3.1 ”% 输出 : 原函数及求得的相应插值多项式的函数的图像以及参数alph 和误差 rx

10、0 = -1:2;y0 = ;n = 3; % n为拟合阶次alph = polyfit(x0, y0, n);y = polyval(alph, x0);r = (y0 -y)*(y0 -y); % 平方误差x = -1:2;y = polyval(alph, x);plot(x, y, k-);xlabel(x); ylabel(y0 * and );hold onplot(x0, y0, *)grid on;disp(平方误差 :, num2str(r)disp(参数 alph:, num2str(alph)实验结果平方误差 :参数 alph: 实验精品资料 - - - 欢迎下载 - -

11、 - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 5 页，共 19 页 - - - - - - - - - - 实验目的：复化求积公式计算定积分. 实验题目：数值计算下列各式右端定积分的近似值. 实验要求：（1）若用复化梯形公式、复化Simpson 公式和复化Gauss-Legendre I 型公式做计算，要求绝对误差限为710*21，分别利用它们的余项对每种算法做出步长的事前估计. （2）分别用复化梯形公式，复化Simpson 公式和复化Gauss-Legendre I 型公式作计算 . （3）将计算结果与精确解做比较，并比较各种算法的计算量

12、. 实验程序：1. 事前估计的Matlab 程序如下：（1） 用复化梯形公式进行事前估计的Matlab 程序format long g x=2:3; f=-4*(3*x.2+1)./(x.2-1).3; %二阶导函数%plot(x,f) %画出二阶导函数图像x=; %计算导函数最大值f=-4*(3*x2+1)/(x2-1)3; h2=*10(-7)*12/f; h=sqrt(abs(h2) %步长n=1/h; n=ceil(1/h)+1 %选取的点数format long g x=0:1; f=8.*(3*x.2-1)./(x.2+1).3; %二阶导函数精品资料 - - - 欢迎下载 - -

13、 - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 6 页，共 19 页 - - - - - - - - - - %plot(x,f) %画出二阶导函数图像x=1; %计算导函数最大值f=8.*(3*x.2-1)./(x.2+1).3; h2=*10(-7)*12/f; h=sqrt(abs(h2) %步长n=1/h n=ceil(1/h)+1 %选取的点数format long g x=0:1; f=log(3).*log(3).*3.x; %二阶导函数%plot(x,f); %画出二阶导函数图像x=1; %计算导函数最大值f=log(3)*l

14、og(3)*3x; h2=*10(-7)*12/f; h=sqrt(abs(h2) %步长n=1/h n=ceil(1/h)+1 %选取的点数format long g x=1:2; f=2.*exp(x)+x.*exp(x);%二阶导函数%plot(x,f) %画出二阶导函数图像x=2; %计算导函数最大值f=2.*exp(x)+x.*exp(x); h2=*10(-7)*12/f; h=sqrt(abs(h2) %步长n=1/h n=ceil(1/h)+1 %选取的点数精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - -

15、 - - -第 7 页，共 19 页 - - - - - - - - - - 估计结果步长h 及结点数n 分别为h = n =1793 h = n =1827 h = n =2458 h = n =7020 （2） 用复化 simpson 公式进行事前估计的Matlab 程序format long g x=2:3; f=-2*(-72*x.2-24).*(x.2-1)-192*x.2.*(x.2+1)./(x.2-1).5;%四阶导函数x=; f=-2*(-72*x2-24)*(x2-1)-192*x2*(x2+1)/(x2-1)5; %计算导函数最大值h4=*10(-7)*180*16/f;

16、 h=sqrt(sqrt(abs(h4) %步长n=1/h; % 求分段区间个数n=2*ceil(1/h)+1 %选取的点数format long g x=0: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=*10(-7)*180*16/f; h=sqrt(sqrt(abs(h4)%步长n=1/h; %求分段区间个数精品资料 - - - 欢迎下载 - - - - - - - - - - -

17、 欢迎下载 名师归纳 - - - - - - - - - -第 8 页，共 19 页 - - - - - - - - - - n=2*ceil(1/h)+1 %选取的点数format long g x=0:1; f=log(3)4*3.x;%四阶导函数x=1; f=log(3)4*3.x;%计算导函数最大值h4=*10(-7)*180*16/f; h=sqrt(sqrt(abs(h4)%步长n=1/h; %求分段区间个数n=2*ceil(1/h)+1 %选取的点数format long g x=1:2; f=4*exp(x)+x.*exp(x);%四阶导函数plot(x,f) %画出原函数x=

18、2; f=4*exp(x)+x.*exp(x); %计算导函数最大值h4=*10(-7)*180*16/f; h=sqrt(sqrt(abs(h4) n=1/h; %求分段区间个数n=2*ceil(1/h)+1 %选取的点数估计结果步长h 及结点数n 分别为h = n =47 h = n =35 h = n =29 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 9 页，共 19 页 - - - - - - - - - - h = n =49 2. 积分计算的Matlab 程序：format long g

19、 promps= 请选择积分公式，若用复化梯形，请输入T，用复化simpson ，输入 S，用复化 Gauss_Legendre，输入 GL ：; result=inputdlg(promps,charpt 4,1,T); Nb=char(result); if(Nb=T&Nb=S&Nb=GL) errordlg(积分公式选择错误); return; end result=inputdlg(请输入积分式题号1-4 ：,实验 ,1,1); Nb_f=str2num(char(result); if(Nb_f4) errordlg(没有该积分式); return; end switch Nb_f

20、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; 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 10 页，共 19 页 - - - - - - - - - - end if(Nb=T)%用复化梯形公式 promps=请输入用复化梯形公式应取的步长：; result=

21、inputdlg(promps,实验 ,1,); 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 公式应取的步长：; result=inputdlg(promps,实验

22、,1,); h=str2num(char(result); if(h=0) errordlg(请输入正确的步长！); return; 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 11 页，共 19 页 - - - - - - - - - - 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

23、)/2; detsum_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,实验 ,1,); h=str2num(char(result); if(h=0) errordlg(请输入正确的步长!); return; end

24、tic; 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 12 页，共 19 页 - - - - - - - - - - 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=) disp(绝对误差： ,num2str(abs(t+); disp(

25、运行时间： ,num2str(time); case 2 disp(精确解： pi=) disp(绝对误差： ,num2str(abs(t-pi); disp(运行时间： ,num2str(time); case 3 disp(精确解： 2/ln3=) disp(绝对误差： ,num2str(abs); disp(运行时间： ,num2str(time); case 4 disp(精确解： e2=) disp(绝对误差： ,num2str(abs); disp(运行时间： ,num2str(time); end精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载

26、 名师归纳 - - - - - - - - - -第 13 页，共 19 页 - - - - - - - - - - 1. 当选用复化梯形公式时：（1）式运行结果为：t = 精确解： ln2-ln3= 绝对误差：运行时间：（2）式运行结果为：t = 精确解： pi= 绝对误差：运行时间：（3）式运行结果为：t = 精确解： 2/ln3= 绝对误差：运行时间：（4）式运行结果为：t = 精确解： e2= 绝对误差：运行时间：2. 当选用复化Simpson 公式进行计算时：（1）式运行结果为：t = 精确解： ln2-ln3= 绝对误差：运行时间：精品资料 - - - 欢迎下载 - - - - -

27、 - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 14 页，共 19 页 - - - - - - - - - - （2）式运行结果为：t = 精确解： pi= 绝对误差： 0 运行时间：（3）式运行结果为：t = 精确解： 2/ln3= 绝对误差：运行时间：（4）式运行结果为：t = 精确解： e2= 绝对误差：运行时间：3. 当选用复化Gauss-Legendre I型公式进行计算时：（1）式运行结果为：t = 精确解： ln2-ln3= 绝对误差：运行时间：（2）式运行结果为：t = 精确解： pi= 绝对误差：运行时间：（3）式运行结果为：精品资料

28、- - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 15 页，共 19 页 - - - - - - - - - - t = 精确解： 2/ln3= 绝对误差：运行时间：（4）式运行结果为：t = 精确解： e2= 绝对误差：运行时间：结果分析：当选用复化梯形公式时，对步长的事前估计所要求的步长很小，选取的节点很多，误差绝对限要达到710*21时，对不同的函数n 的取值需达到1000-10000 之间，计算量是很大。用复化 simpson 公式对步长的事前估计所要求的步长相对大些，选取的节点较少，误差绝对限要达到710*

29、21时， 对不同的函数n 的取值只需在10-100 之间，计算量相对小了很多，可满足用较少的节点达到较高的精度，比复化梯形公式的计算量小了很多。用复化simpson公式计算所得的结果比用复化梯形公式计算所得的结果精度高很多，而且计算量小。当选用 Gauss-Lagrange I型公式进行计算时，选用较少的节点就可以达到很高的精度。实验常微分方程性态和R-K法稳定性试验实验目的：考察下面微分方程右端项中函数y 前面的参数对方程性态的影响 （它可使方程为好条件的或坏条件的）和研究计算步长对R-K 法计算稳定性的影响。实验内容及要求：实验题目：常微分方程初值问题1,01,(0)1,yayaxxy其中

30、，5050a。其精确解为( )axy xex。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 16 页，共 19 页 - - - - - - - - - - 实验要求：本实验题都用4 阶经典 R-K 法计算。（1）对参数 a 分别取 4 个不同的数值：一个大的正值，一个小的正值，一个绝对值小的负值和一个绝对值大的负值。取步长h=，分别用经典的R-K 法计算，将四组计算结果画在同一张图上，进行比较并说明相应初值问题的性态。（2）取参数a 为一个绝对值不大的负值和两个计算步长，一个步长使参数ah在经典 R-K

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

32、间的步长:,实验 ,1,);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);精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 17 页，共 19 页 - - - - - - - - - - k2=func(a,x(n)+h/2,y(n)+k1

33、*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(a*x)+x;plot(x,y0,g+);hold on;plot(x,y,b-);xlabel(x);ylabel(y0=exp(ax)+x: + and R-K(x)-);实验步骤及结果分析：1、对参数 a 分别取 4 个不同的数值 3（题中要取大的正值，但a 大于 5 时其它曲线在图中不好显示） ，1，-8，-40，将四组计算结果画在同一张图上。 1 -8 -40观

34、察图像可知： 4阶经典 R-K算法的绝对稳定区域是0785.2h，本实验中afy。因此，当 h=时，278.50a，该方法才稳定。现有5050a，因此 a 的稳定区域为500a。当 a 处于稳定性范围内时，结果收敛，为好条件。图中显示，四个参数下结果均收敛，是好条件的。2、取参数 a 为一个绝对值不大的负值a=-8，故348.00h时，该方法稳定。取 h=和 h=作比较。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 18 页，共 19 页 - - - - - - - - - - h= h=ix0y 精由表中等距的 10 个点的计算值和图像可以看出，h 值超过稳定区域时结果发散，与精确曲线相差很大， 而 h 值在稳定区域内时， 忽略舍入误差， 计算结果和精确值一样。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 19 页，共 19 页 - - - - - - - - - -