用Java实现 Gauss Seidel 方法
Gauss-Seidel(也称为逐次位移法)是一种数学计算方法,主要用于求线性代数系统的解。听说过雅可比方法吗?好吧,Gauss-Seidel 只不过是 Jacobi 的更好版本。
例子 :
Input:
Enter the number of variables in the equation:
2
Enter the augmented matrix:
1 2 3
6 5 4
Output:
6.0 5.0 4.0
1.0 2.0 3.0
X0 = {0.6666666666666666 1.1666666666666667 }
X1 = {-0.30555555555555564 1.652777777777778 }
X2 = {-0.7106481481481481 1.855324074074074 }
X3 = {-0.8794367283950617 1.9397183641975309 }
X4 = {-0.9497653034979425 1.9748826517489713 }
X5 = {-0.9790688764574759 1.9895344382287379 }
X6 = {-0.9912786985239483 1.9956393492619742 }
X7 = {-0.9963661243849785 1.9981830621924892 }
X8 = {-0.9984858851604077 1.9992429425802039 }
X9 = {-0.9993691188168363 1.999684559408418 }
X10 = {-0.9997371328403484 1.999868566420174 }
X11 = {-0.9998904720168117 1.9999452360084058 }
X12 = {-0.999954363340338 1.999977181670169 }
X13 = {-0.9999809847251406 1.9999904923625702 }
X14 = {-0.9999920769688085 1.9999960384844042 }
X15 = {-0.9999966987370034 1.9999983493685016 }
X16 = {-0.9999986244737512 1.9999993122368755 }
X17 = {-0.9999994268640631 1.9999997134320315 }
X18 = {-0.9999997611933598 1.9999998805966799 }
X19 = {-0.9999999004972331 1.9999999502486165 }
X20 = {-0.9999999585405137 1.9999999792702567 }
X21 = {-0.999999982725214 1.999999991362607 }
X22 = {-0.9999999928021724 1.9999999964010862 }
X23 = {-0.999999997000905 1.9999999985004524 }
X24 = {-0.999999998750377 1.9999999993751885 }
X25 = {-0.9999999994793237 1.9999999997396618 }
X26 = {-0.9999999997830514 1.9999999998915257 }
X27 = {-0.9999999999096048 1.9999999999548024 }
X28 = {-0.9999999999623352 1.9999999999811675 }
X29 = {-0.9999999999843061 1.999999999992153 }
X30 = {-0.9999999999934606 1.9999999999967302 }
X31 = {-0.9999999999972751 1.9999999999986375 }
X32 = {-0.9999999999988646 1.9999999999994322 }
X33 = {-0.9999999999995268 1.9999999999997633 }
X34 = {-0.9999999999998028 1.9999999999999014 }
X35 = {-0.9999999999999176 1.9999999999999587 }
X36 = {-0.9999999999999656 1.9999999999999827 }
X37 = {-0.9999999999999855 1.9999999999999927 }
X38 = {-0.9999999999999938 1.999999999999997 }
X39 = {-0.9999999999999973 1.9999999999999987 }
X40 = {-0.9999999999999988 1.9999999999999993 }
X41 = {-0.9999999999999993 1.9999999999999996 }
Gauss Seidel 的基本数学算法:
给定一组 n 个方程和 n 个未知数:
a11x1+a12x2+a13x3+......+a1nxn=c1
a21x1+a22x2+a23x3+......+a2nxn=c2
. .
. .
. .
a1nx1+a2nx2+a3nx3+......+annxn=cn
如果对角元素不为零,则每个方程都针对相应的未知数进行重写,即第一个方程在左侧用 x 1重写,第二个方程在左侧用 x 2重写,依此类推如下:
x1=(c1-a12x2-a13x3-....-a1nxn)/a11
x2=(c2-a21x1-a23x3-....-a2nxn)/a22
.
.
.
xn=(cn-an1x1-an2x2-....-an.n-1xn-1)/ann
.这些方程可以重写为以下形式:
x1=(c1-nj=1,j!=1∑a1jxj)/a11
x2=(c2-nj=1,j!=2∑a2jxj)/a22
.
.
.
x2=(cn-nj=1,j!=n∑anjxj)/ann
因此,对于任何行 一世 ;
xi=(ci- nj=i,j!=i∑aijxj)/aii , i=1,2,..,n
注:Gauss-Seidel 适用于严格对角占优或对称正定。
简而言之,对于线性方程组,说:
Ax=b
推导公式:
xi(k)=(bi-∑ji a i j x j(k-1))/aii
在矩阵方面,Gauss-Seidel 的定义:
x(k)=(D-L)-1(Ux(k-1)+b)
这里, D表示矩阵A的对角线部分, -L严格表示A的下三角部分, U表示A的严格上三角部分。
编程算法:
- 取方程中的变量数和增广矩阵的值
(通过附加两个给定矩阵的列形成)作为输入。 - 检查给定的增广矩阵是否对角占优。
- 尝试将给定的增广矩阵转换为对角占优矩阵。如果证明不是这样,则向用户返回适当的消息。
- 计算和打印最终结果以及迭代。
在Java中的实现:
Java
// Implementing Gauss Seidel Method in Java
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStreamReader;
import java.io.PrintWriter;
import java.util.Arrays;
import java.util.StringTokenizer;
class GFG {
// we set a max number of iterations to
// prevent an infinite loop
public static final int MAX_ITERATIONS = 100;
private double[][] M;
public GFG(double[][] matrix) { M = matrix; }
public void print() // printing
{
int n = M.length;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n + 1; j++)
System.out.print(M[i][j] + " ");
System.out.println();
}
}
// attempting to change a matrix to dominant
// if proved that it is not
public boolean transformToDominant(int r, boolean[] V,
int[] R)
{
int n = M.length;
if (r == M.length) {
double[][] T = new double[n][n + 1];
for (int i = 0; i < R.length; i++) {
for (int j = 0; j < n + 1; j++)
T[i][j] = M[R[i]][j];
}
M = T;
return true;
}
for (int i = 0; i < n; i++) {
if (V[i])
continue;
double sum = 0;
for (int j = 0; j < n; j++)
sum += Math.abs(M[i][j]);
if (2 * Math.abs(M[i][r]) > sum) {
// diagonally dominant?
V[i] = true;
R[r] = i;
if (transformToDominant(r + 1, V, R))
return true;
V[i] = false;
}
}
return false;
}
// method to check whether matrix is
// diagonally dominant or not
public boolean makeDominant()
{
boolean[] visited = new boolean[M.length];
int[] rows = new int[M.length];
Arrays.fill(visited, false);
return transformToDominant(0, visited, rows);
}
// method to find the solution of the matrix
// after all conditions are satisfied
public void solve()
{
int iterations = 0;
int n = M.length;
double epsilon = 1e-15;
double[] X = new double[n]; // Approximations
double[] P = new double[n]; // Prev
Arrays.fill(X, 0);
while (true) {
for (int i = 0; i < n; i++) {
double sum = M[i][n]; // b_n
for (int j = 0; j < n; j++)
if (j != i)
sum -= M[i][j] * X[j];
// Update xi to use in the next
// row calculation
X[i] = 1 / M[i][i] * sum;
}
System.out.print("X" + iterations + " = {");
for (int i = 0; i < n; i++)
System.out.print(X[i] + " ");
System.out.println("}");
iterations++;
if (iterations == 1)
continue;
boolean stop = true;
for (int i = 0; i < n && stop; i++)
if (Math.abs(X[i] - P[i]) > epsilon)
stop = false;
if (stop || iterations == MAX_ITERATIONS)
break;
P = (double[])X.clone();
}
}
public static void main(String[] args)
throws IOException
{
PrintWriter writer
= new PrintWriter(System.out, true);
int n = 2, k = 1;
double[][] M = new double[n][n + 1];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n + 1; j++)
M[i][j] = k++;
}
GFG gausSeidel = new GFG(M);
if (!gausSeidel.makeDominant()) {
// if it is found that a matrix cannot be
// changed into diagonally dominant then we
// return the message to the user
writer.println(
"The system isn't diagonally dominant: "
+ "The method cannot guarantee convergence.");
}
writer.println();
gausSeidel.print();
gausSeidel.solve();
}
}
输出
The system isn't diagonally dominant: The method cannot guarantee convergence.
1.0 2.0 3.0
4.0 5.0 6.0
X0 = {3.0 -1.2000000000000002 }
X1 = {5.4 -3.1200000000000006 }
X2 = {9.240000000000002 -6.192000000000002 }
X3 = {15.384000000000004 -11.107200000000004 }
X4 = {25.21440000000001 -18.97152000000001 }
X5 = {40.94304000000002 -31.554432000000016 }
X6 = {66.10886400000004 -51.68709120000003 }
X7 = {106.37418240000007 -83.89934592000006 }
X8 = {170.79869184000012 -135.4389534720001 }
X9 = {273.8779069440002 -217.90232555520015 }
X10 = {438.8046511104003 -349.84372088832026 }
X11 = {702.6874417766405 -560.9499534213124 }
X12 = {1124.899906842625 -898.7199254740999 }
X13 = {1800.4398509481998 -1439.1518807585599 }
X14 = {2881.3037615171197 -2303.843009213696 }
X15 = {4610.686018427392 -3687.348814741914 }
X16 = {7377.697629483828 -5900.958103587062 }
X17 = {11804.916207174125 -9442.7329657393 }
X18 = {18888.4659314786 -15109.57274518288 }
X19 = {30222.14549036576 -24176.51639229261 }
X20 = {48356.03278458522 -38683.62622766818 }
X21 = {77370.25245533636 -61895.00196426909 }
X22 = {123793.00392853818 -99033.20314283055 }
X23 = {198069.4062856611 -158454.3250285289 }
X24 = {316911.6500570578 -253528.12004564624 }
X25 = {507059.2400912925 -405646.192073034 }
X26 = {811295.384146068 -649035.1073168544 }
X27 = {1298073.2146337088 -1038457.3717069671 }
X28 = {2076917.7434139343 -1661532.9947311475 }
X29 = {3323068.989462295 -2658453.991569836 }
X30 = {5316910.983139672 -4253527.586511738 }
X31 = {8507058.173023475 -6805645.338418781 }
X32 = {1.3611293676837562E7 -1.088903374147005E7 }
X33 = {2.17780704829401E7 -1.742245518635208E7 }
X34 = {3.484491337270416E7 -2.787592949816333E7 }
X35 = {5.575186199632666E7 -4.460148839706133E7 }
X36 = {8.920297979412267E7 -7.136238263529813E7 }
X37 = {1.4272476827059627E8 -1.1417981341647702E8 }
X38 = {2.2835962983295405E8 -1.8268770266636324E8 }
X39 = {3.653754083327265E8 -2.923003254661812E8 }
X40 = {5.846006539323624E8 -4.6768052194588995E8 }
X41 = {9.353610468917799E8 -7.48288836313424E8 }
X42 = {1.496577675626848E9 -1.1972621393014784E9 }
X43 = {2.394524281602957E9 -1.9156194240823655E9 }
X44 = {3.831238851164731E9 -3.064991079731785E9 }
X45 = {6.12998216246357E9 -4.903985728770856E9 }
X46 = {9.807971460541712E9 -7.84637716723337E9 }
X47 = {1.569275433746674E10 -1.2554203468773392E10 }
X48 = {2.5108406940546783E10 -2.0086725551237427E10 }
X49 = {4.017345110547485E10 -3.2138760883179886E10 }
X50 = {6.427752176935977E10 -5.142201741428782E10 }
X51 = {1.0284403483157564E11 -8.227522786406052E10 }
X52 = {1.6455045573112103E11 -1.3164036458369684E11 }
X53 = {2.6328072917039368E11 -2.1062458333511496E11 }
X54 = {4.212491666732299E11 -3.36999333337384E11 }
X55 = {6.73998666677768E11 -5.391989333410144E11 }
X56 = {1.0783978666850288E12 -8.627182933468231E11 }
X57 = {1.7254365866966462E12 -1.3803492693561172E12 }
X58 = {2.7606985387152344E12 -2.208558830970988E12 }
X59 = {4.417117661944976E12 -3.533694129554781E12 }
X60 = {7.067388259112562E12 -5.65391060728885E12 }
X61 = {1.13078212145807E13 -9.04625697166336E12 }
X62 = {1.809251394332972E13 -1.4474011154662576E13 }
X63 = {2.8948022309328152E13 -2.3158417847461324E13 }
X64 = {4.631683569492565E13 -3.705346855593932E13 }
X65 = {7.410693711188164E13 -5.928554968950412E13 }
X66 = {1.1857109937901123E14 -9.48568795032078E13 }
X67 = {1.897137590064186E14 -1.517710072051337E14 }
X68 = {3.035420144102704E14 -2.4283361152821512E14 }
X69 = {4.8566722305643325E14 -3.8853377844514544E14 }
X70 = {7.770675568902939E14 -6.216540455122339E14 }
X71 = {1.2433080910244708E15 -9.946464728195755E14 }
X72 = {1.989292945639154E15 -1.591434356511322E15 }
X73 = {3.182868713022647E15 -2.5462949704181165E15 }
X74 = {5.092589940836236E15 -4.0740719526689875E15 }
X75 = {8.148143905337978E15 -6.518515124270381E15 }
X76 = {1.3037030248540764E16 -1.042962419883261E16 }
X77 = {2.0859248397665224E16 -1.668739871813218E16 }
X78 = {3.3374797436264364E16 -2.6699837949011492E16 }
X79 = {5.3399675898022984E16 -4.2719740718418392E16 }
X80 = {8.5439481436836784E16 -6.8351585149469432E16 }
X81 = {1.36703170298938864E17 -1.09362536239151104E17 }
X82 = {2.18725072478302208E17 -1.74980057982641792E17 }
X83 = {3.4996011596528358E17 -2.7996809277222688E17 }
X84 = {5.5993618554445376E17 -4.4794894843556301E17 }
X85 = {8.9589789687112602E17 -7.1671831749690086E17 }
X86 = {1.43343663499380173E18 -1.14674930799504141E18 }
X87 = {2.29349861599008282E18 -1.8347988927920663E18 }
X88 = {3.6695977855841326E18 -2.9356782284673065E18 }
X89 = {5.871356456934613E18 -4.697085165547691E18 }
X90 = {9.394170331095382E18 -7.5153362648763064E18 }
X91 = {1.5030672529752613E19 -1.2024538023802092E19 }
X92 = {2.4049076047604183E19 -1.9239260838083346E19 }
X93 = {3.847852167616669E19 -3.0782817340933358E19 }
X94 = {6.1565634681866715E19 -4.925250774549338E19 }
X95 = {9.850501549098675E19 -7.880401239278941E19 }
X96 = {1.5760802478557882E20 -1.2608641982846306E20 }
X97 = {2.5217283965692612E20 -2.017382717255409E20 }
X98 = {4.034765434510818E20 -3.227812347608655E20 }
X99 = {6.45562469521731E20 -5.164499756173848E20 }
时间复杂度: O(N 2 )