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

需要两份东西
1﹑实验程序：输入﹑输出﹑注释
2﹑实验报告：问题描述﹑方法描述﹑方案设计﹑结果分析﹑结论

谢谢，麻烦写的详细些

实验五 解线性方程组的直接方法

实验5.1 （主元的选取与算法的稳定性）
问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
实验内容：考虑线性方程组

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

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;end
[m,n]=size(A);
if m~=n,error('A不是方阵');end
out=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，所以此方程组无解.')
return
end
% 消元过程
for i=1:n-1 ;
i;
if flag==0;
pivot=Ab(i,i);
ri=i;
elseif flag==1;% 按照每列模最小的选取
[pivot,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)));
end
ri=i+min_index-1;
end
if (pivot==0)||(pivot<ptol) ;
warning('系数矩阵奇异！！');
return;
end
if ri~=i; % 交换行
tmp=A(i,:);A(i,:)=A(ri,:);A(ri,:)=tmp;
t=b(i);b(i)=b(ri);b(ri)=t;
end
for kk=i+1:n
L(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);
end
end
if A(n,n)==0
warning('系数矩阵奇异！！');
return;
end
% % 回代求解
x=zeros(n,1)';
for k=n:-1:1
% k
% A(k,k+1:n)
% x(k+1:n)
if k==n
x(n)=b(n)/A(n,n);
else
x(k)=(b(k)-sum( A(k,k+1:n).*x(k+1:n) ) )/(A(k,k));
end
end

out=x';
end

主程序为

N=10时，方程解为
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000

A的条件数如下：
1-条件数 2.5575e+003

2-条件数 1.7276e+003

无穷 条件数2.5575e+003
代码在 里，运行计算即可。Flag=0；，时为自动选取
2） Flag=1，为选取最小模的值为主元，flag=2，为选取最大模为主元，代码是一样的，修改flag值和维数n即可。
3） N=20 ,只需要更改n，重复上述步骤即可。
4）

思考题一：（Vadermonde矩阵）设
，
其中，，
（1）对n=2,5,8，计算A的条件数；随n增大，矩阵性态如何变化？
（2）对n=20，解方程组Ax=b；设A的最后一个元素有扰动10-4，再求解Ax=b
（3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。
（4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？
（5）尝试做出范德蒙矩阵的预优。

‘
答 1） n=2，5，8 ；A的条件数为如下所示，随n增大，条件数越大，矩阵的病态性越严重。

代码为sikaoti_111.m

3) 由此看来，当系数矩阵的条件数越大，则病态性越严重，因而系数矩阵很小相对误差也会让解产生很大的相对误差。
4） 插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因是为了补偿计算中所产生偏差，而且迭代次数越多，偏差就越大。
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) 返回方阵A的特征值和特征向量
norm(A,p) 矩阵或向量的p范数
cond(A,p) 矩阵的条件数
[L，U，P]=lu(A) 选列主元LU分解
R=chol(X) 平方根分解
Hi=hilb(n) 生成n阶Hilbert矩阵

实验六 解线性方程组的迭代法

实验6.1（病态的线性方程组的求解）
问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？
实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵，

这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为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，1]
flag=2为列主元Gauss消去法法，Gauss列主元消去法的结果为：
x=[1，1，1，1，1，1]

Jacobi迭代法程序在jacobi.m，结果为：
[Inf ,Inf,NaN,NaN,NaN,NaN] 。也就是迭代发散，因为迭代阵B的谱半径大于1，为4.3085,所以此方程不适合用Jacobi方法迭
GS迭代法在gauseidel.m，其结果为：
0.9992
1.0125
0.9590
1.0306
1.0267
0.9715
迭代次数为2016
SOR 超级松弛法程序在SOR.m 中，w选择为1.89其结果为：
0.9991
1.0180
0.9136
1.1586
0.8800
1.0304
迭代次数为517次。
主程序如下。
%% 实验6.1
% 实验6.1（病态的线性方程组的求解）
clc;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.m
[x_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 法
omega=1.89;
[x_SOR_5,n5]=SOR(A,b,omega,zeros(size(b)),eps)

2) 将维数增加为100，高斯消元法，结果为
1.0000
1.0000
0.9993
1.0228
0.6754
3.5051
-10.5644
34.3970
-60.7248
77.6571
-77.0817
86.1987
-72.4943
16.5499
23.9588
-7.3192
32.2195
-71.5173
-9.4683
102.1479
-4.3990
-49.2259
-65.4284
77.2985
8.9378
10.7393
-2.8585
-94.7048
18.5097
112.9252
53.2932
-77.7178
-45.6003
0.8528
-20.4096
-56.1519
106.6452
-3.4330
45.7688
8.1846
-62.1443
10.4768
26.8102
-75.7270
76.4337
-55.6616
-22.2777
89.8267
-63.1822
111.9219
-73.8649
28.5538
-22.5285
-44.4054
-20.7347
45.2838
-39.7352
-57.6491
101.3739
57.6956
47.3377
105.3975
-47.4092
-159.0873
-38.1213
-73.9696
18.9134
-69.1624
118.1957
26.6726
141.9736
69.7617
-144.1574
-13.4546
67.8139
-44.4608
-126.1802
-39.0300
192.4826
-67.4864
16.6392
82.7630
-146.8443
-90.7014
164.5643
126.8387
11.1138
-223.7063
-6.5059
67.2569
-139.3375
216.1606
2.0236
-98.0808
114.9628
-3.8427
-55.4881
-100.4703
98.0132
-7.2417
列主元高斯消元法，结果为
[0，…0] 即出现计算错误，弹出，系数矩阵奇异，
jacobi B的谱半径为86.3375，Jacobi迭代法 发散，结果为[NaN，…，NaN]
Gauss-S 方法 0.9984
1.0248
0.9247
1.0273
1.0528
1.0321
1.0024
0.9798
0.9683
0.9662
0.9706
0.9787
0.9882
0.9976
1.0059
1.0127
1.0177
1.0211
1.0228
1.0232
1.0224
1.0207
1.0183
1.0154
1.0122
1.0089
1.0055
1.0022
0.9990
0.9961
0.9934
0.9910
0.9889
0.9872
0.9857
0.9846
0.9838
0.9833
0.9831
0.9832
0.9834
0.9839
0.9846
0.9855
0.9865
0.9876
0.9889
0.9902
0.9916
0.9931
0.9946
0.9961
0.9976
0.9991
1.0006
1.0020
1.0034
1.0048
1.0061
1.0073
1.0084
1.0095
1.0104
1.0113
1.0121
1.0127
1.0133
1.0137
1.0141
1.0143
1.0145
1.0145
1.0144
1.0142
1.0139
1.0134
1.0129
1.0122
1.0115
1.0106
1.0096
1.0085
1.0073
1.0061
1.0047
1.0032
1.0016
0.9999
0.9982
0.9963
0.9944
0.9924
0.9903
0.9881
0.9858
0.9835
0.9811
0.9786
0.9760
0.9734
SOR 迭代法结果为
0.9978
1.0391
0.8492
1.1820
0.9015
1.1599
0.8585
1.1060
0.8456
1.0842
0.8658
1.0830
0.8979
1.0867
0.9281
1.0874
0.9510
1.0831
0.9660
1.0743
0.9744
1.0627
0.9781
1.0499
0.9788
1.0372
0.9778
1.0254
0.9763
1.0152
0.9748
1.0067
0.9738
1.0001
0.9735
0.9953
0.9740
0.9922
0.9753
0.9904
0.9772
0.9900
0.9796
0.9905
0.9826
0.9918
0.9858
0.9936
0.9892
0.9959
0.9927
0.9984
0.9961
1.0010
0.9994
1.0036
1.0025
1.0061
1.0053
1.0083
1.0078
1.0103
1.0100
1.0120
1.0117
1.0133
1.0130
1.0142
1.0138
1.0147
1.0142
1.0147
1.0142
1.0143
1.0137
1.0135
1.0127
1.0123
1.0113
1.0106
1.0094
1.0085
1.0072
1.0060
1.0045
1.0031
1.0014
0.9999
0.9980
0.9962
0.9942
0.9923
0.9901
0.9880
0.9857
0.9834
0.9810
0.9785
0.9759
0.9733

条件数来衡量问题的病态程度。条件数越大，病态可能越严重。

实验6.2 书上P211计算实习题2，其中N=100，第二小问改为用Jacobi迭代、G-S迭代、红黑排序的G-S迭代求解，并比较他们之间的收敛速度；进一步，用BSOR迭代求解，试找出最优松弛因子。习题见图片

解 ： ，由此，可知，x方向上[0,1]上分为了11个小格，y方向也一样。

五点差分方程为

Jacobi 结果为jacobi B的谱半径为0.95949
迭代次数=221

y_jacobi =

1.0139
1.0278
1.0415
1.0550
1.0681
1.0804
1.0914
1.1001
1.1053
1.1049
1.0278
1.0553
1.0825
1.1093
1.1352
1.1596
1.1813
1.1989
1.2101
1.2117
1.0415
1.0825
1.1231
1.1630
1.2017
1.2384
1.2714
1.2988
1.3176
1.3241
1.0550
1.1093
1.1630
1.2160
1.2677
1.3172
1.3625
1.4013
1.4299
1.4441
1.0681
1.1352
1.2017
1.2677
1.3328
1.3958
1.4547
1.5067
1.5477
1.5727
1.0804
1.1596
1.2384
1.3172
1.3958
1.4731
1.5471
1.6147
1.6712
1.7106
1.0914
1.1813
1.2714
1.3625
1.4547
1.5471
1.6379
1.7236
1.7993
1.8577
1.1001
1.1989
1.2988
1.4013
1.5067
1.6147
1.7236
1.8303
1.9296
2.0134
1.1053
1.2101
1.3176
1.4299
1.5477
1.6712
1.7993
1.9296
2.0577
2.1749
1.1049
1.2117
1.3241
1.4441
1.5727
1.7106
1.8577
2.0134
2.1749
2.3363

G-S 迭代结果为
G-S G的谱半径为0.92063
迭代次数=142

y_gs =

1.0139
1.0278
1.0415
1.0551
1.0682
1.0805
1.0914
1.1001
1.1053
1.1049
1.0278
1.0553
1.0826
1.1094
1.1353
1.1597
1.1814
1.1990
1.2101
1.2118
1.0415
1.0826
1.1232
1.1631
1.2019
1.2385
1.2716
1.2990
1.3177
1.3242
1.0551
1.1094
1.1631
1.2162
1.2679
1.3174
1.3627
1.4014
1.4300
1.4441
1.0682
1.1353
1.2019
1.2679
1.3330
1.3960
1.4549
1.5069
1.5479
1.5727
1.0805
1.1597
1.2385
1.3174
1.3960
1.4733
1.5473
1.6149
1.6713
1.7106
1.0914
1.1814
1.2716
1.3627
1.4549
1.5473
1.6381
1.7237
1.7994
1.8578
1.1001
1.1990
1.2990
1.4014
1.5069
1.6149
1.7237
1.8305
1.9297
2.0134
1.1053
1.2101
1.3177
1.4300
1.5479
1.6713
1.7994
1.9297
2.0578
2.1750
1.1049
1.2118
1.3242
1.4441
1.5727
1.7106
1.8578
2.0134
2.1750
2.3363

k_2 =

142

块SOR

CG 迭代