B
Also, ich habe das jetzt nochmal richtig umgebaut und erhalte den gleichen Effekt.
Meine Routinen sehen so aus: (das scheint nun eher Jacobi zu sein als Gauss-Seidel)
double matrix_sum(double** matrix, int index, double lines, double* vin)
{
double vsum = 0;
for (int j=0; j<lines; j++)
if (index != j)
vsum += vin[j] * matrix[index][j];
return vsum;
}
void matrix_solve_gauss_seidel(double** matrix, double* vector_result, double* vector_solution, int lines, double residual)
{
int i;
double residual_tmp_1;
double residual_tmp_2 = 1;
double* vector_old=(double*)malloc(lines * sizeof(double));
for (int i=0; i<lines; i++)
vector_old[i] = 0;
while (residual_tmp_2>residual)
{
residual_tmp_2 = 0;
for(i=0; i<lines; i++)
{
vector_solution[i] =
(vector_result[i] - matrix_sum(matrix, i, lines, vector_old)) / matrix[i][i];
residual_tmp_1 = vector_solution[i] - vector_old[i];
if(fabs(residual_tmp_1)>residual_tmp_2)
residual_tmp_2=fabs(residual_tmp_1);
}
for(i=0; i<lines; i++)
vector_old[i] = vector_solution[i];
}
free(vector_old);
}
Beachte: vector_result bleibt hier stabil der Zielvektor und vector_old wird auf 0|0 gestartet. Endet die Routine, liegt das Resultat in vector_solution.
Damit iteriert er sich bei der ersten 3*3-Matrix wunderschön ran, aber bei der
kleinen haut er ab! Und dann lese ich das:
http://www2.hs-esslingen.de/~koch/numerik/numerik_skript.pdf,
pdf Seite 14, programmiere diesen Mini und verfolge die Entwicklung
int main(void)
{
double xalt[2] = { 0.0, 0.0 };
double xneu[2];
for (int i=0; i<20; i++)
{
// ziel - sum matrix*alt ohne diag diag
xneu[0] = (0.5 -( 0 + 1*xalt[1]))/ 0.01;
xneu[1] = (1 -( 1*xalt[0] + 0 ))/ 1;
xalt[0]=xneu[0];
xalt[1]=xneu[1];
printf("step %d %f %f\n", i, xneu[0], xneu[1]);
}
return 0;
}
Aha! Es gibt da eine schöne Nebenbedingung für die Lösbarkeit:
siehe hier http://de.wikipedia.org/wiki/Diagonaldominante_Matrix
und - genau das ist diese Matrix nicht! Da ist der Fehler
