给出线性方程组 Hn*x = b,其中系数矩阵Hn为希尔伯特矩阵:
假设 x ∗ =(1, 1, . . . , 1)T,b = Hnx ∗。若取 n = 6,8, 10,分别用 Jacobi
迭代法及 SOR迭代(ω = 1, 1:25,1:5)求解,比较计算结果。
MATLAB源码如下,运行Demo_Jacobi_SOR.m文件,附件包含另外两个函数文件,分别为:Jacobi.m与SOR.m。
Demo_Jacobi_SOR.m
clear all
clc
n=[6 8 10];
for i=1:length(n)
H{i}=hilb(n(i));
size_H{i}=size(H{i},1);
x_true{i}=ones(1,size_H{i});
b{i}=x_true{i}*H{i};
x_ini{i}=zeros(1,size_H{i});
%accuracy=5*(10^-6);
accuracy=0.01;
end
%Jacobi
disp('Jacobi');
for i=1:length(n)
fprintf('The dimension is %d\n',n(i))
[x{i},k{i}]=Jacobi(H{i},b{i},x_ini{i},accuracy,x_true{i});
end
%SOR
disp('SOR');
w=1;
disp('SOR w=1');
for i=1:length(n)
fprintf('The dimension is %d\n',n(i))
[q{i},e{i}]=SOR(H{i},b{i},x_ini{i},accuracy,x_true{i},w);
end
w=1.25;
disp('SOR w=1.25');
for i=1:length(n)
fprintf('The dimension is %d\n',n(i))
[z{i},x{i}]=SOR(H{i},b{i},x_ini{i},accuracy,x_true{i},w);
end
w=1.5;
disp('SOR w=1.5');
for i=1:length(n)
fprintf('The dimension is %d\n',n(i))
[v{i},b{i}]=SOR(H{i},b{i},x_ini{i},accuracy,x_true{i},w);
end
%[x,k]=Jacobi(A,b,x_ini,accuracy,x_true);
Jacobi.m
function [x,k]=Jacobi(A,b,x_ini,accuracy,x_true)
% input: A—Left Matrix
% b—Right Matrix
% x_ini—initial value
% accuracy—calculation precision between solution and true value
% x_true—true valuer
% output: x—solution
% k—iteration-1
n=size(A,1);
uint16 k;
k=1;
x{1}=x_ini;
while(norm(x_true-x{k},Inf)>=accuracy)
k=k+1;
%disp(k-1);
for i=1:1:n
temp1=0;
for j=1:1:i-1
temp1=temp1+A(i,j)*x{k-1}(j);
end
temp2=0;
for j=i+1:1:n
temp2=temp2+A(i,j)*x{k-1}(j);
end
x{k}(i)=(b(i)-temp1-temp2)/A(i,i);
end
if (k==200000)
break;
end
end
fprintf('The iteration k=%d\n',k-1);
disp('The solution is');
disp(x{k});
end
SOR.m
function [x,k]=SOR(A,b,x_ini,accuracy,x_true,w)
% input: A—Left Matrix
% b—Right Matrix
% x_ini—initial value
% accuracy—calculation precision between solution and true value
% x_true—true value
% w—relaxation factor
% output: x—solution
% k—iteration-1
n=size(A,1);
uint16 k;
k=1;
x{1}=x_ini;
while(norm(x_true-x{k},Inf)>accuracy)
k=k+1;
%disp(k-1);
for i=1:1:n
temp1=0;
if(i>1)
for j=1:1:i-1
temp1=temp1+A(i,j)*x{k}(j);
end
end
temp2=0;
for j=i:1:n
temp2=temp2+A(i,j)*x{k-1}(j);
end
x{k}(i)=x{k-1}(i)+w*(b(i)-temp1-temp2)/A(i,i);
end
if (k==200000)
break;
end
end
fprintf('The iteration k=%d\n',k-1);
disp('The solution is');
disp(x{k});
end
MATLAB运行结果为:
Jacobi
The dimension is 6
The iteration k=199999
The solution is
Inf Inf Inf Inf Inf Inf
The dimension is 8
The iteration k=199999
The solution is
Inf Inf Inf Inf Inf Inf Inf Inf
The dimension is 10
The iteration k=199999
The solution is
Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
SOR
SOR w=1
The dimension is 6
The iteration k=14508
The solution is
0.9999 1.0016 0.9957 0.9994 1.0100 0.9933
The dimension is 8
The iteration k=42694
The solution is
1.0001 0.9984 1.0076 0.9900 0.9971 1.0073 1.0068 0.9926
The dimension is 10
The iteration k=25508
The solution is
1.0001 0.9980 1.0065 0.9985 0.9912 0.9970 1.0057 1.0091 1.0039 0.9900
SOR w=1.25
The dimension is 6
The iteration k=82840
The solution is
1.0000 0.9995 1.0036 0.9908 1.0100 0.9961
The dimension is 8
The iteration k=50752
The solution is
1.0001 0.9984 1.0073 0.9900 0.9983 1.0064 1.0060 0.9934
The dimension is 10
The iteration k=26267
The solution is
1.0001 0.9983 1.0048 1.0012 0.9900 0.9971 1.0053 1.0087 1.0038 0.9906
SOR w=1.5
The dimension is 6
The iteration k=174587
The solution is
1.0000 0.9994 1.0036 0.9908 1.0100 0.9961
The dimension is 8
The iteration k=52211
The solution is
1.0001 0.9985 1.0071 0.9900 0.9996 1.0049 1.0059 0.9939
The dimension is 10
The iteration k=199999
The solution is
1.0000 1.0007 0.9954 1.0113 0.9909 0.9986 0.9993 1.0049 1.0026 0.9962
运行结果分析:
用雅克比迭代法进行线性方程系数矩阵为希尔伯特矩阵的求解,解是发散的,最终趋于无穷大。
用SOR迭代法进行线性方程系数矩阵为希尔伯特矩阵的求解,解是收敛的,且收敛速度与矩阵的大小有关,但不是单调性的正相关或者负相关关系。