for(i=0, j=0; i<=N; i++, j++)
{
B[i][j] = 1/A[i][j] - A[i+1][j];
B[i][j-1] = -B[i][j] + B[i][j-1];
}
Mam taki kod i mój problem polega na tym żeby przy każdym obrocie pętli wartości w tablicy B odejmowały się maksymalnie od lewej strony. Nie za bardzo wiem jak to osiągnąć. Z góry dzięki za pomoc
Chcę zrobić uproszczonego Gaussa-Jordana, i chcę aby odwrotności pozycji [i][j] z macierzy A odejmowały się od macierzy B (jedynki na diagonali reszta 0). Macierze są NxN (7,8,9) i obecnie tylko diagonala i jedna kolumna w B wygląda poprawnie a pozostała część to dalej 0, a ma powstać macierz dolno-trójkątna.
Dla mnie ten zapis daje podejrzenie 2 wyjsc po za pamięć.... B[i][j-1] dla j = 0, wszystkie dla i = N - dodatkowo jeszcze widze tu i + 1 element (zakładam, że macierz jest N elementowa)
Mogę interować od 1 i wtedy wyliczać dla i-1 ale i tak nie wiem jak wypełniać od lewej strony aż do diagonali.
Nie wiem co kombinujesz (o co chodzi z tą uproszczoną metodą gaussa-jordana), wypisz sobie na sztywno ustalony przypadek np macierze 3x3, tak jakbyś krok po kroku liczył na kartce a potem próbuj to zwinąć w pętlę.
Właśnie cały czas próbuje ale jakoś mi to nie chce działać prawidłowo
Faktycznie mój błąd, chociaż nadal nie wiem jak to zwinąć do pętli. Mam coś takiego:
for(i=0; i<N; i++)
{
for(j=0; j<N; j++)
if(i == j)
{
B[i][j] = 1/A[i][j];
B[i+1][j] = -1/A[i][j] + B[i+1][j];
}
else if(j<i)
{
}
}
Nie mam pomysłu jak zrobić pozostałe elementy pod diagonalą
Tutaj masz zarys algorytmu. Aha, mój kod działa in-place, więc nie da się zastosować bezpośrednio do Twojego przypadku. Relatywnie niewielka modyfikacja powinna Ci też umożliwić rozwiązywanie układu równań.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/* Row-major order */
#define IDX(i, j, n) ((i) * (n) + (j))
#define ERR_MATRIX_SINGULAR 1
#ifdef WITH_PIVOTING
static void
swap_rows(double *matrix, ssize_t n, ssize_t p, ssize_t q)
{
/*
* Swaps row p-th and q-th starting on first non-zero entry.
*/
ssize_t ii = p < q ? p : q;
for (; ii < n; ii++) {
double tmp;
tmp = matrix[IDX(p, ii, n)];
matrix[IDX(p, ii, n)] = matrix[IDX(q, ii, n)];
matrix[IDX(q, ii, n)] = tmp;
}
}
#endif
int
gauss(double *matrix, ssize_t n)
{
ssize_t ii, jj, kk;
double a_inv;
for (ii = 0; ii < n; ii++) {
#ifdef WITH_PIVOTING
ssize_t pivot_row;
double pivot;
pivot_row = ii;
pivot = matrix[IDX(ii, ii, n)];
for (jj = ii + 1; jj < n; jj++) {
if (fabs(matrix[IDX(jj, ii, n)]) > fabs(pivot)) {
pivot_row = jj;
pivot = matrix[IDX(pivot_row, ii, n)];
}
}
/*
* Swap rows pivot_row and ii-th.
*/
if (ii != pivot_row)
swap_rows(matrix, n, ii, pivot_row);
#endif
if (matrix[IDX(ii, ii, n)] == 0.0)
return (ERR_MATRIX_SINGULAR);
a_inv = 1.0 / matrix[IDX(ii, ii, n)];
/* Scale rows below i-th */
for (kk = ii + 1; kk < n; kk++) {
double scale;
scale = a_inv * matrix[IDX(kk, ii, n)];
for (jj = ii; jj < n; jj++) {
matrix[IDX(kk, jj, n)] -= scale *
matrix[IDX(ii, jj, n)];
}
}
/* Scale i-th row. */
for (jj = ii; jj < n; jj++)
matrix[IDX(ii, jj, n)] *= a_inv;
}
return (0);
}
int
main()
{
size_t ii, jj;
size_t test_n = 4;
double test_matrix[] = {
4.0, -2.0, 4.0, -2.0,
3.0, 1.0, 4.0, 2.0,
2.0, 4.0, 2.0, 1.0,
2.0, -2.0, 4.0, 2.0
};
if (gauss(test_matrix, test_n) == ERR_MATRIX_SINGULAR) {
printf("Invalid matrix.\n");
}
for (ii = 0; ii < test_n; ii++) {
for (jj = 0; jj < test_n; jj++) {
printf("%02.5lf ", test_matrix[IDX(ii, jj, test_n)]);
}
printf("\n");
}
return (0);
}
Tu: https://introcs.cs.princeton.edu/java/95linear/GaussianElimination.java.html
Masz implementację w javie, do C++ chyba nie trudno przepisać.