迭代法
这里一共提供了四种迭代法:
+ 雅可比迭代法
+ 高斯赛德迭代法
+ 超松弛迭代法(SOR)
+ 共轭迭代法
随机生成方程组
此处随机生成特征值服从独立同分布的[0,1]间的均匀分布的A矩阵,跟服从独立同分布的正态分布的b向量
算法设计
- A矩阵的实现,首先需要用rand得到一组特征值,将该组特征值通过diag函数生成对角阵Q,之后通过orth函数生成对称矩阵U,再通过UQU’即可得到对称正定矩阵
- b向量的实现同第一题,用randn函数即可
代码实现:
%n代表的是维度
b = randn([n(k) 1]);
Q = diag(rand([1 n(k)]));
U = orth(rand([n(k) n(k)]));
Jacobi 迭代法函数
算法设计
- 传入参数为矩阵A,向量b,以及维度n
- 传出参数为结果与时间(or相对误差数组),具体输出视操作而定
- 通过diag函数,tril函数,triu函数得到A矩阵的对角阵D,上三角矩阵U跟下三角矩阵L
- 通过以下公式得到迭代所需的B跟f
B = D(L+U);
f = D\b;`
- 通过以下公式进行迭代:
x = B*x0 + f;
- 迭代终止的判断,求出前后两次迭代的向量差的范数,直至小于1^-6
- 如果迭代次数超过2000,则视为发散,弹出错误
代码实现:
function [x,y] = Jacobi(n, A ,b)
%y = zeros(1000,1);
eps = 1.0e-6;
D = diag(diag(A));
L = -tril(A, -1);
U = -triu(A, 1);
B = D\(L+U);
f = D\b;
count = 1;
x0 = zeros(n,1);
x = B*x0 + f;
tic;
while norm(x-x0)> eps
x0 = x;
%y(count) = norm(x-A\b);
x = B*x0 + f;
count = count + 1;
if count > 2000
disp('error:该矩阵不收敛');
return;
end
end
toc;
y = toc;
disp(count);
end
Gauss-Seidel 迭代法函数
算法设计
- 传入参数为矩阵A,向量b,以及维度n
- 传出参数为结果与时间(or相对误差数组),具体输出视操作而定
- 通过diag函数,tril函数,triu函数得到A矩阵的对角阵D,上三角矩阵U跟下三角矩阵L
- 通过以下公式得到迭代所需的B跟f
B = (D - L)\U;
f = (D - L)\b;
- 通过以下公式进行迭代:
x = B*x0 + f;
- 迭代终止的判断,求出前后两次迭代的向量差的范数,直至小于1^-6
- 如果迭代次数超过2000,则视为发散,弹出错误
代码实现:
function [x, y]= Gauss_Seidel(n, A ,b)
%y = zeros(1000,1);
eps = 1.0e-6;
D = diag(diag(A));
L = -tril(A, -1);
U = -triu(A, 1);
B = (D - L)\U;
f = (D - L)\b;
count = 1;
x0 = zeros(n,1);
x = B*x0 + f;
tic;
while norm(x-x0) > eps
x0 = x;
% y(count) = norm(x-A\b);
x = B*x0 + f;
count = count + 1;
if count > 2000
disp('error:该矩阵不收敛');
return;
end
end
toc;
y = toc;
disp(count);
end
逐次超松弛迭代法函数
算法设计
- 传入参数为矩阵A,向量b,以及维度n
- 传出参数为结果与时间(or相对误差数组),具体输出视操作而定
- 通过diag函数,tril函数,triu函数得到A矩阵的对角阵D,上三角矩阵U跟下三角矩阵L
- 通过以下公式得到迭代所需的B跟f
B = (D - w*L) \ ((1-w) * D + w*U);
f = w * inv(D - w*L) * b;
- 通过以下公式进行迭代:
x = B*x0 + f;
- 迭代终止的判断,求出前后两次迭代的向量差的范数,直至小于1^-6
- 如果迭代次数超过2000,则视为发散,弹出错误
代码实现:
function [x, y] = SOR(n, A ,b, w)
y = zeros(1000,1);
eps = 1.0e-6;
D = diag(diag(A));
L = -tril(A, -1);
U = -triu(A, 1);
B = (D - w*L) \ ((1-w) * D + w*U);
f = w * inv(D - w*L) * b;
count = 1;
x0 = zeros(n,1);
x = B*x0 + f;
tic;
while norm(x-x0) > eps
x0 = x;
y(count) = norm(x-A\b);
x = B*x0 + f;
count = count + 1;
if count > 2000
disp('error:该矩阵不收敛');
return;
end
end
toc;
y = toc;
disp(count);
end
共轭梯度法函数
算法设计
- 传入参数为矩阵A,向量b,以及维度n
- 传出参数为结果与时间(or相对误差数组),具体输出视操作而定
- 通过以下公式得到迭代所需的r0跟p0
r0 = b - A * x0;
p0 = r0;
- 通过以下公式进行迭代:
a = r0'r0/(p0'A*p0);
x = x0 + a*p0;
r = r0 - a*A*p0;
B = r'*r/(r0'*r0);
p = r + B*p0;
- 迭代终止的判断,求出前后两次迭代的向量差的范数,直至小于1^-6
- 如果迭代次数超过2000,则视为发散,弹出错误
代码实现:
function [x, y] = CG(n, A , b)
y = zeros(1000,1);
eps = 1.0e-6;
x0 = zeros(n,1);
r0 = b - A * x0;
p0 = r0;
count = 1;
tic;
while norm(p0) > eps
a = r0'*r0/(p0'*A*p0);
x = x0 + a*p0;
r = r0 - a*A*p0;
B = r'*r/(r0'*r0);
p = r + B*p0;
p0 = p;
r0 = r;
x0 = x;
% y(count) = norm(x-A\b);
count = count + 1;
if count > 2000
disp('error:该矩阵不收敛');
return;
end
end
toc;
y = toc;
disp(count);
end
相对误差计算
算法设计
- 通过A\b得到近似解
- 通过求出近似解与得到的解的差向量的范数作为相对误差
代码实现
具体的代码实现,已经穿插在四种迭代法中