数值分析实验报告.docx

数值分析实验报告学号： 专业：力学系 姓名： 实验2.1 多项式插值的振荡现象1. 实验目的：考虑在一个固定的区间上用插值逼近一个函数，显然Lagrange插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，Ln(x)是否也更加靠近被逼近的函数。Runge给出的一个例子是极著名并富有启发性的。2. 实验内容：设区间-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插值基函数。3. 实验步骤与结果分析： 3.1 源程序代码function t_charpt2% 数值实验二：“实验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,10);Nd = str2num(char(result);if(Nd1)errordlg(结点输入错误！);return;endswitchNb_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;endx0 = 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)-);%-function y=Lagrange(x0, y0, x);% Lagrange插值n= length(x0); m=length(x);fori=1:m z=x(i); s=0.0;for k=1:n p=1.0;for j=1:nif(j = k) p = p*(z - x0(j)/(x0(k) - x0(j);endend s = s + p*y0(k);endy(i) = s;end3.2 实验结果分析 （1）增大分点n=2，3，时，拉格朗日插值函数曲线如图所示。 n=3 n=4n=7n=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)本身是偶函数的性质相反，因此振荡可能更剧烈。（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=7 g(x), n=8g(x), n=9 g(x), n=10分析两个函数的插值图形，可以看出：随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x)，但是却在两端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)本身是奇函数的性质相反，因此振荡可能更剧烈。实验2.2 样条插值的收敛1. 实验目的：多项式插值是不收敛的，即插值的节点多，效果不一定好。样条函数插值能较好地解决这个问题，通过本实验可以验证这一理论结果。2. 实验内容：通过对实验2.1中的三个函数进行三次样条插值，与Lagrange插值进行对比，验证样条插值的收敛性。3. 实验步骤及内容分析： 3.1源程序代码：function t_charpt2% 数值实验二： “实验2.2：样条插值的收敛”% 输入：函数式选择，插值结点数% 输出：拟合函数及原函数的图形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,10);Nd = str2num(char(result);if(Nd1)errordlg(结点输入错误！);return;endswitchNb_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;endx0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b;cs = spline(x0, y0); y = ppval(cs, x);plot(x0, y0, o); hold on; plot(x, y, k-);xlabel(x); ylabel(y = f(x) o and y = Spline(x)-); 3.2实验结果分析：f(x), n=5 f(x), n=10f(x), n=20 f(x), n=30h(x), n=5 h(x), n=10h(x), n=20 h(x), n=30g(x), n=5 g(x), n=10g(x), n=20 g(x), n=30从图中可以看出，三次样条插值随着插值结点的增加，由于其采用了分段三次多项式拟合的方法，因此并没有出现振荡现象。样条插值思想最早产生于工业部门，如表，某汽车制造商用三次样条插值设计车门曲线，其中一段数据如表所示。xk012345678910yk0.00.791.532.192.713.033.272.893.063.193.29yk0.80.2functionCarDoorDesignx0=0:10;y0=0.8 ,0.0, 0.79, 1.53, 2.19, 2.71, 3.03, 3.27, 2.89, 3.06, 3.19, 3.29,0.2;x=0:0.1:10;cs=spline(x0,y0);y = ppval(cs, x);plot(x0, y0(2:12), o); hold on;plot(x, y, k-);xlabel(x); ylabel(y = f(x) o and y = Spline(x)-);设计的车门曲线如图所示.实验3.1 多项式最小二乘拟合1. 实验目的：编制以函数xkk=0,n；为基的多项式最小二乘拟合程序。2. 实验内容：对表中的数据作三次多项式最小二乘拟合。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的拟合函数的图形。3. 实验步骤与结果分析：3.1 实验源程序:function charpt3% 数值实验三:“实验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(x); ylabel(y0 * and polyfit.y-);hold onplot(x0, y0, *)grid on;disp(平方误差:, num2str(r)disp(参数alph:, num2str(alph)3.2 实验结果分析:平方误差:2.1762e-005参数alph:1.9991 -2.9977 -3.9683e-005 0.54912实验3.2 正交化多项式的最小二乘拟合1. 实验目的：编制正交化多项式最小二乘拟合程序。2. 实验内容：用正交化多项式最小二乘拟合求解3.1中的三次多项式最小二乘拟合问题，作拟合曲线的图形，计算平方误差。3. 实验步骤及结果分析：3.1 实验源程序:function charpt3% 数值实验三:“实验3.2”% 子函数调用:dlsa% 输出:原函数及求得的相应插值多项式的函数的图像以及参数alph和误差rx0 = -1:0.5:2;y0 = -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552; result = inputdlg(请输入权矢量w:,charpt_3, 1, 1 1 1 1 1 1 1); w = str2num(char(result); n = 3; % n为拟合阶次 a, b, c, alph, r = dlsa(x0, y0, w, n);disp(平方误差:, num2str(r)disp(参数alph:, num2str(alph)%-function a, b, c, alph, r = dlsa(x, y, w, n)% 功能：用正交化方法对离散数据作多项式最小二乘拟合m = length(x) - 1;if(n = m)errordlg(错误: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; end s3 = (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% 求平方误差r=y.*y*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+1fori=3:k+1 T(k,i)=T(k-1,i)-a(k-1)*T(k-1,i-1)-b(k-1)*T(k-2,i-2);endendendfori=1:n+1for k=i:n+1alph(n+2-i)=alph(n+2-i)+c(k)*T(k,k+2-i);endend% 用秦九韶方法计算s(t)的输出序列(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+1 s=s.*t+alph(k);end% 输出点列x-y和拟合曲线t-s的图形plot(x,y,*,t,s,-);title(离散数据的多项式拟合);xlabel(x); ylabel(y);grid on;3.2 实验结果分析权矢量取1 1 1 1 1 1 1平方误差:2.1762e-005参数alph:1.9991 -2.9977 -3.9683e-005 0.54912与3.1相比，拟合系数和误差完全相同，拟合曲线也完全一致。通过比较可以看出三种方法得到的结果十分的相似，差别很小（表中舍去了一些位数，所以看起来没有任何差别，实际上还是有十分微小的差别，几乎可以忽略不计）。说明这三种方法均是对多项式进行最小二乘拟合的有效方法。表1-2 最小二乘拟合得到的拟合曲线中各参数的值平方误差多项式拟合函数法1.9991111-2.9976666-0.00003960.54911902.1761904E-005求解系数矩阵法1.9991111-2.9976666-0.00003960.54911902.1761904E-005正交多项式法1.9991111-2.9976666-0.00003960.54911902.1761904E-005注：由于表的大小所限，各类参数均取到小数点后七位。实验4.1复化求积公式计算定积分1 实验目的：复化求积公式计算定积分2 实验内容：数值计算下列各式右端定积分的近似值（1）（2）（3）（4）3 实验要求：（1） 若用复化梯形公式，复化Simpson公式和复化Gauss-Legendre I型公式做计算，要求绝对误差限为，分别利用它们的余项对每种算法做出步长的事前估计。（2） 分别用复化梯形公式，复化Simpson公式和复化Gauss-Legendre I型公式做计算。（3） 将计算结果与精确解做比较，并比较各种算法的计算量。4 源程序代码：function charpter4_1%数值实验4 .1 复化求积公式计算定积分clcformat longresult=inputdlg(请选择积分法，若用复化梯形，输入T，用复化Simpson法，输入S，用复化Gauss_Legendre，输入GL:,charpt 4,1,S);Nb=char(result);if(Nb=T&Nb=S&Nb=GL)errordlg(积分公式选择错误！);return;endresult=inputdlg(请输入积分式题号14：,实验4.1,1,1);Nb_f=str2num(char(result);switchNb_fcase 1fun=inline(-2./(x.2-1);a=2;b=3;case 2fun=inline(4./(x.2+1);a=0;b=1;case 3fun=inline(3.x);a=0;b=1;case 4fun=inline(x.*exp(x);a=1;b=2;otherwisedisp(没有该积分式，请选择1,2,3,4);return;end if(Nb=T) %复化梯形求积公式其中a，b分别为积分上下限，h为步长，n为分步计数器，y为积分值N=0;tol=0.5e-7;h=(b-a)/2N;y0=(fun(a)+fun(b)*(b-a)/2;y=0;while(abs(y0-y)=3*tol)N=N+1;h=(b-a)/2N;y=y0;y0=y/2+h*sum(fun(a+h:2*h:b-h);endy=y0n=2Nelseif(Nb=S) %复化辛普森求积公式其中a，b分别为积分上下限，h为步长，n为分步计数器，y为积分值 N=0;tol=0.5e-7;h=(b-a)/2N;y0=(fun(a)+fun(b)+4*fun(a+b)/2)*(b-a)/6;y=0;while(abs(y0-y)=15*tol)N=N+1;h=(b-a)/2N;y=y0;s1=(fun(a)+fun(b)*h/6;s2=sum(fun(a+h/2.0:h:b-h/2.0)*4.0*h/6.0;s3=sum(fun(a+h:h:b-h)*2.0*h/6.0;y0=s1+s2+s3;endy=y0n=2Nelseif(Nb=GL) %复化高斯勒让德求积公式其中a，b分别为积分上下限，h为步长，n为分步计数器，y为积分值 N=0;tol=0.5e-7;h=(b-a)/2N;b1=-1/sqrt(3);b2=1/sqrt(3);a1=1;a2=1;y0=(a1*fun(a+b)/2.0+(b-a)/2.0*b1)+a2*fun(a+b)/2.0+(b-a)/2.0*b2)*(b-a)/2;y=0;while(abs(y0-y)=tol)N=N+1;h=(b-a)/2N;y=y0;sum2=0;fori=1:2Ns=0.5*h*(a1*fun(2*a+(2*i-1)*h)/2.0+h/2.0*b1)+a2*fun(2*a+(2*i-1)*h)/2.0+h/2.0*b2);sum2=sum2+s;endy0=sum2;endy=y0n=2Nelseenddisp(实验题4.1 (,num2str(Nb_f),)的计算结果：,num2str(y);switchNb_fcase 1disp( 精确解：ln2-ln3=-0.40546510810816)disp( 绝对误差：,num2str(abs(y-log(2)+log(3);case 2disp( 精确解：pi=3.14159265358979)disp( 绝对误差：,num2str(abs(y-pi);case 3disp( 精确解：2/ln3=1.82047845325368)disp( 绝对误差：,num2str(abs(y-2/log(3);case 4disp( 精确解：e2=7.38905609893065)disp( 绝对误差：,num2str(abs(y-exp(2);endend实验数据及分析：（1） 第一问要求可知复化梯形公式的余项表达式为：（1）复化辛普森公式的余项表达式为：（2）单区间的Gauss-Legendre I型公式的余项表达式为：（3）再将函数1-4分别求二阶导函数和和四阶导函数，在在响应区间内取其绝对值的最大值，最后带入（1）式和（2）式，求得满足条件下复化梯形公式和复化辛普森公式的最大步长（由于不清楚复化Gauss-Legendre I型公式的整体余项表达式，（3）式先不考虑）（1） 函数1：的二阶导函数为： （4） 四阶导函数为：（5）（2） 函数2：的二阶导函数为：（6）四阶导函数为：（7）（3） 函数3：的二阶导函数为：（8）四阶导函数为：（9）（4） 函数4：的二阶导函数为：（10）四阶导函数为：（11）表2-1 复化梯形法的步长事前估计函数1函数2函数3函数4二阶导函数的最大值1.92598.00003.620829.5562最大步长0.5582e-30.2739e-30.4071e-3 0.1425e-3最小区间数1792365224577019表2-2 复化Simpson法的步长事前估计函数1函数2函数3函数4四阶导函数的最大值23.901296.00004.370244.3343最大步长0.04950.03450.07580.0425最小区间数21291324（2） 程序计算及结果分析表2-3 定积分为的数值解份数步长数值解精确解绝对误差复化梯形法20480.4883e-3-0.40546512204351-0.405465108108161.3935e-008复化Simpson320.0313-0.40546513756984-0.405465108108162.9462e-008复化G-L I型320.0313-0.40546510687839-0.405465108108161.2298e-009表2-4 定积分为的数值解份数步长数值解精确解绝对误差复化梯形法20480.4883e-33.141592613853363.141592653589793.9736e-008复化Simpson160.06253.141592651224823.141592653589792.365e-009复化G-L I型160.06253.141592653616073.141592653589792.6281e-011表2-5 定积分为的数值解份数步长数值解精确解绝对误差复化梯形法20480.4883e-31.820478496908611.820478453253684.3655e-008复化Simpson160.06251.820478467302181.820478453253681.4049e-008复化G-L I型320.03131.820478452668261.820478453253685.8541e-010表2-6 定积分为的数值解份数步长数值解精确解绝对误差复化梯形法81920.1221e-37.389056119706107.389056098930652.0775e-008复化Simpson320.03137.389056107563767.389056098930658.6331e-009复化G-L I型640.01567.389056098570937.389056098930653.5972e-010实验分析1. 从事前步长估计和事后实际步长的比较可以看出步长的事前估计一般都要大于实际计算时满足条件的临界步长，说明程序的编写和运行是合理的。（由于程序是利用不断二分区间的形式减小步长，所以最后的得到的区间均为2的倍数，与事先估计的区间数有一点差别，但一般计算精度也比0.5e-7要高）。2. 比较复化梯形法和复化Simpson法的步长和绝度误差可以看出，在绝对误差基本相同的时候，复化Simpson法的程序收敛速度要快很多（具体表现在步长较大，比复化梯形法一般大两个量级），选择该算法可以减少不少计算量。3. 利用复化Gauss-Legendre I型公式进行计算，计算量和复化Simpson法基本相同，但通过绝对误差的比较可以看出，Gauss-Legendre I型公式法的精度较高，因此，若在精度相同的情况下，复化Gauss-Legendre I型公式的计算量应该是最小的，复化Simpson法次之，复化梯形法计算量最大。4. 由于不知道复化Gauss-Legendre I型公式的余项递推式，所以程序中未能有效的将其绝对误差逼近绝对误差限，导致精度的增加，步长的减小，造成不能直接与复化Simpson法进行比较，这方面仍需要改进。实验5.1常微分方程性态和R_K法稳定性试验1 实验目的：常微分方程的初值问题其中，其精确解为。2 实验要求：本实验题都用4阶经典R_K法计算。（1）对参数，分别取四个不同的值：一个大的正值，一个小的正值，一个绝对值小的负值和一个绝对值大的负值。步长取为h=0.01，分别用经典R_ K法计算，将四组计算结果画在同一张图上，进行比较并说明相应初值问题的性态。（2）取参数为一个绝对值不大的负数和两个计算步长，一个步长使参数在经典R-K法的稳定区域里，另一个步长在经典R-K法的稳定域外，分别用经典R-K法计算并比较计算结果。取全域等距的10个点上的计算值，列表说明。3源程序代码：function charpt5_1%数值实验5.1：常微分方程性态和R-K法稳定性试验%输入：参数a，步长h%输出：精确解和数值解图形对比clcresult=inputdlg(请输入-50 50间的参数a:,实验5.1,1,-1);a=str2num(char(result);if(a50)errordlg(请输入正确的参数a！);return;endresult=inputdlg(请输入（0,1）之间的份数（步长的倒数）：,实验5.1,1,100);n=str2num(char(result);h=1.0/n;if(h1)errordlg(请输入正确的步长h！);return;endfun=inline(1+(y-x).*a);x=0:h:1;y=x;y(1)=1;fori=1:n k1=fun(a,x(i),y(i); k2=fun(a,x(i)+h/2,y(i)+k1*h/2); k3=fun(a,x(i)+h/2,y(i)+k2*h/2); k4=fun(a,x(i)+h,y(i)+k3*h);y(i+1)=y(i)+h*(k1+2*k2+3*k3+k4)/6;endy0=exp(a*x)+x;plot(x,y0,r*);hold on;plot(x,y,kd);hold on;xlabel(x);ylabel(y0=exp(ax)+x: + and R_K(x)-);实验分析：（3） 第一问中对于参数的讨论，本实验中分别取=-20，=-1，=1和=20来进行计算作图和讨论：图中线条代表的含义：=-20的函数曲线用红色的线条表示=-1 的函数曲线用绿色的线条表示=-1 的函数曲线用蓝色的线条表示=20的函数曲线用黑色的线条表示两条线之中*点线表示的是精确解点线表示的是R-K法数值解 图3-1 =-20，=-1，=1三条曲线的比较由于=20的函数曲线在坐标轴上的跨度远远大于其他三组结果，四线放在同一幅图上会严重影响其他组数据的可视化比较，所以仅将其与=1曲线单独比较，从纵坐标跨度上我们可以清晰的了解到其发散程度远远超过其他三组数据。因此图1为=-20，=-1，=1三条曲线的比较，图2为=1，=20两条曲线的比较。图3-2 =20，=1，两条曲线的比较分析结论：1.从图1的结果来看当0的时候，该算法是不稳定的，也是不收敛的。2.从图2的对比结果来看当0