Risolvere un sistema lineare 4x4

[Eugenio]

Dovrei risolvere due sistemi lineari di 4 equazioni in 4 incognite,
di cui ho la matrice dei coefficienti (float) e il vettore dei termini
noti. Una matrice è del tipo

a 0 c 0
0 f g 0
i j k l
m n o p

e devo risolverla solo una volta. L'altra è del tipo

a b c 0
e f g 0
i j 0 l
m n 0 p

e devo risolverla decine di volte al secondo. Cosa potrei usare che
non mi complichi troppo la vita?

Ipotesi 1: appoggiarmi ad una libreria specializzata tipo Lapack++
(che, per quello che devo fare, mi sembra come prendere un aereo per
andare a comprare il pane, e le cui firme dei metodi non mi sembrano
molto user-friendly)

Ipotesi 2: usare il metodo di Cramer, cosa che però implica il calcolo
del determinante della matrice; e il determinante di una matrice 4x4 è
dispendioso da calcolare senza usare programmazione dinamica (a meno
che non debba implementare un metodo veloce... ma allora tanto vale
implementare una fattorizzazione LU con cui risolvere tutto il
sistema!)

Ipotesi 3: sviluppare a manina le formule (es. per sostituzione) e
scriverle direttamente...

Ci sono altre librerie C++ non commerciali che lo fanno in modo
"semplice"? Mi servirebbe un singolo metodo che risolve un sistema 4x4
(o che calcola il determinante in modo veloce). Se può servire,
sviluppo su Linux...

Consigli?

Grazie...
--
Narcolessico
http://binaryunit.blogspot.com

[Jonas Willard (AKA "Archeopteryx")]

Post by Eugenio
Consigli?

Mh, immagino che non ti aiuterà molto ma io avevo un
integrale orripilante che però sapevo essere calcolabile
analiticamente. Ho fatto digerire la formula al tool di
calcolo simbolico di matlab e il programma, ovviamente, va
che è una scheggia rispetto alle prime prove col calcolo
numerico anche se sul sorgente la formula occupa sei o
sette righe :) Penso che potresti fare lo stesso,
oltretutto, matlab restituisce una stringa col nome delle
variabili che gli hai dato tu, e fare copincolla
nell'editor del sorgente è cosa breve. Però magari non hai
matlab a portata di mano, e non so se tra le cose che fa
c'è la soluzione dei sistemi lineari.

ciao!

Apx.

[Eugenio]

On Mar 23, 6:03 pm, "Jonas Willard (AKA \"Archeopteryx\")"

Post by Jonas Willard (AKA "Archeopteryx")
Mh, immagino che non ti aiuterà molto ma io avevo un
integrale orripilante che però sapevo essere calcolabile
analiticamente. Ho fatto digerire la formula al tool di
calcolo simbolico di matlab e il programma, ovviamente, va
che è una scheggia rispetto alle prime prove col calcolo
numerico anche se sul sorgente la formula occupa sei o
sette righe :) Penso che potresti fare lo stesso,
oltretutto, matlab restituisce una stringa col nome delle
variabili che gli hai dato tu, e fare copincolla
nell'editor del sorgente è cosa breve. Però magari non hai
matlab a portata di mano, e non so se tra le cose che fa
c'è la soluzione dei sistemi lineari.

Grazie, se in poco trovo un equivalente in Octave seguirò il tuo
consiglio. Ti ricordi mica il nome esatto della funzione/comando?

PS: perchè chiedi l'annullamento automatico del messaggio dopo un tot
di giorni?

--
Narcolessico
http://binaryunit.blogspot.com

[pan]

In article

Post by Eugenio
On Mar 23, 6:03 pm, "Jonas Willard (AKA \"Archeopteryx\")"

Post by Jonas Willard (AKA "Archeopteryx")
Mh, immagino che non ti aiuterà molto ma io avevo un
integrale orripilante che però sapevo essere calcolabile
analiticamente. Ho fatto digerire la formula al tool di
calcolo simbolico di matlab e il programma, ovviamente, va
che è una scheggia rispetto alle prime prove col calcolo
numerico anche se sul sorgente la formula occupa sei o
sette righe :) Penso che potresti fare lo stesso,
oltretutto, matlab restituisce una stringa col nome delle
variabili che gli hai dato tu, e fare copincolla
nell'editor del sorgente è cosa breve. Però magari non hai
matlab a portata di mano, e non so se tra le cose che fa
c'è la soluzione dei sistemi lineari.

Grazie, se in poco trovo un equivalente in Octave seguirò il
tuoconsiglio. Ti ricordi mica il nome esatto della funzione/comando?

Per il calcolo simbolico, matlab e'... lasciamo perdere le volgarita',
e diciamo semplicemente "scomodo". Ti consiglio di provare
maple/mathematica (se hai un universita' che fornisce uno di questi,
come capitava a me qualche anno fa) oppure macsyma (se non sbaglio lo
spelling, dovrebbe essere l'equivalente opensource). Sono fatti
esplicitamente per il calcolo simbolico.

Ciao,
Marco

--
I'm using an evaluation license of nemo since 218 days.
You should really try it!
http://www.malcom-mac.com/nemo

[Eugenio]

Ho provato Maxima... è fantastico! Il determinante di una matrice
4x4 mi è venuto così:

a*(f*(k*p-l*o)-g*(j*p-l*n)+h*(j*o-k*n))+
-b*(e*(k*p-l*o)-g*(i*p-l*m)+h*(i*o-k*m))+
+c*(e*(j*p-l*n)-f*(i*p-l*m)+h*(i*n-j*m))+
-d*(e*(j*o-k*n)-f*(i*o-k*m)+g*(i*n-j*m))

Non è cortissima, ma devo ancora sistemare gli zeri al loro posto. Ho
chiesto la soluzione del sistema e ho ottenuto delle formule simili a
questa; tanto vale evitare Cramer e calcolare direttamente le
soluzioni. Per ottimizzare, userò formule specifiche per ogni matrice.

@?manu*: quali altri ng frequentiamo io e te? :)

Grazie a tutti per i consigli.

--
Narcolessico
http://binaryunit.blogspot.com

[pan]

In article

Post by Eugenio
Ho provato Maxima... è fantastico! Il determinante di una
a*(f*(k*p-l*o)-g*(j*p-l*n)+h*(j*o-k*n))+
-b*(e*(k*p-l*o)-g*(i*p-l*m)+h*(i*o-k*m))+
+c*(e*(j*p-l*n)-f*(i*p-l*m)+h*(i*n-j*m))+
-d*(e*(j*o-k*n)-f*(i*o-k*m)+g*(i*n-j*m))
Non è cortissima, ma devo ancora sistemare gli zeri al loro posto.
Ho chiesto la soluzione del sistema e ho ottenuto delle formule
simili a questa; tanto vale evitare Cramer e calcolare direttamente
le soluzioni. Per ottimizzare, userò formule specifiche per ogni
matrice.

Metti gli zeri nella forma simbolica (fai due versioni della cos), e
usa la funzione di semplificazione per semplificare.Poi fatti un po'
di raccoglimenti a mano, e il gioco e' fatto :)

E non usare Cramer, usa la sua funzione di soluzione dei sistemi.

--*PaN!*

--
I'm using an evaluation license of nemo since 219 days.
You should really try it!
http://www.malcom-mac.com/nemo

[pan]

Post by pan
In article

Post by Eugenio
Ho provato Maxima... è fantastico! Il determinante di una
a*(f*(k*p-l*o)-g*(j*p-l*n)+h*(j*o-k*n))+
-b*(e*(k*p-l*o)-g*(i*p-l*m)+h*(i*o-k*m))+
+c*(e*(j*p-l*n)-f*(i*p-l*m)+h*(i*n-j*m))+
-d*(e*(j*o-k*n)-f*(i*o-k*m)+g*(i*n-j*m))
Non è cortissima, ma devo ancora sistemare gli zeri al loro
posto. Ho chiesto la soluzione del sistema e ho ottenuto delle
formule simili a questa; tanto vale evitare Cramer e calcolare
direttamente le soluzioni. Per ottimizzare, userò formule specifiche
per ogni
matrice.

Metti gli zeri nella forma simbolica (fai due versioni della cos), e
usa la funzione di semplificazione per semplificare.Poi fatti un
po'di raccoglimenti a mano, e il gioco e' fatto :)
E non usare Cramer, usa la sua funzione di soluzione dei sistemi.

Ovviamente non ho letto la seconda frase, in cui dici esattamente
quello che ti ho scritto :)

--*PaN!*


--
I'm using an evaluation license of nemo since 219 days.
You should really try it!
http://www.malcom-mac.com/nemo

[Max]

Eugenio ha scritto:

Ciao,

Solo una osservazione. Il calcolo simbolico e' potente ma deve essere
usato con attenzie. La formula che hai trovato contiene molti termini
del tipo

(a*b - c*d)

e termini del tipo

a *(c*d -e*f) - g *(h*i-l*m)

questi termini possono avere un problema detto "cancellazione numerica",
la perdita di precisione nel sottrarre quantita' "vicine". La
letteratura e' piena di esempi e consigli. In double e 4x4 avrai
problemi solo per matrici "patologiche" ma salendo di ordine scoprirai
presto che e' il metodo errato. Volendo potrei crearti una matrice "mal
condizionata", ad esempio una matric di Hilbert (nell'elemento (i,j)
inserisci il numero razionale 1 / (i + j - 1) ) che ti darebbe seri
problemi di precisione anche per un sistema 4x4. Insomma la soluzione
proposta non e' "robusta" anche se "normalmente ti funzionera'". Le
matrici di Hilbert le incontri nei problemi di "minimi quadrati", che e'
un tipo di calcolo piuttosto comune.

Il numero di condizionamento (prodotto della norma della matrice per la
norma del'inversa della matrice) di una matrice di Hilbert 4x4 e' 1550.
Una piccola variazione nei dati può causare una variazione di circa 1500
volte nel risultato finale. Per matrici mal condizionate il metodo
eficace e' spesso riformulare il problema. Questo per dirti che una
"ricetta numerica" che vada sempre bene non esiste.

Dipende dal problema che stai affrontando. Se e' un esercizio OK ma se
stai controllando un elettromedicale non vorrei mai finirci "sotto" :-))

Consulta un testo di analisi numerica od usa librerie scritte da chi ha
fatto "studi regolari". Se non vuoi usare Lapack prova GSL.

Post by Eugenio
Ho provato Maxima... è fantastico! Il determinante di una matrice
a*(f*(k*p-l*o)-g*(j*p-l*n)+h*(j*o-k*n))+
-b*(e*(k*p-l*o)-g*(i*p-l*m)+h*(i*o-k*m))+
+c*(e*(j*p-l*n)-f*(i*p-l*m)+h*(i*n-j*m))+
-d*(e*(j*o-k*n)-f*(i*o-k*m)+g*(i*n-j*m))
Non è cortissima, ma devo ancora sistemare gli zeri al loro posto. Ho
chiesto la soluzione del sistema e ho ottenuto delle formule simili a
questa; tanto vale evitare Cramer e calcolare direttamente le
soluzioni. Per ottimizzare, userò formule specifiche per ogni matrice.
@?manu*: quali altri ng frequentiamo io e te? :)
Grazie a tutti per i consigli.
--
Narcolessico
http://binaryunit.blogspot.com

[pan]

Post by Max
Ciao,
Solo una osservazione. Il calcolo simbolico e' potente ma deve
essere usato con attenzie. La formula che hai trovato contiene molti
termini del tipo
(a*b - c*d)
e termini del tipo
a *(c*d -e*f) - g *(h*i-l*m)
questi termini possono avere un problema detto "cancellazione
numerica", la perdita di precisione nel sottrarre quantita' "vicine".

Sottoscrivo.
Anche se mi par difficile trovare queste condizioni nelle matrici
citate dall'op, e' una cosa di cui tenere conto.

--*PaN!*

--
I'm using an evaluation license of nemo since 219 days.
You should really try it!
http://www.malcom-mac.com/nemo

[Jonas Willard (AKA "Archeopteryx")]

Post by pan
Sottoscrivo.

Beh, pure io. Avevo suggerito il calcolo simbolico dando
per scontato che ce ne fossero le condizioni. Esiste un
tipo di decomposizione (mi pare si chiami "singular value"
ma devo rivedere gli antichi ricordi) per sistemi
malissimo condizionati e a quel punto, se quello è il
caso, si può sempre correre ai ripari con tecniche
numeriche adeguate.

ciao ciao

Apx.

[Aka]

In effetti non ho specificato nulla sull'applicazione. Nella prima
matrice i coefficenti sono tutti numeri naturali minori di 1300; i
valori delle incognite nella prima matrice diventano i coefficienti
del secondo sistema. Per il primo, sistemando gli zeri, Maxima ha
prodotto:

x = (c1*(h*k*n-f*k*p)+b*(c2*k*p-c3*g*p+h*(c3*o-c4*k))+d*(f*(c4*k-c3*o)-
c2*k*n+c3*g*n)) /
(b*(e*k*p-g*i*p+h*(i*o-k*m))+d*(f*(k*m-i*o)-e*k*n+g*i*n))
y = (...)

c'è una sola divisione tra numeri interi, quindi quello che mi
preoccupa sono solo eventuali overflow. Devo ancora controllare il
range effettivo dei numeri in input, ma credo che con un long con
segno dovrei farcela con ampio margine di sicurezza (nel caso peggiore
moltiplicherò 4 numeri minori di 1300 tra loro).

Per la seconda matrice, ottengo:

x=(c*(c2*(j*p-l*n)+f*(c4*l-c3*p))+c1*g*(l*n-j*p)+b*g*(c3*p-c4*l)) /
(c*(e*(j*p-l*n)+f*(l*m-i*p))+a*g*(l*n-j*p)+b*g*(i*p-l*m))
y = (...)

anche qui, una sola divisione. Il caso di parametri uguali a zero è
rarissimo (e, comunque, non ci possono essere più di k zeri
contemporaneamente, credo con k=2), quindi sia denominatore che
numeratore non dovrebbero (de)crescere fino a dare problemi. Infine
devo arrotondare ad un intero, quindi la precisione che mi interessa è
solo al primo decimale.

Mi rendo conto che il mio ragionamento è zeppo di condizionali e
controlli "a occhio", ma al momento non serve una grande precisione né
accuratezza. Al momento mi serve un sistema appena funzionante, ma
ottenuto quello riprenderò il mio bel libro di analisi numerica e
andrò a calcolare il numero di condizionamento delle matrici e la
precisione del calcolo in base al range di input, che è sottoposto a
diversi vincoli.

PS: Max può stare tranquillo: niente elettromedicali per ora, solo
finger tracking ;)

--
Narcolessico
http://binaryunit.blogspot.com

[pan]

In article

Post by Aka
In effetti non ho specificato nulla sull'applicazione. Nella prima
matrice i coefficenti sono tutti numeri naturali minori di 1300; i
valori delle incognite nella prima matrice diventano i
coefficientidel secondo sistema. Per il primo, sistemando gli zeri,

[...]

Post by Aka
c'è una sola divisione tra numeri interi, quindi quello che mi
preoccupa sono solo eventuali overflow. Devo ancora controllare il
range effettivo dei numeri in input, ma credo che con un long con
segno dovrei farcela con ampio margine di sicurezza (nel caso
peggioremoltiplicherò 4 numeri minori di 1300 tra loro).

Semplice allora.
per rappresentare un numero fino a 1300 ci vogliono 11 bit, 11*4=44
bit. Te la cavi con numeri a 64 bit fino alla divisione, o se
rappresenti il tutto in forma frazionaria avrai il risultato _esatto_.
L'unico caso cattivo e' quello in cui il determinante sia 0, ma in
quel caso non c'e' nulla da fare.

Se riuscissi a tenere in forma razionale anche il secondo sistema,
potresti giungere ad un risultato razionale esatto (controllando
sempre il numero di bit necessari ad ogni passaggio).

Post by Aka
Il caso di parametri uguali a zero è
rarissimo (e, comunque, non ci possono essere più di k zeri
contemporaneamente, credo con k=2), quindi sia denominatore che
numeratore non dovrebbero (de)crescere fino a dare problemi.

Il numero di zeri purtroppo non influisce sul numero di
condizionamento (anzi, se non rende il sistema indeterminato, uno zero
in piu' e' facile che migliori le cose).

Post by Aka
Mi rendo conto che il mio ragionamento è zeppo di condizionali e
controlli "a occhio", ma al momento non serve una grande precisione
né accuratezza. Al momento mi serve un sistema appena funzionante, ma
ottenuto quello riprenderò il mio bel libro di analisi numerica e
andrò a calcolare il numero di condizionamento delle matrici e la
precisione del calcolo in base al range di input, che è sottoposto
adiversi vincoli.

Se riesci a lavorare solo con naturali, non ti servira' nessun numero
di condizionamento. Vedi cosa riesci a fare :)

Post by Aka
PS: Max può stare tranquillo: niente elettromedicali per ora,
solofinger tracking ;)

Potrebbe essere sempre usato in futuro per trackare la mano di un
chirurgo ed eseguire in modo robotizzato un'operazione a cuore aperto
:)

--*PaN!*

--
I'm using an evaluation license of nemo since 219 days.
You should really try it!
http://www.malcom-mac.com/nemo

[Eugenio]

Mannaggia... Ho postato dall'account gmail l'ultima volta... Speriamo
di non essere sommerso di spam!

Post by pan
L'unico caso cattivo e' quello in cui il determinante sia 0, ma in
quel caso non c'e' nulla da fare.

Il sistema ha sempre una e una sola soluzione, quindi da questo punto
di vista non c'è problema.

Post by pan

Post by Aka
Il caso di parametri uguali a zero è
rarissimo (e, comunque, non ci possono essere più di k zeri
contemporaneamente, credo con k=2), quindi sia denominatore che
numeratore non dovrebbero (de)crescere fino a dare problemi.

Il numero di zeri purtroppo non influisce sul numero di
condizionamento (anzi, se non rende il sistema indeterminato, uno zero
in piu' e' facile che migliori le cose).

Però influisce sulle frazioni: se b=d=0, ad esempio, nella prima
formula ho un denominatore nullo (e ci sono altri casi, più complessi
per cui si azzera). Ma se b e d sono nulli, per la prima equazione
anche c1 sarà nullo; quindi anche il numeratore si annulla. Dato che
so che il sistema ha sempre soluzione, il denominatore non dovrebbe
mai annullarsi se il numeratore non è nullo. A meno di perdite di
precisione che, nel primo sistema, non ho perchè lavoro con i
naturali. Sul secondo, invece, qualche controllino non guasterebbe...

Post by pan

Post by Aka
PS: Max può stare tranquillo: niente elettromedicali per ora,
solofinger tracking ;)

Potrebbe essere sempre usato in futuro per trackare la mano di un
chirurgo ed eseguire in modo robotizzato un'operazione a cuore aperto
:)

Beh, prima o poi di qualcosa si muore :D

--
Narcolessico
http://binaryunit.blogspot.com

[pan]

In article

Post by Eugenio
Mannaggia... Ho postato dall'account gmail l'ultima volta...
Speriamodi non essere sommerso di spam!

Post by pan
L'unico caso cattivo e' quello in cui il determinante sia 0, ma in
quel caso non c'e' nulla da fare.

Il sistema ha sempre una e una sola soluzione, quindi da questo
puntodi vista non c'è problema.

[...]

Post by Eugenio
se b=d=0, ad esempio, nella prima
formula ho un denominatore nullo (e ci sono altri casi, più
complessi per cui si azzera). Ma se b e d sono nulli, per la prima
equazione anche c1 sarà nullo; quindi anche il numeratore si annulla.
Dato che so che il sistema ha sempre soluzione, il denominatore non
dovrebbe mai annullarsi se il numeratore non è nullo. A meno di
perdite di precisione che, nel primo sistema, non ho perchè lavoro
con inaturali. Sul secondo, invece, qualche controllino non
guasterebbe...

Se il denominatore s'annulla e il numeratore pure significhera' anche
che il sistema ha soluzione, ma non e' unica, il che e' in contrasto
con quanto hai detto prima (una e una sola soluzione). Inoltre il
fatto che c1 sia nullo se b e d sono nulli, mi fa sospettare che ci
sia qualcosa in piu' rispetto ad un sistema 4x4...

Non hai in realta' 5 incognite e un vincolo non lineare (chesso'
quadratico) vero?!? :)

--*PaN!*

--
I'm using an evaluation license of nemo since 219 days.
You should really try it!
http://www.malcom-mac.com/nemo

[Eugenio]

Post by pan
Se il denominatore s'annulla e il numeratore pure significhera' anche
che il sistema ha soluzione, ma non e' unica, il che e' in contrasto
con quanto hai detto prima (una e una sola soluzione).

Sì, hai ragione, ma quello che temo è che numeratore e denominatore si
annullino per cancellazione numerica (dato che nel secondo sistema i
coefficienti non sono naturali come nel primo). Ma la cosa, come
accennava Max, dovrebbe essere abbastanza improbabile,

Post by pan
Inoltre il
fatto che c1 sia nullo se b e d sono nulli, mi fa sospettare che ci
sia qualcosa in piu' rispetto ad un sistema 4x4...

No, è che ho copiato la prima matrice invertendo le prime due righe
(l'avrò presa da altri appunti). E mi sono appena accorto di un
errore: la seconda riga della prima matrice dovrebbe essere "0 f 0 h"
e non "0 f g 0". Invertendo le prime due righe si ha:

0 b 0 d
e 0 g 0
i j k l
m n o p

da cui si vede subito che b=d=0 => c1=0

--
Narcolessico
http://binaryunit.blogspot.com

[pan]

In article

Post by Eugenio
No,

Bene cosi', non voglio ricominciare a svegliarmi nel cuore della notte
urlando "PLUCKER!!!".

Post by Eugenio
è che ho copiato la prima matrice invertendo le prime due righe
(l'avrò presa da altri appunti). E mi sono appena accorto di un
errore: la seconda riga della prima matrice dovrebbe essere "0 f 0
0 b 0 d
e 0 g 0
i j k l
m n o p
da cui si vede subito che b=d=0 => c1=0

Ok, mi mancava il fatto che c1 fosse la soluzioe del primo sistema.
Ma quindi son sistemi omogenei?

--*PaN!*

--
I'm using an evaluation license of nemo since 219 days.
You should really try it!
http://www.malcom-mac.com/nemo

[Eugenio]

Post by pan
Ok, mi mancava il fatto che c1 fosse la soluzioe del primo sistema.

Scusa, sono stato impreciso in molti passaggi.

Post by pan
Ma quindi son sistemi omogenei?

No... Volevo vedere quando il denominatore della frazione si
azzerasse, e uno dei casi in cui si azzera è quando b e d sono nulli.
Ma se b e d sono nulli, allora anche c1 lo è (perchè la prima
equzione del sistema è b*y+d*w = c1, con y e w incognite); ovvero, la
matrice ha rango inferiore a 4 e ci sono infinite soluzioni. Ma
parlavo dell'ipotesi in cui b=d=0, non di una cosa che avviene sempre.

L'ultimo passaggio era: so per certo che i valori della matrice
producono sempre un sistema determinato, con una ed una sola
soluzione; per cui b=d=0 (o, più in genere, una frazione indeterminata
o una divisione per zero) non dovrebbero avvenire mai, se non per
cancellazione numerica (ovvero per errori di precisione).

--
Narcolessico
http://binaryunit.blogspot.com

[pan]

In article

Post by Eugenio

Post by pan
Ok, mi mancava il fatto che c1 fosse la soluzioe del primo sistema.

Scusa, sono stato impreciso in molti passaggi.

Post by pan
Ma quindi son sistemi omogenei?

No...

Ah ecco, mi tornano piu' cose... altrimenti se avessero avuto sempre
una e una sola soluzione, questa sarebbe sempre stata nulla, ma a
questo punto mi pareva troppo strano.

--
Marco

--
I'm using an evaluation license of nemo since 220 days.
You should really try it!
http://www.malcom-mac.com/nemo

[Manlio Perillo]

Il Sun, 23 Mar 2008 18:03:04 +0100, Jonas Willard (AKA \"Archeopteryx\")

Post by Eugenio
Consigli?

Mh, immagino che non ti aiuterà molto ma io avevo un integrale
orripilante che però sapevo essere calcolabile analiticamente. Ho fatto
digerire la formula al tool di calcolo simbolico di matlab e il
programma, ovviamente, va che è una scheggia rispetto alle prime prove
col calcolo numerico anche se sul sorgente la formula occupa sei o sette
righe :) Penso che potresti fare lo stesso, oltretutto, matlab
restituisce una stringa col nome delle variabili che gli hai dato tu, e
fare copincolla nell'editor del sorgente è cosa breve.

Attenzione che una soluzione analitica *non* è detto che sia l'ideale
quando poi la applichi a numeri in virgola mobile.

Andrebbe comunque controllata.

Però magari non
hai matlab a portata di mano, e non so se tra le cose che fa c'è la
soluzione dei sistemi lineari.

Per il calcolo simbolico c'è Maxima, che è Open Source.

Se invece lo vuoi fare numericamente, installati Atlas
http://math-atlas.sourceforge.net/
e usa questo esempio:
http://seehuhn.de/media/data/testlapack.c
preso da:
http://seehuhn.de/pages/linear#lapack-easy


Atlas implementa BLAS con in più alcune funzioni di Lapack, quindi non
devi fare cose strane per interfacciarti con la libreria in Fortran.


P.S.: incredibile come, in questi tutti anni, non si sia ancora pensato a
standardizzare LAPACK e altre API per il calcolo numerico per il C.

ciao!
Apx.

Saluti ed auguri di buona Pasqua
Manlio Perillo

[?manu*]

Post by Eugenio
Dovrei risolvere due sistemi lineari di 4 equazioni in 4 incognite,
di cui ho la matrice dei coefficienti (float) e il vettore dei termini
noti. Una matrice è del tipo
a 0 c 0
0 f g 0
i j k l
m n o p
e devo risolverla solo una volta. L'altra è del tipo
a b c 0
e f g 0
i j 0 l
m n 0 p
e devo risolverla decine di volte al secondo. Cosa potrei usare che
non mi complichi troppo la vita?
Ipotesi 2: usare il metodo di Cramer, cosa che però implica il calcolo
del determinante della matrice; e il determinante di una matrice 4x4 è
dispendioso da calcolare senza usare programmazione dinamica (a meno
che non debba implementare un metodo veloce... ma allora tanto vale
implementare una fattorizzazione LU con cui risolvere tutto il
sistema!)
Ipotesi 3: sviluppare a manina le formule (es. per sostituzione) e
scriverle direttamente...

Io utilizzerei la 2, (o se preferisci la 3, che è lo stesso). Infatti un
algoritmo generico ti creerebbe comunque un ciclo, che in questo caso,
penso, ti farebbe perdere tempo. Inoltre la presenza di alcuni zeri
potrebbe semplificarti un po' la formula.

E.

[HappyCactus]

Post by Eugenio
Dovrei risolvere due sistemi lineari di 4 equazioni in 4 incognite,
di cui ho la matrice dei coefficienti (float) e il vettore dei termini
noti. Una matrice è del tipo
a 0 c 0
0 f g 0
i j k l
m n o p
e devo risolverla solo una volta. L'altra è del tipo
a b c 0
e f g 0
i j 0 l
m n 0 p

In entrambi i casi, utilizzerei il metodo specifico, visto che queste
forme sono note e fisse, quindi andrei con il metodo 3. Di sicuro il
metodo è molto ottimizzabile, quindi ti garantisce il massimo delle
performance.
In alternativa utilizzerei la riduzione LU che viene fatta in due tempi;
se nel secondo sistema a cambiare sono i soli termini noti, si fa il
primo passo offline e il secondo passo online, risparmiando una metà del
tempo.