Good morning,
I'm developing a code to perform the solution of a linear system in C using iterative methods (Gauss-Jacobi, Gauss-Seidel), the resolution method is correct, but at the end of code execution the value shown in the console is wrong (variable error, it shows this value -1. # QNAN), how can I fix this?
Below the code ...
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
//n: é a ordem do sistema.
//metodo: guarda o método a ser utilizado para a resulução: 1 = jacobi, 2 = gauss seidel.
//matriz a[][]: é a matriz dos coeficientes.
//matriz b[]: é a matriz dos termos independentes.
//matriz x[]: guarda os valores das variaveis após a resolução.
//e: é o valor de tolerancia de erro aceitavel no resultado, para ser utilizado como criterio de parada.
//k_max: é o numero de iterações maxima permitida. k_max=0 significa ilimitado (mas na realidade faz 100000 iterações, por segurança).
//O criterio de parada utilizado é o erro relativo: dr = max|xi(k) - xi(k-1)| / max|xi(k)|
//retorna 0 caso haja sucesso, 1 para erros de entrada de argumentos invalidos, 2 para erros de divisão por zero.
int iterativos_jacobi_gauss_seidel (int n, int metodo, float a[][n], float b[], float x[], float e, int k_max)
{
int k, i, j;
float s, dr, x_max;
float x_aux[n]; //variavel auxiliar que guarda os valores das variaveis xi's da iteracao anterior
if (n < 2 || e < 0.0 || k_max < 0) //testa se existe argumentos invalidos
{return 1;}
if (k_max == 0) //simula um numero de iteracoes "ilimitado" (maximo de 100000 iteracoes)
{k_max = 100000;}
switch (metodo) //aplica o metodo iterativo escolhido
{
//Jacobi ===============================================================
case 1:
{
for (k=1; k<=k_max; k++)
{
for (i=0; i<n; i++)
{x_aux[i] = x[i];} //guarda os valores de x da iteracao anterior
for (i=0; i<n; i++)
{
if (a[i][i] != 0.0) //Evita erros de divisao por zero
{
s = a[i][i]*x_aux[i]; //começa com o termo que nao deveria entrar nas subtrações...para anular sua participação
for (j=0; j<n; j++)
{s = s - a[i][j]*x_aux[j];}
s = b[i] + s;
x[i] = s / a[i][i];
}
else
{return 2;} //Houve erro de divisão por zero
}
//tolerancia de erro - busca o valor maximo entre |xi(k) - xi(k-1)|/|xi(k)|, testando se é menor do que erro aceitavel
dr = fabs(x[0] - x_aux[0]);
x_max = fabs(x[0]);
for (i=1; i<n; i++)
{
if (fabs(x[i] - x_aux[i]) > dr)
{dr = fabs(x[i] - x_aux[i]);}
if (fabs(x[i]) > x_max)
{x_max = fabs(x[i]);}
}
dr = dr / x_max;
if (dr < e) //se dr (erro relativo) for menor do que a tolerancia, a solução é aceitavel e termina a execucao do metodo
{k++; break;}
}
}
break;
//Gauss Seidel =========================================================
case 2:
{
for (k=1; k<=k_max; k++)
{
for (i=0; i<n; i++)
{x_aux[i] = x[i];} //guarda os valores de x da iteracao anterior
for (i=0; i<n; i++)
{
if (a[i][i] != 0.0) //Evita erros de divisao por zero
{
s = a[i][i]*x[i]; //começa com o termo que nao deveria entrar nas subtrações...para anular sua participação
for (j=0; j<n; j++)
{s = s - a[i][j]*x[j];}
s = b[i] + s;
x[i] = s / a[i][i];
}
else
{return 2;} //Houve erro de divisão por zero
}
//tolerancia de erro - busca o valor maximo entre |xi(k) - xi(k-1)|/|xi(k)|, testando se é menor do que erro aceitavel
dr = fabs(x[0] - x_aux[0]);
x_max = fabs(x[0]);
for (i=1; i<n; i++)
{
if (fabs(x[i] - x_aux[i]) > dr)
{dr = fabs(x[i] - x_aux[i]);}
if (fabs(x[i]) > x_max)
{x_max = fabs(x[i]);}
}
dr = dr / x_max;
if (dr < e) //se dr (erro relativo) for menor do que a tolerancia, a solução é aceitavel e termina a execucao do metodo
{k++; break;}
}
}
break;
//Caso não exista a opção de método passada para a função, implica em erro. Retorna 1, indicando esse erro.
default:
return 1;
} //fim do switch-case com os metodos iterativos
printf ("Foram realizadas %d iteracoes\nErro relativo = max |xi(k) - xi(k-1)| / max |xi(k)| = %f\n\n", k-1, dr);
return 0;
}
//==========================================================================================================
//funcao que implementa o criterio de convergencia das linhas.
//a[][] é a matriz dos coeficientes.
//retorna o valor de alfa. Se alfa < 1 o método de Jacobi ou Gauss Seidel gera uma sequencia convergente
//retona -1 caso haja elementos iguais a zero na diagonal principal, o que inviabiliza os metodos iterativos de Jacobi e Gauss Seidel
float criterio_linhas (int n, float a[][n])
{
float alfa, aux;
int i, j;
alfa=0.0;
for (i=0; i<n; i++)
{
if (a[i][i] != 0.0) //evita divisões por zero
{
aux = - fabs(a[i][i]); //inicia alfa com o negativo do termo que nao deve entrar no somatorio, para anular sua participacao
for (j=0; j<n; j++)
{
aux = aux + fabs(a[i][j]);
}
aux = aux/fabs(a[i][i]);
if (aux > alfa)
{alfa = aux;}
}
else
{return -1;} //Se houver zeros na diagonal principal retorna -1 indicando esse erro
}
return alfa;
}
//==========================================================================================================
//funcao que implementa o criterio de convergencia de Sassenfeld (somente para o metodo de Gauss Seidel).
//a[][] é a matriz dos coeficientes.
//retorna o valor de beta. Se beta < 1 o método de Gauss Seidel gera uma sequencia convergente
//retona -1 caso haja elementos iguais a zero na diagonal principal, o que inviabiliza o metodo iterativo de Gauss Seidel
float criterio_sassenfeld (int n, float a[][n])
{
float beta[n], beta_max, aux;
int i, j;
beta_max=0;
for (i=0; i<n; i++)
{beta[i] = 1.0;} //inicializa todos os betas com 1, elemento neutro da multiplicacao
for (i=0; i<n; i++)
{
if (a[i][i] != 0.0) //evita divisões por zero
{
aux = - fabs(a[i][i]); //inicia alfa com o negativo do termo que nao deve entrar no somatorio, para anular sua participacao
for (j=0; j<n; j++)
{
aux = aux + fabs(a[i][j])*beta[j];
}
beta[i] = aux/fabs(a[i][i]);
if (beta[i] > beta_max)
{beta_max = beta[i];}
}
else
{return -1;} //Se houver zeros na diagonal principal retorna -1 indicando esse erro
}
return beta_max;
}
//==========================================================================================================
//Função principal
int main()
{
int n, i, j, k_max, metodo, teste;
float e, alfa, beta;
char op, nome_arq[260];
FILE *fp;
do
{
printf ("Resolucao de sistemas lineares por metodos iterativos de Jacobi e Gauss Seidel\n\n");
do
{
printf ("Digite a ordem do sistema (n), no maximo n=100: ");
fflush(stdin); scanf ("%d",&n);
}
while (n < 2 || n > 100);
printf ("\nDigite a opcao do metodo iterativo:\n\n[1] Jacobi\n[2] Gauss Seidel\n\n");
do
{
printf ("Opcao: ");
fflush(stdin); scanf ("%d",&metodo);
}
while (metodo < 1 || metodo > 2);
do
{
printf ("\nDigite o numero maximo de iteracoes permitidas (digite 0 para ilimitado): ");
fflush(stdin); scanf ("%d",&k_max);
}
while (k_max < 0);
do
{
printf ("\nDigite a tolerancia de erro aceitavel: ");
fflush(stdin); scanf ("%f",&e);
}
while (e < 0.0);
//Cria as matrizes
float a[n][n], x[n], b[n]; //cria as matrizes dos coeficientes (A), das variaveis (X) e dos termos independentes (B)
//carrega do arquivo
do
{
printf ("\nDigite o caminho do arquivo de texto (com a extensao) que contem os dados:\n");
fflush(stdin);
scanf ("%259[^\n]",nome_arq);
fp=fopen(nome_arq,"r");
if (!fp)
{
printf ("\nErro na abertura do arquivo! Pressione [ENTER] para continuar\n");
fflush(stdin); op = getchar();
}
}
while(!fp);
//inicializa e lê os coeficientes das matrizes A e B e a solução inicial
for (i=0; i<n; i++)
{
for (j=0; j<n; j++)
{
a[i][j]=0.0;
fscanf (fp, "%f", &a[i][j]); //matriz A
}
b[i]=0.0;
fscanf (fp, "%f", &b[i]); //matriz B
}
for (i=0; i<n; i++)
{fscanf (fp, "%f", &x[i]);} //solucao inicial
fclose(fp);
//criterio de convergencia das linhas
system ("cls");
printf ("Deseja verificar a convergencia pelo criterio das linhas?\nDigite [S] para sim e [N] para nao: ");
fflush(stdin); op = getchar();
if (op == 's' || op == 'S')
{
alfa = criterio_linhas (n, a);
if (alfa >= 0)
{
if (alfa < 1)
{printf ("\nalfa = %f. alfa < 1, convergencia garantida!\n", alfa);}
else
{printf ("\nalfa = %f. alfa >= 1, nao existe convergencia garantida!\n", alfa);}
}
else
{
printf ("\nHa elementos iguais a zero na diagonal principal, isso inviabiliza\nos metodos iterativos!");
printf (" Tente reagrupar as equacoes para que nao existam zeros\nna diagonal principal!\n\nPressione [ENTER] para sair\n");
fflush(stdin); op = getchar(); return 0;
}
}
//criterio de convergencia de Sassenfeld (somente para o metodo de Gauss Seidel)
if (metodo == 2) //metodo de Gauss Seidel
{
printf ("\nDeseja verificar a convergencia pelo criterio de Sassenfeld?\nDigite [S] para sim e [N] para nao: ");
fflush(stdin); op = getchar();
if (op == 's' || op == 'S')
{
beta = criterio_sassenfeld (n, a);
if (beta >= 0)
{
if (beta < 1)
{printf ("\nbeta = %f. beta < 1, convergencia garantida!\n", beta);}
else
{printf ("\nbeta = %f. beta >= 1, nao existe convergencia garantida!\n", beta);}
}
else
{
printf ("\nHa elementos iguais a zero na diagonal principal, isso inviabiliza o metodo\nde Gauss Seidel!");
printf (" Tente reagrupar as equacoes para que nao existam\nzeros na diagonal principal!\n\nPressione [ENTER] para sair\n");
fflush(stdin); op = getchar(); return 0;
}
}
}
printf ("\nDados carregados! Pressione [ENTER] para continuar e resolver o sistema linear\n");
fflush(stdin); op = getchar(); system ("cls");
//Chama a função que aplica o metodo iterativo - retorna zero para teste caso haja sucesso
teste = iterativos_jacobi_gauss_seidel (n, metodo, a, b, x, e, k_max);
if (teste == 1) //se a função retornar 1, argumentos invalidos foram passados para a função
{
printf ("Erro! Argumentos invalidos foram digitados!\n\nPressione [ENTER] para sair\n");
fflush(stdin); op = getchar(); return 0;
}
if (teste == 2) //se a função retornar 2 houve erro de divisão por zero durante o calculo
{
printf ("Ocorreu erro de divisao por zero!\n");
printf ("Tente reagrupar as equacoes para que nao existam zeros na diagonal principal!\n\nPressione [ENTER] para sair\n");
fflush(stdin); op = getchar(); return 0;
}
//exibe o resultado na tela
printf ("Resultados:\n\n");
for (i=0; i<n; i++)
{printf ("x%d = %f\n", i+1, x[i]);}
printf ("\nDeseja reiniciar o programa? (digite [S] para sim, ou [N] para nao): ");
fflush(stdin); op = getchar(); system ("cls");
}
while (op == 's' || op == 'S');
return 0;
}
PS: The array to be loaded into the text file must be formatted as follows:
If this array is:
10 x1 + 2 x2 + 1 x3 = 7
1 x1 + 5 x2 + 1 x3 = -8
2 x1 + 3 x2 + 10 x3 = 6
It should be as follows in the file:
10 2 1 7
1 5 1 -8
2 3 10 6