数值分析实验报告

数值分析实验报告数值分析实验报告学院：学生姓名：学号：2016年11月页共25页实验一插值法本人选择Newton插值法进行编程。选择的例题为：教材（P25，例2.7）进行验证：题目：设，并已知的数据如表1-1试用2次Newton插值多项式计算计算的近似值，并讨论其误差。表1-1x2.02.12.21.4142141.4491381.483240一实验内容根据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中CommandWindow输出该插值多项式。(3)在教材中随机选取某一例题来验证所编程序的正确性。三所编程序清单所编程序对应的文件为Newton.m%Newton插值法function[y]=Newton(Xi,Yi,x,n)%%%计算均差表%根据教材，均差表一般列举到三阶均差即可%一阶均差f1=zeros(1,n);%一阶均差数值保存在f1中fori=1:nf1(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中fori=1:(n-1)f2(1,i)=(f1(1,i+1)-f1(1,i))/(Xi(1,i+2)-Xi(1,i));enddisp('二阶均差数值为：');disp(f2);%三阶均差ifn>2f3=zeros(1,n-2);%三阶均差数值保存在f3中fori=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次插值多项式的值ifn==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));elseifn==2y=Yi(1,1)+f1(1,1)*(x-Xi(1,1))+f2(1,1)*(x-Xi(1,1))*(x-Xi(1,2));elseifn==1y=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.02.12.2];Yi=[1.4142141.4491381.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）。计算克莱姆矩阵(2)根据下式计算系数j=0,1,...,n(3)将代入(1)即可求得待求拟合模型多项式。4.根据3所求的多项式画出函数曲线三难点分析克莱姆矩阵的求解，需要使用嵌套一次的循环，容易出错。四所编算法清单所编算法对应的文件为duoxiangshinihe.mfunction[ai]=duoxiangshinihe(i,Xi,Yi,N)wi=ones(1,i);%定义权值矩阵，并不防设权值均为1ai=zeros(N+1,1);Cramer=zeros(N+1,N+1);forii=1:(N+1)%列foriii=1:(N+1)%行Cramer(ii,iii)=0;%克莱姆矩阵元素初始化foriiii=1:iCramer(ii,iii)=Cramer(ii,iii)+Xi(1,iiii)^(ii+iii-2);endendendYFAI=zeros(N+1,1);%求(y,φ)矩阵foriiiii=1:(N+1)foriiiiii=1:iYFAI(iiiii,1)=YFAI(iiiii,1)+Yi(1,iiiiii)*Xi(1,iiiiii)^(iiiii-1);endendai=inv(Cramer)*YFAI;%由(a)=inv(φ)*(y,φ)求系数矩阵Ai=ai';Aii=rot90(Ai);%以上两行的作用是将（a0,a1,a2,a3)改成(a3,a2,a1,a0)x1=0:1;axis([0102]);plot(x1,polyval(Aii,x1),'r');holdonplot(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.000.250.500.751.00];%离散点横坐标yi=[0.100.350.811.091.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)求解积分。通过梯形公式的误差公式(2)求解求积的误差并判断，当误差不超过允许值时即可认为求积结束。找教材上的例题进行验证，验证所编程序的正确性。二实验步骤首先要对区间[a,b]进行n等分(n>1),并且对都要对(2)进行判断误差是否已经不超过允许值，如果超过了，则对区间[a,b]进行(n+1)等分(记作n=n+1)，然后再根据(2)判断误差是否不超过允许值，依次类推，直至找到能使误差不超过允许值时的n.在(1)的基础上，根据公式(1)求积，得出计算结果并输出。找教材上例题进行验证。三难点分析因为误差需要和给定的误差要求比较大小，所以的取值范围是算法编程过程的一个难点。本人的思路是先求的极值点（假设在给定区间上只有一个极值点），然后求边界点的函数值以及极值点处的函数值分别为，通过比较这三个值的大小就可以找出函数的取值范围。在确定步长的基础上，通过公式(*)求积。通过分析可知，公式(*)可分为常量和累加量两部分，我的编程思路是先给赋初值，然后再采用循环加上累加量。四所编算法清单所编算法对应的文件为Complex_Quadrature.mfunction[Tnt]=Complex_Quadrature(f,a,b,deta)%T复化梯形求积结果，nt为运算次数%f待积分函数表达式，[a,b]待积分区间，nt表示复化梯形求积公式满足误差要求时需要将[a,b]区间nt等分，ns表示复化Simpson求积公式满足%误差要求时需要将[a,b]区间ns等分,deta表示积分时误差要求%首先复化梯形求积公式%symsx;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;ifdeta1>=0&&deta2>=0tbiggest=ta;tsmallest=tx0;elseifdeta1>=0&&deta2<=0&&deta3>=0tbiggest=ta;tsmallest=tb;elseifdeta1>=0&&deta2<=0&&deta3<=0tbiggest=tx0;tsmallest=tb;endifdeta1<=0&&deta2<=0tbiggest=tx0;tsmallest=ta;elseifdeta1<=0&&deta2>=0&&deta3>=0tbiggest=tb;tsmallest=tx0;elseifdeta1<=0&&deta2>=0&&deta3<=0tbiggest=tb;tsmallest=ta;endfori=2:1000%从2等分、3等分、4等分……n等分依次逐个验证误差精度h=(b-a)/i;%计算步长RTmax=((-1)*(b-a)/12)*h^2*tsmallest;ifRTmax<=detant=i;T=(h/2)*(f(a)+f(b));forii=1:iT=T+h*f(a+h*(ii-1));endbreakelsecontinueendendns=0;S=0;end五例题验证选教材题目（P79，例4.6）作为验证例题，题目如下：题目（P79，例4.6）分别用复化梯形公式和复化Simpson公式计算时，要使误差不超过，问各取多少个点。根据所选例题进行编程，如下：(对应文件为t3.m)symsx;f=inline('sin(x)','x');a=0;b=pi;deta=2*10^(-5);[TntnsS]=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.找教材上习题进行验证所编程序的正确性。二实验步骤通过公式（1）和公式（2）将系数矩阵化为上三角阵，增广矩阵的格式如式（3）所示

(3)通过反向逆推法，依次求得，其求解公式如下：(4)i=n-1,n-2,...1.(5)用教材例题（P113，例5.3）对程序进行验证，以验证程序的正确性。三难点分析1.式（5）里面含有累加过程，在设置循环时容易出错。四所编算法清单所编算法对应的文件为Gauss_dc.m%线性方程直接接法之高斯消去法function[ATX]=Gauss_dc(A)[a,b]=size(A);%首先找出增广矩阵的行数和列数fork=1:(a-1)%总共需要（a-1）次消元fori=(k+1):aL=A(i,k)/A(k,k);forj=1:bA(i,j)=A(i,j)-L*A(k,j);endendendX=zeros(1,a);X(1,a)=A(a,b)/A(a,b-1);form=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);forn=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表示形式如下：构造迭代公式如下：k=0,1,2,...(1)(2)(3)验证下式是否成立，(4)选择教材上例题进行验证，以判断所编程序是否正确。三难点分析如何把式(4)的判别与迭代循环结合在一起是本程序编程过程中的一个难点。本人的解决方法是设默认循环次数1000次（1000次是一个比较大的数值了，教材上的例题很少有需要循环1000次才能得出满足误差要求的近似解的），并且每循环一次对式(4)进行一次判断，只要式(4)成立，则直接跳出循环即可。四所编算法序清单所编的算法对应的文件为Jacobi.m%线性方程组的迭代解法之Jacobi迭代法function[Xn]=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矩阵fori=1:aforj=1:aifi<jU(i,j)=-A(i,j);elseifi>jL(i,j)=-A(i,j);elseifi==jD(i,j)=A(i,j);endendendBJ=inv(D)*(L+U);fJ=inv(D)*B;x0=[0,0,0]';x1=BJ*x0+fJ;D=x0-x1;n=1;fori=1:1000ifnorm(D,inf)>detat=x1;x1=BJ*x1+fJ;D=x1-t;n=n+1;endifnorm(D,inf)<=detabreak;endendX=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)进行验证，具体题目如下;方程的根是二重根，用三种方法求解。说明：只做该题的部分题目，并且给所给答案的表格进行验证。一实验内容使用Newton迭代法求解非线性方程。使用教材上的例题验证所编程序的正确性。二实验步骤先求解的一阶导数；根据下式进行迭代：k=0,1,...(1)与误差限定值进行比对，如果误差满足需要则终止迭代，否者则一直持续下去。找教材上的例题对所编程序进行验证，验证程序的正确性。三实验难点因为需要向matlab自定义函数中代入函数解析式，所以需要用到inline()函数。向微分后的函数解析式df代入自变量时需要使用subs()函数。四所编程序详细清单所编程序对应的文件为Newton_DD.mfunction[R]=Newton_DD(f,x0)symsx;y=f(x);dy=diff(y,x,1);%求f(x)微分R=zeros(2,11);R(2,1)=x0;fori=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)函数symsx;f=inline('x^4-4*x^2+4','x');%取迭代初值x0x0=1.5;R=Newton_DD(f,x0);%经过Newton迭代法求得的结果矩阵disp('计算结果为：（第一行为迭代次数，第二行为对应迭代次数的计算结果）')disp(R);实验结果如下：经过与课本例题所给答案进行比对可以发现，所编程序完全正确。实验七矩阵特征值问题的数值计算之幂法本人选择幂法进行矩阵特征值问题的数值计算所选题目（P212，例9.3）用幂法求矩阵的主特征值和主特征向量。一实验内容使用幂法求取矩阵的主特征值和主特征向量。使用课本例题验证所编程序的正确性。二实验步骤根据以下迭代公式求取主特征向量，主特征值。(1)(2)，k=1,2,...(3)因为主特征向量(4)主特征值(5)采用迭代算法时，可以这样理解：每一次迭代求得的和近似为本次迭代所得到的主特征向量和主特征值。将所选教材例题（P212，例9.3）进行编程，用来验证所编程序的正确性。三难点分析因为每迭代一次就要有一次输出，所以输出结果显得比较凌乱。四具体所编算法清单所编算法对应的文件为MiMethod.mfunction[lanbudaukn]=MiMethod(A,u0)v0=u0;n=0;fori=1:20v1=A*u0;[a,b]=size(v1);vmax=v1(1,1);%先设max(v1)初值为v1第一个元素forj=1:aifvmax<v1(j,1)vmax=v1(j,1);elsecontinue;endendu0=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方法解一实验内容使用改进的Euler编写程序。选择教材上例题验证所编程序的正确性。二实验步骤使用下式进行迭代求解：(1)代入所选例题进行验证。三实验难点1.编程过程中设中间变量有一定的困难。四具体所编算法本人所编算法对应的文件为EulerX.m。%改进的Euler公式function[EX]=EulerX(h,dy,y0)symsxy;f=dy;disp(f);EX=zeros(3,8);fori=1:9EX(1,i)=0.1*(i-1);endEX(2,1)=y0;fori=1:8yba=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%为方便比较误差，下面求不同自变量下的精确解fori=1:9EX(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);经过与教材例题进行比对，完全正确。目录第一章总论 11.1项目概要 11.2可行性研究报告编制依据 11.3项目区简介及建设单位概况 21.4可研报告研究内容 61.5可研报告研究结论、问题及建议 7第二章项目背景及建设的必要性和可行性 92.1建设背景 92.2项目建设的必要性 132.3项目建设的可行性 15第三章项目建设内容及规模 173.1项目建设内容 173.2建设规模 17第四章场址选择及建设条件 194.1场址现状 194.2建设条件 19第五章工程方案 245.1方案设计指导思想及原则 245.2建筑设计 245.3结构设计 295.4给排水设计 335.5暖通设计 365.6电气及弱电设计 395.7消防设计 42第六章能源和资源节约措施 456.1能源节约措施 456.2建筑节能具体措施 466.3给排水资源节约措施 466.4电气节能具体措施 466.5供热系统节能技术措施 47第七章环境影响评价 487.1环境影响 487.2保护措施 507.3安全保护措施 517.4环境影响评价结论 52第八章组织机构 538.1管理机构 538.2项目组织管理 53第九章工程管理及实施计划 549.1项目建设管理原则 549.2工程管理 549.3项目实施步骤 55HYPERLINK\l"