Calcoli e approssimazione dei dati in MATLAB

Risoluzione di sistemi di equazioni algebriche lineari e problemi spettrali

La risoluzione di sistemi di equazioni algebriche lineari è uno dei principali problemi computazionali, poiché ad esso si riduce un numero enorme di problemi che sorgono in un'ampia varietà di aree di applicazione. Metodi numerici utilizzati per la soluzione approssimativa di problemi meccanici, per il calcolo del flusso di liquidi e gas, altri processi fisici in mezzi continui, metodi per il calcolo di circuiti elettrici, metodi per l'approssimazione dei dati portano alla necessità di risolvere sistemi di equazioni algebriche lineari (e questo non è un elenco completo delle fonti di occorrenza dei sistemi algebrici lineari). I problemi spettrali sorgono, ad esempio, nell'analisi in frequenza delle strutture, nello studio della stabilità dei processi nei circuiti elettrici e nei sistemi di controllo e in molti altri settori applicativi.

Panoramica delle capacità di MATLAB per la risoluzione di sistemi di equazioni algebriche lineari

L'ambiente MATLAB è stato originariamente sviluppato per lavorare con le matrici (MATLAB è l'abbreviazione di Matrix Laboratory), quindi l'arsenale di strumenti MATLAB per risolvere sistemi di equazioni algebriche lineari è piuttosto ricco e comprende:

• sistemi risolutivi con matrici quadrate e rettangolari;
• sistemi di risoluzione mediante metodi diretti e iterativi (compresa la possibilità di precondizionamento);
• decomposizioni di matrice;
• memorizzazione di grandi matrici sparse in una forma compatta e algoritmi speciali per la risoluzione di sistemi con tali matrici.

Risolvere sistemi utilizzando le barre rovesciate

Il modo più semplice per risolvere i sistemi è utilizzare il segno barra rovesciata. Supponiamo di dover risolvere il sistema

Per fare ciò, compila la matrice e il vettore colonna sul lato destro (il lato destro deve essere una colonna, altrimenti verrà visualizzato un errore relativo alla mancata corrispondenza delle dimensioni)

e usa il segno della barra rovesciata

X = A\fx = 1 1 1

Invece del segno della barra rovesciata, potresti chiamare la funzione mldivide

X = mldividi(A, f)

Il risultato sarà lo stesso

Calcolando la discrepanza ci assicuriamo che la soluzione sia stata trovata correttamente (in generale, piccoli valori delle componenti della discrepanza non sempre indicano una soluzione trovata correttamente)

F - A*x ans = 0 0 0

È molto importante non confondere A e f, perché quando si esegue un'operazione

non ci sarà alcun errore, ma verrà visualizzato

X = 0,2707 0,3439 0,3854

che non ha nulla a che vedere con la soluzione del sistema in esame. Vediamo perché questo accade scrivendo un sistema con una “matrice” fe un “lato destro” A. In questo caso ci sarà una sola incognita, perché f dimensione 3 per 1:

Se c'è una matrice sul lato destro, verranno risolti tanti sistemi quante sono le colonne della matrice e ciascuna colonna è un vettore del lato destro del sistema corrispondente (ovvero, l'operatore barra rovesciata consente di risolvere diversi sistemi contemporaneamente). Quindi i sistemi sono risolti

A rigor di termini, nessuno di questi sistemi ha una soluzione. Tuttavia, se la matrice del sistema è rettangolare e il numero delle sue righe è maggiore del numero di colonne (cioè il numero di equazioni è maggiore del numero di incognite), allora tale sistema è chiamato sistemi sovradeterminati e la sua approssimazione la soluzione si cerca minimizzando la norma euclidea del residuo. Quindi per il primo sistema la soluzione sarà X, che fornisce un minimo alla seguente espressione

(4-7x 1) 2 +(3-11x 1) 2 +(2-12x 1) 2

È facile assicurarsi che il minimo fornisca solo . Usiamo, ad esempio, il Symbolic Math Toolbox

Simboli x R = (7*x-4)^2 + (11*x-3)^2 + (12*x-2)^2 dRdx = diff(R) x = risolvere(dRdx) x = 85/314 double(x) ans = 0,2707

Quindi il vettore ottenuto scrivendo x = f\A contiene la soluzione di tre sistemi sovradeterminati, il membro destro di ciascuno dei quali è la corrispondente colonna della matrice A.

Naturalmente la matrice del sistema sovrascritto non contiene necessariamente una sola colonna. Ad esempio, consideriamo il problema della selezione dei parametri UN E B modello lineare secondo i dati tabellari forniti

K	1	2	3	4	5
xk	1.0	1.5	2.0	2.5	3.0
sì k	2.99	2.81	2.89	3.03	3.21

Questo problema si riduce ad un sistema sovradeterminato rispetto a UN E B, se pretendiamo l’uguaglianza y(xk)=sì k Per K=1,2,...,5:

Componiamo la matrice del sistema e il vettore del lato destro

X = (1:0.5:3)" A = ; f = ; e risolvendolo c = A\f otteniamo c = 2.9903 2.0100 cioè il grafico dei dati e delle approssimazioni conferma la correttezza dei risultati ottenuti figura fun = @(x ) c(1)/x + c(2)*log(x); fplot(fun, ) tieni premuto plot(x, y, "o")

Se un sistema ha più equazioni che incognite, allora tale sistema è chiamato sistema sottodeterminato e può essere risolto in MATLAB utilizzando il segno barra rovesciata. Avrà infinite soluzioni e si trova una soluzione che contiene quanti più valori zero possibili. Consideriamo questo esempio

UN = ; f = ; x = A\f x = 1.3333 1.0000 0 1.6667

Oltre al carattere barra rovesciata viene utilizzato anche il carattere barra normale. Sono collegati dalla seguente regola

B/A è equivalente a (A"\B")"
dove l'apostrofo significa trasposizione. Invece della barra normale, puoi anche utilizzare la funzione mrdivide.

Per ora non approfondiremo gli algoritmi per risolvere i sistemi incorporati nell'operazione backslash. Questa sezione è dedicata a questo. In particolare, la scomposizione LU viene utilizzata per risolvere sistemi con una matrice quadrata generale. Inoltre, viene verificata la condizionalità della matrice. Un classico esempio di matrice mal condizionata è la matrice di Hilbert, i cui elementi sono determinati dalla formula. Per creare una matrice di Hilbert di una determinata dimensione in MATLAB, utilizzare la funzione hilb. Risolviamo, ad esempio, un sistema del 13° ordine, il cui membro destro è tale che la sua soluzione deve essere composta da tutte le unità ( K L'esimo elemento del secondo membro è la somma degli elementi K-esima riga della matrice). A = impugnatura(13); f = somma(A, 2); x = A\f

Di conseguenza, riceviamo un messaggio di cattivo condizionamento

Attenzione: la matrice è vicina al singolare o mal scalata. I risultati potrebbero essere imprecisi. RCCOND = 2.339949e-018.

(RCOND è una stima per l'inverso del numero di condizione) e una soluzione errata

X = 1,0000 1,0000 1,0004 0,9929 1,0636 0,6552 2,2023 -1,7860 5,3352 -3,4773 3,9431 -0,1145 1,1851

I pericoli generati da un cattivo condizionamento della matrice sono descritti nella sezione. Passiamo ora ad una panoramica delle scomposizioni di matrice disponibili in MATLAB.

Scomposizioni di matrici di Cholesky, LU e QR

MATLAB dispone di funzioni per le seguenti espansioni:

• Espansione di Cholesky - funzione chol;
• Espansione LU - funzione lu;
• Scomposizioni QR - funzione QR.

Le espansioni aiutano a risolvere in modo efficiente un intero insieme di sistemi di equazioni algebriche lineari con la stessa matrice e diversi vettori a destra. Consideriamo il seguente insieme di sistemi:

Ascia=f(k), K=1, 2,...,N.

Supponiamo che la matrice UN si presenta come un prodotto di due matrici e sistemi con matrici A=BC e sistemi con matrici B E C risolto molto più velocemente che con una matrice UN. Quindi le soluzioni dei sistemi originali si ottengono come segue. Perché il Ascia=f(k) E A=BC, Quello BCx=f (K) e risolvere in sequenza Per=f(k) E Cx=y troviamo una soluzione K sistema. Per le matrici piene, la scomposizione viene eseguita come segue: O(N 3) operazioni dove N- la dimensione della matrice (cioè contemporaneamente alla risoluzione del sistema) e la risoluzione del sistema con ciascun fattore di espansione in O(N 2) operazioni. Quindi la soluzione N sistemi c con tempi di espansione preliminari O(N 3)+O(N 2)N. Sarebbe necessario risolvere ciascun sistema senza decomposizione della matrice O(N 3)N operazioni aritmetiche.

Decomposizione di Cholesky per una data matrice UN consiste nel trovare tale matrice triangolare superiore R con elementi diagonali positivi, che A=R T R. È noto che se la matrice UN simmetrico e definito positivo (cioè per qualsiasi vettore Xè vero: x T Ax0, oppure, che è lo stesso, tutti gli autovalori della matrice sono positivi), allora la decomposizione di Cholesky esiste ed è unica.

Consideriamo come esempio la scomposizione di Cholesky della matrice

UN = ; R = col(A) R = 2,0000 0,5000 0,5000 0 1,9365 0,3873 0 0 1,8974

Un semplice controllo garantisce che la scomposizione venga eseguita correttamente (entro la precisione degli errori derivanti durante le operazioni sui numeri reali)

A - R"*R ans = 1.0e-015 * 0 0 0 0 0 0 0 0 0.4441

Se la matrice non è definita positiva, ciò risulterà chiaro durante il processo di decomposizione, perché sarà necessario estrarre la radice di un numero negativo e verrà visualizzato un avviso corrispondente:

UN = ; R = col(A)??? Errore nell'utilizzo di ==> chol La matrice deve essere definita positiva.

La funzione chol non controlla la simmetria della matrice. I calcoli utilizzano elementi della matrice originale situati sulla diagonale e sopra:

UN = ; R = chol(A) R = 2.0000 0.5000 1.0000 0 1.9365 0.2582 0 0 1.7127 A - R"*R ans = 0 0 0 0 0 0 -1 0 0

Per una matrice quadrata arbitraria si può fare L.U.-decomposizione, cioè trovare la matrice triangolare inferiore l e matrice triangolare superiore U tale che A=LU. Un'affermazione più rigorosa è formulata come segue.

Se Ak- principale minore di una matrice quadrata UN misurare N, compilato dal primo K righe e colonne (ad esempio A(1:k, 1:k)) e det( A k)0 per K=1,2,..., N-1, allora esiste un'unica matrice triangolare inferiore l, la cui diagonale è composta da unità, e l'unica matrice triangolare superiore U tale che A=LU.

Durante il calcolo L.U. Un'espansione LU potrebbe richiedere la permutazione delle righe per garantire la stabilità numerica del processo di espansione, quindi la funzione lu può restituire una matrice l, che è triangolare inferiore fino alla permutazione delle righe, ad esempio:

UN = ; = lu(A) L = 0 1.0000 0 1.0000 0 0 0.5000 0.5000 1.0000 U = 2 3 1 0 1 1 0 0 4

Lavoro L.U. Dove l- triangolare fino alla permutazione delle righe e sarà uguale alla matrice UN(in generale, tenendo conto dell'errore nelle operazioni con numeri reali).

A-L*U ans = 0 0 0 0 0 0 0 0 0

Se per la matrice UN Se viene eseguita la scomposizione di Cholesky o la scomposizione LU, allora la soluzione del sistema con la matrice UN, come accennato in precedenza, si riduce a risolvere due sistemi con matrici triangolari. Questa soluzione può essere realizzata utilizzando l'operazione barra rovesciata, poiché l'algoritmo dietro di essa identifica le matrici triangolari e applica un modo efficiente per risolvere il sistema con esse, richiedendo O(N 2) operazioni aritmetiche. Consideriamo un esempio di risoluzione del sistema

con preliminare L.U.-decomposizione della matrice.

UN = ; f = ;

Eseguiamo L.U.-decomposizione

Lu(A);

e risolvere in sequenza due sistemi con matrici triangolari, prima con L, poi con U

Y = L\f; x = U\y x = 1 1 1

La soluzione di due sistemi può essere scritta in un'unica espressione

X = U\(L\f)

e il risultato sarà lo stesso. Occorre prestare attenzione all'importanza dell'uso delle parentesi per determinare l'ordine di soluzione dei sistemi con matrici triangolari. Espressione senza parentesi

X = U\L\f

porterà a un risultato completamente diverso, poiché viene eseguito per primo U\L, che equivale (come discusso sopra) a risolvere sistemi con matrice U e colonne L come vettori del lato destro. Il risultato è una matrice, ciascuna colonna della quale è una soluzione del sistema corrispondente, quindi con questa matrice e il vettore f di destra il sistema viene risolto. Questo processo ovviamente non ha nulla a che vedere con la soluzione del sistema originario.

Consideriamo l'ultima scomposizione della matrice - QR-decomposizione eseguita dalla funzione qr. QR-la scomposizione può essere fatta per matrici rettangolari, precisamente se UN matrice dimensionale M SU N, allora esiste una matrice ortogonale Q misurare M SU M(cioè tale che Q -1 =Q T) e la matrice R misurare M SU N con tutti gli elementi zero sotto la diagonale principale, che A=QR.

QR-l'espansione ha una buona stabilità computazionale e viene utilizzata, in particolare, quando si approssimano i dati utilizzando il metodo dei minimi quadrati. Ad esempio, considera di risolvere il sistema precedente utilizzando la scomposizione QR:

Qr(A) Q = -0,7428 0,6000 -0,2971 -0,5571 -0,8000 -0,2228 -0,3714 0,0000 0,9285 R = -5,3852 -5,3852 -5,0138 0 -5,0000 0,4000 0 0 6,6108

Una volta fatto QR-decomposizione, soluzione del sistema Ascia=f può essere trovato come soluzione al sistema di matrice triangolare R, Perché perché il QR=f, Quello Rx=Q T f:

>> x = R\(Q"*f) x = 1.0000 1.0000 1.0000

Maggiori dettagli sulle scomposizioni di matrici, le corrispondenti funzioni MATLAB chol, lu e qr e i loro argomenti sono scritti nella sezione Scomposizioni di matrici.

Algoritmo per la risoluzione di un sistema di equazioni lineari utilizzando il segno della barra rovesciata

Quando si risolve un sistema di equazioni algebriche lineari in MATLAB utilizzando il segno della barra rovesciata
x = A\b
funziona un algoritmo che, a seconda del tipo di matrice, risolve il sistema utilizzando il metodo più adatto, implementato in una delle procedure del pacchetto LAPACK o UMFPACK. Qui b può anche essere una matrice, il cui numero di righe coincide con il numero di righe della matrice A. Quindi viene restituita la matrice x, ciascuna i-esima colonna contiene la soluzione del sistema Ax (i) = b (i ) , i=1,2,... ,k, dove b (i) è la i-esima colonna della matrice b = [ b (1) | b(2) ...| b(k)].

L'algoritmo implementato nell'operazione \ consiste dei seguenti passaggi:

• 1. Nel caso banale, se A è sparsa e diagonale (la sezione Matrici sparse è dedicata alle matrici sparse in MATLAB), la soluzione si trova utilizzando la semplice formula x k = b k /a kk , dove k=1,2,. ..N.

2. Se A è una matrice quadrata, sparsa e a bande, viene utilizzato un risolutore di matrici a bande. Potrebbero esserci due opzioni:

UN. Se A è una matrice tridiagonale e b è una singola colonna di numeri reali, allora il sistema viene risolto mediante eliminazione gaussiana (le sue operazioni vengono eseguite solo sugli elementi sulle diagonali). Se durante il processo di eliminazione è necessario eseguire permutazioni di riga per mantenere la stabilità, o se la matrice A non è tridiagonale, allora funziona il punto seguente.
B. Se A è un nastro con una densità di elementi diversa da zero maggiore del parametro fasciato, per default pari a 0,5, oppure non sono soddisfatte le condizioni del paragrafo precedente, quindi, a seconda della tipologia A e b ( Doppio O separare) vengono chiamate le seguenti procedure della libreria LAPACK:

A e b sono veri doppi: procedure DGBTRF, DGBTRS
Doppie di tipo complesso Aeb - procedure ZGBTRF, ZGBTRS
A o b tipo reale singolo - procedure SGBTRF, SGBTRS
Aob tipo complesso singolo - procedure CGBTRF, CGBTRS

La densità degli elementi diversi da zero è intesa come il rapporto tra il numero di elementi diversi da zero nel nastro e il numero di tutti gli elementi nel nastro della matrice, ad esempio, per una matrice 7 per 7, il cui modello è dato sotto, il numero di elementi diversi da zero

nz = 1 (sullo schema N.-2) + 6 (sullo schema N.-1) + 7 (sullo schema N. 0) + 6 (sullo schema N. 1) + 1 (sullo schema N. 2) + 1 (sul diagramma n. 3) = 22,

e il numero di tutti gli elementi nel nastro

fascia = 5 (sullo schema n.-2) + 6 (sullo schema n.-1) + 7 (sullo schema n. 0) + 6 (sullo schema n. 1) + 5 (sullo schema n. 2) + 4 (sul diagramma n. 3) = 33

e la densità degli elementi diversi da zero sarà 2/3. Poiché 2/3 > 0.5 verrà utilizzato un solutore per matrici a strisce
Parametro fasciato, come altri parametri che controllano gli algoritmi di matrice sparsa in MATLAB, viene impostato utilizzando la funzione spparms.



Il significato di questo passaggio dell'algoritmo è che la risoluzione di un sistema con una matrice a nastro non richiede azioni al di fuori del nastro. Se il feed è abbastanza pieno, non ci saranno molte azioni con zero elementi. Al contrario, se il nastro è piuttosto scarso, gli approcci forniti nei passaggi successivi dell'algoritmo possono avere un effetto maggiore.

• 3. Se A è una matrice triangolare superiore o inferiore, viene utilizzato il metodo di sostituzione inversa o, di conseguenza, il metodo di sostituzione diretta, in cui il componente della soluzione viene trovato dall'ultima (o prima equazione) e quindi i componenti trovati sono sostituito nelle seguenti equazioni per trovare i successivi sistemi di componenti della soluzione di equazioni.
• 4. Se A - può essere ridotto a una forma triangolare mediante permutazioni, viene eseguita la riduzione corrispondente e quindi il sistema di equazioni viene risolto come nel passaggio 3.

5. Se A è una matrice simmetrica (o, nel caso complesso, hermitiana), e la sua diagonale contiene solo elementi reali positivi, allora a seconda che A sia sparso o meno, vale il punto a) o il punto b).

UN. Se A non è sparso, si tenta di eseguire la decomposizione di Cholesky A = R T R seguita dalla risoluzione di sistemi con matrici triangolari R T e R: R T y = b e Rx = y. In questo caso vengono richiamate le seguenti procedure della libreria LAPACK:

A e b sono reali e doppi: DLANGE, DPOTRF, DPOTRS, DPOCON
A e b sono complessi e doppi: ZLANGE, ZPOTRF, ZPOTRS, ZPOCON
A e b sono reali e di tipo singolo: SLANGE, SPOTRF, SPOTRS, SDPOCON
A e b tipo complesso e singolo: CLANGE, CPOTRF, CPOTRS, CPOCON

La decomposizione di Cholesky fallirà se la matrice non è definita positiva. Ma non esiste un controllo preliminare per la definitezza positiva, ed è determinata proprio durante il processo di decomposizione, poiché molto spesso è sufficiente eseguire diversi passaggi della decomposizione di Cholesky per chiarire il fatto che la matrice non è definita positiva (nel processo di calcolo diventa necessario estrarre la radice quadrata di un numero negativo). Se la decomposizione di Cholesky non è soddisfatta si passa al punto successivo.

B. Se A è sparso, le permutazioni simmetriche di righe e colonne vengono prima eseguite utilizzando l'algoritmo di grado minimo simmetrico (funzione symmmd) per ridurre il riempimento del fattore di decomposizione di Cholesky, ovvero ridurre il numero di nuovi elementi diversi da zero che si verificano durante il processo di riempimento:

p = symmmd(A) - il vettore p contiene permutazioni

R = chol(A(p, p));

e vengono risolti due sistemi con moltiplicatori di Cholesky, il primo con un moltiplicatore trasposto e corrispondenti permutazioni nel vettore a destra

e la seconda, con un fattore di espansione e con l'inserimento delle componenti della soluzione nelle corrispondenti posizioni del vettore x
x(p, :) = R\y

• 6. Se A è una matrice di Hessenberg, allora viene ridotta a una matrice triangolare superiore (con opportuna modifica del membro di destra) e quindi il sistema viene risolto mediante sostituzioni

7. Se A è una matrice quadrata che non soddisfa le condizioni dei paragrafi 1-6, a seconda che sia sparsa o meno, vengono eseguite le seguenti azioni

UN. Se A non è una matrice sparsa, allora utilizzando l'eliminazione gaussiana con la scelta di un elemento guida (per garantire la stabilità della scomposizione), viene eseguita una scomposizione LU della matrice A = LU, dove

U - matrice triangolare superiore
L è una matrice ridotta a triangolare per permutazioni di righe

e la soluzione del sistema Ax = b si trova risolvendo sequenzialmente i sistemi con matrici triangolari Ly = b, Ux = y.
Per eseguire la scomposizione della LU, vengono chiamate le seguenti procedure della libreria LAPACK:

A e b sono reali e doppi: DLANGE, DGESV, DGECON
A e b sono complessi e doppi: ZLANGE, ZGESV, ZGECON
A e b sono reali e di tipo singolo: SLANGE, SGESV, SGECON
A e b sono complessi e di tipo singolo: CLANGE, CGESV, CGECON

B. Se A è una matrice sparsa, le colonne vengono riorganizzate per ridurre il riempimento dei fattori L e U nel processo di ricerca della scomposizione. Inoltre, riorganizzando le righe nel processo di scomposizione LU, viene garantita la stabilità computazionale, dopodiché i sistemi con matrici triangolari vengono nuovamente risolti. Per eseguire la scomposizione LU di una matrice sparsa vengono utilizzate le procedure della libreria UMFPACK

8. Se i punti precedenti dell'algoritmo non hanno funzionato e, quindi, A non è una matrice quadrata, allora tutto è determinato dalla sua scarsità e uno dei punti seguenti funziona:

UN. Se A non è una matrice sparsa, viene eseguita la scomposizione QR AP = QR, dove P è la matrice di permutazione delle colonne, Q è la matrice ortogonale (Q T Q = I) e R è triangolare superiore. Se A è la dimensione m per n, allora Q è la dimensione m per m e R è la dimensione m per n. Successivamente, la soluzione del sistema si trova come segue:

x = P*(R \ (Q" * b))

poiché da Ax = b e AP = QR segue: QRx = bP e Rx = Q T bP.

B. Nel caso di una matrice A sparsa e rettangolare, non ha senso calcolare la matrice Q, poiché molto probabilmente sarà riempita. Pertanto, l'algoritmo di decomposizione QR calcola c = Q T b (cioè il lato destro modificato) e risolve il sistema Rx = c con una matrice triangolare per ottenere una soluzione del sistema originale Ax = b.

In aggiunta a tutto quanto sopra, l'algoritmo di cui sopra valuta la condizionalità della matrice e visualizza un avviso su un possibile errore nella soluzione se la matrice è scarsamente condizionata (vedere la sezione L'influenza della condizionalità di una matrice sull'accuratezza di risolvere un sistema con esso). L'algoritmo sopra riportato contiene vari controlli che potrebbero richiedere molto tempo aggiuntivo. La sezione successiva, La funzione linsolve per la risoluzione di sistemi di equazioni lineari, introduce la funzione linsolve, che consente di specificare il tipo di matrice, riducendo così il tempo di calcolo.

funzione linsolve per risolvere sistemi di equazioni lineari

Invece di risolvere un sistema (o sistemi con più membri destri definiti in una matrice) di equazioni algebriche lineari utilizzando il segno della barra rovesciata, è possibile utilizzare la funzione linsolve, che non esegue tutti i controlli della matrice inerenti all'algoritmo dell'operazione \ ( vedere la sezione precedente).

La funzione linsolve richiamata nella sua forma più semplice con due argomenti di input e un argomento di output
x = linsolve(A, b)
risolve il sistema Ax = b in uno dei modi, a seconda che la matrice sia quadrata o meno.

• Se A è una matrice quadrata, viene prima calcolata la sua scomposizione LU e poi vengono risolti due sistemi con matrici triangolari L e U.
• Se A è una matrice rettangolare, viene prima calcolata la sua scomposizione QR e poi viene risolto il sistema con matrici triangolari.

(Maggiori dettagli sulle azioni con i risultanti fattori di espansione sono scritti nella sezione). Inoltre, se la matrice A è mal condizionata o A è una matrice di rango incompleto, viene visualizzato un avviso corrispondente, ad esempio:

UN = ; b = ; x = linsolve(A,b) Avvertenza: rango carente, rango = 1, tol = 4.4686e-015. x = 0 0 2.0000

Per evitare che questo messaggio venga stampato nella finestra di comando, chiamare la funzione linsolve con due argomenti di output
= linsolve(A, b)
quindi la soluzione risultante verrà scritta in x e in r - il rango della matrice (se A è una matrice rettangolare) o il reciproco della stima per il suo numero di condizione (vedi sezione), ad esempio

UN = ; b = ; = linsolve(A,b) x = -1.0000e+000 9.9999e-001 3.3307e-006 r = 6.9444e-013

I principali vantaggi della funzione linsolve rispetto all'operazione \ si manifestano quando si specifica il tipo di matrice del sistema di equazioni. Per fare ciò, prima di chiamare la funzione linsolve, è necessario compilare una struttura con informazioni sulla matrice con i seguenti campi

• SYM - simmetrico;
• LT - triangolare inferiore;
• UT - triangolare superiore;
• UHESS - Hessenberg;
• POSDEF - simmetrico, definito positivo;
• RETTO - rettangolare;
• TRANSA - se è necessario risolvere un sistema con una matrice trasposta ad una data.

Ogni campo può essere vero o falso. Naturalmente non tutte le combinazioni dei valori dei campi sono valide; ad esempio, una matrice non può essere triangolare e allo stesso tempo simmetrica e definita positiva. Le corrette combinazioni dei valori dei campi sono raccolte nella tabella seguente

LT	UT	UHESS	SIMBOLO	POSDEF	RETT	TRANSAZIONE
VERO	falso	falso	falso	falso	vero falso	vero falso
falso	VERO	falso	falso	falso	vero falso	vero falso
falso	falso	VERO	falso	falso	falso	vero falso
falso	falso	falso	VERO	VERO	falso	vero falso
falso	falso	falso	falso	falso	vero falso	vero falso

Se la matrice del sistema è definita positiva, è necessario tenerne conto durante la risoluzione, poiché per le matrici definite positive la soluzione si basa sulla decomposizione di Cholesky, che richiede meno operazioni rispetto alla decomposizione LU utilizzata quando si risolvono sistemi con matrici quadrate generali . Questo è facile da verificare con l'aiuto del seguente esempio, in cui viene creata una matrice definita positiva simmetrica (una matrice di numeri casuali viene aggiunta con una trasposizione ad essa e grandi numeri vengono aggiunti alla diagonale) e un sistema di Le equazioni con questa matrice vengono prima risolte come un sistema con una matrice di forma generale (opts.SYM = false e opts.POSDEF = false), e poi come con una matrice simmetrica e definita positiva (opts.SYM = true e opts .POSDEF = vero).

% imposta tutti i campi della struttura ops, eccetto SYM e POSDEF opts.TRANSA = false; opts.UHESS = falso; opts.RECT ​​= false; opts.UT = falso; opts.LT = falso; % crea un vettore di dimensioni della matrice N = 2.^(8:12); % crea array vuoti per registrare i tempi di soluzione TimeGeneral = ; TempoPosSym = ; % in un ciclo creiamo matrici e confrontiamo i tempi di soluzione per n = N % creiamo una matrice simmetrica e definita positiva % e il vettore del membro destro A = rand(n); A = A + A" + 2*n*eye(n); b = sum(A, 2); % risolve il sistema come un sistema con una matrice generale opts.SYM = false; opts.POSDEF = false; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimeGeneral = ; % risolve il sistema come un sistema con una matrice di simmetria e posizione opts.SYM = true; opts.POSDEF = true; Tstart = cputime; x = linsolve( A,b, opts); Tend = cputime; TimePosSym = ; end % visualizza i grafici temporali figura loglog(N, TimeGeneral, N, TimePosSym) legend("TimeGeneral", "TimePosSym")

I costi computazionali risultanti da questi metodi di risoluzione del sistema sono mostrati nel grafico sottostante

Naturalmente, se la matrice è triangolare, allora è molto importante indicarlo, poiché la risoluzione di un sistema con una matrice triangolare viene eseguita in operazioni O(n 2) e se un algoritmo di soluzione progettato per una matrice di si applichi una forma generale ad un sistema a matrice triangolare, allora occorreranno circa O(n 3) operazioni.

La funzione linsolve non controlla se la matrice soddisfa le proprietà specificate, il che potrebbe causare un errore. Ad esempio, se, quando si risolve un sistema con la matrice seguente e il membro di destra

UN = ; b = ; imposta true per il campo UT (e imposta tutti gli altri su false) opts.UT = true; opts.TRANSA = false; opts.LT = falso; opts.UHESS = falso; opts.SYM = falso; opts.POSDEF = falso; opts.RECT ​​= false; quindi, quando risolverà il sistema, la funzione linsolve considererà la matrice come una matrice triangolare superiore e selezionerà gli elementi di cui ha bisogno dal triangolo superiore A x = linsolve(A,b, opts) Ciò risulterà in una soluzione x = 1 1 1 corrispondente alla matrice A = ;

L'influenza della condizionalità della matrice sull'accuratezza della risoluzione di un sistema con essa

In questa sezione vedremo come gli errori negli elementi del vettore a destra e nella matrice di un sistema di equazioni lineari Ax = b possono influenzare la soluzione di quel sistema. Passiamo prima al caso di introduzione di disturbi nel vettore del lato destro b. Quindi risolviamo due sistemi

Inoltre, si presuppone di risolvere esattamente ciascuno dei sistemi, e la domanda principale che sorge è quanto sarà diversa la soluzione x del sistema (1) dalla soluzione del sistema (2) con una perturbazione sul lato destro δb. Questa è una questione piuttosto importante, poiché gli elementi sulla parte destra possono essere ottenuti come risultato di alcune misurazioni, contenenti rispettivamente l'errore δb k in ciascuna componente b k. Vorrei che misurando la stessa quantità (ogni volta con i suoi piccoli errori), anche le soluzioni corrispondenti di sistemi con lati destri leggermente diversi non differiscano molto l'una dall'altra. Sfortunatamente, non è sempre così. La ragione di ciò è una caratteristica della matrice chiamata its condizionalità. Questo sarà discusso ulteriormente.

Innanzitutto, è necessario introdurre una misura della prossimità dei vettori e x, ad es. vettore di errore. La misura della grandezza di un vettore è la norma (può essere definita in diversi modi). Per ora, prendiamo come norma la solita norma euclidea di un vettore (la radice quadrata della somma dei quadrati delle sue componenti), cioè

Di conseguenza

dove n è il numero di incognite ed equazioni nel sistema. Oltre alla norma vettoriale per stimare la grandezza del vettore, abbiamo bisogno anche di una norma matriciale per stimare la grandezza della matrice. Prendiamo la norma della matrice, che è definita come segue (è chiamata norma spettrale):

quelli. la norma della matrice spettrale è uguale alla radice quadrata dell'autovalore massimo della matrice A T A. MATLAB ha una funzione norma per il calcolo delle norme di matrici e vettori, che, in particolare, può calcolare le norme di cui sopra. Il motivo per cui abbiamo scelto queste particolari norme di vettori e matrici è spiegato in dettaglio nella sezione, dove vengono forniti alcuni calcoli e definizioni. Questo è legato alla stima che utilizzeremo per l'errore nella risoluzione del sistema (la derivazione di questa disuguaglianza è data anche nella sezione citata):

Qui indichiamo il numero di condizione della matrice, che è definito come segue:

Dalla disuguaglianza di cui sopra ne consegue che maggiore è il numero di condizione della matrice di sistema, maggiore può essere l'errore relativo nella soluzione con un piccolo errore sul lato destro.

Consideriamo un classico esempio di matrice mal condizionata - la matrice di Hilbert - e scopriamo come l'errore sul lato destro del sistema influenza l'errore nella soluzione. La matrice di Hilbert è definita come segue

Per creare una matrice di Hilbert, MATLAB fornisce la funzione hilb, il cui argomento di input specifica la dimensione della matrice. Prendiamo una piccola matrice 6 per 6:

N = 6; H = hilb(n) H = 1,0000 0,5000 0,3333 0,2500 0,2000 0,1667 0,5000 0,3333 0,2500 0,2000 0,1667 0,1429 0,3333 0,2500 0,2000 0,1667 0 0,1429 0,1250 0,2500 0,2000 0,1667 0,1429 0,1250 0,1111 0,2000 0,1667 0,1429 0,1250 0,1111 0,1000 0,1667 0,1429 0,1250 0,1111 0,1000 0,0909

X = unità(n, 1); b = H*x b = 2,4500 1,5929 1,2179 0,9956 0,8456 0,7365

Vediamo che non c'è nulla di “sospetto” né nella matrice né nel vettore sul lato destro; tutti i numeri non sono molto diversi tra loro. Formiamo ora il membro destro perturbato b + δb sommando al vettore b piccoli numeri dell'ordine di 10 -5 e risolviamo il sistema con il membro destro perturbato per ottenere il vettore .

Delta_b = 1e-5*(1:n)"; x_tilda = H\(b + delta_b) x_tilda = 0,9978 1,0735 0,4288 2,6632 -1,0160 1,8593

Si può vedere che la soluzione risultante è lontana da quella esatta in cui dovrebbero trovarsi tutte le unità. Calcoliamo l'errore relativo nella soluzione e sul lato destro (la funzione norma calcola per impostazione predefinita la norma euclidea del vettore):

Delta_x = x - x_tilda; SINISTRA = norm(delta_x)/norm(x) SINISTRA = 1.1475 DESTRA = norm(delta_b)/norm(b) DESTRA = 2.7231e-005

Quindi, l'errore nella soluzione è dell'ordine dell'unità, sebbene i cambiamenti sul lato destro fossero dell'ordine di 10 -5. Ciò si adatta perfettamente alla disuguaglianza di errore di cui sopra. Infatti, calcoliamo il numero di condizione cond(H) utilizzando la funzione MATLAB chiamata cond e per impostazione predefinita calcola per la norma spettrale della matrice

C = cond(H) c = 1.4951e+007

SINISTRA ≈ 1.1475 meno

C* GIUSTO? 1.4951e+07 * 2.7231e-05 ≈ 407

e la disuguaglianza è soddisfatta (anche con un certo margine).

All'aumentare della dimensione della matrice di Hilbert, l'errore nella soluzione non farà altro che aumentare (questo è facile da verificare ponendo n = 7, 8, ...). Inoltre per n = 12 verrà visualizzato un messaggio che informa che la matrice è mal condizionata e la soluzione potrebbe essere errata

Attenzione: la matrice è vicina al singolare o mal scalata.
I risultati potrebbero essere imprecisi. RCOND = 2.409320e-017.

Come misura di condizionalità, il valore scelto qui è RCOND uguale a uno diviso per la stima del numero di condizione (il numero di condizione viene stimato utilizzando un algoritmo abbastanza veloce che funziona molto più velocemente di cond, che è descritto più dettagliatamente nella guida sulla seconda funzione). Se il valore di RCOND è piccolo, la matrice è considerata mal condizionata.

L'aumento dell'errore nella soluzione all'aumentare della dimensione della matrice di Hilbert è dovuto al fatto che il numero di condizione della matrice di Hilbert cresce molto rapidamente con la sua dimensione. Questo è facile da verificare utilizzando un semplice ciclo e la funzione semilogica (la scala lungo l'asse delle ordinate è logaritmica):

N = 1:20; C = zeri(1, 20); per n = N H = hilb(n); C(n) = cond(H); fine semilogia(N, C) griglia su, titolo("cond(H)"), xlabel("n")

A proposito, poiché la funzione cond trova il numero di condizione utilizzando un metodo numerico (vale a dire, espansione singolare per trovare numeri singolari), il numero di condizione dopo n = 12 non viene più calcolato correttamente; infatti, dovrebbe crescere ulteriormente, il che può essere verificati mediante calcoli simbolici in MATLAB e operazioni con un dato numero di cifre significative

N = 1:20; C = zeri(1, N); cifre(60) per n = N H = vpa(sym(hilb(n))); Calcolo % della matrice di Hilbert accurato alla 60a cifra sigma = svd(H); % trovare i valori singolari della matrice di Hilbert % calcolare il numero di condizione della matrice di Hilbert C(n) = max(double(sigma))/min(double(sigma)); fine semilogia(N, C) griglia su, titolo("cond(H)"), xlabel("n")

Consideriamo ora tre punti importanti riguardo a questo esempio.

Il primo di essi - risolvere il sistema Hx = b (con il vettore di destra corrispondente alla soluzione in unità) utilizzando l'operatore barra rovesciata non fornirà unità esatte, l'errore sarà già nella decima cifra (sebbene in MATLAB tutti i calcoli vengono eseguiti con doppia precisione per impostazione predefinita)

Formato lungo H\b ans = 0,999999999999916 1,00000000002391 0,99999999983793 1,00000000042209 0,99999999953366 1,00000000018389

Ciò accade perché nel calcolo del vettore b moltiplicando la matrice H per un vettore di tutti uno, abbiamo già incluso qualche errore in esso. Inoltre, anche gli errori di arrotondamento nel processo di risoluzione del sistema hanno avuto un ruolo e la scarsa condizionalità della matrice (anche di piccole dimensioni) ha portato a tali errori nella soluzione. Infatti, per matrici piccole con un numero di condizione piccolo, tale effetto non sarà osservato. Prendiamo una matrice 6 per 6 di numeri casuali e calcoliamo il suo numero di condizione

A = Rand(n); cond(A) ans = 57,35245279907571

Quindi formiamo il lato destro corrispondente alla soluzione esatta di tutte le unità

X = unità(n, 1); b = A*x;

e risolvere il sistema, ottenendo una buona precisione

>> A\b ans = 1.00000000000000 1.0000000000000000 1.000000000000000 1.00000000000000 1.0000000000000000 1.00000000000000

Il secondo punto importante è legato al determinante della matrice. Calcoliamo il determinante della matrice di Hilbert 6 x 6 con cui abbiamo risolto il sistema

>> det(hilb(6)) ans = 5.3673e-018

Il determinante risulta essere sufficientemente piccolo, da cui si può trarre la conclusione errata che la fonte del problema di un grande errore nella soluzione con una piccola perturbazione del lato destro del sistema è la piccolezza del determinante. In realtà non è così; se il determinante del sistema è piccolo, ciò non significa che la matrice sia mal condizionata e, quindi, ci sia un grosso errore nella soluzione. Infatti, prendiamo la matrice

UN = ;

Il suo determinante è molto piccolo, circa 10 -121

Det(A) ans = 9,9970e-121

(la parola “molto” è ovviamente condizionale, ma è più piccola del determinante di una matrice di Hilbert 6 x 6, con cui abbiamo avuto problemi durante la risoluzione del sistema). Tuttavia, la piccolezza del determinante non influenzerà in alcun modo l'errore nella soluzione quando è perturbato il lato destro del sistema, cosa facile da dimostrare formando un sistema con una soluzione nota, introducendo disturbi nel lato destro e risolvendo il sistema:

X = ; b = A*x; delta_b = 1e-5*; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; SINISTRA = norma(delta_x)/norm(x) DESTRA = norma(delta_b)/norm(b) SINISTRA = 2.1272e-005 DESTRA = 2.1179e-005

Quindi, l'errore relativo nella soluzione è quasi uguale all'errore relativo sul lato destro, poiché il numero di condizione della matrice A sopra è leggermente maggiore di 1, vale a dire:

C = cond(A) c = 1.0303

E infine, consideriamo la terza domanda riguardante il raggiungimento del segno uguale nella disuguaglianza per un errore nella soluzione

Cioè, in altre parole, scopriamo: può esserci una situazione in cui facciamo una piccola perturbazione relativa del lato destro del sistema, diciamo 10 -5, il numero di condizione della matrice del sistema è uguale a 10 10 , e nella soluzione si ottiene un errore relativo di 10 5. Puoi fare molti esempi, provare matrici diverse, non raggiungere l'uguaglianza e concludere erroneamente che questa è solo una sovrastima dall'alto per un errore nella soluzione. Ma non è così, come ci convince il seguente esempio.

in cui la perturbazione relativa del lato destro è 10 -5

DESTRA = norma(delta_b)/norma(b) DESTRA = 1.0000e-005

L'errore relativo nella risoluzione del sistema risulta essere 10 5

X = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; SINISTRA = norma(delta_x)/norma(x) SINISTRA = 1.0000e+005

E questo accade perché
1) in questo esempio il numero di condizione della matrice A è 10 10;

C = cond(A) c = 1.0000e+010

2) la disuguaglianza per l'errore nella soluzione si trasforma in uguaglianza.

Se impostiamo il formato su long, vediamo che LEFT, RIGHT e il numero della condizione non sono esattamente 10 5 , 10 -5 e 10 10 rispettivamente, ma ciò è dovuto a errori di arrotondamento. Se risolta in aritmetica esatta, l'uguaglianza sarebbe esatta, cosa che può essere mostrata analiticamente (vedi sezione).

Negli esempi precedenti, abbiamo assunto che la matrice non contenga errori e che gli errori possano trovarsi solo sul lato destro del sistema di equazioni lineari. In pratica può succedere che gli elementi della matrice di sistema siano specificati con errori, ad es. invece del sistema

|| Ascia || V = || b || V⇒ || A || M || x || V≥ || b || V

dove || || M è una norma di matrice. Per fare ciò, la norma della matrice || || M e norma vettoriale || || V deve soddisfare la seguente disuguaglianza

|| A || M || x || V≥ || Ascia || V

per qualsiasi matrice A e vettore x. Se questa disuguaglianza vale, allora la norma della matrice || || Si chiama M concordato con norma vettoriale || || V. È noto, ad esempio, che la norma spettrale

(radice quadrata dell'autovalore massimo della matrice A T) è coerente con la norma vettoriale euclidea

Ecco perché in questa sezione abbiamo effettuato tutti gli esperimenti utilizzando queste norme.

Avendo ricevuto la disuguaglianza || A || M || x || V≥ || b || V notiamo inoltre che da Ax = b segue che . Poiché la matrice A non è singolare, ne consegue che δx = A -1 δb e || δx || V = || A -1δb || V. Ancora una volta usiamo la coerenza delle norme e otteniamo la disuguaglianza || δx || V ≤ || A -1 || M || δb || V. Inoltre, nelle due disuguaglianze ottenute

|| A || M || x || V≥ || b || V e || δx || V ≤ || A -1 || M || δb || V

dividere rispettivamente il valore minore di una delle disuguaglianze per il valore maggiore dell'altra disuguaglianza e il valore maggiore per quella minore

e con una semplice trasformazione otteniamo finalmente la disuguaglianza richiesta

dove cond(A) = || A || M* || A -1 || M.

Si noti che nel derivare la disuguaglianza per l'errore, abbiamo utilizzato solo il fatto che la norma della matrice è coerente con la norma del vettore. Ciò è vero non solo per la norma spettrale di una matrice e la norma euclidea di un vettore, ma anche per altre norme. Quindi, ad esempio, la norma massima della matrice della colonna, calcolata dalla formula

coerente con la prima norma vettoriale

e la norma massima della matrice di riga, calcolata dalla formula

coerente con la norma vettoriale infinita

La funzione norma calcola non solo la norma euclidea di un vettore e la norma spettrale di una matrice, ma anche le norme dei vettori e delle matrici sopra elencate. Per fare ciò, devi chiamarlo con un secondo argomento aggiuntivo:

• q = norm(A, 1) - norma massima della colonna della matrice A
• q = norm(A, inf) - norma massima di riga della matrice A
• q = norma(a, 1) - prima norma vettoriale di a
• q = norma(a, inf) - norma vettoriale infinita a

Numero di condizione della matrice cond(A) = || A || M* || A -1 || M rispetto a varie norme di matrice può essere calcolato utilizzando la funzione cond. Se cond viene chiamato con un argomento di input (una matrice), il numero di condizione viene calcolato rispetto alla norma della matrice spettrale. Chiamare la funzione cond con un argomento aggiuntivo restituisce i numeri di condizione relativi alla norma della matrice specificata:

• c = cond(A, 1) - numero di condizione relativo alla norma massima della colonna della matrice
• ñ = cond(A, inf) - numero di condizione relativo alla norma di riga massima della matrice

Come esempio che dimostra l'importanza della coerenza della norma della matrice con la norma del vettore nella disuguaglianza dell'errore, diamo un esempio con la matrice A, il vettore del lato destro b e il vettore dell'errore sul lato destro delta_b.

UN = ; b = [ -5.7373057243726833e-001 -1.5400413072907607e-001 -5.3347697688693385e-001 -6.0209490373259589e-001]; delta_b = [-0,71685839091451e-5 0,54786619630709e-5 0,37746931527138e-5 0,20850322383081e-5];

Calcoliamo le parti destra e sinistra della stima utilizzando la prima norma vettoriale e il numero di condizione della matrice in relazione alla norma della matrice spettrale

DESTRA = norma(delta_b,1)/norma(b,1); c = cond(A); x = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; SINISTRA = norma(delta_x,1)/norma(x,1); formato breve e disp()

Otteniamo i seguenti valori per i lati sinistro e destro della disuguaglianza
9.9484e+004 9.9323e+004

Il lato sinistro non dovrebbe superare il destro (come dimostrato sopra), ma si è rivelato più grande a causa della scelta di norme incoerenti.

Esaminiamo ora perché il numero di condizione di una matrice non può essere inferiore a uno. Il numero di condizione della matrice A è definito come cond(A) = || A || M* || A -1 || M , dove || A || M è una norma di matrice di A. La norma di matrice (cioè la regola secondo cui un numero è associato a ciascuna matrice) non può essere arbitraria; deve soddisfare i seguenti quattro assiomi:

A1. || A || ≥ 0 per qualsiasi matrice A e || A || = 0 se e solo se A è una matrice nulla.

A2. || αA || = | α | * || A || per qualsiasi matrice A e numero α.

A3. || A + B || ≤ || A || + || B || per qualsiasi matrice A e B

A4. || AB || ≤ || A || * || B || per qualsiasi matrice A e B.

Dall'ultimo assioma è chiaro che la norma è definita solo per matrici quadrate (sebbene nelle formule sopra riportate per il calcolo di varie norme, in linea di principio, non vi sia tale restrizione). Inoltre, dall'ultimo assioma segue che qualsiasi norma della matrice identità I non è inferiore a uno, anzi

|| Io || = || Io*Io || ≤ || Io || 2 ⇒ || Io || ≥ 1.

Quindi, sempre utilizzando il quarto assioma, troviamo che il numero di condizione della matrice è sempre maggiore di uno (vero per il numero di condizione della matrice rispetto ad una norma matriciale arbitraria)

1 ≤ || Io || = || AA-1 || ≤ || A || * || A -1 || = cond(A).

Completiamo questa sezione esaminando il motivo per cui appare l'esatta uguaglianza nella stima dell'errore nella soluzione quando è perturbato il lato destro:

e costruire esempi corrispondenti in MATLAB. (La discussione che segue è contenuta, ad esempio, nel libro di J. Forsythe, K. Mohler. Soluzione numerica di sistemi di equazioni algebriche lineari. M: “Mir”, 1969.)

Un ruolo essenziale nel ragionamento è giocato dal teorema sulla scomposizione singolare di una matrice, secondo il quale per ogni matrice reale di dimensione n esistono due matrici ortogonali U e V di dimensione n per n (U T U=UU T e V T V = VV T) tale che il prodotto D = U T AV sia una matrice diagonale, e possiamo scegliere U e V in modo che

dove μ 1 ≥ μ 2 ≥…≥μ r ≥ μ r+1 =…=μ n =0,

e r è il rango della matrice A. I numeri μ k sono chiamati numeri spettrali della matrice A. Per una matrice A non singolare, vale quanto segue:

μ1 ≥ μ2 ≥ … ≥μn ≥ 0.

Il successivo fatto importante è che la moltiplicazione per una matrice ortogonale non cambia la norma euclidea del vettore, cioè per qualsiasi vettore x con n elementi e qualsiasi matrice ortogonale U di dimensione n per n l'uguaglianza è vera

|| Ux || = || x ||.

Poiché la moltiplicazione per una matrice ortogonale non cambia la norma spettrale, allora

pertanto, per il numero di condizione della matrice vale la seguente uguaglianza:

Abbiamo due sistemi: Ax = b (con soluzione esatta x) e (con soluzione esatta ). Ovviamente, l'errore δx è una soluzione di un sistema il cui membro destro è una perturbazione δb, cioè sistemi Aδx = δb. Sia D = U T AV la scomposizione in valori singolari della matrice A, allora UDV T = A per il fatto che U e V sono matrici ortogonali. Ulteriore

Ax = b ⇒ UDV T x = b.

Introduciamo la notazione

x" = V T x, b" = U T b,

allora i seguenti sistemi sono equivalenti

Ax = b ⇔ Dx" = b"

In modo del tutto simile, consideriamo il sistema rispetto all'errore

Aδx = δb ⇒ UDV T δx = δb

Introduciamo la notazione

δx" = V T δx, δb" = UT δb,

per cui i seguenti sistemi risultano equivalenti

Aδx = δb ⇔ Dδx" = δb"

Quelli. abbiamo ottenuto sistemi equivalenti semplici con matrice diagonale, sulla cui diagonale si trovano i numeri spettrali della matrice A.

Scegliamo ora in modo particolare i membri destri di questi sistemi, cioè poniamo

dove β > 0, allora sarà la corrispondente soluzione del sistema Dx" = b".

dove μ 1 è il massimo valore singolare della matrice A. Faremo lo stesso con il sistema Dδx" = δ b", ovvero sia

dove γ > 0, allora la corrispondente soluzione del sistema Dδx" = δb" sarà

Dalla conservazione della norma di un vettore quando viene moltiplicato per una matrice ortogonale ne consegue che

β/μ1 = || x" || = || V T x || = || x || e γ/μ n = || δx" || = || V T δx || = ||δx ||.

Esattamente nello stesso modo otteniamo le uguaglianze:

β = || b" || = || U T b || = || b || e γ = || δb" || = || U Tδb || = || δb ||.

e da allora

poi finalmente otteniamo

Quindi, per ogni matrice A è possibile costruire un vettore sul lato destro del sistema e la sua perturbazione tale che l'errore nella soluzione sarà il prodotto del numero di condizione della matrice per l'errore sul lato destro. In MATLAB, la costruzione corrispondente può essere eseguita senza utilizzare la scomposizione di valori singolari (sebbene possa essere trovata utilizzando la funzione svd).

Per prima cosa impostiamo n e otteniamo due matrici ortogonali U e V eseguendo una scomposizione QR di n per n matrici da numeri casuali:

N = 4; = qr(rand(n)); = qr(rand(n));

D = diag();

Successivamente, costruiamo la matrice A utilizzando la matrice diagonale D e le matrici ortogonali U e V

A = U*D*V";

Il numero di condizione della matrice A coinciderà con il numero di condizione della matrice D, che nel nostro esempio è pari a 10 10

Beta = 1; gamma = 1e-5;

e costruisci i vettori b" e δb"

B1 = "; db1 = ";

da cui troviamo i vettori b e δb

X = A\b; x_tilda = A\(b+db);

calcolare i lati sinistro e destro della disuguaglianza

dx = x - x_tilda; DESTRA = norma(db)/norma(b); SINISTRA = norma(dx)/norma(x);

e portarli fuori

Formato breve e disp()

Otteniamo l'uguaglianza

1.Finestra di comando(Finestra di comando).

Le espressioni matematiche vengono scritte sulla riga di comando dopo il prompt " >> ". Per esempio,

Per eseguire l'azione, premere il tasto "Invio".

Per impostazione predefinita, il programma scrive il risultato nella variabile speciale ans.

Per salvare il risultato con un altro nome, utilizzare ad esempio il nome della variabile e il segno di assegnazione “="

>> z = 1,25/3,11

Modifica dentro Finestra di comando Puoi utilizzare solo la riga di comando corrente. Per modificare un comando immesso in precedenza, è necessario posizionare il cursore nella riga di input e utilizzare i tasti “ ” o “ ”.

Se un comando termina con ";", il risultato della sua azione non viene visualizzato sulla riga di comando.

La finestra di comando può essere chiusa con il pulsante “ ”, e il pulsante “ ” serve a separare la finestra dall'interfaccia di sistema. È possibile tornare al modulo della finestra standard utilizzando il menu:

Menu principale → Desktop → Disposizione della scrivania → Predefinito.

È possibile cancellare la finestra di comando utilizzando il menu:

Menu principale → Modificare → Cancella finestra di comando.

Menù principale del sistema MatLab.

Menu principale MatLab contiene i seguenti elementi:

File(File) – lavora con i file;

Modificare(Modifica) – modifica;

Visualizzazione(Visualizza) – gestione delle finestre;

ragnatela– comunicazione con la società sviluppatrice tramite Internet;

Finestra(Window) – connessione con le finestre di sistema;

Aiuto(Aiuto) – collegamento al sistema di aiuto MatLab.

Barra degli strumenti del sistema MATLAB.

I pulsanti della barra degli strumenti hanno i seguenti scopi:

Nuovo file(Nuovo) – visualizza le finestre dell'editor di file;

Apri il file(Apri) – apre le finestre di download dei file;

Taglio(Taglia) – taglia il frammento selezionato e lo inserisce negli appunti;

copia(Copia) – copia il frammento selezionato negli appunti;

Impasto(Incolla) – trasferisce il frammento selezionato dagli appunti alla riga di input;

Disfare(Annulla) – annulla il risultato dell'operazione precedente;

Rifare(Ripeti) – ripete il risultato dell'operazione precedente;

Simulink– crea un modello Simulink (estensioni MatLab);

Finestra di aiuto(Aiuto) – apre le finestre della guida.

4. Formato di output dei risultati del calcolo .

Quando si inseriscono numeri reali, viene utilizzato un punto per separare la parte frazionaria!

>> s = 0,3467819

Il risultato del calcolo viene visualizzato nel formato corto(notazione breve per un numero), definita per impostazione predefinita. Puoi cambiare il formato in lungo(notazione di numeri lunghi).

>> formato lungo

0.34678190000000

Sulla lista Formato numerico formati numerici disponibili

corto– una breve annotazione del numero;

lungo– notazione lunga dei numeri;

breve e– breve notazione di un numero in formato a virgola mobile;

lungo e– registrazione lunga di un numero in formato virgola mobile;

corto g– la seconda forma di notazione breve di un numero;

lungo G– la seconda forma della notazione lunga di un numero;

Il formato per la visualizzazione dei dati numerici può essere impostato nel menu File(file) elemento Preferenze(preferenze). Vai alla scheda Finestra di comando(finestra di comando). In opzione Visualizzazione del testo(visualizzazione testo), nell'elenco Formato numerico(formato numerico) impostato corto g, in opzione Visualizzazione numerica(visualizzare i numeri) impostato compatto. Questi formati di output generano numeri in una forma universale di cinque cifre significative e sopprimono gli spazi bianchi tra le righe.

Nozioni di base sui calcoli in MatLab.

Per eseguire semplici operazioni aritmetiche in MatLab si utilizzano gli operatori:

· addizione e sottrazione +, – ;

· moltiplicazione e divisione *, / ;

· esponenziazione ^ .

Alcune variabili speciali:

ans – il risultato dell'ultima operazione senza segno di assegnazione;

eps – errore relativo nei calcoli in virgola mobile;

pi greco – numero;

i o j – unità immaginaria;

Inf – infinito;

NaN – valore indefinito.

Alcune funzioni elementari integrateMatLab:

exp(x) – esponente del numero x;

log(x) – logaritmo naturale del numero x;

sqrt(x) – radice quadrata del numero x;

abs(x) – modulo del numero x;

sin(x), cos(x), tan(x), cot(x) – seno, coseno, tangente, cotangente del numero x;

asin(x), acos(x), atan(x), acot(x) – arcoseno, arcocoseno, arcotangente, arcocotangente del numero x;

sec(x), csc(x) – secante, cosecante del numero x;

round(x) – arrotondando il numero x all'intero più vicino;

mod(x,y) – resto della divisione intera di x per y, tenendo conto del segno;

sign(x) – restituisce il segno del numero x.

Calcoliamo il valore dell'espressione

>> exp(–2.5)*log(11.3)^0.3 – sqrt((sin(2.45*pi)+cos(3.78*pi))/tan(3.3))

Se l'operatore non può essere posizionato su una riga, è possibile continuare a inserirlo nella riga successiva indicando il segno di continuazione "..." alla fine della prima riga, ad esempio:

>> exp(–2.5)*log(11.3)^0.3 – ...

sqrt((sin(2,45*pi)+cos(3,78*pi))/tan(3,3))

Funzioni per lavorare con numeri complessi:

Immissione di un numero complesso

>> z = 3 + 4i

3.0000 + 4.0000i

Le funzioni abs(z), angolo(z) restituiscono il modulo e l'argomento di un numero complesso, dove , ;

complesso(a,b) costruisce un numero complesso dalle sue parti reale e immaginaria:

conj(z)restituisce il numero complesso coniugato;

imag(z), real(z) restituisce la parte immaginaria e reale del numero complesso z.

6. Vettori.

Inserire, aggiungere, sottrarre, moltiplicare per un numero.

Vettore dentro MatLab formato utilizzando l'operatore parentesi quadre. In questo caso, gli elementi di un vettore colonna sono separati da un punto e virgola “;”, e gli elementi di un vettore riga sono separati da uno spazio “” o una virgola “,”.

Introduciamo un vettore colonna.

>> x =

Introduciamo un vettore riga .

>> y =

Per trasporre un vettore, utilizzare l’apostrofo “’”:

>> z = y’

Per trovare la somma e la differenza dei vettori, utilizzare i segni “+” e “–”:

>> c = x+z

La moltiplicazione di un vettore per un numero viene eseguita sia a destra che a sinistra utilizzando il segno "*".

I vettori possono essere argomenti per funzioni integrate, ad es.

>> d = peccato(c)

Per riferirsi agli elementi dei vettori si usano le parentesi (), ad esempio,

>> x_2 = x(2)

L'ultimo elemento del vettore può essere selezionato digitando il comando

>> X_fine = x(fine)

Da diversi vettori puoi crearne uno, ad esempio

>> r =

1.3 5.4 6.9 7.1 3.5 8.2

Il carattere due punti ":" viene utilizzato, ad esempio, per separare più elementi da un vettore

>> w = r(3:5)

I due punti " : " consentono anche di sostituire elementi di un vettore, ad esempio,

>> r(3:5)= 0

1.3 5.4 0 0 0 8.2

Il simbolo “:” può essere utilizzato anche per costruire un vettore, ogni elemento del quale differisce dal precedente per un numero costante, cioè passo, ad esempio

>> h =

1 1.2 1.4 1.6 1.8 2

Il passo può essere negativo (in questo caso il numero iniziale deve essere maggiore del numero finale).

Un passo uguale a uno può essere omesso

>> k =

Funzioni di base per lavorare con i vettori.

• length(x) – determinazione della lunghezza del vettore x;
• prod(x) – moltiplicazione di tutti gli elementi del vettore x;
• sum(x) – somma di tutti gli elementi del vettore x;
• max(x) – trovare l'elemento massimo del vettore x;
• min(x) – trovare l'elemento minimo del vettore x.

Se chiami la funzione min o max con due argomenti di output = min(x),

quindi alla prima variabile viene assegnato il valore dell'elemento minimo (massimo) e alla seconda variabile viene assegnato il numero di questo elemento.

7 Matrici.

Vari metodi di input della matrice.

1. Una matrice può essere inserita come un vettore colonna costituito da due elementi, ciascuno dei quali è un vettore riga ed è separato da un punto e virgola. Ad esempio, introduciamo la matrice

>> A =

2. È possibile inserire la matrice riga per riga eseguendo la sequenza di comandi:

>> A = )