Miglior approccio per risolvere matrice pentadiagonale

[Dânêl]

Ciao,
mi ritrovo a dover scrivere un programma in grado di risolvere una matrice penta o eptadiagonale ed ero alla ricerca di un metodo di risoluzione.
Ovviamente come requisiti richiesti sono di avere un basso costo computazionale, bassa occupazione di memoria e limitare per quanto possibile la propagazione degli errori

Cercando in rete ho trovato poco per non dire nulla, a parte qualche subrourtine fortran poco commentate....

Prendendo esempio dalla risoluzione delle matrici tridiagonali è un buon approccio utilizzare una versione modificata dell'algoritmo di Thomas?
Ossia applicando una forward substitution->back substitution (algoritmo di Thomas) e aggiungendo una nuova forward substitution?

Ci sono altre vie per ottenere lo stesso/migliore risultato?

[^TiGeRShArK^]

ehmm...
che intendi con risoluzione di una matrice?
la risoluzione delle equazioni lineari associate o qualcos'altro?

[Dânêl]

si scusate ho scritto una cosa per un altra...
ovviamente si tratta di risolverei un sistema di equazioni lineari associate ad una matrice

[^TiGeRShArK^]

credo che meglio di quest'algoritmo non puoi ottenere dato che ha complessità O(n).
Al massimo puoi migliorare l'occupazione di memoria nel caso tu non l'abbia già fatto memorizzando le diagonali in un oggetto che mantiene in memoria un vettore per ogni diagonale e fornisce le opportune trasformazioni di coordinate quando accedi ad un elemento della matrice.....

[Dânêl]

grazie della risposta.
Si le diagonali sono già memorizzate in vettori, anche perchè viene richiesto nella traccia

Quello che mi ha fatto sorgere qualche dubbio su questo metodo è che non viene utilizzato...
Se devo essere sincero non l'ho provato ancora non avendo a portata di mano un sistema pentadiagonale con relative soluzioni per fare il confronto. Però pensandoci non vedo il motivo per il quale non debba funzionare

Ho girato parecchio sulla rete prima di scrivere qui. Quello che ho trovato è stato un sito dove prima di applicare Thomas procedeva all'eliminazione opportuna di alcuni termini e una libreria (su netlib) che applica, anche se non vorrei sbagliarmi, il metodo di eliminazione di gauss (ovviamente senza pivot)

Dato i vantaggi di questo metodo (semplicità e velocità) non vorrei che applicando gauss ci fossero dei vantaggi tali da giustificare un O[ (n^3)/3+n^2-n/3)], come per esempio una limitata propagazione degli errori

PS ma il costo di 3 sostituzioni (backward forward) non dovrebbe essere O( 3n(n+1)/2)?

[^TiGeRShArK^]

Originariamente inviato da Dânêl

grazie della risposta.
Si le diagonali sono già memorizzate in vettori, anche perchè viene richiesto nella traccia

Quello che mi ha fatto sorgere qualche dubbio su questo metodo è che non viene utilizzato...
Se devo essere sincero non l'ho provato ancora non avendo a portata di mano un sistema pentadiagonale con relative soluzioni per fare il confronto. Però pensandoci non vedo il motivo per il quale non debba funzionare

Ho girato parecchio sulla rete prima di scrivere qui. Quello che ho trovato è stato un sito dove prima di applicare Thomas procedeva all'eliminazione opportuna di alcuni termini e una libreria (su netlib) che applica, anche se non vorrei sbagliarmi, il metodo di eliminazione di gauss (ovviamente senza pivot)

Dato i vantaggi di questo metodo (semplicità e velocità) non vorrei che applicando gauss ci fossero dei vantaggi tali da giustificare un O[ (n^3)/3+n^2-n/3)], come per esempio una limitata propagazione degli errori

PS ma il costo di 3 sostituzioni (backward forward) non dovrebbe essere O( 3n(n+1)/2)?

mmm...no.
Perchè dovrebbe essere 3n(n+1)?
L'implementazione che ho visto su wikipedia di thomas è O(n), se fosse come dici tu sarebbe quadratico...
invece la complessità di gauss che hai riportato è corretta dato che è O(n^3).
Ovviamente la propagazione potrebbe essere un problema, ma in generale non si può sapere dato che dipende dal problema specifico (in particolare dalla dimensione dei dati) e dai requisiti che hai....

[Dânêl]

quel valore l'avevo tirato fuori tenendo conto del costo computazionale delle sostituzioni che venivano ripetute 3 volte.
Evidentemente non è il metodo corretto per calcolarlo

Vorrà dire che ci dovrò riflettere meglio....

[^TiGeRShArK^]

Originariamente inviato da Dânêl

quel valore l'avevo tirato fuori tenendo conto del costo computazionale delle sostituzioni che venivano ripetute 3 volte.
Evidentemente non è il metodo corretto per calcolarlo

Vorrà dire che ci dovrò riflettere meglio....

appunto, ma è 3(n +1), ovvero O(n) per tre sostituzioni, non O(N^2) che ti verrebbe fuori da 3n(n + 1)...

[Dânêl]

Qualcosa non quadra

quell' n(n+1)/2 me lo sono ritrovato su alcune dispense.
Viene calcolato il costo computazionale del metodo di eliminazione di gauss e al costo della trasformazione del sistema lineare in triangolare superiore viene sommato quel valore, riportato come costo delle sostituzioni all'indietro.

Trattandosi di un sistema "particolare" risulta un valore diverso?

[^TiGeRShArK^]

non lo so a cosa si riferiva quel valore delle tue dispense...
ma vedendo da wikipedia questo codice per le tridiagonali:

void TridiagonalSolve(const double *a, const double *b, double *c, double *d, double *x, unsigned int n){
	int i;
 
	/* Modify the coefficients. */
	c[0] /= b[0];				/* Division by zero risk. */
	d[0] /= b[0];				/* Division by zero would imply a singular matrix. */
	for(i = 1; i < n; i++){
		double id = (b[i] - c[i-1] * a[i]);	/* Division by zero risk. */
		c[i] /= id;				/* Last value calculated is redundant. */
		d[i] = (d[i] - d[i-1] * a[i])/id;
	}
 
	/* Now back substitute. */
	x[n - 1] = d[n - 1];
	for(i = n - 2; i >= 0; i--)
		x[i] = d[i] - c[i] * x[i + 1];
}

Si hanno due cicli separati e non annidati, e dunque la complessità è lineare e non quadratica.....

[Dânêl]

Originariamente inviato da ^TiGeRShArK^

non lo so a cosa si riferiva quel valore delle tue dispense...

esattamente a quello che avevo detto in precedenza .
Ti riporto il testo

Inoltre il costo delle sostituzioni al'indietro è dato da n(n+1)/2 moltiplicazioni o divisioni

Non credo dipenda dal fatto che si stesse analizzando il metodo di Gauss perchè le sostituzioni all'indietro sono indipendenti da qualsiasi altro ciclo....e non credo di interpretare male il testo...

[^TiGeRShArK^]

si, ma nel caso dell'algoritmo preso da wikipedia la back-sostituion è O(n), con n tra l'altro uguale a 5 nel caso di matrice pentadiagonale....

[Dânêl]

si non dico che hai torto....bho
comunque dato che questo sembra essere il metodo migliore per risolvere un sistema con matrice pentadiagonale, sbattere la testa su questa discrepanza di costo computazionale sarebbe solo questione di principio

La risposta fondamentale è già stata data