Il metodo iterativo di Gauss-Seidel si applica se i coefficienti delle incognite soddisfano il seguente requisito: i coefficienti sulla diagonale principale devono essere dominanti rispetto agli altri coefficienti sulla stessa riga. Cioè il valore assoluto dell'elemento sulla diagonale principale deve essere maggiore della somma dei valori assoluti degli altri elementi della riga.
Il metodo inizia assegnando valori di tentativo alle incognite e valutando ciascuna equazione. Per illustrare il procedimento, si consideri un sistema di 3 equazioni:
A11·X1 + A12·X2 + A13·X3 = B1
A21·X1 + A22·X2 + A23·X3 = B2
A31·X1 + A32·X2 + A33·X3 = B3
Si scelgano dei valori, ad esempio tutti 0, e li si assegnino alle incognite Xi (i = 1, 2, 3).
Così si può risolvere la prima equazione, ottenendo un nuovo valore dell'incognita, indicato con X'1:
X'1 = (B1 - A12·X2 - A13·X3 ) / A11
Usando questo nuovo valore dell'incognita X1 e il valore di X3, dalla seconda equazione, si calcola un nuovo valore di X2:
X'2 = (B2 - A21·X'1 - A23·X3 ) / A22
Con i due valori calcolati, X'1 e X'2, dalla terza equazione, si calcola un nuovo valore di X3:
X'3 = (B3 - A31·X'1 - A32·X'2 ) / A33
In questo modo sono stati ottenuti dei nuovi valori approssimati delle incognite. Iterando il procedimento si deve giungere a trovare valori delle incognite che differiscono da quelli calcolati nell'iterazione precedente di un'errore prefissato. In altri termini, il procedimento termina quando tra due iterazioni successive si osserva che i valori delle incognite sono sufficientemente vicini tra loro.
Bisogna definire come si può misurare la vicinanza tra due insiemi di valori delle incognite. In questo caso si propone di calcolare la somma dei valori assoluti della differenza tra i nuovi valori delle incognite e i precedenti valori delle incognite:
norma =
Al termine di una iterazione si deve stabilire se questa quantità è minore dell'approssimazione fissata. Se lo è il procedimento termina: le incognite sono state determinate con una accuratezza sufficiente, altrimenti si prosegue con un'altra iterazione.
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;
int main(int argc, char *argv[]) {
dichiarazioni delle variabili.
norma: di tipo doppia precisione, mantiene la somma dei valori assoluti delle differenze tra un nuovo valore dell'incognita e il precedente;
approssimazione: di tipo double, è un valore costante che potrebbe essere scelto in input, per fissare l'accuratezza dei risultati;
somma: di tipo double, mantiene il risultato della somma dei prodotti tra i coefficienti e le incognite (Aij·Xj)
tmp: di tipo double, una variabile per scambiare gli elementi durante la fase di ordinamento della matrice dei coefficienti, allo scopo di portare l'elemento dominante sulla diagonale principale;
b[100]: vettore dei termini noti. Viene scelto un valore massimo di 100 equazioni;
x[100]: vettore delle incognite;
Coeff[101][100]: la matrice dei coefficienti, con al massimo 100 equazioni, più la colonna dei termini noti;
N: di tipo intero, numero di equazioni.
double norma, approssimazione, somma, tmp, b[100], x[100];
double Coeff[100][101], max;
int N, pos;
FILE *pf;
approssimazione =0.000001;
Si assegnino i valori alle variabili che saranno le incognite:
x[0] = 2.0;
x[1] = 3.0;
x[2] = 4.0;
Creare, scegliendo numeri a caso, la matrice dei coefficienti, evitando che i coefficienti di una riga siano multipli dei coefficienti di un'altra riga:
Coeff[0][0] = -3.0;
Coeff[0][1] = +12.0;
Coeff[0][2] = +5.0;
Coeff[1][0] = 4.0;
Coeff[1][1] = -3.0;
Coeff[1][2] = -5.0;
Coeff[2][0] = 15.0;
Coeff[2][1] = +2.0;
Coeff[2][2] = -8.0;
Si scrivano 3 equazioni e si assegnino i risultati ai termini noti:
b[0] = Coeff[0][0]*x[0] + Coeff[0][1]*x[1] + Coeff[0][2]*x[2];
b[1] = Coeff[1][0]*x[0] + Coeff[1][1]*x[1] + Coeff[1][2]*x[2];
b[2] = Coeff[2][0]*x[0] + Coeff[2][1]*x[1] + Coeff[2][2]*x[2];
NOTA: l'elemento dominante di ciascuna riga non è stato posto sulla diagonale, allo scopo di proporre un metodo di ordinamento delle equazioni che porti l'elemento dominante sulla diagonale principale.
Scrivere le istruzioni per salvare in un file di testo la matrice dei coefficienti e il vettore dei termini noti, secondo la seguente organizzazione: La prima riga contiene un intero N che rappresenta il numero di equazioni, le restanti N righe del file contengono gli N coefficienti e il termine noto di ciascuna equazione. Ecco un esempio:
3 -3.0 +12.0 +5.0 50.0 +4.0 -3.0 -5.0 -21.0 +15.0 +2.0 -8.0 4.0
Assegnare il valore a N e creare il file:
N = 3;
pf = fopen("SisEq.txt","w");
if (pf==NULL) cout << "non sono riuscito a creare il file" << endl;
Scrivere dapprima il numero di righe del sistema di equazioni:
fprintf(pf, "%d\n", N);
Scrivere nel file di testo le righe dei coefficienti e, al termine di ciascuna riga, il rispettivo termine noto:
for (int i=0; i<N; i++) {
for (int j=0; j<N; j++)
fprintf(pf, "%lf\t", Coeff[i][j]);
fprintf(pf, "%lf\n", b[i]);
}
fclose(pf);
Dopo aver salvato i coefficienti e i termini noti del sistema di equazioni, il file viene chiuso, per riaprirlo in lettura. Le operazioni svolte fino a questo punto hanno lo scopo di collaudare il programma, in seguito potranno essere eliminate. In fase operativa il programma inizierà da questo punto, cioè con la lettura dei coefficienti dal file:
pf = fopen("SisEq.txt","r");
fscanf(pf, "%d", &N);
for (int i=0; i<N; i++)
for (int j=0; j<N+1; j++)
fscanf(pf, "%lf", &Coeff[i][j]);
fclose(pf);
Si noti che la lettura di un double dal file di testo avviene con lo specificatore di formato lf (long float, dove il tipo float specifica un numero decimale in semplice precisione).
La colonna dei termini noti è stata memorizzata nella (N+1)-ma colonna della matrice dei coefficienti.
Può iniziare la stampa di verifica del sistema di equazioni, per accertarsi che i dati siano stati acquisiti correttamente. Anche in questo caso, in fase operativa questa stampa potrà essere eliminata, ma sarà sempre opportuno prevedere delle verifiche sulla corretta acquisizione dei dati dal file.
cout << endl << endl << "Sistema di equazioni: " << endl;;
for (int i=0; i<N; i++) {
for (int j=0; j<N; j++) {
cout << Coeff[i][j] << "x" << j+1 << " \t";
}
cout << " = " << Coeff[i][N] << endl;
}
Si deve ricordare che il metodo iterativo di Gauss-Seidel converge alla soluzione se l'elemento dominante di ogni riga si trova sulla diagonale della matrice dei coefficienti. Una tale condizione non dovrebbe essere imposta all'utilizzatore del programma. Si propone un metodo di ordinamento delle righe della matrice dei coefficienti che potrebbe anche essere inefficace, in tal caso il programma deve comunicare che il sistema di equazioni non può essere risolto.
Ordinamento delle righe, per portare sulla diagonale principale l'elemento dominante di ogni riga:
for (int k=0; k<N; k++) { // contatore di riga
max = abs(Coeff[k][k]);
pos = k;
for (int j=k+1; j<N; j++) { // ricerca dell'elemento dominante nella colonna k
if (abs(Coeff[j][k]) > max) {
max = abs(Coeff[j][k]);
pos = j;
}
}
// scambia riga k con riga pos
for (int i=0; i<N+1; i++) {
tmp = Coeff[k][i];
Coeff[k][i] = Coeff[pos][i];
Coeff[pos][i] = tmp;
}
}
cout << endl;
cout << "Sistema di equazioni con elemento dominante sulla diagonale" << endl;
for (int i=0; i<N; i++) {
for (int j=0; j<N; j++) {
cout << Coeff[i][j] << "x" << j+1 << " \t";
}
cout << " = " << Coeff[i][N] << endl;
}
Il primo passo del metodo iterativo di Gauss consiste nell'assegnare un valore di tentativo a ognuna delle incognite Xi. Nell'istruzione seguente si svolge una doppia assegnazione dei termini noti sia al vettore x delle incognite sia al vettore b dei termini noti (quest'ultima assegnazione non è necessaria, ma si è preferito mettere in evidenza la posizione del termine noto all'interno dell'espressione di calcolo del nuovo valore dell'incognita).
Provare ad indovinare la soluzione:
for (int i=0; i<N; i++) x[i] = b[i] = Coeff[i][N];
Iterazioni per migliorare l'approssimazione:
do {
norma=0.0;
for (int i=0; i<N; i++) {
somma=0.0;
for (int j=0; j<N; j++) {
if (i==j) continue;
somma += Coeff[i][j]*x[j];
}
tmp = (b[i] - somma)/Coeff[i][i];
norma += abs(tmp - x[i]);
x[i] = tmp;
}
} while (norma > approssimazione);
cout << endl << "Soluzioni:";
for (int i=0; i<N; i++) cout << endl << "x[" << i << "] = " << x[i];
cout << endl;
system("PAUSE");
return EXIT_SUCCESS;
}
Individuare un controllo per riconoscere i casi in cui il sistema non può essere risolto; (Un elemento sulla diagonale è uguale a zero - il metodo non converge)
Modificare il programma affinchè l'approssimazione venga specificata in input dallo stesso file da cui si acquisiscono i coefficienti;
|
Equazioni di Kirchoff: |
Ridisponendo i termini delle equazioni alle maglie, il sistema di equazioni assume la seguente forma:
| V1 = | (R1 + R2 + R3)·i1 | -R2·i2 | + 0 | + 0 | + 0 | -R3·i6 |
| V2 = | -R2·i1 | (R2 + R4 + R5)·i2 | -R4·i3 | + 0 | -R5·i5 | + 0 |
| V3 = | 0 | -R4·i2 | (R4 + R6 + R7)·i3 | -R7·i4 | + 0 | + 0 |
| V4 = | 0 | + 0 | -R7·i3 | (R7 + R8 + R9)·i4 | -R9·i5 | + 0 |
| 0 = | 0 | + 0 | + 0 | -R9·i4 | (R9 + R10 + R11)·i5 | -R11·i6 |
| -V6 = | -R3·i1 | + 0 | + 0 | + 0 | -R11·i5 | (R3 + R11 + R12)·i6 |
Osservazioni:
La matrice dei coefficienti del sistema di equazioni viene ottenuta disponendo in colonna le correnti. Valgono le seguenti proprietà:
I termini noti coincidono con le tensioni dei generatori della maglia;
I coefficienti sulla diagonale principale sono dati dalla somma delle resistenze della maglia;
I termini non nulli della matrice dei coefficienti in coordinate (i, j) sono dati dal valore della resistenza comune alle due maglie i e j, cambiato di segno.
Determinare le correnti di maglia della rete di figura, usando il metodo di Gauss-Seidel per risolvere il sistema di equazioni.
Usare i seguenti valori:
| V1 = 10V | R1 = 100Ω | R7 = 75Ω | ||
| V2 = 8V | R2 = 180Ω | R8 = 150Ω | ||
| V3 = 12V | R3 = 500Ω | R9 = 400Ω | ||
| V4 = 15V | R4 = 120Ω | R10 = 250Ω | ||
| V6 = 5V | R5 = 220Ω | R11 = 80Ω | ||
| R6 = 50Ω | R12 = 40Ω |