数值分析实验报告.doc

数值分析实验报告数值分析实验报告学 院：学生姓名：学号：2016年11月目录实验一 插值法1一 实验内容1二 实验步骤1三 所编程序清单2四 例题验证3实验二 数据拟合5一 实验内容5二 实验步骤5三 难点分析6四 所编算法清单6五 例题验证6实验三 数值积分8一 实验内容8二 实验步骤8三 难点分析8四 所编算法清单9五 例题验证10实验四 线性方程组的直接求法之高斯消去法12一 实验内容12二 实验步骤12三 难点分析13四 所编算法清单13五 例题验证13实验五 线性方程组的迭代解法之雅可比迭代法15一 实验内容15二 实验步骤15三 难点分析16四 所编算法序清单16四 例题验证17实验六 非线性方程的数值解法之Newton迭代法19一 实验内容19二 实验步骤19三 实验难点19四 所编程序详细清单19五 例题验证20实验七 矩阵特征值问题的数值计算之幂法21一 实验内容21二 实验步骤21三 难点分析21四 具体所编算法清单21五 例题验证22实验八 常微分方程的数值解法之改进的Euler法24一 实验内容24二 实验步骤24三 实验难点24四 具体所编算法24五 例题验证25II实验一 插值法本人选择Newton插值法进行编程。选择的例题为：教材（P25，例2.7）进行验证：题目：设，并已知的数据如表1-1试用2次Newton插值多项式计算计算的近似值，并讨论其误差。表1-1x2.02.12.21.4142141.4491381.483240一 实验内容1. 根据Newton插值多项式 (*)使用Matlab编程实现Newton插值算法。2.使用教材上的例题验证所编算法的正确性。二 实验步骤Newton插值法的关键是求出的k阶均差，然后将k阶均差代入上式(*)，求出插值多项式，最后将所选例题的条件代入Newton插值多项式，与正确答案进行比对，进而得出结论。1 难点分析在编写程序时需要考虑的难点之一是如何确定所选用的Newton插值法的次数，本人拟选用直接从Matlab中m文件中直接人工输入的方法，设一个未知量n作为选用Newton插值法的次数确定选用几次的Newton插值多项式。在编写程序时的另一个难点是经过Newton插值公式计算之后，输出是什么？是直接列写多项式，还是以系数矩阵的形式输出？本人考虑到编程的简便性，拟直接以均差表形式输出。2 具体步骤详细步骤如下：(1)已知k阶均差的求解公式如下： k=1,2,.; (1) 通过上式分别求取一阶均差、二阶均差、.、k阶均差，并保存在与函数自变量x同阶的同一个大矩阵中。在大矩阵中，第一列用来存放自变量x，第二列用来存放因变量，第三列用来存放一阶均差，.，第（k+2）列用来存放k阶均差。具体矩阵表示如下表1-2所示表1-2一阶均差二阶均差三阶均差.注：表1-2未标示出的元素都以0填充。(2)根据公式(*)求解n次(n为人为设定的Newton插值多项式系数)Newton插值多项式，并且在Matlab中Command Window输出该插值多项式。(3)在教材中随机选取某一例题来验证所编程序的正确性。三 所编程序清单所编程序对应的文件为Newton.m%Newton插值法function y=Newton(Xi,Yi,x,n)%-%计算均差表%根据教材，均差表一般列举到三阶均差即可%一阶均差f1=zeros(1,n);%一阶均差数值保存在f1中for i=1:n f1(1,i)=(Yi(1,i+1)-Yi(1,i)/(Xi(1,i+1)-Xi(1,i);enddisp(一阶均差数值为：);disp(f1);%二阶均差 f2=zeros(1,n-1);%二阶均差数值保存在f2中 for i=1:(n-1) f2(1,i)=(f1(1,i+1)-f1(1,i)/(Xi(1,i+2)-Xi(1,i); enddisp(二阶均差数值为：);disp(f2);%三阶均差if n2f3=zeros(1,n-2);%三阶均差数值保存在f3中for i=1:(n-2) f3(1,i)=(f2(1,i+1)-f2(1,i)/(Xi(1,i+3)-Xi(1,i);enddisp(三阶均差数值为：);disp(f3);end%-%求3次插值多项式的值if n=3y=Yi(1,1)+f1(1,1)*(x-Xi(1,1)+f2(1,1)*(x-Xi(1,1)*(x-Xi(1,2)+f3(1,1)*(x-Xi(1,1)*(x-Xi(1,2)*(x-Xi(1,3);elseif n=2 y=Yi(1,1)+f1(1,1)*(x-Xi(1,1)+f2(1,1)*(x-Xi(1,1)*(x-Xi(1,2);elseif n=1 y=Yi(1,1)+f1(1,1)*(x-Xi(1,1);endend四 例题验证选教材（P25，例2.7进行验证）题目：设，并已知的数据如表1-3.试用2次Newton插值多项式计算计算的近似值，并讨论其误差。表1-3x2.02.12.21.4142141.4491381.483240编写m文件如下：(对应文件为t1.m)clear;clc;Xi=2.0 2.1 2.2;Yi=1.414214 1.449138 1.483240;x=2.15;n=2;result=Newton(Xi,Yi,x,n);disp(所求近似值为：)disp(result);输出结果为经与教材例题所给答案比对，完全正确。实验二 数据拟合本人选择多项式拟合进行编程。选教材（p62，例3.10）作为验证例题，题目如下：题目（P62,例3.10）用多项式拟合表2-1中的离散数据表2-1i01234xi0.000.250.500.751.00yi0.100.350.811.091.96一 实验内容 1.所谓多项式拟合是根据多项式： (1)做曲线拟合，其中为基函数。2.找教材上例题进行验证，根据程序仿真结果判断所编算法的正确性。二 实验步骤1. 设定子空间的基函数为，设定权数（为了使运算简便，不妨都设为1）。2. 计算克莱姆矩阵 (2)3. 根据下式计算系数 j=0,1,.,n (3)将代入(1)即可求得待求拟合模型多项式。4.根据3所求的多项式画出函数曲线三 难点分析1. 克莱姆矩阵的求解，需要使用嵌套一次的循环，容易出错。四 所编算法清单所编算法对应的文件为duoxiangshinihe.mfunction ai=duoxiangshinihe(i,Xi,Yi,N)wi=ones(1,i);%定义权值矩阵，并不防设权值均为1ai=zeros(N+1,1);Cramer=zeros(N+1,N+1);for ii=1:(N+1) %列 for iii=1:(N+1)%行 Cramer(ii,iii)=0;%克莱姆矩阵元素初始化 for iiii=1:i Cramer(ii,iii)=Cramer(ii,iii)+Xi(1,iiii)(ii+iii-2); end endendYFAI=zeros(N+1,1); %求(y,)矩阵for iiiii=1:(N+1) for iiiiii=1:i YFAI(iiiii,1)=YFAI(iiiii,1)+Yi(1,iiiiii)*Xi(1,iiiiii)(iiiii-1); endend ai=inv(Cramer)*YFAI;%由(a)=inv()*(y,)求系数矩阵 Ai=ai; Aii=rot90(Ai);%以上两行的作用是将（a0,a1,a2,a3)改成(a3,a2,a1,a0) x1=0:1; axis(0 1 0 2); plot(x1,polyval(Aii,x1),r); hold on plot(Xi,Yi,*); end五 例题验证选教材（p62，例3.10）作为验证例题，题目如下：题目（P62,例3.10）用多项式拟合表2-2中的离散数据表2-2i01234xi0.000.250.500.751.00yi0.100.350.811.091.96编写例题m文件如下：（即t2.m文件）I=5;%离散点个数xi=0.00 0.25 0.50 0.75 1.00;%离散点横坐标yi=0.10 0.35 0.81 1.09 1.96;%离散点纵坐标N=2;%选用2次多项式拟合AI=duoxiangshinihe(I,xi,yi,N);disp(采用多项式拟合出来的多项式系数为：（从低次到高次）);disp(AI);仿真结果为：图2-1 拟合曲线经与教材例题答案比对可知，所编算法完全正确。实验三 数值积分本人选择复化梯形求积公式（Complex_Quadrature）。选教材题目（P79，例4.6）作为验证例题，题目如下：题目（P79，例4.6）分别用复化梯形公式和复化Simpson公式计算时，要使误差不超过，问各取多少个点。一 实验内容1. 通过复化梯形求积公式 (1)求解积分。2. 通过梯形公式的误差公式 (2)求解求积的误差并判断，当误差不超过允许值时即可认为求积结束。3. 找教材上的例题进行验证，验证所编程序的正确性。二 实验步骤1. 首先要对区间a,b进行n等分(n1),并且对都要对(2)进行判断误差是否已经不超过允许值，如果超过了，则对区间a,b进行(n+1)等分(记作n=n+1)，然后再根据(2)判断误差是否不超过允许值，依次类推，直至找到能使误差不超过允许值时的n.2. 在(1)的基础上，根据公式(1)求积，得出计算结果并输出。3. 找教材上例题进行验证。三 难点分析1. 因为误差需要和给定的误差要求比较大小，所以的取值范围是算法编程过程的一个难点。本人的思路是先求的极值点（假设在给定区间上只有一个极值点），然后求边界点的函数值以及极值点处的函数值分别为，通过比较这三个值的大小就可以找出函数的取值范围。2. 在确定步长的基础上，通过公式(*)求积。通过分析可知，公式(*)可分为常量和累加量两部分，我的编程思路是先给赋初值，然后再采用循环加上累加量。四 所编算法清单所编算法对应的文件为Complex_Quadrature.mfunction T nt=Complex_Quadrature(f,a,b,deta)%T复化梯形求积结果，nt为运算次数%f待积分函数表达式，a,b待积分区间，nt表示复化梯形求积公式满足误差要求时需要将a,b区间nt等分，ns表示复化Simpson求积公式满足%误差要求时需要将a,b区间ns等分,deta表示积分时误差要求%首先复化梯形求积公式%-syms x;y=f(x);ddy=diff(y,x,2);dddy=diff(y,x,3);x0=solve(dddy,x,0);%求取函数f(x)的二阶导数极点 （为便于简便，不妨假设f(x)二阶导数只有一个极点）t=double(x0);ddfa=subs(ddy,x,a);%f(a)ddfb=subs(ddy,x,b);%f(b)ddfx0=subs(ddy,x,t);%f(x0)ta=double(ddfa);tb=double(ddfb);%转换为常量tx0=double(ddfx0);%判断f(a),f(b),f(x0)值的最大值和最小值deta1=ta-tb;deta2=tb-tx0;deta3=ta-tx0; if deta1=0&deta2=0 tbiggest=ta; tsmallest=tx0;elseif deta1=0&deta2=0 tbiggest=ta; tsmallest=tb;elseif deta1=0&deta2=0&deta3=0 tbiggest=tx0; tsmallest=tb;end if deta1=0&deta2=0 tbiggest=tx0; tsmallest=ta;elseif deta1=0&deta3=0 tbiggest=tb; tsmallest=tx0;elseif deta1=0&deta3=0 tbiggest=tb; tsmallest=ta;end for i=2:1000 %从2等分、3等分、4等分n等分依次逐个验证误差精度 h=(b-a)/i;%计算步长 RTmax=(-1)*(b-a)/12)*h2*tsmallest; if RTmax=deta nt=i; T=(h/2)*(f(a)+f(b); for ii=1:i T=T+h*f(a+h*(ii-1); end break else continue endend ns=0;S=0; end 五 例题验证选教材题目（P79，例4.6）作为验证例题，题目如下：题目（P79，例4.6）分别用复化梯形公式和复化Simpson公式计算时，要使误差不超过，问各取多少个点。根据所选例题进行编程，如下：(对应文件为t3.m)syms x;f=inline(sin(x),x);a=0;b=pi;deta=2*10(-5);T nt ns S=Complex_Quadrature(f,a,b,deta);disp(T);disp(nt);仿真结果： 经过与教材例题所给答案进行比对，结果一致，所以所编程序完全正确。实验四 线性方程组的直接求法之高斯消去法本人选择高斯消去法。选取教材上例题（P113，例5.3）进行验证，例题如下：用Gauss消去法解方程组 一 实验内容 1.应用Guass消去法编写求解线性方程的直接求法，其公式推导过程如下： (1) (2) i=k+1,k+2,.n; j=k+1,k+2,.,n+1.2. 找教材上习题进行验证所编程序的正确性。二 实验步骤1. 通过公式（1）和公式（2）将系数矩阵化为上三角阵，增广矩阵的格式如式（3）所示 (3)2. 通过反向逆推法，依次求得，其求解公式如下： (4) i=n-1,n-2,.1. (5)3. 用教材例题（P113，例5.3）对程序进行验证，以验证程序的正确性。三 难点分析1.式（5）里面含有累加过程，在设置循环时容易出错。四 所编算法清单所编算法对应的文件为Gauss_dc.m%线性方程直接接法之高斯消去法function AT X=Gauss_dc(A) a,b=size(A);%首先找出增广矩阵的行数和列数for k=1:(a-1)%总共需要（a-1）次消元for i=(k+1):a L=A(i,k)/A(k,k); for j=1:b A(i,j)=A(i,j)-L*A(k,j); endendendX=zeros(1,a);X(1,a)=A(a,b)/A(a,b-1);for m=1:(a-1)%X(1,a-1)=(A(a-1,b)-A(a-1,b-1)*X(1,a)/A(a-1,b-2);%X(1,a-2)=(A(a-2,b)-A(a-2,b-1)*X(1,a)-A(a-2,b-2)*X(1,a-1)/A(a-2,b-3);X(1,a-m)=A(a-m,b);for n=1:mX(1,a-m)=X(1,a-m)-A(a-m,b-n)*X(1,a-n+1);endX(1,a-m)=X(1,a-m)/A(a-m,b-m-1);endAT=A;end五 例题验证选取教材上例题（P113，例5.3）进行验证，例题如下：用Gauss消去法解方程组 根据例题进行编程，程序如下：clear;clc;%题目（P113，例5.3）%给出增广矩阵A=1,2/3,1/3,2; 9/20,1,11/20,2; 2/3,1/3,1,2;AT,X=Gauss_dc(A);disp(消元后的增广矩阵为：);disp(AT);disp(求得的方程组的解为：);disp(X);仿真实验结果为： 经过与教材所给例题答案进行比对，完全正确，由此可知，所编程序完全正确。实验五 线性方程组的迭代解法之雅可比迭代法本人选择线性方程组的迭代解法为雅可比迭代法，本人选择教材上的例题为（P153，例6.5），其具体题目如下：证明用J法和GS法解下列方程组 必收敛，并比较满足的迭代次数.说明：验证时只挑选此题的一小部分，并不完全做。一 实验内容1.编写Jacobi迭代法程序，并解线性方程组的近似解。2.使用教材上例题对所编程序的正确性进行验证。二 实验步骤1.由系数矩阵A分别求出L、U、D矩阵，A、L、U、D表示形式如下： 2. 构造迭代公式如下： k=0,1,2,. (1) (2) (3)3. 验证下式是否成立， (4)4. 选择教材上例题进行验证，以判断所编程序是否正确。三 难点分析1. 如何把式(4)的判别与迭代循环结合在一起是本程序编程过程中的一个难点。本人的解决方法是设默认循环次数1000次（1000次是一个比较大的数值了，教材上的例题很少有需要循环1000次才能得出满足误差要求的近似解的），并且每循环一次对式(4)进行一次判断，只要式(4)成立，则直接跳出循环即可。四 所编算法序清单所编的算法对应的文件为Jacobi.m%线性方程组的迭代解法之Jacobi迭代法function X n=Jacobi(A,B,deta)%X为最后求得的近似解；n为迭代次数；A为系数增广矩阵；deta为允许误差a,b=size(A);L=zeros(a,a);U=zeros(a,a);D=zeros(a,a);%下面内容是分别求L、U和D矩阵for i=1:a for j=1:a if ij L(i,j)=-A(i,j); elseif i=j D(i,j)=A(i,j); end endendBJ=inv(D)*(L+U);fJ=inv(D)*B;x0=0,0,0;x1=BJ*x0+fJ;D=x0-x1;n=1;for i=1:1000 if norm(D,inf)deta t=x1; x1=BJ*x1+fJ; D=x1-t; n=n+1; end if norm(D,inf)=deta break; endend X=x1;end 四 例题验证本人选择教材上的例题为（P153，例6.5），其具体题目如下：证明用J法和GS法解下列方程组 必收敛，并比较满足的迭代次数.说明：验证时只挑选此题的一小部分，并不完全做。根据例题编程如下：(对应的文件为t5.m)clear;clc;%验证例题 （P153，例6.5）%系数矩阵A=10,-2,-2; -2,10,-1; -1,-2,3;%方程右端矩阵B=1;0.5;1;a,b=size(A);X=zeros(a,1);wx=10(-5);%允许误差X,n=Jacobi(A,B,wx);disp(通过Jacobi迭代法求得的近似解为：);disp(X);disp(使用Jacobi方法迭代次数为：);disp(n);实验结果如下：经过与教材例题所给答案比对，完全正确，由此可知，所编程序正确。实验六 非线性方程的数值解法之Newton迭代法本人所选择的非线性方程组的数值解法为Newton迭代法。本人选择课本例题(P183，例7.11)进行验证，具体题目如下;方程的根是二重根，用三种方法求解。说明：只做该题的部分题目，并且给所给答案的表格进行验证。一 实验内容1. 使用Newton迭代法求解非线性方程。2. 使用教材上的例题验证所编程序的正确性。二 实验步骤1. 先求解的一阶导数；2. 根据下式进行迭代： k=0,1,. (1)3. 与误差限定值进行比对，如果误差满足需要则终止迭代，否者则一直持续下去。4. 找教材上的例题对所编程序进行验证，验证程序的正确性。三 实验难点1. 因为需要向matlab自定义函数中代入函数解析式，所以需要用到inline()函数。2. 向微分后的函数解析式df代入自变量时需要使用subs()函数。四 所编程序详细清单所编程序对应的文件为Newton_DD.mfunction R=Newton_DD(f,x0)syms x;y=f(x);dy=diff(y,x,1);%求f(x)微分R=zeros(2,11);R(2,1)=x0;for i=2:11%设定迭代10次,第一列表征初始值 R(1,i)=i-1;%R矩阵第一行表示第几次迭代 R(2,i)=R(2,i-1)-subs(y,x,R(2,i-1)/subs(dy,x,R(2,i-1);end 五 例题验证本人选择课本例题(P183，例7.11)进行验证，具体题目如下;方程的根是二重根，用三种方法求解。说明：只做该题的部分题目，并且给所给答案的表格进行验证。根据题目编写程序如下：(对应文件为t6.m)%(P183，例7.11)clear;clc;%设f(x)=0中的f(x)函数syms x;f=inline(x4-4*x2+4,x);%取迭代初值x0x0=1.5;R=Newton_DD(f,x0);%经过Newton迭代法求得的结果矩阵disp(计算结果为：（第一行为迭代次数，第二行为对应迭代次数的计算结果）)disp(R);实验结果如下：经过与课本例题所给答案进行比对可以发现，所编程序完全正确。实验七 矩阵特征值问题的数值计算之幂法本人选择幂法进行矩阵特征值问题的数值计算所选题目（P212，例9.3）用幂法求矩阵 的主特征值和主特征向量。一 实验内容1. 使用幂法求取矩阵的主特征值和主特征向量。2. 使用课本例题验证所编程序的正确性。二 实验步骤1. 根据以下迭代公式求取主特征向量，主特征值。 (1) (2) ，k=1,2,. (3)因为 主特征向量 (4) 主特征值 (5)采用迭代算法时，可以这样理解：每一次迭代求得的和近似为本次迭代所得到的主特征向量和主特征值。2. 将所选教材例题（P212，例9.3）进行编程，用来验证所编程序的正确性。三 难点分析1. 因为每迭代一次就要有一次输出，所以输出结果显得比较凌乱。四 具体所编算法清单所编算法对应的文件为MiMethod.mfunction lanbuda uk n=MiMethod(A,u0)v0=u0;n=0;for i=1:20 v1=A*u0; a,b=size(v1); vmax=v1(1,1);%先设max(v1)初值为v1第一个元素 for j=1:a if vmaxv1(j,1) vmax=v1(j,1); else continue; end end u0=v1/vmax; n=n+1; disp(第); disp(n); disp(次运算结果输出的主特征值为：); disp(vmax); disp(第); disp(n); disp(次运算结果输出的主特征向量为：); disp(u0);endlanbuda=vmax;uk=u0;end五 例题验证所选题目（P212，例9.3）用幂法求矩阵 的主特征值和主特征向量。将所选例题进行验证，对例题所编程序代码如下：(对应文件为t7.m)clear;clc;%教材P212 例9.3 求A矩阵的住特征值和主特征限量A=1,1,0.5; 1,1,0.25; 0.5,0.25,2;u0=1,1,1;lanbuda,uk,n=MiMethod(A,u0);%lanbuda为主特征值；uk为最后计算的主特征向量；n为第几次运算%disp(主特征向量为：);%disp(uk);%disp(主特征值为：);%disp(lanbuda);输出结果为：注：本人对程序设定的输出是每迭代一次就输出一次，以上截图仅取最后一次输出结果。实验八 常微分方程的数值解法之改进的Euler法本人选择常微分方程的数值解法为改进的Euler法，所选的验证题目为教材（P238,例10.5）。具体题目如下：取h=0.1，用改进的的Euler方法解 一 实验内容1. 使用改进的Euler编写程序。2. 选择教材上例题验证所编程序的正确性。二 实验步骤1. 使用下式进行迭代求解： (1)2. 代入所选例题进行验证。三 实验难点1.编程过程中设中间变量有一定的困难。四 具体所编算法本人所编算法对应的文件为EulerX.m。%改进的Euler公式function EX=EulerX(h,dy,y0)syms x y;f=dy;disp(f);EX=zeros(3,8);for i=1:9 EX(1,i)=0.1*(i-1);end EX(2,1)=y0;for i=1:8 yba=EX(2,i)+h*(f(EX(1,i),EX(2,i);%预估 EX(2,i+1)=EX(2,i)+(h/2)*(f(EX(1,i),EX(2,i)+f(EX(1,i+1),yba);%校正end%为方便比较误差，下面求不同自变量下的精确解for i=1:9 EX(3,i)=sqrt(1+2*EX(1,i);end五 例题验证所选的验证题目为教材（P238,例10.5）。具体题目如下：取h=0.1，用改进的的Euler方法解 根据例题编写程序如下：(对应文件为t8.m)。%教材例题 (P236，例10.5)h=0.1;%设定步长dy=inline(y-2*x/y,x,y);y0=1;EX=EulerX(h,dy,y0);disp(改进的Euler方法：（第一行为x，第二行为近似解,第三行为精确解）)disp(EX); 经过与教材例题进行比对，完全正确。 1. 基于C8051F单片机直流电动机反馈控制系统的设计与研究2. 基于单片机的嵌入式Web服务器的研究 3. MOTOROLA单片机MC68HC（8）05PV8/A内嵌EEPROM的工艺和制程方法及对良率的影响研究 4. 基于模糊控制的电阻钎焊单片机温度控制系统的研制 5. 基于MCS-51系列单片机的通用控制模块的研究 6. 基于单片机实现的供暖系统最佳启停自校正（STR）调节器7. 单片机控制的二级倒立摆系统的研究8. 基于增强型51系列单片机的TCP/IP协议栈的实现 9. 基于单片机的蓄电池自动监测系统 10. 基于32位嵌入式单片机系统的图像采集与处理技术的研究11. 基于单片机的作物营养诊断专家系统的研究 12. 基于单片机的交流伺服电机运动控制系统研究与开发 13. 基于单片机的泵管内壁硬度测试仪的研制 14. 基于单片机的自动找平控制系统研究 15. 基于C8051F040单片机的嵌入式系统开发 16. 基于单片机的液压动力系统状态监测仪开发 17. 模糊Smith智能控制方法的研究及其单片机实现 18. 一种基于单片机的轴快流CO,2激光器的手持控制面板的研制 19. 基于双单片机冲床数控系统的研究 20. 基于CYGNAL单片机的在线间歇式浊度仪的研制 21. 基于单片机的喷油泵试验台控制器的研制 22. 基于单片机的软起动器的研究和设计 23. 基于单片机控制的高速快走丝电火花线切割机床短循环走丝方式研究 24. 基于单片机的机电产品控制系统开发 25. 基于PIC单片机的智能手机充电器 26. 基于单片机的实时内核设计及其应用研究 27. 基于单片机的远程抄表系统的设计与研究 28. 基于单片机的烟气二氧化硫浓度检测仪的研制 29. 基于微型光谱仪的单片机系统 30. 单片机系统软件构件开发的技术研究 31. 基于单片机的液体点滴速度自动