数值分析-实验报告

数值分析上机实验报告学院：专业：班级：学号：学生姓名：指导教师：实验五解线性方程组的直接方法实验5.1〔主元的选取与算法的稳定性〕问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：〔1〕取矩阵，那么方程有解。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？〔2〕现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。假设每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。〔3〕取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。〔4〕选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。实验过程：程序：建立M文件：functionx=gauss(n,r)n=input('请输入矩阵A的阶数:n=')A=diag(6*ones(1,n))+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)b=A*ones(n,1)p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)pause[m,n]=size(A);nb=n+1;Ab=[Ab]r=input('请输入是否为手动，手动输入1，自动输入0：r=')fori=1:n-1ifr==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);disp(Ab);pauseendendifr==1i=iip=input('输入i列所选元素所处的行数：ip=');Ab([iip],:)=Ab([ipi],:);disp(Ab);pauseendpivot=Ab(i,i);fork=i+1:nAb(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end数值实验结果及分析：⑴取矩阵A的阶数:n=10，自动选取主元：>>formatlong>>gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=1p=1pp=2.557500000000000e+003请输入是否为手动，手动输入1，自动输入0：r=0r=0⑵取矩阵A的阶数:n=10，手动选取主元：①选取绝对值最大的元素为主元：>>gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=2p=2pp=1.727556024913903e+003请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=1111111111②选取绝对值最小的元素为主元：>>gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=2p=2pp=1.727556024913903e+003请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999999991.000000000000010.999999999999981.00000000000003⑶取矩阵A的阶数:n=20，手动选取主元：选取绝对值最大的元素为主元：>>gauss请输入矩阵A的阶数:n=20条件数对应的范数是p-范数：p=1p=1ans=11111111111111111111选取绝对值最小的元素为主元：>>gauss请输入矩阵A的阶数:n=20.n=20条件数对应的范数是p-范数：p=2p=2pp=1.789670565881683e+006请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000010.999999999999971.000000000000060.999999999999891.000000000000230.999999999999551.000000000000900.999999999998211.000000000003520.999999999993181.000000000012730.999999999978171.00000000002910⑷将M文件中的第三行：A=diag(6*ones(1,n))+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)改为：A=hilb(n)①>>gauss请输入矩阵A的阶数:n=7n=7条件数对应的范数是p-范数：p=1p=1请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000510.999999999972511.000000000313540.999999998641331.000000002688050.999999997541811.00000000084337②>>gauss请输入矩阵A的阶数:n=7n=7条件数对应的范数是p-范数：p=2p=2pp=4.753673569067072e+008请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=0.999999999998691.000000000043370.999999999642991.000000001211430.999999998030381.000000001528250.99999999954491该问题在主元选取与算出结果有着很大的关系，取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比拟准确，即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大。讨论：在gauss消去法解线性方程组时，主元的选择与算法的稳定性有密切的联系，选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。条件数越大对用gauss消去法解线性方程组时，对结果产生的误差就越大。实验总结：对用gauss消去法解线性方程组时，主元的选取与算法的稳定性有密切的联系，选取适当的主元有利于得出稳定的算法，在算法的过程中，选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定，选取绝对值最大的元素作为主元时，得出的结果相对较准确较稳定。条件数越小，对用这种方法得出的结果更准确。在算除法的过程中要尽量防止使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。实验5.2〔线性代数方程组的性态与条件数的估计〕问题提出：理论上，线性代数方程组的摄动满足矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算通常要比求解方程还困难。实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动，使得充分小。实验要求：〔1〕假设方程Ax=b的解为x，求解方程，以1-范数，给出的计算结果。〔2〕选择一系列维数递增的矩阵〔可以是随机生成的〕，比拟函数“condest”所需机器时间的差异.考虑假设干逆是的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比拟。〔3〕利用“condest”给出矩阵A条件数的估计，针对〔1〕中的结果给出的理论估计，并将它与〔1〕给出的计算结果进行比拟，分析所得结果。注意，如果给出了cond(A)和的估计，马上就可以给出的估计。〔4〕估计著名的Hilbert矩阵的条件数。实验过程：程序：⑴n=input('pleaseinputn:n=')%输入矩阵的阶数a=fix(100*rand(n))+1%随机生成一个矩阵ax=ones(n,1)%假设知道方程组的解全为1b=a*x%用矩阵a和以知解得出矩阵bdata=rand(n)*0.00001%随即生成扰动矩阵datadatb=rand(n,1)*0.00001%随即生成扰动矩阵datbA=a+dataB=b+datbxx=geshow(A,B)%解扰动后的解x0=norm(xx-x,1)/norm(x,1)%得出的理论结果保存为：fanshu.mfunctionx=geshow(A,B)%用高斯消去法解方程组[m,n]=size(A);nb=n+1;AB=[AB];fori=1:n-1pivot=AB(i,i);fork=i+1:nAB(k,i:nb)=AB(k,i:nb)-(AB(k,i)/pivot)*AB(i,i:nb);endendx=zeros(n,1);x(n)=AB(n,nb)/AB(n,n);fori=n-1:-1:1x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n))/AB(i,i);end保存为：geshow.m⑵functioncond2(A)%自定义求二阶条件数B=A'*A;[V1,D1]=eig(B);[V2,D2]=eig(B^(-1));cond2A=sqrt(max(max(D1)))*sqrt(max(max(D2)))end保存为：cond2.mformatlongforn=10:10:100n=n%n为矩阵的阶A=fix(100*randn(n));%随机生成矩阵AcondestA=condest(A)%用condest求条件数cond2(A)%用自定义的求条件数condA2=cond(A,2)%用cond求条件数pause%运行一次暂停end保存为：shiyan52.m⑶n=input('pleaseinputn:n=')%输入矩阵的阶数a=fix(100*rand(n))+1;%随机生成一个矩阵ax=ones(n,1);%假设知道方程组的解全为1b=a*x;%用矩阵a和以知解得出矩阵bdata=rand(n)*0.00001;%随即生成扰动矩阵datadatb=rand(n,1)*0.00001;%随即生成扰动矩阵datbA=a+data;B=b+datb;xx=geshow(A,B);%利用第一小问的geshow.m求出解阵x0=norm(xx-x,1)/norm(x,1)%得出的理论结果x00=cond(A)/(1-norm(inv(A))*norm(xx-x))*(norm((xx-x))/(norm(A))+norm(datb)/norm(B))%得出的估计值datx=abs(x0-x00)%求两者之间的误差保存为：sy5_2.m⑷formatlongforn=4:11n=n%n为矩阵的阶数Hi=hilb(n);%生成Hilbert矩阵cond1Hi=cond(Hi,1)%求Hilbert矩阵得三种条件数cond2Hi=cond(Hi,2)condinfHi=cond(Hi,inf)pauseend数值实验结果及分析：⑴>>fanshupleaseinputn:n=6n=6a=142516881989329385489260144088501316235219292324010100737241437227701x=111111b=251410221157218187data=1.0e-005*datb=1.0e-005*0.335632943352170.275099821466210.04452752039203A=1.0e+002*0.320000064986810.930000023756510.850000054592420.480000097047240.920000035801710.60000002215793B=1.0e+002*xx=0.999998307797201.000000225695551.000000193415550.999999093880730.999999968940211.00000066032794x0=的计算结果为：〔2〕NcondestAcond2AcondA210e+00232.8905456307542132.89054563075420203.470959631940668e+002306.050503865112835e+002e+002e+002403.549487892582470e+00261.3753756968344861.3753756968336550e+002601.082004656409367e+0041.704830815154781e+0031.704830815108527e+003703.234679145192132e+0033.878481155980936e+0023.878481155978439e+00280e+002902.063634143407935e+003e+002e+0021001.536592818758897e+003e+002e+002⑶>>sy5_2pleaseinputn:n=8n=8x0=1.095033343195828e-006x00=1.705456352162135e-005datx=1.595953017842553e-005给出对的估计是：1.705456352162135e-005的理论结果是：1.095033343195828e-006结果相差：1.595953017842553e-005(4)ncond1Hicond2HicondinfHi42.837499999999738e+004e+0042.837499999999739e+00459.436559999999364e+0054.766072502414135e+0059.436559999999336e+00562.907027900294878e+007e+0072.907027900294064e+0077e+0084.753673565864586e+008e+00883.387279082022742e+0101.525757545841988e+0103.387279081949470e+01091.099650993366047e+012e+0111.099650991701052e+012103.535372424347474e+0131.602528637652488e+0133.535372455375642e+013111.230369955362001e+0155.223946340715823e+0141.230369938308720e+015讨论：线性代数方程组的性态与条件数有着很重要的关系，既矩阵的条件数是刻画矩阵性质的一个重要的依据，条件数越大，矩阵“病态”性越严重，在解线性代数方程组的过程中较容易产生比拟大的误差，那么在实际问题的操作过程中，我们必须要减少对条件数来求解，把条件数较大的矩阵化成条件数较小的矩阵来进行求解。实验总结：在本次实验中，使我们知道了矩阵条件数对线性代数方程组求解的影响，条件数越大，对最后解的影响的越大，hilbert矩阵是一个很”病态”的矩阵，他的条件数随着阶数的增加而增大，每增加一阶，条件数就增大一个数量级，在求解的过程中要尽量防止hilbert矩阵思考题一：〔Vadermonde矩阵〕设，其中，，〔1〕对n=2,5,8，计算A的条件数；随n增大，矩阵性态如何变化？〔2〕对n=5，解方程组Ax=b；设A的最后一个元素有扰动10-4，再求解Ax=b〔3〕计算〔2〕扰动相对误差与解的相对偏差，分析它们与条件数的关系。〔4〕你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？实验六解线性方程组的迭代法实验6.1〔病态的线性方程组的求解〕问题提出：理论的分析说明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵，这是一个著名的病态问题。通过首先给定解〔例如取为各个分量均为1〕再计算出右端b的方法给出确定的问题。实验要求：〔1〕选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比拟，结论如何？〔2〕逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？〔3〕讨论病态问题求解的算法第〔1〕问：选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比拟，结论如何？源程序：Gauss消去法functionx=gauss(n)formatshortdisp('请输入当前你要计算的Hilbert矩阵的阶数:')n=input('');disp('构造出的Hilbert矩阵为:')A=hilb(n)d=diag((diag(A))')*eye(n,n);b=A*ones(n,1);[m,n]=size(A);n1=n+1;disp('该线性方程组的增广矩阵为:')Ab=[Ab]fori=1:n-1[zhuyuan,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);endzhuyuan=Ab(i,i);fork=i+1:nAb(k,i:n1)=Ab(k,i:n1)-(Ab(k,i)/zhuyuan)*Ab(i,i:n1);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,n1)/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,n1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endJ迭代法formatshortdisp('请输入当前你要计算的Hilbert矩阵的阶数:')n=input('');disp('构造出的Hilbert矩阵为:')A=hilb(n)d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);x=[0;0;0;0;0;0];fori=1:1000x1=inv(d)*(l+u)*x+inv(d)*b;x=x1;enddisp('采用J迭代法结果为:')x1GS迭代法formatshortdisp('请输入当前你要计算的Hilbert矩阵的阶数:')n=input('');disp('构造出的Hilbert矩阵为:')A=hilb(n)d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);x=[0;0;0;0;0;0];fori=1:1000x1=inv(d-l)*u*x+inv(d-l)*b;x=x1;enddisp('采用GS迭代法结果为:')x1SOR迭代法formatshortdisp('请输入你当前要计算的Hilbert矩阵的阶数:')n=input('');disp('构造出的Hilbert矩阵为:')A=hilb(n)d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);disp(请输入松弛因子w:')w=input('');x=[0;0;0;0;0;0];fori=1:1000x1=inv(d-w*l)*((1-w)*d+w*u)*x+inv(d-w*l)*w*b;x=x1;enddisp('采用SOR迭代法结果为:')x1实验结果：Gauss消去法请输入当前你要计算的Hilbert矩阵的阶数:6构造出的Hilbert矩阵为:A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909该线性方程组的增广矩阵为:Ab=1.00000.50000.33330.25000.20000.16672.45000.50000.33330.25000.20000.16670.14291.59290.33330.25000.20000.16670.14290.12501.21790.25000.20000.16670.14290.12500.11110.99560.20000.16670.14290.12500.11110.10000.84560.16670.14290.12500.11110.10000.09090.73651.00000.50000.33330.25000.20000.16672.450000.08330.08330.07500.06670.05950.367900.08330.08890.08330.07620.06940.401200.07500.08330.08040.07500.06940.383100.06670.07620.07500.07110.06670.355600.05950.06940.06940.06670.06310.32821.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.401200-0.0056-0.0083-0.0095-0.0099-0.0333000.00330.00540.00640.00690.0221000.00510.00830.01020.01110.0347000.00600.00990.01220.01350.04161.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.0416000-0.0002-0.0004-0.0006-0.0013000-0.0001-0.0003-0.0004-0.00090000.00090.00190.00270.00551.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.04160000.00090.00190.00270.00550000-0.0000-0.0000-0.00010000-0.0000-0.0001-0.00011.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.04160000.00090.00190.00270.00550000-0.0000-0.0001-0.000100000-0.0000-0.0000ans=1.00001.00001.00001.00001.00001.0000J迭代法请输入当前你要计算的Hilbert矩阵的阶数:6构造出的Hilbert矩阵为:A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909采用J迭代法结果为:x1=NaNNaNNaNNaNNaNNaNGS迭代法disp('采用GS迭代法结果为:')x1请输入当前你要计算的Hilbert矩阵的阶数:6构造出的Hilbert矩阵为:A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909采用GS迭代法结果为:x1=0.99931.01310.95371.03741.02960.9662SOR迭代法请输入你当前要计算的Hilbert矩阵的阶数:6构造出的Hilbert矩阵为:A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909请输入松弛因子w:1.5采用SOR迭代法结果为:x1=0.99851.02510.90761.10480.99180.9713结果分析整理以上采用四种方法求出的希尔伯特矩阵的解，如下表所示：GUESS迭代法J迭代法GS迭代法SOR迭代法1NAN0.99930.99851NAN1.01311.02511NAN0.95370.90761NAN1.03741.10481NAN1.02960.99181NAN0.96620.9713从上表可以看到，在希尔伯特矩阵为六阶的时候。采用GUESS迭代法可以得到精确解。在J迭代法时候迭代法发散，不能得出解，改良GUESS的GS迭代法得到了近似解，在采用松弛因子为1.5时，用SOR迭代法也得到了近似解。但是我们必不能就此得出GUESS迭代法最好等结论。为了得到更加合理的结果，接下来继续对不同阶矩阵进行讨论。第〔2〕问：逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？源程序：GUESS迭代法依次改变高斯迭代法中矩阵的阶数，分别对N等于7到12运行程序，得到解。J迭代法formatshortforn=7:12A=hilb(n);d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);w=1.2;x=zeros(n,1);fori=1:1000x1=inv(d)*(l+u)*x+inv(d)*b;x=x1;endndisp('采用J迭代法结果为:')x1endGS迭代法formatshortforn=7:12A=hilb(n);d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);w=1.2;x=zeros(n,1);fori=1:1000x1=inv(d-l)*u*x+inv(d-l)*b;x=x1;endndisp('采用GS迭代法结果为:')x1endSOR迭代法formatshortforn=7:12A=hilb(n);d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);w=1.2;x=zeros(n,1);fori=1:1000x1=inv(d-w*l)*((1-w)*d+w*u)*x+inv(d-w*l)*w*b;x=x1;endndisp('采用SOR迭代法结果为:')x1end实验结果GUESS迭代法采用GUESS迭代法，分别计算阶数为7到12阶时，对应的解如下：7891011121.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00000.99991.00001.00001.00001.00020.99981.00051.00001.00001.00000.99961.00070.99811.00001.00001.00001.00070.99841.00431.00001.00000.99931.00220.99401.00001.00040.99821.00520.99991.00080.99740.99981.00061.0000J迭代法采用J迭代法，分别计算阶数为7到12阶时，对应的解如下：789101112NANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANGS迭代法采用GS迭代法，分别计算阶数为7到12阶时，对应的解如下：7891011120.99810.99720.99670.99680.99730.99821.02521.03291.03441.02971.02021.00780.93210.92720.93820.96120.99031.02051.02381.00240.98130.96530.95620.95361.04851.04581.02991.00850.98710.96881.01651.04111.04581.03741.02211.00450.95451.00351.03081.04091.03931.03080.94790.99451.02261.03621.03970.94580.9891.01661.03230.94590.98531.01190.94650.98230.9467SOR迭代法采用SOR迭代法，分别计算阶数为7到12阶时，对应的解如下：7891011120.99770.99680.99650.99680.99760.99881.03171.03871.03781.02951.01621.00060.9110.90920.92860.96251.00281.04281.04461.01830.98860.96270.94380.93251.04761.04631.03131.01030.98910.97061.01441.04011.04541.03741.02271.00560.95161.00161.02951.04041.03981.03240.94690.99381.0221.03621.04050.9460.98891.01651.03250.94660.98551.01190.94710.98230.9467结果分析在以上四种运算中，我们可以看到最不稳定的是采用J迭代法，在第一问中我们看到，当希尔伯特矩阵为六阶的时候没有得到近似解，随着矩阵阶数的继续增加，依然没有得到解。故此在解高度病态矩阵的时候我们应该尽量防止使用J迭代法。在采用GUESS迭代法的时候我们看到，针对特定的希尔伯特矩阵求解，一直到九阶以下都能够得到精确解，但是随着矩阵阶数的提高，可以看到解出现了不稳定的情况，在1的附近出现左右震荡的情况。且从上表可以看出，随着矩阵阶数的增加，有震荡越来越剧烈的倾向。故此在求解阶数较小的矩阵的时候，可以优先选择高斯迭代法，但是在阶数较大的矩阵的求解过程，应该尽量选择更有的数值解法。GS迭代法和SOR迭代法在阶数逐渐增加的时候，解依然在一附近波动，并且没有出现误差增大的情况。我们可以认为这是两种比拟优秀的解高阶病态矩阵的方法。特别值得提出的就是SOR迭代法的收敛性还和迭代因子的选择有密切的关系。在实际计算过程中，需要通过简要计算选择较好的迭代因子，或者是在MATLAB中编程，在1~2之间选择微小的步长，通过编程比拟求得其最正确的步长因子，使得其迭代收敛并且收敛速度较快。总之，在求解不同的高解线性方程的时候，我们应该根据具体情况选择适宜的迭代法求解矩阵。而迭代矩阵并不仅限于上面提到的四种，我们可以自己构造更加恰当的矩阵，使到达事半功倍效果。第〔3〕问：讨论病态问题求解的算法关于病态问题的求解，主要的方法是对原方程作写预处理，以降低系数矩阵的条件数。例如选择非奇异矩阵，使得方程化为等价方程组的，原方程的解。原那么上应使PAQ的条件数比A有所改善。一般P和Q可选择为三角形矩阵或者对角矩阵。理论上最好选择对角阵P和Q，并且应当满足：源代码：cleark=500；i=i:k；forn=1:k；H=hilb(n)；K(n)=cond(H)；D1=diag(sqrt(diag(H)))；D=inv(D1)；Hy=D*H*D；Ky1(n)=cond(Hy)；Ky(n)=log10(Ky1(n)/K(n))；endplot(i,Ky)xlabel(‘n’,’fontsize’,20)ylabel(‘lg(cond(hn))’,’fontsize’,20)实验结果：实验分析:上图中的〔Dn为Hn的对角线元素开放构成的对角矩阵〕，N为希尔伯特矩阵的阶数。可以看出，随着矩阵阶数的增大，函数在之间波动。还可以看出，对于上图，大多数的函数值都是小于或者等于零的，这说明经过处理以后的矩阵的条件数较小，在一定程度上改变了矩阵的性质。实验总结：对于希尔伯特矩阵病态的性质，对于原来希尔伯特矩阵进行据处理后的新的希尔伯特矩阵在一定区间里面呈现波动的变化规律。从整体上观察，对于大多说的N，进行预处理后的希尔伯特矩阵的条件数有了明显的降低，这就说明预处理改善了大多数阶数的希尔伯特矩阵的性质。实验七非线性方程求根实验7.1〔迭代法、初始值与收敛性〕实验目的：初步认识非线性问题的迭代法与线性问题迭代法的差异，探讨迭代法及初始值与迭代收敛性的关系。问题提出：迭代法是求解非线性方程的根本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。实验内容：考虑一个简单的代数方程针对上述方程，可以构造多种迭代法，如在实轴上取初始值x0，请分别用迭代〔7.1〕-〔7.3〕作实验，记录各算法的迭代过程。实验要求：〔1〕取定某个初始值，分别计算〔7.1〕-〔7.3〕迭代结果，它们的收敛性如何？重复选取不同的初始值，反复实验。请自选设计一种比拟形象的记录方式〔如利用MATLAB的图形功能〕，分析三种迭代法的收敛性与初值选取的关系。〔2〕对三个迭代法中的某个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？〔3〕线性方程组迭代法的收敛性是不依赖初始值选取的。比拟线性与非线性问题迭代的差异，有何结论和问题。实验过程：程序：clearclcs=input('请输入要运行的方程，运行第几个输入几s=');clfifs==1%决定坐标轴的范围和初始值a=-1.5;b=2.5;y00