数值分析实验报告(包括高斯消去、二分法、牛顿迭代法)

1开课学院、实验室：实验时间：2014年1月1日课程名称数值分析基础性实验实验项目名称数值计算算法及实现实验项目类型验证演示综合设计其他指导教师成绩√实验报告一：Gauss列主元消去法求解方程组Axb一、算法（1）选主元，对1,2,,1kn确定ki，使,maxkikikkinaa，若,0kika，停机.（2）若kik转（3），否则换行,,()kkkjijkiaabbkjn（3）ikikkkaaa,(1,2,,)ijijikkjiiikkaaaabbabikkn（4）若0nna，输出无解信息，停机.（5）回代求解1,1,2,,3,2,1nnnnniijjjniiibbababbinna（6）输出解12(,,,)Tnbbb.二、Matlab源程序functionx=gausscolumn(A,b)%解Ax=b的高斯列主元消去法[m,n]=size(A);ifm~=nerror('输入错误，系数矩阵只能是方阵')endifn~=length(b)error('输入错误,常数项的个数应与方程的个数相同')end2x=zeros(1,n);fork=1:(n-1)%选主元[m,r]=max(abs(A(k:n,k)));r=k+r-1;%判断A是否为奇异阵ifm==0disp('A是奇异阵，请重新输入非奇异阵');returnend%交换增广矩阵的行ifrktemp=A(k,:);A(k,:)=A(r,:);A(r,:)=temp;c=b(k);b(k)=b(r);b(r)=c;end%进行初等行变换使主元下面的元素为0fori=(k+1):nforj=k+1:nA(i,j)=A(i,j)-A(i,k)/A(k,k)*A(k,j);endb(i)=b(i)-A(i,k)/A(k,k)*b(k);A(i,k)=0;endend%解上三角方程组ifA(n,n)==0disp('A是奇异阵，请重新输入非奇异阵');returnelsex(n)=b(n)/A(n,n);endforp=(n-1):-1:1forj=(p+1):nx(p)=(b(p)-A(p,j)*x(j))/A(p,p);endend三、结果讨论及分析用课本P13的例2.2对上述程序加以验证：求解线性方程组1231231230.501.13.16.02.04.50.360.0205.00.966.50.96xxxxxxxxx用matlab调用gausscolumn.m求解方程组，界面显示如下：A=[0.501.13.1;2.04.50.36;5.00.966.5];3b=[6.00.0200.96]';gausscolumn(A,b)ans=-2.60001.00002.0000(inv(A)*b)'ans=-2.60001.00002.0000在此例中，利用inv(A)*b可得到精确解，且通过程序得到的解和精确解一致，验证了程序的正确性。结果分析：列主元消去法相对于Gauss顺序消去法多了选主元这个方法，避免了零主元和小主元的出现，从而保证Gauss消去法能进行下去或保证方程组解的数值稳定性。从上例中对比课本给出的顺序消去法的结果即可得出此结论。实验报告二：二分法求解非线性方程()0fx的根一、算法输入函数()fx、根的存在区间[,]ab、误差限0、最大的迭代次数N（1）对1,2,kN做到第5步；（2）2abx；（3）2ebax；（4）若ex或()fx，则输出,(),xfxk，停机；否则（5）若()()0fafx，则bx，否则ax；（6）输出超过最大迭代次数N的信息，停机。二、Matlab源程序function[k,x,fx]=dichotomy(f,a,b,e,N)%输入函数f，根的存在区间[a,b]，误差限e，最大迭代次数Nf=input('pleaseenterafunction:f(x)=');a=input('pleaseentertheinitialvalue:a=');b=input('b=');e=input('pleaseentererror:e=');N=input('pleaseenterthelargestnumberofiterations:N=');fork=1:Nx=(a+b)/2;fx=feval(f,x);fa=feval(f,a);ifabs((b-a)/2)e||abs(fx)edisp('thenumberofiterationsis');kdisp('theapproximatesolutionis');x4disp('f(x)is');fxreturnendiffa*fx0b=x;elsea=x;endenddisp('Themaximumnumberofiterationsisreached,stopcalculation');end三、结果分析及讨论用课本P75的例5.4对上述程序加以验证：求解sin()0.5xx在0.7附近的根，误差限710，取根的存在区间为[0.5,1.0]。解：令()sin()0.5fxxx，即求()0fx的根。Matlab界面显示如下：dichotomypleaseenterafunction:f(x)=@(x)x*sin(x)-0.5pleaseentertheinitialvalue:a=0.5b=1pleaseentererror:e=10^(-7)pleaseenterthelargestnumberofiterations:N=10^4结果为：求得的近似解为0.7408，对应的函数值为85.280910，迭代次数17k。近似解和课本上的近似解0.7408410相近，程序可得验证。实验报告三：Newton迭代法求()0fx的根一、算法输入函数()fx，初值0x，误差限0，最大迭代次数N（1）对1,2,,kN，做到第6步；（2）计算0'()fx；（3）若0'()0fx，停止计算；否则（4）000()'()fxxxfx；5（5）若0xx或()fx，则输出,(),xfxk，停机；否则（6）0xx（7）若kN，输出超过最大迭代次数的信息，停机。二、Matlab源程序function[x,k,fx]=newton(f,x0,e,N)f=input('pleaseenterafunction:f(x)=');x0=input('pleaseentertheinitialvalue:x0=');e=input('pleaseentererror:e=');N=input('pleaseenterthelargestnumberofiterations:N=');g=diff(sym(f));fork=1:Nx=x0;f0=eval(f);g0=eval(g);ifg0==0disp('pleasechooseanotherinitialvalue');returnelsex=x0-f0/g0;fx=eval(f);ifabs(x-x0)e||abs(fx)edisp('theapproximatesolutionis');xdisp('f(x)is');fxdisp('thenumberofiterationsis');kreturnelsex0=x;endendenddisp('Themaximumnumberofiterationsisreached,stopcalculation');end三、结果讨论及分析和二分法一样，用课本P75的例5.4对上述程序加以验证：用Newton法求解sin()0.5xx在0.7附近的根，误差限710。解：令()sin()0.5fxxx，即求()0fx的根。Matlab界面显示如下：newtonpleaseenterafunction:f(x)='x*sin(x)-0.5'6pleaseentertheinitialvalue:x0=0.7pleaseentererror:e=10^(-7)pleaseenterthelargestnumberofiterations:N=10^4结果为：近似解为0.7408，对应的函数值为142.309310，迭代次数3k。结果分析：求解上述同一个方程的根，二分法迭代17次，Newton法迭代3次，但两种方法求得的近似解相同，从而得到以下结论。对于二分法，只要能够保证在给定的区间内有根，是能够收敛的，收敛速度和给定的区间有关，而且收敛速度较慢，为线性收敛。而Newton法收敛速度比二分法快，但是其最终收敛结果与初值的选取有关。学生签名年月日备注：1、此表表头必须按此格式制作。2、表头以下的栏目和内容，可根据实验内容的具体需要和要求确定，表中所列内容仅供参考。该栏可以根据需要加页。3、一门课程有多个实验项目的，应每一个实验项目一份，并将该课程所有实验项目内页与封面一起装订成册。