Gauss-Seidel method convergence algorithm
3
$begingroup$
From Wikipedia:
The convergence properties of the Gauss–Seidel method
are dependent on the matrix A. Namely, the procedure
is known to converge if either:
1. A is symmetric positive-definite, or
2. A is strictly or irreducibly diagonally dominant.
Is there algorithm for transforming matrix to meet the convergence criteria or must I do it manually?
I am implementing Gauss-Seidel method in C++:
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
int n;
double eps = 0.001;
//yygnanma sherti
bool converge(double *xk, double *xkp)
{
double norm = 0;
for (int i = 0; i < n; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm) >= eps)
return false;
return true;
}
int main()
{
//olcheg
printf("Olchegi giriz: ");
scanf("%d", &n);
//koeffisientler uchin
double a[n][n];
for(int i=0; i<n; i++)
{
for(int j=0; j<n; j++)
{
double tmp;
printf("a[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
a[i][j] = tmp;
}
}
//azat agzalar uchin
double b[n];
for(int i=0; i<n; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i] = tmp;
}
double x[n];
double p[n];
for(int i=0; i<n; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for (int i = 0; i < n; i++)
p[i] = x[i];
for (int i = 0; i < n; i++)
{
double var = 0;
for (int j = 0; j < i; j++)
var += (a[i][j] * x[j]);
for (int j = i; j < n; j++)
var += (a[i][j] * p[j]);
x[i] = (b[i] - var) / a[i][i];
}
}while (!converge(x, p));
return 0;
}
Now I want to implement matrix transformation algorithm (if there is one) to meet convergence criteria.
linear-algebra numerical-methods numerical-linear-algebra
asked May 16 '14 at 18:56
torayefftorayeff
1164
$endgroup$
add a comment |
3
$begingroup$
From Wikipedia:
The convergence properties of the Gauss–Seidel method
are dependent on the matrix A. Namely, the procedure
is known to converge if either:
1. A is symmetric positive-definite, or
2. A is strictly or irreducibly diagonally dominant.
Is there algorithm for transforming matrix to meet the convergence criteria or must I do it manually?
I am implementing Gauss-Seidel method in C++:
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
int n;
double eps = 0.001;
//yygnanma sherti
bool converge(double *xk, double *xkp)
{
double norm = 0;
for (int i = 0; i < n; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm) >= eps)
return false;
return true;
}
int main()
{
//olcheg
printf("Olchegi giriz: ");
scanf("%d", &n);
//koeffisientler uchin
double a[n][n];
for(int i=0; i<n; i++)
{
for(int j=0; j<n; j++)
{
double tmp;
printf("a[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
a[i][j] = tmp;
}
}
//azat agzalar uchin
double b[n];
for(int i=0; i<n; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i] = tmp;
}
double x[n];
double p[n];
for(int i=0; i<n; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for (int i = 0; i < n; i++)
p[i] = x[i];
for (int i = 0; i < n; i++)
{
double var = 0;
for (int j = 0; j < i; j++)
var += (a[i][j] * x[j]);
for (int j = i; j < n; j++)
var += (a[i][j] * p[j]);
x[i] = (b[i] - var) / a[i][i];
}
}while (!converge(x, p));
return 0;
}
Now I want to implement matrix transformation algorithm (if there is one) to meet convergence criteria.
linear-algebra numerical-methods numerical-linear-algebra
asked May 16 '14 at 18:56
torayefftorayeff
1164
$endgroup$
$begingroup$
Where does your matrix A come from? Is it symmetric? Clearly arbitrarily "transforming" the matrix to meet some criteria may change the problem you are solving -- so if you tell me more about the problem I can make some suggestions.
$endgroup$
– RRL
May 16 '14 at 19:21
$begingroup$
A is not symmetric matrix. I am going to implement for abitrary matrix
$endgroup$
– torayeff
May 16 '14 at 19:38
1
$begingroup$
@torayeff you left out the fact that it will often converge without those conditions. Usually, if the system was small I would manually work it till there was as much diagonal dominance as possible.
$endgroup$
– bobbym
May 16 '14 at 19:51
$begingroup$
Why are you implementing this algorithm for a dense matrix?
$endgroup$
– Pawel Kowal
Jul 14 '16 at 10:01
add a comment |
3
3
3
$begingroup$
From Wikipedia:
The convergence properties of the Gauss–Seidel method
are dependent on the matrix A. Namely, the procedure
is known to converge if either:
1. A is symmetric positive-definite, or
2. A is strictly or irreducibly diagonally dominant.
Is there algorithm for transforming matrix to meet the convergence criteria or must I do it manually?
I am implementing Gauss-Seidel method in C++:
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
int n;
double eps = 0.001;
//yygnanma sherti
bool converge(double *xk, double *xkp)
{
double norm = 0;
for (int i = 0; i < n; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm) >= eps)
return false;
return true;
}
int main()
{
//olcheg
printf("Olchegi giriz: ");
scanf("%d", &n);
//koeffisientler uchin
double a[n][n];
for(int i=0; i<n; i++)
{
for(int j=0; j<n; j++)
{
double tmp;
printf("a[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
a[i][j] = tmp;
}
}
//azat agzalar uchin
double b[n];
for(int i=0; i<n; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i] = tmp;
}
double x[n];
double p[n];
for(int i=0; i<n; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for (int i = 0; i < n; i++)
p[i] = x[i];
for (int i = 0; i < n; i++)
{
double var = 0;
for (int j = 0; j < i; j++)
var += (a[i][j] * x[j]);
for (int j = i; j < n; j++)
var += (a[i][j] * p[j]);
x[i] = (b[i] - var) / a[i][i];
}
}while (!converge(x, p));
return 0;
}
Now I want to implement matrix transformation algorithm (if there is one) to meet convergence criteria.
linear-algebra numerical-methods numerical-linear-algebra
asked May 16 '14 at 18:56
torayefftorayeff
1164
$endgroup$
From Wikipedia:
The convergence properties of the Gauss–Seidel method
are dependent on the matrix A. Namely, the procedure
is known to converge if either:
1. A is symmetric positive-definite, or
2. A is strictly or irreducibly diagonally dominant.
Is there algorithm for transforming matrix to meet the convergence criteria or must I do it manually?
I am implementing Gauss-Seidel method in C++:
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
int n;
double eps = 0.001;
//yygnanma sherti
bool converge(double *xk, double *xkp)
{
double norm = 0;
for (int i = 0; i < n; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm) >= eps)
return false;
return true;
}
int main()
{
//olcheg
printf("Olchegi giriz: ");
scanf("%d", &n);
//koeffisientler uchin
double a[n][n];
for(int i=0; i<n; i++)
{
for(int j=0; j<n; j++)
{
double tmp;
printf("a[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
a[i][j] = tmp;
}
}
//azat agzalar uchin
double b[n];
for(int i=0; i<n; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i] = tmp;
}
double x[n];
double p[n];
for(int i=0; i<n; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for (int i = 0; i < n; i++)
p[i] = x[i];
for (int i = 0; i < n; i++)
{
double var = 0;
for (int j = 0; j < i; j++)
var += (a[i][j] * x[j]);
for (int j = i; j < n; j++)
var += (a[i][j] * p[j]);
x[i] = (b[i] - var) / a[i][i];
}
}while (!converge(x, p));
return 0;
}
Now I want to implement matrix transformation algorithm (if there is one) to meet convergence criteria.
linear-algebra numerical-methods numerical-linear-algebra
linear-algebra numerical-methods numerical-linear-algebra
asked May 16 '14 at 18:56
torayefftorayeff
1164
asked May 16 '14 at 18:56
torayefftorayeff
1164
asked May 16 '14 at 18:56
torayefftorayeff
1164
asked May 16 '14 at 18:56
torayefftorayeff
1164
asked May 16 '14 at 18:56
torayefftorayeff
1164
1164
$begingroup$
Where does your matrix A come from? Is it symmetric? Clearly arbitrarily "transforming" the matrix to meet some criteria may change the problem you are solving -- so if you tell me more about the problem I can make some suggestions.
$endgroup$
– RRL
May 16 '14 at 19:21
$begingroup$
A is not symmetric matrix. I am going to implement for abitrary matrix
$endgroup$
– torayeff
May 16 '14 at 19:38
1
$begingroup$
@torayeff you left out the fact that it will often converge without those conditions. Usually, if the system was small I would manually work it till there was as much diagonal dominance as possible.
$endgroup$
– bobbym
May 16 '14 at 19:51
$begingroup$
Why are you implementing this algorithm for a dense matrix?
$endgroup$
– Pawel Kowal
Jul 14 '16 at 10:01
add a comment |
$begingroup$
Where does your matrix A come from? Is it symmetric? Clearly arbitrarily "transforming" the matrix to meet some criteria may change the problem you are solving -- so if you tell me more about the problem I can make some suggestions.
$endgroup$
– RRL
May 16 '14 at 19:21
$begingroup$
A is not symmetric matrix. I am going to implement for abitrary matrix
$endgroup$
– torayeff
May 16 '14 at 19:38
1
$begingroup$
@torayeff you left out the fact that it will often converge without those conditions. Usually, if the system was small I would manually work it till there was as much diagonal dominance as possible.
$endgroup$
– bobbym
May 16 '14 at 19:51
$begingroup$
Why are you implementing this algorithm for a dense matrix?
$endgroup$
– Pawel Kowal
Jul 14 '16 at 10:01
$begingroup$
Where does your matrix A come from? Is it symmetric? Clearly arbitrarily "transforming" the matrix to meet some criteria may change the problem you are solving -- so if you tell me more about the problem I can make some suggestions.
$endgroup$
– RRL
May 16 '14 at 19:21
$begingroup$
Where does your matrix A come from? Is it symmetric? Clearly arbitrarily "transforming" the matrix to meet some criteria may change the problem you are solving -- so if you tell me more about the problem I can make some suggestions.
$endgroup$
– RRL
May 16 '14 at 19:21
$begingroup$
A is not symmetric matrix. I am going to implement for abitrary matrix
$endgroup$
– torayeff
May 16 '14 at 19:38
$begingroup$
A is not symmetric matrix. I am going to implement for abitrary matrix
$endgroup$
– torayeff
May 16 '14 at 19:38
1
1
$begingroup$
@torayeff you left out the fact that it will often converge without those conditions. Usually, if the system was small I would manually work it till there was as much diagonal dominance as possible.
$endgroup$
– bobbym
May 16 '14 at 19:51
$begingroup$
@torayeff you left out the fact that it will often converge without those conditions. Usually, if the system was small I would manually work it till there was as much diagonal dominance as possible.
$endgroup$
– bobbym
May 16 '14 at 19:51
$begingroup$
Why are you implementing this algorithm for a dense matrix?
$endgroup$
– Pawel Kowal
Jul 14 '16 at 10:01
$begingroup$
Why are you implementing this algorithm for a dense matrix?
$endgroup$
– Pawel Kowal
Jul 14 '16 at 10:01
add a comment |
1 Answer
1
active
oldest
votes
0
$begingroup$
I have found solution for my question and implemented it in C++:
/*
* www.twitter.com/torayeff
*
* Gauss-Seidel method:
* http://math.semestr.ru/optim/methodzeidel.php
* http://en.wikipedia.org/wiki/Gauss–Seidel_method
* http://ru.wikipedia.org/wiki/Метод_Гаусса_—_Зейделя
*
* To solve A*x = b transform to
* A(t)*A*x = A(t)*b <---> C*x = d
* (A(t) transpose of matrix A)
* (in this state, Gauss-Seidel method always congerges)
*/
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
#define MAXN 100
int N; //matrix dimension
double eps = 0.001;
double err = eps;
double A[MAXN][MAXN];
double b[MAXN][MAXN];
double At[MAXN][MAXN];
double C[MAXN][MAXN];
double d[MAXN][MAXN];
//convergence check
bool converge(double *xk, double *xkp)
{
double norm = 0;
for(int i=0; i<N; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm)>=eps)
return false;
err = eps;
return true;
}
int main()
{
printf("Input dimension N: ");
scanf("%d", &N);
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
double tmp;
printf("A[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
A[i][j] = tmp;
}
}
for(int i=0; i<N; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i][0] = tmp;
}
//transpose A
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
At[i][j] = A[j][i];
}
}
//multiply At*A = C
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
C[i][j] = 0;
for(int k=0; k<N; k++)
{
C[i][j] = C[i][j] + At[i][k]*A[k][j];
}
}
}
//multiply At*b = d
for(int i=0; i<N; i++)
{
for(int j=0; j<1; j++)
{
d[i][j] = 0;
for(int k=0; k<N; k++)
{
d[i][j] = d[i][j] + At[i][k]*b[k][j];
}
}
}
//Solve Gauss-Seidel method for C*x = d
double x[N]; //current values
double p[N]; //previous values
for(int i=0; i<N; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for(int i = 0; i<N; i++)
p[i] = x[i];
for(int i=0; i<N; i++)
{
double v = 0.0;
for(int j=0; j<i; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*x[j];
}
for(int j=i; j<N; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*p[j];
}
v += 1.0*d[i][0]/C[i][i];
x[i] = v;
}
}while (!converge(x, p));
//print solution
printf("Err. val.: %lfn", err);
for(int i=0; i<N; i++)
{
printf("%lf ", x[i]);
}
return 0;
}
answered May 17 '14 at 6:35
torayefftorayeff
1164
$endgroup$
add a comment |
Your Answer
StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
});
});
}, "mathjax-editing");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "69"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
noCode: true, onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f798300%2fgauss-seidel-method-convergence-algorithm%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
0
$begingroup$
I have found solution for my question and implemented it in C++:
/*
* www.twitter.com/torayeff
*
* Gauss-Seidel method:
* http://math.semestr.ru/optim/methodzeidel.php
* http://en.wikipedia.org/wiki/Gauss–Seidel_method
* http://ru.wikipedia.org/wiki/Метод_Гаусса_—_Зейделя
*
* To solve A*x = b transform to
* A(t)*A*x = A(t)*b <---> C*x = d
* (A(t) transpose of matrix A)
* (in this state, Gauss-Seidel method always congerges)
*/
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
#define MAXN 100
int N; //matrix dimension
double eps = 0.001;
double err = eps;
double A[MAXN][MAXN];
double b[MAXN][MAXN];
double At[MAXN][MAXN];
double C[MAXN][MAXN];
double d[MAXN][MAXN];
//convergence check
bool converge(double *xk, double *xkp)
{
double norm = 0;
for(int i=0; i<N; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm)>=eps)
return false;
err = eps;
return true;
}
int main()
{
printf("Input dimension N: ");
scanf("%d", &N);
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
double tmp;
printf("A[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
A[i][j] = tmp;
}
}
for(int i=0; i<N; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i][0] = tmp;
}
//transpose A
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
At[i][j] = A[j][i];
}
}
//multiply At*A = C
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
C[i][j] = 0;
for(int k=0; k<N; k++)
{
C[i][j] = C[i][j] + At[i][k]*A[k][j];
}
}
}
//multiply At*b = d
for(int i=0; i<N; i++)
{
for(int j=0; j<1; j++)
{
d[i][j] = 0;
for(int k=0; k<N; k++)
{
d[i][j] = d[i][j] + At[i][k]*b[k][j];
}
}
}
//Solve Gauss-Seidel method for C*x = d
double x[N]; //current values
double p[N]; //previous values
for(int i=0; i<N; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for(int i = 0; i<N; i++)
p[i] = x[i];
for(int i=0; i<N; i++)
{
double v = 0.0;
for(int j=0; j<i; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*x[j];
}
for(int j=i; j<N; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*p[j];
}
v += 1.0*d[i][0]/C[i][i];
x[i] = v;
}
}while (!converge(x, p));
//print solution
printf("Err. val.: %lfn", err);
for(int i=0; i<N; i++)
{
printf("%lf ", x[i]);
}
return 0;
}
answered May 17 '14 at 6:35
torayefftorayeff
1164
$endgroup$
add a comment |
0
$begingroup$
I have found solution for my question and implemented it in C++:
/*
* www.twitter.com/torayeff
*
* Gauss-Seidel method:
* http://math.semestr.ru/optim/methodzeidel.php
* http://en.wikipedia.org/wiki/Gauss–Seidel_method
* http://ru.wikipedia.org/wiki/Метод_Гаусса_—_Зейделя
*
* To solve A*x = b transform to
* A(t)*A*x = A(t)*b <---> C*x = d
* (A(t) transpose of matrix A)
* (in this state, Gauss-Seidel method always congerges)
*/
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
#define MAXN 100
int N; //matrix dimension
double eps = 0.001;
double err = eps;
double A[MAXN][MAXN];
double b[MAXN][MAXN];
double At[MAXN][MAXN];
double C[MAXN][MAXN];
double d[MAXN][MAXN];
//convergence check
bool converge(double *xk, double *xkp)
{
double norm = 0;
for(int i=0; i<N; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm)>=eps)
return false;
err = eps;
return true;
}
int main()
{
printf("Input dimension N: ");
scanf("%d", &N);
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
double tmp;
printf("A[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
A[i][j] = tmp;
}
}
for(int i=0; i<N; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i][0] = tmp;
}
//transpose A
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
At[i][j] = A[j][i];
}
}
//multiply At*A = C
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
C[i][j] = 0;
for(int k=0; k<N; k++)
{
C[i][j] = C[i][j] + At[i][k]*A[k][j];
}
}
}
//multiply At*b = d
for(int i=0; i<N; i++)
{
for(int j=0; j<1; j++)
{
d[i][j] = 0;
for(int k=0; k<N; k++)
{
d[i][j] = d[i][j] + At[i][k]*b[k][j];
}
}
}
//Solve Gauss-Seidel method for C*x = d
double x[N]; //current values
double p[N]; //previous values
for(int i=0; i<N; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for(int i = 0; i<N; i++)
p[i] = x[i];
for(int i=0; i<N; i++)
{
double v = 0.0;
for(int j=0; j<i; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*x[j];
}
for(int j=i; j<N; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*p[j];
}
v += 1.0*d[i][0]/C[i][i];
x[i] = v;
}
}while (!converge(x, p));
//print solution
printf("Err. val.: %lfn", err);
for(int i=0; i<N; i++)
{
printf("%lf ", x[i]);
}
return 0;
}
answered May 17 '14 at 6:35
torayefftorayeff
1164
$endgroup$
add a comment |
0
0
0
$begingroup$
I have found solution for my question and implemented it in C++:
/*
* www.twitter.com/torayeff
*
* Gauss-Seidel method:
* http://math.semestr.ru/optim/methodzeidel.php
* http://en.wikipedia.org/wiki/Gauss–Seidel_method
* http://ru.wikipedia.org/wiki/Метод_Гаусса_—_Зейделя
*
* To solve A*x = b transform to
* A(t)*A*x = A(t)*b <---> C*x = d
* (A(t) transpose of matrix A)
* (in this state, Gauss-Seidel method always congerges)
*/
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
#define MAXN 100
int N; //matrix dimension
double eps = 0.001;
double err = eps;
double A[MAXN][MAXN];
double b[MAXN][MAXN];
double At[MAXN][MAXN];
double C[MAXN][MAXN];
double d[MAXN][MAXN];
//convergence check
bool converge(double *xk, double *xkp)
{
double norm = 0;
for(int i=0; i<N; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm)>=eps)
return false;
err = eps;
return true;
}
int main()
{
printf("Input dimension N: ");
scanf("%d", &N);
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
double tmp;
printf("A[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
A[i][j] = tmp;
}
}
for(int i=0; i<N; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i][0] = tmp;
}
//transpose A
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
At[i][j] = A[j][i];
}
}
//multiply At*A = C
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
C[i][j] = 0;
for(int k=0; k<N; k++)
{
C[i][j] = C[i][j] + At[i][k]*A[k][j];
}
}
}
//multiply At*b = d
for(int i=0; i<N; i++)
{
for(int j=0; j<1; j++)
{
d[i][j] = 0;
for(int k=0; k<N; k++)
{
d[i][j] = d[i][j] + At[i][k]*b[k][j];
}
}
}
//Solve Gauss-Seidel method for C*x = d
double x[N]; //current values
double p[N]; //previous values
for(int i=0; i<N; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for(int i = 0; i<N; i++)
p[i] = x[i];
for(int i=0; i<N; i++)
{
double v = 0.0;
for(int j=0; j<i; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*x[j];
}
for(int j=i; j<N; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*p[j];
}
v += 1.0*d[i][0]/C[i][i];
x[i] = v;
}
}while (!converge(x, p));
//print solution
printf("Err. val.: %lfn", err);
for(int i=0; i<N; i++)
{
printf("%lf ", x[i]);
}
return 0;
}
answered May 17 '14 at 6:35
torayefftorayeff
1164
$endgroup$
I have found solution for my question and implemented it in C++:
/*
* www.twitter.com/torayeff
*
* Gauss-Seidel method:
* http://math.semestr.ru/optim/methodzeidel.php
* http://en.wikipedia.org/wiki/Gauss–Seidel_method
* http://ru.wikipedia.org/wiki/Метод_Гаусса_—_Зейделя
*
* To solve A*x = b transform to
* A(t)*A*x = A(t)*b <---> C*x = d
* (A(t) transpose of matrix A)
* (in this state, Gauss-Seidel method always congerges)
*/
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
#define MAXN 100
int N; //matrix dimension
double eps = 0.001;
double err = eps;
double A[MAXN][MAXN];
double b[MAXN][MAXN];
double At[MAXN][MAXN];
double C[MAXN][MAXN];
double d[MAXN][MAXN];
//convergence check
bool converge(double *xk, double *xkp)
{
double norm = 0;
for(int i=0; i<N; i++)
{
norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
}
if(sqrt(norm)>=eps)
return false;
err = eps;
return true;
}
int main()
{
printf("Input dimension N: ");
scanf("%d", &N);
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
double tmp;
printf("A[%d][%d] = ", i+1, j+1);
scanf("%lf", &tmp);
A[i][j] = tmp;
}
}
for(int i=0; i<N; i++)
{
double tmp;
printf("b[%d] = ", i+1);
scanf("%lf", &tmp);
b[i][0] = tmp;
}
//transpose A
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
At[i][j] = A[j][i];
}
}
//multiply At*A = C
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
C[i][j] = 0;
for(int k=0; k<N; k++)
{
C[i][j] = C[i][j] + At[i][k]*A[k][j];
}
}
}
//multiply At*b = d
for(int i=0; i<N; i++)
{
for(int j=0; j<1; j++)
{
d[i][j] = 0;
for(int k=0; k<N; k++)
{
d[i][j] = d[i][j] + At[i][k]*b[k][j];
}
}
}
//Solve Gauss-Seidel method for C*x = d
double x[N]; //current values
double p[N]; //previous values
for(int i=0; i<N; i++)
{
x[i] = 0;
p[i] = 0;
}
do
{
for(int i = 0; i<N; i++)
p[i] = x[i];
for(int i=0; i<N; i++)
{
double v = 0.0;
for(int j=0; j<i; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*x[j];
}
for(int j=i; j<N; j++)
{
double cij;
if(i==j) cij = 0;
else cij = -1.0*C[i][j]/C[i][i];
v += cij*p[j];
}
v += 1.0*d[i][0]/C[i][i];
x[i] = v;
}
}while (!converge(x, p));
//print solution
printf("Err. val.: %lfn", err);
for(int i=0; i<N; i++)
{
printf("%lf ", x[i]);
}
return 0;
}
answered May 17 '14 at 6:35
torayefftorayeff
1164
answered May 17 '14 at 6:35
torayefftorayeff
1164
answered May 17 '14 at 6:35
torayefftorayeff
1164
answered May 17 '14 at 6:35
torayefftorayeff
1164
1164
add a comment |
add a comment |
draft saved
draft discarded
Thanks for contributing an answer to Mathematics Stack Exchange!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f798300%2fgauss-seidel-method-convergence-algorithm%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
$begingroup$
Where does your matrix A come from? Is it symmetric? Clearly arbitrarily "transforming" the matrix to meet some criteria may change the problem you are solving -- so if you tell me more about the problem I can make some suggestions.
$endgroup$
– RRL
May 16 '14 at 19:21
$begingroup$
A is not symmetric matrix. I am going to implement for abitrary matrix
$endgroup$
– torayeff
May 16 '14 at 19:38
1
$begingroup$
@torayeff you left out the fact that it will often converge without those conditions. Usually, if the system was small I would manually work it till there was as much diagonal dominance as possible.
$endgroup$
– bobbym
May 16 '14 at 19:51
$begingroup$
Why are you implementing this algorithm for a dense matrix?
$endgroup$
– Pawel Kowal
Jul 14 '16 at 10:01