Ministero dell'Istruzione della Repubblica di Bielorussia

Istituto d'Istruzione

· contorno3(X, Y, Z) – costruisce linee di contorno per la superficie ottenute sovrapponendo una figura tridimensionale da un numero di piani secanti posti parallelamente al piano di riferimento della figura.

L'esempio seguente mostra come utilizzare i comandi descritti per tracciare la superficie Z = sin(X) cos(X).

>> =griglia(-3:0.1:3,-3:0.1:3); Z=peccato(X).*cos(X);

>> subplot(3,2,1), plot3(X, Y, Z) % Figura 4.3 (a)

>> sottotrama(3,2,2), mesh(X, Y,Z) % Figura 4.3 (b)

>> sottotrama(3,2,3), surf(X, Y, Z) % Figura 4.3 (c)

>> sottotrama(3,2,4), surfc(X, Y, Z) % Figura 4.3(d)

>> sottotrama(3,2,5),meshz(X, Y,Z) % Figura 4.3(e)

>> sottotrama(3,2,6),contorno3(X, Y, Z) % Figura 4.3(e)

4.3. Formattazione del grafico

Il sistema Matlab offre la possibilità di regolare e correggere le proprietà dei grafici sia tramite l'interfaccia della finestra grafica che impostando gli opportuni comandi e parametri grafici. Tabella 4.3.1. vengono fornite alcune semplici tecniche di formattazione dei grafici.

Tabella 4.2 - Formattazione del grafico

Azione

Strumenti della GUI

Passa alla modalità di modifica.

Fare clic sul pulsante ModificareComplotto nella barra degli strumenti della finestra grafica.

Linee di formattazione e marcatori punti di riferimento grafici.

In modalità di modifica, fare doppio clic sulla linea del grafico con il tasto sinistro del mouse. Nella finestra che appare proprietàEditore–linea impostare tutti i parametri di linea necessari (spessore, stile, colore, ecc.).

traccia(X, Y, S), traccia3(X, Y, Z, S)

(La descrizione dei comandi è data al paragrafo 4.1 e al paragrafo 4.2)

Formattazione degli assi del grafico.

In modalità di modifica, fare doppio clic sull'asse del grafico. Nella finestra che appare proprietàEditore–assi impostare tutti i parametri necessari degli assi.

Il titolo del grafico e le etichette degli assi possono essere utilizzati anche con i comandi InserireTitolo, InserireEtichetta X,InserireEtichetta Y il menu principale della finestra grafica.

assi– controlla le proprietà degli assi.

griglia- abilita e disabilita la griglia di coordinate.

etichetta x(S),etichetta y(S),etichetta (S)– imposta le etichette vicino agli assi. Qui S è una costante o variabile di stringa.

titolo(S)- visualizza il titolo del grafico

Disegnare iscrizioni direttamente sul grafico.

Fare clic sul pulsante InserireTesto, fissa la posizione dell'iscrizione con un clic del mouse e inserisci il testo desiderato.

testo(X,Y,S)– aggiunge il testo specificato dalla stringa S al grafico bidimensionale, in modo che l'inizio del testo si trovi nel punto con coordinate (X, Y).

testo(X,Y,Z,S)- aggiunge testo a un grafico 3D.

Tracciare linee e linee con le frecce direttamente sul grafico

Fare clic su uno dei pulsanti InserireFreccia o Inserirelinea. Posiziona il puntatore del mouse nella posizione desiderata sulla carta e, tenendo premuto il pulsante sinistro del mouse, traccia una linea.

Costruire una leggenda

inserire, e poi il comando leggenda.

leggenda(S1,S2,S3,…)– aggiunge una legenda al grafico corrente con spiegazioni sotto forma di stringhe specificate nell'elenco dei parametri.

Uscita scala colore

Nel menu principale della finestra grafica, selezionare il comando inserire, e poi il comando barra dei colori.

barra dei colori('vert'),barra dei colori('orizzontale')– Visualizza una scala di colori verticale o orizzontale.

Rotazione del grafico

Fare clic sul pulsante Ruota 3D e ruotare il grafico con il mouse (può essere utilizzato anche per grafici bidimensionali).

ruotare3d– imposta la rotazione di una figura tridimensionale.

con condizioni al contorno y(t 0 , t finire, p) = y, dove t finire, t 0 punti di inizio e fine degli intervalli. Parametro t(variabile indipendente) non significa necessariamente tempo, sebbene il più delle volte la soluzione DE sia ricercata nel dominio del tempo. Il sistema DE nella forma Cauchy è scritto in modo simile a (1.1), ma sotto y in questo caso è implicito un vettore colonna di variabili dipendenti. Vettore p stabilisce le condizioni iniziali.

Risolvere il DE del secondo e ordine superiore devono essere ridotti a un sistema DE di prim'ordine.

Esistono possibili funzioni di controllo non consentite rispetto alla derivata:

F(t, y, dio/dt) = 0.

(1.2)

Le equazioni (1.2) di solito non possono essere ridotte analiticamente alla forma (1.1). Tuttavia, la soluzione numerica di difficoltà speciali non causa abbastanza per determinare f(y, t) risolvere (1.2) numericamente rispetto alla derivata di dato y e t.

Risolutori di ODE

Per risolvere i sistemi ODE in MATLAB, vengono implementati vari metodi numerici. Le loro implementazioni sono denominate risolutori ODE.

In questa sezione, il nome generico solutore indica uno dei possibili metodi numerici per risolvere le ODE: ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb, bvp4c o pdepe.

I solutori implementano i seguenti metodi per risolvere i sistemi DE:

Ode45 metodi Runge-Kutta espliciti in un passaggio di 4° e 5° ordine modificati da Dormand e Prinz. Questo è il metodo classico consigliato per una prima prova della soluzione. In molti casi dà buoni risultati se il sistema di equazioni da risolvere non è rigido.

Ode23 metodi Runge-Kutta espliciti in un passaggio del 2° e 4° ordine nella modifica di Bogacki e Champin. Con una moderata rigidità del sistema ODE e requisiti di precisione bassi, questo metodo può dare un guadagno nella velocità della soluzione.

Ode113 Metodo multi-step Adams-Bashworth-Multon della classe del predittore-correttore di ordine variabile. È un metodo adattivo in grado di fornire un'elevata precisione della soluzione.

Metodo dell'ordine variabile multi-step di Ode15 (da 1 a 5, predefinito 5) utilizzando formule numeriche di "differenziazione all'indietro". Questo è un metodo adattivo e dovrebbe essere utilizzato se il solutore ode45 non fornisce una soluzione e il sistema di controllo è rigido.

Ode23s è un metodo one-step che utilizza una formula Rosenbrock modificata del 2° ordine. Può fornire un'elevata velocità di calcolo con una bassa precisione nella risoluzione di un sistema DE rigido.

Ode23t metodo trapezoidale implicito con interpolazione. Questo metodo fornisce buoni risultati quando si risolvono problemi che descrivono sistemi oscillatori con un segnale di uscita quasi armonico. Per i sistemi moderatamente difficili, il controllo remoto può fornire un'elevata precisione della soluzione.

Ode23tb Metodo implicito di Runge Kutta all'inizio della soluzione e un metodo che utilizza formule di "differenziazione all'indietro" del 2° ordine in seguito. Nonostante la precisione relativamente bassa, questo metodo potrebbe essere più efficiente di ode15s.

Bvp4c viene utilizzato per il problema dei valori al contorno dei sistemi di controllo della forma y′ = f(t, y), F(y(un), y(b), p) = 0 (forma lunga sistema di equazioni di Cauchy). I problemi che risolve sono chiamati problemi al contorno a due punti, poiché la soluzione viene cercata impostando condizioni al contorno sia all'inizio che alla fine dell'intervallo di soluzione.

Tutti i risolutori possono risolvere sistemi di equazioni esplicite y′ = F(t, y), e per risolvere sistemi di equazioni rigidi, si consiglia di utilizzare solo solutori speciali ode15s, ode23s, ode23t, ode23tb.

Utilizzo di risolutori ODE

tspan vettore che definisce l'intervallo di integrazione [ t 0 t finale]. Per ottenere soluzioni in momenti specifici t 0 , t 1 , …, t finale(ordinato in ordine decrescente o crescente) dovrebbe essere utilizzato tspan = [t 0 t 1 … t finale];

y 0 vettore delle condizioni iniziali;

Argomento Options generato dalla funzione odeset (un'altra funzione odeget o bvpget (solo per bvp4c) permette di visualizzare le opzioni impostate di default o usando la funzione odeset/bvpset);

p 1, p 2,… parametri arbitrari passati alla funzione F;

T, Y matrice decisionale Y, dove ogni riga corrisponde al tempo restituito nel vettore colonna T.

Passiamo alla descrizione della sintassi delle funzioni per la risoluzione dei sistemi di controllo remoto (il nome solutore indica una qualsiasi delle funzioni presentate sopra).

[T,Y]=risolutore(@ F,tspan,y 0) integra un sistema di controllo remoto del form y′ = F(t, y) sull'intervallo tspan con condizioni iniziali y 0 . @F descrittore di una funzione ODE (puoi anche specificare una funzione nella forma " F"). Ogni riga nell'array di soluzioni Y corrisponde al valore temporale restituito nel vettore colonna T.

[T,Y]=risolutore(@ F,tspan,y 0 ,opzioni) fornisce una soluzione simile a quella sopra, ma con opzioni determinate dai valori dell'argomento opzioni creato dalla funzione odeset. I parametri comunemente usati includono l'errore relativo consentito RelTol (predefinito 1e3) e il vettore valori consentiti errore assoluto AbsTol (tutti i componenti sono 1e6 per impostazione predefinita).

[T,Y]=risolutore(@ F,tspan,y 0 ,opzioni, p 1 ,p 2 ...) fornisce una soluzione simile a quella sopra descritta, passando Opzioni extra p 1 , p 2, … a m-file F ogni volta che viene chiamato. Usa opzioni= se non vengono fornite opzioni.

Soluzione ODE del primo ordine

PROCEDURA DI LAVORO

· frontespizio;

dati iniziali della variante;

· la soluzione del problema;

i risultati della risoluzione del problema.

Esempio

Trova una soluzione all'equazione differenziale sull'intervallo per cui a(1,7) = 5,3.

Crea una funzione utente nella finestra dei comandi

[email protetta](x,y);

Nella sintassi della funzione @(x,y) X variabile indipendente, y variabile dipendente, X-cos( y/pi) il lato destro del DE.

Il processo risolutivo viene eseguito chiamando il risolutore nella finestra di comando con la seguente istruzione:

Ode23(g,,);

La costruzione di un grafico con una griglia viene effettuata dai seguenti operatori:

Il risultato è mostrato in fig. 1.1

Riso. 1.2.1. Visualizzazione della soluzione numerica

ESERCIZIO

1. Trova soluzioni DE del primo ordine , soddisfacendo le condizioni iniziali y(x 0 ) = y 0 sull'intervallo [ un,b].

2. Costruire grafici della funzione.

Opzioni attività.

numero di opzione		y(x 0 )=y 0	[un,b]
		y 0 (1,8)=2,6
		y 0 (0,6)=0,8
		y 0 (2,1)=2,5
		y 0 (0,5)=0,6
		y 0 (1,4)=2,2
		y 0 (1,7)=5,3
		y 0 (1,4)=2,5
		y 0 (1,6)=4,6
		y 0 (1,8)=2,6
		y 0 (1,7)=5,3
		y 0 (0,4)=0,8
		y 0 (1,2)=1,4

Laboratorio n. 2

Soluzione di sistemi ODE

SCOPO DEL LAVORO

Per formare le idee degli studenti sull'uso dei sistemi di controllo remoto in vari campi; per infondere la capacità di risolvere il problema di Cauchy per i sistemi di controllo remoto.

PROCEDURA DI LAVORO

1. Studia la parte teorica. Completa i compiti corrispondenti al numero della tua opzione e mostrali all'insegnante.

2. Preparare un rapporto di laboratorio, che dovrebbe contenere:

· frontespizio;

dati iniziali della variante;

· la soluzione del problema;

i risultati della risoluzione del problema.

Esempio

Risolvi il sistema

utilizzando il risolutore ode23().

Soluzione:

1. Creare nell'editor un file m della funzione per il calcolo delle parti giuste del DE.

Lascia che il nome nell'editor di file sia sisdu.m, quindi la funzione potrebbe assomigliare a questa:

funzione z=sisdu(t,y)

z1=-3*y(2)+cos(t)-exp(t);

z2=4*y(2)-cos(t)+2*exp(t);

>> t0=0;tf=5;y0=[-3/17,4/17];

>> =ode23("sisdu",,y0);

>>traccia(t,y)

Riso. 1.3.1. Visualizzazione della soluzione numerica ottenuta utilizzando la funzione ode23.

1. Cosa significa risolvere il problema di Cauchy per un sistema di controllo?

2. Quali metodi esistono per risolvere i sistemi di controllo remoto?

ESERCIZIO

1. Trova una soluzione di controllo remoto

soddisfare le condizioni iniziali sull'intervallo;

2. Costruisci grafici di funzioni.

Ad esempio, la funzione di decisione dell'8a opzione è data:

funzione z=ssisdu(t,y)

% opzione 8

z1=-a*y(1)+a*y(2);

z2=a*y(1)-(a-m)*y(2)+2*m*y(3);

z3=a*y(2)-(a-m)*y(3)+3*m*y(4);

z4=a*y(3)-3*m*y(4);

>> =ode23("ssdu",,);

>>traccia(t,100*y)

Riso. 1.3.2. Visualizzazione della soluzione numerica ottenuta utilizzando la funzione ode23.

Opzioni attività.

numero di opzione	Compiti
un	m
	0,1	1,2
	0,2	1,5
	0,3	1,7
	0,4	1,9
	0,5
	0,6	1,9
	0,7	2,3
	0,8	2,7
	0,9
	0,1	1,5
	0,2	1,1
	0,3

Laboratorio n. 3

1.4 Soluzione dell'ODE n-esimo ordine

SCOPO DEL LAVORO

Formare le idee degli studenti sull'uso di sistemi di controllo di ordine superiore in vari campi; per infondere la capacità di risolvere il problema di Cauchy per l'utilizzo di DE di ordine superiore programmi applicativi; sviluppare le competenze per verificare i risultati ottenuti.

PROCEDURA DI LAVORO

1. Studia la parte teorica. Completa i compiti corrispondenti al numero della tua opzione e mostrali all'insegnante.

2. Preparare un rapporto di laboratorio, che dovrebbe contenere:

· frontespizio;

dati iniziali della variante;

· la soluzione del problema;

i risultati della risoluzione del problema.

Esempio 1

Risolvi DE del secondo ordine in determinate condizioni iniziali .

Soluzione:

Innanzitutto, portiamo il DE al sistema:

1. Creare un file m della funzione per calcolare le parti giuste del DE.

Lascia che il nome del file sia sisdu_3.m, quindi la funzione potrebbe assomigliare a questa:

funzione z=sisdu_3(x,y)

z2=6*x*exp(x)+2*y(2)+y(1);

2. Procedi come segue:

>> x0=0;xf=10;y0=;

>> =ode23("sisdu_3",,y0);

>>traccia(x,y(:,1))

Riso. 1.4.1. Visualizzazione della soluzione numerica ottenuta utilizzando la funzione ode23.

ESEMPIO DOMANDE SULLA DIFESA LAVORATIVA

1. Cosa significa risolvere il problema di Cauchy per DE di ordine superiore?

2. Come portare DU m l'ordine al sistema di controllo remoto?

ESERCIZIO

1. Trovare una soluzione DE che soddisfi le condizioni iniziali sull'intervallo.

2. Costruisci grafici di funzioni.

Opzioni attività.

numero di opzione	Compiti
Equazioni	Condizioni iniziali

Laboratorio n. 4 – 5

Sistemi dinamici (DS)

SCOPO DEL LAVORO

Conoscenza degli studenti con i concetti di base del DS, la loro classificazione, lo spazio delle fasi del DS, l'interpretazione cinematica del sistema DE, l'evoluzione del DS. L'equazione del moto del pendolo. Dinamica dell'oscillatore Van der Pol.

2. Sistema dinamico (DS) un oggetto matematico corrispondente a sistemi reali (fisici, chimici, biologici, ecc.), la cui evoluzione è determinata unicamente dallo stato iniziale. Il DS è determinato da un sistema di equazioni (differenziale, differenza, integrale, ecc.) che consentono l'esistenza di un'unica soluzione per ogni condizione iniziale su un intervallo di tempo infinito.

Lo stato del DS è descritto da un insieme di variabili scelte per ragioni di naturalezza della loro interpretazione, semplicità di descrizione, simmetria, ecc. L'insieme degli stati DS forma uno spazio delle fasi, ogni stato corrisponde a un punto in esso e l'evoluzione è rappresentata da traiettorie (fase). Per determinare la prossimità degli stati, il concetto di distanza viene introdotto nello spazio delle fasi di DS. L'insieme degli stati in un punto fisso nel tempo è caratterizzato da un volume di fase.

Anche la descrizione del DS nel senso di fissare la legge dell'evoluzione consente una grande varietà: viene effettuata con l'ausilio di equazioni differenziali, mappature discrete, con l'ausilio della teoria dei grafi, della teoria delle catene di Markov, ecc. La scelta di uno dei metodi descrittivi imposta un tipo specifico di modello matematico del corrispondente DS.

Modello matematico di DSè considerato dato se vengono introdotte variabili dinamiche (coordinate) del sistema che ne determinano univocamente lo stato, e viene specificata la legge di evoluzione dello stato nel tempo.

A seconda del grado di approssimazione, allo stesso sistema possono essere assegnati diversi modelli matematici. Lo studio dei sistemi reali segue il percorso dello studio dei corrispondenti modelli matematici, il cui miglioramento e sviluppo è determinato dall'analisi dei risultati sperimentali e teorici quando vengono confrontati. A questo proposito, per sistema dinamico si intende proprio il suo modello matematico. Indagando lo stesso DS (ad esempio il movimento di un pendolo), a seconda del grado di considerazione di vari fattori, otterremo diversi modelli matematici.

Il linguaggio di programmazione interpretativo del sistema MATLAB è progettato in modo tale che qualsiasi calcolo (a volte molto complesso) possa essere eseguito in modalità di calcolo diretto, cioè senza che l'utente prepari il programma. Allo stesso tempo, MATLAB svolge le funzioni di un supercalcolatore e lavora in modalità riga di comando.

Il lavoro con il sistema è di natura interattiva e avviene secondo la regola "ha posto una domanda - ha ricevuto una risposta". L'utente digita un'espressione calcolata sulla tastiera, la modifica (se necessario) sulla riga di comando e completa l'immissione premendo il tasto INVIO. A titolo di esempio, la figura mostra i calcoli più semplici e abbastanza ovvi.

Anche da esempi così semplici si possono trarre alcune conclusioni istruttive:

* il simbolo >> è utilizzato per indicare l'inserimento dei dati iniziali;

* i dati vengono inseriti utilizzando un semplice editor di riga;

* per bloccare l'output del risultato dei calcoli di qualche espressione successiva, è necessario impostare un segno; (punto e virgola);

* se non viene specificata alcuna variabile per il valore del risultato del calcolo, MATLAB assegna a tale variabile il nome ans;

* il segno di assegnazione è il segno di uguale familiare ai matematici =, e non il segno combinato: =, come in molti altri linguaggi di programmazione e sistemi matematici;

* le funzioni integrate (ad esempio sin) sono scritte in lettere minuscole e i loro argomenti sono specificati tra parentesi;

* il risultato dei calcoli viene visualizzato nelle righe di output (senza il segno >>);

* il dialogo si svolge nello stile di "ha posto una domanda - ha ricevuto una risposta".

Gli esempi seguenti illustrano come MATLAB può essere utilizzato per eseguire alcune operazioni vettoriali più semplici. La figura mostra anche la finestra del browser file system, disponibile nella scheda Directory corrente. In modalità comando, è più comodo richiamare la finestra del browser del file system dalla barra degli strumenti attivando il pulsante dopo l'elenco delle directory di sistema MATLAB. Potrebbero verificarsi casi di rifiuto dei calcoli quando la directory corrente è impostata in modo errato, se non vengono trovati gli m-file necessari per i calcoli.

Nella maggior parte dei sistemi matematici, il calcolo di sin(V) o exp(V), dove V è un vettore, risulterebbe in un errore, perché sin ed exp devono avere un argomento scalare. Tuttavia, MATLAB è un sistema di matrici e un vettore è un tipo di matrice con una dimensione di 1 × n o n × 1. Pertanto, nel nostro caso, il risultato del calcolo sarà un vettore della stessa dimensione dell'argomento V, ma gli elementi del vettore restituito saranno seni o esponenti degli elementi del vettore V.

La matrice è specificata come una serie di vettori che rappresentano le sue righe e racchiusi tra parentesi quadre. Uno spazio o una virgola viene utilizzato per separare gli elementi dei vettori e un punto e virgola viene utilizzato per separare un vettore da un altro. Per selezionare un singolo elemento della matrice M, viene utilizzata un'espressione della forma M(j,i), dove M è il nome della matrice, j è il numero di riga e i è il numero di colonna.

È conveniente utilizzare il browser Workspace per visualizzare il contenuto degli array. Ogni vettore e matrice in esso contenuto è rappresentato come un quadrato con celle, a destra del quale è indicata la dimensione dell'array. Facendo doppio clic sul quadrato con il mouse viene visualizzata la finestra dell'editor di array. Lavorare con l'editor di array è abbastanza ovvio: non solo puoi visualizzare gli elementi dell'array, ma anche modificarli e sostituirli.

Come si può vedere dagli esempi precedenti, l'inserimento delle espressioni iniziali per i calcoli nel sistema MATLAB viene effettuato nel formato di testo più comune. I risultati dei calcoli vengono visualizzati nello stesso formato, ad eccezione di quelli grafici. Di seguito sono riportati esempi di registrazione dei calcoli eseguiti dal sistema MATLAB sulla riga di comando:

Lavorare con l'editor di array

Per iniziare, seleziona "MATLAB Help" dal menu Help.

>> digita peccato

sin è una funzione incorporata.

>> aiuta il peccato

SIN(X) è il seno degli elementi di X.

Metodi sovraccaricati

>>V=

0.8415 0.9093 0.1411 -0.7568

Errore utilizzando ==> ^

La matrice deve essere quadrata.

È possibile prestare attenzione alla forma delle risposte quando si eseguono operazioni semplici senza specificare una variabile a cui è assegnato il risultato. In questi casi, MATLAB stesso assegna la variabile ans a cui è assegnato il risultato e il cui valore viene poi visualizzato.

Modulo di output e interruzione di riga nella sessione

Da notare le caratteristiche dell'output nel sistema MATLAB. L'output inizia su una nuova riga, con i dati numerici rientrati e i dati di testo non rientrati. Per risparmiare spazio in questo libro, la derivazione sarà data in quanto segue senza traduzione in nuova linea. Ad esempio, l'output del vettore di riga

sarà dato nella forma:

Un'eccezione è l'output dei vettori di colonna e delle matrici: qui verrà mantenuto il modulo di output MATLAB più descrittivo e predefinito.

In alcuni casi, l'espressione matematica inserita può essere così lunga che una riga non è sufficiente per essa. Quindi parte dell'espressione può essere spostata su una nuova riga utilizzando il segno di ellissi "..." (3 o più punti), ad esempio:

s = 1 - 1/2 + 1/3 - 1/4 + 1/5 - 1/6 + 1/7 ...

1/8 + 1/9 - 1/10 + 1/11 - 1/12;

Il numero massimo di caratteri in una riga della modalità di comando è 4096 e in un file m non è limitato, ma è scomodo lavorare con righe così lunghe. A prime versioni non c'erano più di 256 caratteri in una riga.

Esecuzione di esempi di applicazioni MATLAB dalla riga di comando

MATLAB ha molti esempi di applicazioni, alcuni dei quali possono essere eseguiti direttamente dalla riga di comando. Ad esempio, il comando

esegue il m-file bench.m della demo di test del sistema.

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 vari campi applicativi. I metodi numerici utilizzati per la soluzione approssimativa di problemi di meccanica, per il calcolo dei flussi di liquidi e gas, altri processi fisici in mezzi continui, metodi per il calcolo dei circuiti elettrici, metodi per l'approssimazione dei dati portano alla necessità di risolvere sistemi di equazioni algebriche lineari ( e questo non è un elenco completo di fonti per l'occorrenza di sistemi algebrici lineari). I problemi spettrali sorgono, ad esempio, nell'analisi in frequenza delle strutture, nello studio della stabilità dei processi in schemi elettrici e sistemi di controllo e in molti altri ambiti applicativi.

Una 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 la risoluzione di sistemi di equazioni algebriche lineari è piuttosto ricco e include:

• soluzione di sistemi con matrici quadrate e rettangolari;
• risolvere i sistemi con metodi diretti e iterativi (compresi quelli con possibilità di precondizionamento);
• espansioni matriciali;
• memorizzazione di grandi matrici sparse in una forma compatta e algoritmi speciali per la risoluzione di sistemi con tali matrici.

Risolvere i sistemi utilizzando il segno della barra rovesciata

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

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

e usa la barra rovesciata

X = A\f x = 1 1 1

Invece di una barra rovesciata, puoi chiamare la funzione mldivide

X = mldivide(A,f)

Il risultato sarà lo stesso

Calcolando il residuo, ci assicuriamo che la soluzione sia trovata correttamente (in genere non sempre valori piccoli dei componenti residui 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 output

X = 0,2707 0,3439 0,3854

che non ha nulla a che vedere con la soluzione del sistema in esame. Analizziamo perché ciò avvenga scrivendo il sistema con la "matrice" f e il "lato destro" A. In questo caso ci sarà una sola incognita, perché f taglia 3 per 1:

Se c'è una matrice sul lato destro, allora tanti sistemi saranno risolti quante sono le colonne nella matrice, e ogni colonna è un vettore del lato destro del sistema corrispondente (cioè, l'operatore backslash ti permette 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 sovradeterminato (sistemi sovradeterminati) e la sua soluzione approssimativa si cerca minimizzando la norma euclidea del residuo. Quindi per il primo sistema la soluzione è X, che fornisce il minimo all'espressione seguente

(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

Syms 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 lato destro di ciascuno dei quali è la colonna corrispondente della matrice A.

Naturalmente, la matrice del sistema ridefinito non contiene necessariamente una sola colonna. A titolo di 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
y 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 richiediamo l'uguaglianza y(x k)=y k per K=1,2,...,5:

Componiamo la matrice del sistema e il vettore del lato destro

X = (1:0.5:3)" A = ;f = ; c(1)/x + c(2)*log(x); fplot(fun, ) tieni premuto plot(x, y, "or")

Se il sistema ha più equazioni che incognite, allora un tale sistema è chiamato un sistema sottodeterminato e viene risolto anche in MATLAB usando il segno della barra rovesciata. Avrà infinite soluzioni e troverà una soluzione contenente quanti più valori zero possibili. Considera questo esempio

A = ; 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 correlati dalla seguente regola

B/A equivale a (A"\B")"
dove l'apostrofo significa trasporre. Puoi anche usare la funzione mrdivide invece del normale segno di barra.

Per ora, non approfondiremo gli algoritmi per la risoluzione dei sistemi incorporati nell'operazione di backslash. Questa sezione è dedicata a questo. In particolare, la scomposizione LU viene utilizzata per risolvere sistemi a 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 data dimensione in MATLAB, usa la funzione hilb. Risolviamo, ad esempio, un sistema del 13° ordine, il cui lato destro è tale che tutte le unità devono essere la sua soluzione ( K-esimo elemento del lato destro è la somma degli elementi K riga della matrice). A = hilb(13); f = somma(A, 2); x = A\f

Di conseguenza, riceviamo il messaggio di malcondizionamento

Attenzione: Matrix è vicino al singolare o mal ridimensionato. I risultati potrebbero essere imprecisi. RCOND = 2.339949e-018.

(RCOND è una stima per il reciproco del numero di condizione) e una decisione 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

A proposito dei pericoli generati dalla cattiva condizionalità della matrice, è scritto nella sezione. Passiamo ora a una panoramica delle espansioni di matrici disponibili in MATLAB.

Espansioni della matrice di Cholesky, LU e QR

MATLAB ha funzioni per le seguenti espansioni:

• Espansioni di Cholesky - funzione chol;
• Decomposizioni LU - funzione lu;
• Espansioni 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 di destra. Considera il seguente insieme di sistemi:

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

Assumiamo che la matrice UNè rappresentato 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 la matrice UN. Quindi le soluzioni dei sistemi originali si ottengono come segue. Perché il Ax=f(k) e A=BC, poi BCx=f (K) e risolvendo in sequenza Per=f(k) e Cx=y trovare una soluzione K esimo sistema. Per le matrici riempite, la scomposizione viene eseguita in o(n 3) operazioni, dove n- la dimensione della matrice (cioè, nello stesso tempo della soluzione del sistema), e la soluzione del sistema con ogni fattore di espansione in o(n 2) operazioni. Pertanto la decisione N sistemi c con predecomposizione prende o(n 3)+o(n 2)N. La soluzione di ogni sistema senza scomposizione della matrice richiederebbe o(n 3)N operazioni aritmetiche.

Decomposizione di Cholesky per una data matrice UNè trovare una tale matrice triangolare superiore R con elementi diagonali positivi, che è A=R T R. È noto che se la matrice UNè simmetrica e definita positiva (cioè, per qualsiasi vettore Xè vero: x T Ax0, ovvero, che è lo stesso, tutti gli autovalori della matrice sono positivi), allora la scomposizione di Cholesky esiste e, inoltre, è unica.

Si consideri, ad esempio, la scomposizione della matrice di Cholesky

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

Un semplice controllo assicura che l'espansione sia corretta (entro gli errori che si verificano 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, questo diventerà chiaro nel processo di scomposizione, perché sarà necessario estrarre la radice da numero negativo e viene visualizzato il relativo avviso:

A = ; R = col(A) ??? Errore utilizzando ==> chol Matrix deve essere definito positivo.

La funzione chol non verifica la simmetria della matrice. Durante il calcolo, vengono utilizzati gli elementi della matrice originale situati sulla diagonale e sopra:

A = ; R = col(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 LU-decomposizione, cioè trova la matrice triangolare inferiore l e matrice triangolare superiore u tale che A=LU. Un'affermazione più rigorosa è formulata come segue.

Se una Akè il minore principale della matrice quadrata UN taglia n, composto dal primo K righe e colonne (cioè A(1:k, 1:k)) e det( Ak)0 per K=1,2,..., n-1, allora c'è un'unica matrice triangolare inferiore l, la cui diagonale è costituita da uno, e l'unica matrice triangolare superiore u tale che A=LU.

Quando si calcola LU La decomposizione LU potrebbe richiedere la permutazione delle righe per garantire la stabilità numerica del processo di decomposizione, quindi la funzione lu può restituire una matrice l, che è triangolare inferiore fino alla permutazione della linea, ad esempio:

A = ; = 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

Opera LU 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 viene eseguita la decomposizione Cholesky o LU, quindi 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 eseguita utilizzando l'operazione backslash, poiché l'algoritmo incorporato in essa determina le matrici triangolari e applica un modo efficiente per risolvere il sistema con esse, richiedendo o(n 2) operazioni aritmetiche. Considera un esempio di risoluzione del sistema

con preliminare LU-scomposizione della matrice.

A = ; f = ;

Effettuiamo LU-decomposizione

Lu(A);

e risolvere successivamente 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 di utilizzare le parentesi per determinare l'ordine di risoluzione dei sistemi con matrici triangolari. Espressione senza parentesi

X = U\L\f

porterà a un risultato molto diverso, perché prima il U/L, che è equivalente (come discusso sopra) alla risoluzione di sistemi con matrice U e colonne L come vettori di destra. Di conseguenza, si ottiene una matrice, ciascuna colonna della quale è una soluzione del sistema corrispondente, quindi il sistema viene risolto con questa matrice e il vettore del lato destro f. Questo processo ovviamente non ha nulla a che fare con la soluzione del sistema originale.

Considera l'ultima scomposizione della matrice - QR-scomposizione eseguita dalla funzione qr. QR-la decomposizione può essere fatta per matrici rettangolari, precisamente se UN matrice di dimensioni m sul n, allora c'è una matrice ortogonale Q taglia m sul m(cioè tale che Q -1 =Q T) e la matrice R taglia m sul n con tutti gli elementi zero sotto la diagonale principale, che è A=QR.

QR-l'espansione è caratterizzata da una buona stabilità computazionale e viene utilizzata, in particolare, per approssimare i dati con il metodo dei minimi quadrati. A titolo di esempio, consideriamo la soluzione del 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

Dopo che è stato fatto QR-decomposizione, soluzione del sistema Ax=f può essere trovato come una soluzione al sistema a matrice triangolare R, perché perché il QR=f, poi Rx=Q T f:

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

Per ulteriori informazioni sulle espansioni di matrici corrispondenti alle funzioni MATLAB chol, lu e qr e ai loro argomenti, vedere la sezione Scomposizioni di matrici.

Algoritmo per risolvere 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ù appropriato implementato in una delle procedure del pacchetto LAPACK o UMFPACK. Qui b può anche essere una matrice, il cui numero di righe è uguale al numero di righe della matrice A. Quindi viene restituita la matrice x, ciascuna i-esima delle sue colonne 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 nei seguenti passaggi:

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

2. Se A è una matrice quadrata, sparsa ea bande, viene utilizzato il risolutore della matrice di bande. Ci possono essere due opzioni per questo:

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

A e b tipo reale doppie - procedure DGBTRF, DGBTRS
A e b procedure doppie di tipo complesso ZGBTRF, ZGBTRS
A o b procedure singole di tipo reale SGBTRF, SGBTRS
A o b procedure singole di tipo complesso CGBTRF, CGBTRS

La densità di 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 a matrice, ad esempio, per una matrice 7 per 7, il cui modello è fornito sotto, il numero di elementi diversi da zero

nz = 1 (su diag.#-2) + 6 (su diag.#-1) + 7 (su diag.#0) + 6 (su diag.#1) + 1 (su diag.#2) + 1 (nella Fig. n. 3) = 22,

e il numero di tutti gli elementi nel nastro

banda = 5 (su diag.#-2) + 6 (su diag.#-1) + 7 (su diag.#0) + 6 (su diag.#1) + 5 (su diag.#2) + 4 (nella Fig. n. 3) = 33

e la densità degli elementi diversi da zero sarà 2/3. Poiché 2/3 > 0,5, verrà applicato un risolutore a matrice di bande
Parametro bandito, come altri parametri che controllano gli algoritmi per matrici sparse in MATLAB, viene impostato utilizzando la funzione spparms.



Il significato di questo passaggio dell'algoritmo è che la soluzione di un sistema con una matrice di bande non richiede azioni al di fuori della banda. Se il nastro è sufficientemente pieno, non ci saranno molte azioni con zero elementi. Al contrario, se il nastro è sufficientemente sparso, gli approcci forniti nelle fasi successive dell'algoritmo possono portare un effetto maggiore.

• 3. Se A è una matrice triangolare superiore o inferiore, viene applicato il metodo di sostituzione all'indietro, o, di conseguenza, il metodo di sostituzione diretta, in cui la componente soluzione è trovata dall'ultima (o prima equazione) e quindi le componenti trovate sono sostituito nelle seguenti equazioni per trovare i sistemi di equazioni dei componenti della soluzione successiva.
• 4. Se A - può essere ridotto per permutazioni a una forma triangolare, si esegue la corrispondente riduzione e si risolve il sistema di equazioni come al paragrafo 3.

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

un. Se A non è sparso, si tenta di eseguire la scomposizione di Cholesky A = R T R seguita dalla soluzione di sistemi con matrici triangolari R T e R: R T y = b e Rx = y. Questo chiama le seguenti procedure della libreria LAPACK:

Aeb sono reali e doppi - DLANGE, DPOTRF, DPOTRS, DPOCON
A e b sono complessi e doppi: ZLANGE, ZPOTRF, ZPOTRS, ZPOCON
Aeb 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 viene eseguito un controllo preliminare per la definitività positiva, ed è determinato proprio nel processo di decomposizione, poiché abbastanza spesso è sufficiente eseguire diversi passaggi della decomposizione di Cholesky per scoprire il fatto che la matrice non è definita positiva ( durante il calcolo diventa necessario estrarre la radice quadrata di un numero negativo). Se la decomposizione di Cholesky fallisce, si procede al passaggio successivo.

b. Se A è sparso, le permutazioni simmetriche di righe e colonne vengono preliminarmente eseguite secondo l'algoritmo simmetrico di grado minimo (mediante la funzione symmmd) per ridurre il riempimento del fattore di espansione di Cholesky, cioè per 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 = col(A(p, p));

e si risolvono due sistemi con fattori di Cholesky, il primo con il fattore trasposto e le corrispondenti permutazioni nel vettore di destra

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

• 6. Se A è una matrice di Hessenberg, allora si riduce a una matrice triangolare superiore (con l'opportuna modifica del lato destro) e quindi si risolve il sistema per sostituzioni

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

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

U - matrice triangolare superiore
L - una matrice che può essere ridotta a una triangolare mediante permutazioni di righe

e la soluzione del sistema Ax = b si trova per soluzione successiva di sistemi con matrici triangolari Ly = b, Ux = y.
Per eseguire la scomposizione LU, vengono chiamate le seguenti routine della libreria LAPACK:

Aeb sono reali e doppi - DLANGE, DGESV, DGECON
A e b sono complessi e di tipo double - ZLANGE, ZGESV, ZGECON
Aeb sono reali e di tipo singolo - SLANGE, SGESV, SGECON
A e b tipo complesso e unico - CLANGE, CGESV, CGECON

b. Se A è una matrice sparsa, le colonne vengono permutate per ridurre il riempimento dei fattori L e U nel processo di scomposizione. Inoltre, permutando le righe nel processo di decomposizione LU, viene assicurata 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 di libreria UMFPACK

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

un. Se A non è una matrice sparsa, viene eseguita una scomposizione QR AP = QR, dove P è una matrice di permutazione di colonna, Q è una matrice ortogonale (Q T Q = I) e R è triangolare superiore. Se A è m per n, allora Q è m per m e R è m per n. Inoltre, la soluzione del sistema è la seguente:

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 scomposizione 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.

Oltre 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 è mal condizionata (vedere la sezione Effetto della condizionalità della matrice sull'accuratezza della risoluzione un sistema con esso). L'algoritmo di cui sopra contiene tutti i tipi di controlli che possono richiedere una notevole quantità di tempo extra. La sezione seguente, La funzione linsolve per la risoluzione di sistemi di equazioni lineari, fornisce la funzione linsolve, che consente di specificare il tipo di matrice, riducendo così i tempi di calcolo.

funzione linsolve per la risoluzione di sistemi di equazioni lineari

Invece di risolvere un sistema (o sistemi con più lati destri specificati in una matrice) di equazioni algebriche lineari usando il segno della barra rovesciata, puoi usare la funzione linsolve, che non esegue tutti i controlli matriciali inerenti all'algoritmo del \ operazione (vedere la sezione precedente).

La funzione linsolve, chiamata 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, la sua scomposizione LU viene precalcolata e quindi vengono risolti due sistemi con matrici triangolari L e U.
• Se A è una matrice rettangolare, viene precalcolata la sua scomposizione QR e quindi risolto il sistema con matrici triangolari.

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

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

Per eliminare l'output di questo messaggio nella finestra di comando, chiamare la funzione linsolve con due argomenti di output
= linsolve(A, b)
quindi la soluzione ottenuta sarà scritta in x, e il rango della matrice (se A è una matrice rettangolare) o il reciproco della stima per il suo numero di condizione (vedi sezione ), per esempio

A = ; 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 \ compaiono quando si specifica il tipo di matrice del sistema di equazioni. Per fare ciò, prima di chiamare la funzione linsolve, è necessario compilare la struttura con le informazioni sulla matrice con i seguenti campi

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

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

LT	UT	UHESS	SIM	POSDEF	RETT	TRANSA
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, questo deve essere preso in considerazione durante la risoluzione, poiché per matrici definite positive la soluzione si basa sulla decomposizione di Cholesky, che richiede meno operazioni rispetto alla scomposizione LU utilizzata quando si risolvono sistemi con matrici quadrate generali . Questo è facile da verificare utilizzando l'esempio seguente, in cui viene creata una matrice simmetrica definita positiva (a quella trasposta viene aggiunta una matrice di numeri casuali e alla diagonale vengono aggiunti numeri grandi) e il sistema di equazioni con questa matrice si risolve prima come un sistema con una matrice 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 delle operazioni, ad eccezione di SYM e POSDEF opts.TRANSA = false; opts.UHESS = falso; opts.RECT ​​= falso; opts.UT = falso; opts.LT = falso; % crea un vettore di dimensioni di matrice N = 2.^(8:12); % crea array vuoti per registrare i tempi di soluzione TimeGeneral = ; TimePosSym = ; % crea matrici in un ciclo e confronta i tempi di soluzione per n = N % crea una matrice simmetrica e definita positiva % e il vettore del lato 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 simmetria e posit matrix opts.SYM = true; opts.POSDEF = true; Tstart = cputime; x = linsolve( A,b, opts); Tend = cputime; TimePosSym = ; end

I costi computazionali risultanti di questi modi di risolvere il sistema sono mostrati nel grafico sottostante.

Naturalmente, se la matrice è triangolare, allora è molto importante indicarlo, poiché la soluzione di un sistema con una matrice triangolare viene eseguita in O(n 2) operazioni e se un algoritmo di soluzione progettato per un la matrice generale viene applicata a un sistema con una matrice triangolare, quindi l'ordine di O(n 3) operazioni.

La funzione linsolve non verifica se la matrice soddisfa le proprietà specificate, il che potrebbe causare un errore. Ad esempio, se durante la risoluzione di un sistema con la seguente matrice e il lato destro

A = ; b = ; imposta il campo UT su true (e imposta tutti gli altri su false) opts.UT = true; opts.TRANSA = falso; opts.LT = falso; opts.UHESS = falso; opts.SYM=falso; opts.POSDEF = falso; opts.RECT ​​= falso; quindi quando risolve il sistema, la funzione linsolve considererà la matrice come un 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 = ;

Influenza della condizionalità di una matrice sull'accuratezza della risoluzione di un sistema con essa

In questo paragrafo consideriamo come gli errori negli elementi del vettore di destra e nella matrice del sistema di equazioni lineari Ax = b possono influenzare la soluzione di questo sistema. Consideriamo innanzitutto il caso di introdurre perturbazioni nel vettore di destra b. Quindi risolviamo due sistemi

inoltre, si assume di risolvere esattamente ciascuno dei sistemi, e la domanda principale che si pone è quanto la soluzione x del sistema (1) differirà dalla soluzione del sistema (2) con una perturbazione sul lato destro δb. Questa è una domanda piuttosto importante, poiché gli elementi del membro destro possono essere ottenuti come risultato di alcune misurazioni, contenenti rispettivamente l'errore δb k in ciascuna componente b k . Sarebbe auspicabile che quando si misura la stessa quantità (ogni volta con i propri piccoli errori), anche le soluzioni corrispondenti di sistemi con lati destri leggermente diversi non differissero molto l'una dall'altra. Sfortunatamente, questo non è sempre il caso. La ragione di ciò è la caratteristica della matrice, chiamata sua condizionata. Questo sarà discusso ulteriormente.

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

Rispettivamente

dove n è il numero di incognite ed equazioni nel sistema. Oltre alla norma vettoriale per stimare la dimensione di un vettore, abbiamo anche bisogno di una norma matriciale per stimare la dimensione di una matrice. Prendi 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 di norma per il calcolo delle norme di matrici e vettori, che, in particolare, può calcolare le norme sopra riportate. Il motivo per cui abbiamo scelto tali norme di vettori e matrici è spiegato in dettaglio nella sezione, dove vengono forniti alcuni calcoli e definizioni. Questo è legato alla stima che useremo per l'errore nella risoluzione del sistema (la derivazione di questa disuguaglianza è data anche nella sezione menzionata):

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

Dalla disuguaglianza di cui sopra segue 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 di .

Considera un classico esempio di matrice mal condizionata - la matrice di Hilbert - e scopri come l'errore sul lato destro del sistema influisca sull'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.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 = uno(n, 1); b = H*x b = 2,4500 1,5929 1,2179 0,9956 0,8456 0,7365

Vediamo che né nella matrice né nel vettore del lato destro c'è qualcosa di "sospetto", tutti i numeri non differiscono molto tra loro. Ora formiamo il lato destro perturbato b+δb sommando piccoli numeri dell'ordine di 10-5 al vettore b e risolviamo il sistema con il lato 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, dove dovrebbero trovarsi tutte le unità. Calcoliamo l'errore relativo nella soluzione e sul lato destro (la funzione norma calcola di default la norma euclidea del vettore):

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

Quindi, l'errore è nella soluzione dell'ordine di uno, sebbene le modifiche sul lato destro fossero dell'ordine di 10 -5 . Questo rientra pienamente nella disuguaglianza di cui sopra per l'errore. Infatti, calcoliamo il numero della condizione cond(H) usando la funzione MATLAB chiamata cond e per default calcola per la norma spettrale della matrice

C = cond(H) c = 1,4951e+007

SINISTRA 1.1475 meno

C*DESTRO? 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 aumenterà solo (questo è facile da verificare impostando n = 7, 8, ...). Inoltre, per n = 12, verrà visualizzato un messaggio in cui si afferma che la matrice è mal condizionata e la soluzione potrebbe risultare errata

Attenzione: Matrix è vicino al singolare o mal ridimensionato.
I risultati potrebbero essere imprecisi. RCOND = 2.409320e-017.

RCOND viene scelto qui come misura di condizionalità. uguale a uno diviso per la stima del numero di condizione (il numero di condizione viene stimato utilizzando un algoritmo abbastanza veloce che è molto più veloce di cond, come descritto più dettagliatamente nella guida alla funzione rcond). 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 usando un semplice ciclo e la funzione di semilogia (la scala lungo l'asse y è logaritmica):

N = 1:20; C = zeri(1, 20); per n = N H = hilb(n); C(n) = cond(H); end semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

A proposito, poiché la funzione cond trova il numero della condizione utilizzando un metodo numerico (vale a dire, la scomposizione del valore singolare per trovare valori singolari), il numero della condizione dopo n = 12 è già calcolato in modo errato, infatti, dovrebbe crescere ulteriormente, come può essere visto utilizzando calcoli simbolici in MATLAB e operazioni con un dato numero di cifre significative

N = 1:20; C = zeri(1, N); digits(60) per n = N H = vpa(sym(hilb(n))); % calcolo della matrice di Hilbert fino alla 60a cifra sigma = svd(H); % trova i valori singolari della matrice di Hilbert % calcola il numero di condizione della matrice di Hilbert C(n) = max(double(sigma))/min(double(sigma)); end semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

Consideriamo ora tre punti importanti relativi a questo esempio.

Il primo di questi - risolvendo il sistema Hx = b (con il vettore sul lato destro corrispondente alla soluzione delle unità) usando l'operatore backslash non darà unità esatte, l'errore sarà già nella decima cifra (sebbene in MATLAB tutti i calcoli siano eseguito di default in doppia precisione)

Formato lungo H\b ans = 0.999999999999916 1.00000000002391 0.999999999983793 1.00000000042209 0.99999999953366 1.00000000018389

Ciò è dovuto al fatto che quando si calcola il vettore b moltiplicando la matrice H per un vettore di tutte le unità, abbiamo già inserito qualche errore. Inoltre, anche gli errori di arrotondamento nel processo di risoluzione del sistema hanno giocato un ruolo e la cattiva condizionalità della matrice (anche di piccole dimensioni) ha portato a tali errori nella soluzione. Infatti, per piccole matrici con un piccolo numero di condizione, questo effetto non sarà osservato. Prendi una matrice 6 per 6 di numeri casuali, calcola il suo numero di condizione

A = rand(n); cond(A) ans = 57.35245279907571

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

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

e risolvere il sistema, ottenendo una buona precisione

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

Secondo punto importante associato al determinante matriciale. Calcola il determinante della matrice di Hilbert 6 per 6 con cui abbiamo risolto il sistema

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

Il determinante risulta essere abbastanza piccolo, da cui si può trarre una 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 effetti, questo non è il caso, se il determinante del sistema è piccolo, ciò non significa che la matrice sia mal condizionata e, di conseguenza, c'è un grande errore nella soluzione. Infatti, prendi la matrice

A = ;

Il suo determinante è molto piccolo, dell'ordine di 10 -121

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

(la parola "molto" è ovviamente arbitraria, ma è inferiore al determinante della matrice di Hilbert 6 per 6, quando si risolve il sistema con cui abbiamo avuto problemi). Tuttavia, la piccolezza del determinante non influirà sull'errore nella soluzione quando il lato destro del sistema è perturbato, il che è facile da mostrare formando un sistema con una soluzione nota, introducendo perturbazioni sul 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)/norma(x) DESTRA = norma(delta_b)/norma(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, ovvero:

C = cond(A) c = 1,0303

E, infine, si consideri la terza questione, relativa al raggiungimento del segno di uguale nella disuguaglianza per l'errore nella soluzione

Cioè, in altre parole, scopriamo: può esserci una situazione del genere quando facciamo una piccola perturbazione relativa del lato destro del sistema, diciamo 10 -5 , il numero di condizione della matrice del sistema è ​​10 10, e nella soluzione si ottiene un errore relativo di 10 5. Puoi passare attraverso molti esempi, iterando su diverse matrici, non ottenendo uguali e concludendo erroneamente che questa è solo una sovrastima dell'errore nella soluzione. Tuttavia, questo non è il caso, poiché l'esempio seguente ci convince.

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

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

L'errore relativo nella soluzione 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, possiamo vedere 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. Una volta risolta in aritmetica esatta, l'uguaglianza risulterebbe esatta, che può essere mostrata analiticamente (vedi sezione ).

Negli esempi precedenti, abbiamo assunto che la matrice non contenga errori e che gli errori possono trovarsi solo nella parte destra del sistema di equazioni lineari. In pratica, può risultare che gli elementi della matrice di sistema sono specificati con errori, ad es. invece del sistema

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

dove || || M - qualche norma matriciale. Per poterlo fare, la norma matrice || || M e la 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 matriciale || || Si chiama M concordato con norma vettoriale || || V. È noto, ad esempio, che la norma spettrale

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

Ecco perché nella sezione abbiamo condotto tutti gli esperimenti che coinvolgono queste norme.

Avendo ottenuto la disuguaglianza || A || M || x || V ≥ || b || V nota inoltre che da Ax = b e ne consegue che . Poiché la matrice A non è singolare, ne consegue che δx = A -1 δb e || δx || V = || A -1 δb || V. Usiamo ancora la consistenza 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

dividiamo il valore minore di una delle disuguaglianze per il valore maggiore dell'altra disuguaglianza, e quello maggiore, rispettivamente, per quello minore

e con una semplice trasformazione otteniamo infine 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 matriciale è coerente con quella vettoriale. Questo è 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 di colonne, calcolata dalla formula

coerente con la prima norma vettoriale

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

coerente con la norma del vettore infinito

La funzione norma calcola non solo la norma euclidea di un vettore e la norma spettrale di una matrice, ma anche le norme vettoriali e matriciali sopra elencate. Per fare questo, devi chiamarlo con secondo aggiuntivo discussione:

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

Il numero della condizione della matrice cond(A) = || A || M* || A -1 || M rispetto a varie norme matriciali può essere calcolato usando la funzione cond. Se cond viene chiamato con un argomento di input (una matrice), il numero della condizione viene calcolato rispetto alla norma della matrice spettrale. Quando si chiama la funzione cond con un argomento aggiuntivo, i numeri di condizione vengono restituiti rispetto alla norma matrice specificata:

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

A titolo di esempio per dimostrare l'importanza della coerenza del foro della matrice con la norma del vettore nella disuguaglianza di errore, faremo un esempio con la matrice A, il vettore di destra b e il vettore di errore di destra delta_b .

A = ; 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,208503222383081e-5];

Calcoliamo il diritto e lato sinistro stime utilizzando la prima norma vettoriale e il numero di condizione della matrice rispetto 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 le parti sinistra e destra della disuguaglianza
9.9484e+004 9.9323e+004

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

Analizziamo 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 - qualche norma matriciale A. La norma matriciale (cioè la regola con cui a ciascuna matrice viene assegnato un numero) 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 zero.

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.

Si può vedere dall'ultimo assioma che la norma è definita solo per matrici quadrate (sebbene, in linea di principio, non vi sia tale restrizione nelle formule precedenti per il calcolo di varie norme). Inoltre, dall'ultimo assioma consegue che ogni norma della matrice identitaria I non è inferiore a una, anzi

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

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

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

Completiamo questa sezione studiando il motivo della comparsa dell'esatta uguaglianza nella stima dell'errore nella soluzione sotto la perturbazione del lato destro:

e costruire gli esempi corrispondenti in MATLAB. (Il ragionamento riportato di seguito è contenuto, 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 è svolto dal teorema sulla scomposizione al valore singolare di una matrice, secondo il quale, per ogni matrice reale di dimensione on, 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 detti numeri spettrali della matrice A. Per una matrice A non singolare vale quanto segue:

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

Il prossimo fatto importante è che la moltiplicazione per una matrice ortogonale non cambia la norma euclidea di un vettore, cioè per qualsiasi vettore x con n elementi e qualsiasi matrice ortogonale U di dimensione n per n, abbiamo

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

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

pertanto, per il numero di condizione della matrice, l'uguaglianza è vera

Abbiamo due sistemi: Ax = b (con soluzione esatta x) e (con soluzione esatta). È ovvio che l'errore δx è la soluzione del sistema, il cui lato destro è la perturbazione δb, cioè sistemi Aδx = δb. Sia D = U T AV la scomposizione al valore singolare 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"

Esattamente allo stesso modo, consideriamo il sistema rispetto all'errore

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

Introduciamo la notazione

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

in cui i seguenti sistemi sono equivalenti

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

Quelli. abbiamo ottenuto sistemi equivalenti semplici con una matrice diagonale la cui diagonale contiene i numeri spettrali della matrice A.

Scegliamo ora i lati destri di questi sistemi in modo speciale, cioè lett

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

dove μ 1 è il massimo valore singolare della matrice A. Procederemo in modo analogo con il sistema Dδx" = δ b" , cioè sia

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

Poiché la norma di un vettore viene preservata quando viene moltiplicata per una matrice ortogonale, ne consegue che

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

Allo stesso modo otteniamo le uguaglianze:

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

e da allora

poi finalmente arriviamo

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

Innanzitutto, specifichiamo 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 è 10 10

beta=1; gamma = 1e-5;

e costruisci i vettori b" e δb"

B1="; db1=";

per cui troviamo i vettori b e δb

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

calcola le parti sinistra e destra della disuguaglianza

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

e tirali fuori

formato breve e disp()

Otteniamo l'uguaglianza

1. Puoi avviare il programma utilizzando l'icona sul desktop o tramite il pulsante Start (nell'angolo in basso a sinistra dello schermo). Si aprirà l'area di lavoro del pacchetto MatLab.

2. Fare clic nella casella della finestra dei comandi per attivarla.

Digitare 1+2 nella riga con l'icona » e un cursore verticale lampeggiante e premere Invio. Di conseguenza, nella finestra di comando di MatLab viene visualizzato quanto segue:

Il risultato del calcolo della somma 1+2 viene scritto nella variabile speciale ans e il suo valore, pari a 3, viene visualizzato nella finestra di comando. Sotto la risposta c'è una riga di comando con un cursore lampeggiante, che indica che MatLab è pronto per ulteriori calcoli. Puoi digitare nuove espressioni sulla riga di comando e trovarne i valori.

3. Per continuare a lavorare con l'espressione precedente, ad esempio per calcolare (1 + 2) / 4.5, è consigliabile utilizzare il risultato già esistente, che è memorizzato nella variabile ans.

Digita ans/4.5 (un punto viene utilizzato quando si inseriscono i decimali) e premi Invio, risulta:

4. Il record "ans = 0,6667" è chiamato eco.

L'esecuzione di ogni comando in MatLab è accompagnata da un'eco simile, che spesso rende difficile percepire il funzionamento del programma.

Per disattivare l'eco, metti un carattere dopo il comando; (punto e virgola). Per esempio:

Qui, il risultato della moltiplicazione di ans per 3, che è intermedio, non viene visualizzato. Viene mostrata solo la risposta finale.

5.Salvare i valori delle variabili. Per questo:

- selezionare Salva spazio di lavoro con nome dal menu File;

— nella finestra di dialogo Salva variabili dell'area di lavoro visualizzata, specificare la directory e il nome del file (per impostazione predefinita, si suggerisce di salvare il file nella sottodirectory di lavoro della directory principale di MatLab). I risultati del lavoro verranno salvati in un file con estensione mat.

6. Chiudi MatLab.

7.Avvia nuovamente MatLab. Per ripristinare i valori delle variabili della sessione precedente, aprire il file salvato utilizzando la sottovoce Apri del menu File. Le variabili salvate possono essere utilizzate nei comandi appena inseriti.

8. Annota i comandi che esegui e i risultati file di testo, da cui è quindi possibile leggere o stampare editor di testo. Per questo:

- inserire il comando diario;

- specificare il nome del file in cui verrà archiviato il registro di lavoro come argomento del comando diary.

Un esempio è fornito nella clausola 1.3.

9. Per terminare il lavoro nel sistema MatLab, digitare il comando quit.

1. Studia la parte teorica.

2. Ottieni un'opzione per l'attività.

3. Fare l'esempio di cui al paragrafo 2.

4. Esegui i calcoli in base alla tua opzione.

5. Invia un rapporto in in formato elettronico.

6. Proteggi lavoro di laboratorio rispondendo alle domande dell'insegnante.

Opzioni

Articoli da leggere:

Nozioni di base sull'utilizzo di Mathcad. I calcoli più semplici. Lezione 4