计算方法矩阵直接法与迭代法

1、需要两份东西1、实验程序：输入、输出、注释2、实验报告：问题描述、方法描述、方案设计、结 果分析、结论谢谢，麻烦写的详细些实验五解线性方程组的直接方法实验5.1 (主元的选取与算法的稳定性)问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为 数值算法的稳定性呢？ Gauss消去法从理论算法到数值算法，其关键是主元的选 择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组Ax =b, A Rn n,b Rn编制一个能自动选取主元，又能手动选取主元的求解线性方程组

2、的Gauss消去过程。实验要求：一715,b =:,则方程有解11561.14TX =（1,11，1）T。一6 18 61(1)取矩阵A二|三三|86I8取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？(2)现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或 按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按 模最大的元素作为主元，结果又如何？分析实验的结果。(3)取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（ 4）选取其他你感兴趣的问题或者随机生成

3、矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。1）首先编写gauss消元法代码；输入为矩阵A,向量b,精度ptol,输出为方程 组的解x。程序如下function out = gausle（ A,b,flag,ptol ）%UNTITLED 高斯消元法可选择的主元消去法% A n x n 矩阵% b n x 1% flag 标志主元为自动选取还是手动选取，%0，自动（对角为主元）1，选取最小模为主元，%2，选取模最大，（误差最小）3 ，次最小的模为主元% ptol 精度% out 输出值，为nx1 的解 Ax=b 的解%if nargin<4,ptol=50*eps;endm

4、,n=size（A）;if m=n,error（'A 不是方阵'）;endout=zeros（n,1）;% 预先设定解的维数nb=n+1;Ab=A,b;%扩维矩阵% disp('开始用扩维阵计算)；% disp(Ab);RA=rank(A);RB=rank(Ab);zhica=RB-RA;if zhica>0,disp('请注意：因为RA=RB,所以此方程组无解.') returnend% 消元过程for i=1:n-1 ;i;if flag=0;pivot=Ab(i,i);ri=i;elseif flag=1;% 按照每列模最小的选取pivot,

5、min_index=min(abs(A(i:n,i);ri=i+min_index-1;elseif flag=2;% 按照模最大的选取pivot,max_index=max(abs(A(i:n,i);ri=i+max_index-1;elseif flag=3%方程最小非0 数tA=A;pivot,min_index=min(abs(tA(i:n,i);while pivot=0;tA(min_index+i-1,i)=inf;pivot,min_index=min(abs(tA(i:n,i);endri=i+min_index-1;endif (pivot=0)|(pivot<pto

6、l) ;warning('系数矩阵奇异！ ！)；return;endif ri=i; % 交换行tmp=A(i,:);A(i,:)=A(ri,:);A(ri,:)=tmp;t=b(i);b(i)=b(ri);b(ri)=t;endfor kk=i+1:nL(kk,i)=A(kk,i)/A(i,i);A(kk,i+1:n)=A(kk,i+1:n)-L(kk,i)*A(i,i+1:n);b(kk)=b(kk)-L(kk,i)*b(i);endendif A(n,n)=0warning('系数矩阵奇异！ ！)；return;end% % 回代求解x=zeros(n,1)'fo

7、r k=n:-1:1%k%A(k,k+1:n)%x(k+1:n)if k=nx(n)=b(n)/A(n,n);elsex(k)=(b(k)-sum( A(k,k+1:n).*x(k+1:n) ) )/(A(k,k); endendout=x'end主程序为N=10时，方程解为1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000A的条件数如下：1-条件数 2.5575e+0032-条件数 1.7276e+003无穷条件数2.5575e+003代码在 里，运行计算即可。Flag=0；,时为自动选取2) Flag=1,为选取

8、最小模的值为主元，flag=2,为选取最大模为主元，代码是一样的，修改flag值和维数n即可。3) N=20，只需要更改n,重复上述步骤即可。4)-nIiXon 1i=0XonnXizx1i zz0n,X2，b =n3zX2i=0anxn 1n1zx：_i=0思考题一 ：(Vadermond即阵)设1XoX2d21X1X1A = 1x2x2a鼻鼻JXnX2其中，xk =1 +0.1k,k =0,1,，n ,(1)对n=2,5,8,计算A的条件数；随n增大，矩阵性态如何变化？(2)对n=20,解方程组Ax=b;设A的最后一个元素有扰动10-4,再求解Ax=b(3)计算(2)扰动相对误差与解的相对

9、偏差，分析它们与条件数的关系。(4)你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？(5)尝试做出范德蒙矩阵的预优。答1) n=2, 5, 8 ; A的条件数为如下所示，随n增大，条件数越大，矩阵的 病态性越严重。n =2,cond1(A) =16060;cond2(A)=11044; cond 二(A) =16016;n =5,cond1(A) =1.0522 107; cond2(A)=0.5594 107;cond:;( A) =1.0668 107;n =8;cond1(A) =5.0565 1010; cond2(A)=2.2739 1010;

10、cond 二(A) =4.45788 1010;代码为 sikaoti_111.m2)无扰动时x =1，,1T,x为20父1的向量；有扰动时，x = 1.0e+006 *0.0027-0.02700.1250-0.35610.6956-0.98381.0344-0.81320.4692-0.18470.03570.0093-0.01040.0043-0.00100.00010.0000-0.00000.00000.00000.0000,代码文件为sikaoti _112.3)由此看来，当系数矩阵的条件数越大，则病态性越严重，因而系数矩阵很小 相对误差也会让解产生很大的相对误差。4)插值函数存在

11、定理直接求插值函数而要用拉格朗日或牛顿插值法的原因是为 了补偿计算中所产生偏差，而且迭代次数越多，偏差就越大。5)预优，相关MATLAB函数提示：zeros(m,n) 生成m行，n列的零矩阵ones(m,n) 生成m行，n列的元素全为1的矩阵eye(n)生成n阶单位矩阵rand(m,n) 生成m行,n列(0,1)上均匀分布的随机矩阵diag(x)返回由向量x的元素构成的对角矩阵tril(A)提取矩阵A的下三角部分生成下三角矩阵triu(A)提取矩阵A的上三角部分生成上三角矩阵rank(A)返回货!阵A的秩det(A)返回方阵A的行列式inv(A)返回可逆方阵A的逆矩阵V , D=eig(A)

12、返回方阵A的特征值和特征向量norm(A,p)矩阵或向量的p范数cond(A,p)矩阵的条件数L , U ,巳=lu(A)选列主元LU分解R=chol(X)平方根分解Hi=hilb(n) 生成 n 阶 Hilbert 矩阵实验六 解线性方程组的迭代法实验6.1 (病态的线性方程组的求解)问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况 是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵，.1. . _H =(hi,j)nn, hi,j, i, j =1,2, ,ni j -1这是一个著名的病态问题。通过首先给定解(例如取为各个分

13、量均为1)再 计算出右端b的办法给出确定的问题。实验要求：(1)选择问题的维数为6,分别用Gauss消去法、列主元Gauss消去法、J 迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果 与问题的解比较，结论如何？(2)逐步增大问题的维数(至少到100),仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（ 3）讨论病态问题求解的算法答：1) Guass消元法的程序在gausle.m程序里，如下所示。Flag=0,是自 动消元，Guass消去法的结果为：x=1 ，1，1，1，1，1flag=2为列主元Gauss消去法法，Gauss列主元消去法的结果为：x=1

14、 ，1，1，1，1，1Jacobi迭代法程序在jacobi.m,结果为：Inf ,Inf,NaN,NaN,NaN,NaN 。 也就是迭代发散，因为迭代阵B 的谱半径大于1,为4.3085所以此方程不适合用Jacobi方法迭GS迭代法在gauseidel.m,其结果为：0.99921.01250.95901.03061.02670.9715迭代次数为2016SOR 超级松弛法程序在SOR.m 中， w 选择为 1.89其结果为：0.99911.01800.91361.15860.88001.0304迭代次数为517 次。主程序如下。% 实验 6.1% 实验 6.1(病态的线性方程组的求解)clc

15、;clear;close all;n=6; % 矩阵维数A=hilb(n);eps=1e-5;R=max(abs(eig(A) ;%disp('A 的谱半径为',num2str(R);b=sum(A,2);x_gauss=gausle(A,b,0)x_gauss_2=gausle(A,b,2)% 编写 jacobi.m 和 gauseidel.mx_jaco_3,n3=jacobi(A,b,zeros(size(b),eps ) %( A,b,omega,x0,eps,N)x_gs_4,n4=gauseidel(A,b,zeros(size(b),eps )% SOR 法ome

16、ga=1.89;x_SOR_5,n5=SOR(A,b,omega,zeros(size(b),eps)2) 将维数增加为100，高斯消元法，结果为1.00001.00000.99931.02280.67543.5051-10.564434.3970-60.724877.6571-77.081786.1987-72.494316.549923.9588-7.319232.2195-71.5173-9.4683102.1479-4.3990-49.2259-65.428477.29858.937810.7393-2.8585-94.704818.5097112.925253.2932-77.717

17、8-45.60030.8528-20.4096-56.1519106.6452-3.433045.7688 8.1846-62.144310.476826.8102-75.727076.4337-55.6616-22.277789.8267-63.1822 111.9219 -73.864928.5538-22.5285-44.4054-20.734745.2838-39.7352-57.6491 101.373957.6956 47.3377105.3975-47.4092-159.0873-38.1213-73.969618.9134-69.1624 118.195726.6726141.

18、973669.7617-144.1574-13.454667.8139-44.4608-126.1802-39.0300 192.4826 -67.486416.6392 82.7630-146.8443-90.7014164.5643126.838711.1138-223.7063-6.505967.2569-139.3375216.16062.0236-98.0808114.9628-3.8427-55.4881-100.470398.0132-7.2417列主元高斯消元法，结果为0,即出现计算错误，弹出，系数矩阵奇异，jacobi B的谱半径为86.3375, Jacobi迭代法发散，结

19、果为NaN ,，NaNGauss-S 方法 0.99841.02480.92471.02731.05281.03211.00240.97980.96830.96620.97060.97870.98820.99761.00591.01271.01771.02111.02281.02321.02241.02071.01831.01541.01221.00891.00551.00220.99900.99610.99340.99100.98890.98720.98570.98460.98380.98330.98310.98320.98340.98390.98460.98550.98650.98760.

20、98890.99020.99160.99310.99460.99610.99760.99911.00061.00201.00341.00481.00611.00731.00841.00951.01041.01131.01211.01271.01331.01371.01411.01431.01451.01451.01441.01421.01391.01341.01291.01221.01151.01061.00961.00851.00731.00611.00471.00321.00160.99990.99820.99630.99440.99240.99030.98810.98580.98350.

21、98110.97860.97600.9734SOR 迭代法结果为0.99781.03910.84921.18200.90151.15990.85851.10600.84561.08420.8658 1.0830 0.89791.08670.92811.0874 0.9510 1.08310.96601.0743 0.97441.06270.97811.0499 0.97881.0372 0.9778 1.02540.97631.0152 0.97481.00670.97381.00010.9735 0.9953 0.9740 0.9922 0.97530.99040.9772 0.9900 0

22、.97960.99050.98260.99180.98580.9936 0.9892 0.99590.9927 0.9984 0.99611.00100.99941.00361.00251.00611.00531.00831.00781.01031.01001.01201.01171.01331.01301.01421.01381.01471.01421.01471.01421.01431.01371.01351.01271.01231.01131.01061.00941.00851.00721.00601.00451.00311.00140.99990.99800.99620.99420.9

23、9230.99010.98800.98570.98340.98100.97850.97590.9733条件数来衡量问题的病态程度。条件数越大，病态可能越严重实验6.2书上P211计算实习题2,其中N=100,第二小问改为用Jacobi迭 代、G-S迭代、红黑排序的G-S迭代求解，并比较他们之间的收敛速度；进一步， 用BSOR迭代求解，试找出最优松弛因子。习题见图片- x=Xi=ih,y=yj = jh,解：j ,由此，可知，x万向上0,1上分为了 11个小格，y万向也i, j =0,1，,N,N 1一样。五点差分方程为11rrui.j 1 " "TTui j,j hh114

24、 一、了ui1,j 一声ui,j 1 产ui,j = f(Xi,yj),即-121ui,j 1 121ui121ui,j 1 - 121ui 1,j484ui,j = f (xi, yj)Jacobi结果为jacobi B的谱半径为 0.95949 迭代次数=221y_jacobi =1.01391.02781.04151.05501.06811.08041.09141.10011.10531.10491.02781.05531.08251.10931.13521.15961.18131.19891.21011.21171.04151.08251.12311.16301.20171.23841.27141.29881.31761.32411.05501.10931.16301.21601.26771.31721.36251.40131.42991.44411.06811.13521.20171.26771.33281.39581.45471.50671.54771.57271.08041.15961.23841.31721.39581.47311.54711.61471.67121.71061.09141.18131.27141.36251.45471.54711.63791.72361.79931.85771.10011.1