Gauss Seidel求解线性方程组的Matlab程序
在matlab file exchange找了一圈，没找到合适的GS代码，于是自己写了一个, forward and backward
这里没有提及GS的理论，have fun!
- 2021.12.2-测试，该程序计算大规模稀疏矩阵时极其耗时，今日已优化，待程序整理完毕后上传
- 2021.12.3-updating, have optimized the GS algorithm to solve the large-scale sparsity problems. An example is given in .m file.
2021.12.6-计划，编写Blocked GS程序\ [已放弃]- 2022.6.14-测试，该GS程序应该利用矩阵分解和mldivide，可以实现较好的迭代效率，但更高效的迭代方式还得使用C/C++，该部分代码暂时未上传
- 2022.6.15- mexGS代码已上传
- 2022.6.17- 添加可并行的red-black GS，mexRBGS代码已上传
- 2022.6.18- 预告，Chebyshev semi iteration [懒了，先用matlab版凑合一下]
%% ********************************************************************
% % Matlab implementation of Gauss seidel to solve Ax = b iteratively
% % >>> Gauss seidel idea is derived here:
% % A*x = b -> (L + D + U)*x = b
% % where L is the lower triangular elements, U is the upper triangular
% % elements, D is the diagonal elements.
% % >>> then:
% % (L + D)*x_new = b - U*x_old
% % >>> obviously (Forward):
% % x_new = (L + D)^{-1} * [b - U*x_old]
% % >>> similarly (backward):
% % x_new = (D + U)^{-1} * [b - L*x_old]
% % *********************************************2021.11.30*By-LitBro*
function x = Gauss_seidel(A,x,b,tol,maxit,row_start,row_stop,row_step)
% NEED: N*N: A - (matrix, diag dominant)
% N*1: x - (initial choice, could be the start of GS)
% N*1: b - (right hand)
% float: tol - (tolerence)
% int: maxit - (max iteration number)
% int: row_start - (beginning of the sweep)
% int: row_stop - (end of the sweep)
% int: row_step - (stride used during the sweep, may be negative)
%
% OUT: N*1: x_new (approximate solution, sweep x after maxit or satisfy tol)
%
% ==== Example: Run the Code or Press F5 ==============================
A = [10 0 1 -5;1 8 -3 0;3 2 8 1;1 -2 2 7];
x = [1;1;1;1]; % initial random choice
b = [-7;11;23;17];
tol = 1d-3;
maxit = 1000; % max iteration number
row_start = 1;
row_stop = 4;
row_step = 1; % forward
%
% show: it = 8
% x = [ 0.3156; 2.0625; 1.9386; 2.4189]
% A*x = [-6.9993; 10.9998; 22.9999; 17.0000] \approx b
% =========================================================================
I = diag(A);
A = A - diag(diag(A)); % remove diag A
err = 1; it = 1; % start
while (err >= tol) ~= 0 && it <= maxit
x_old = x; % store x_old
for ii = row_start:row_step:row_stop
B = b(ii) - A(ii,:) * x;
x(ii) = B / I(ii);
end
err = sum(abs(x - x_old));
it = it+1;
end
end
最后更新于 2022年6月18日 --- 最初发表于 2021年11月30日
原创作者：LitBro
关于作者：百度、CSDN开禁用复制，**#
本文链接: [https://www.cnblogs.com/LitBro/p/15626226.html]
版权声明：本文采用 BY-NC-SA协议，转载或引用请注明出处!
关于后续：碍于学业不精，如有描述不当，还请见谅并非常感谢指出
