dispensamatlab

Annamaria Mazzia

Esercitazioni di CalcoloNumerico in MATLAB e

Octave

Dipartimento di Ingegneria Civile Edile e AmbientaleUniversit degli Studi di Padova

Creative Commons Attribuzione- Non commerciale -Non opere derivate 3.0 Italia License

a.a. 2012/2013

Indice

Indice ii

1 Introduzione 11.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Avvio di MATLAB e Octave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2.1 Le variabili . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.2 Prima introduzione (breve) su matrici e vettori in MATLAB . . . . . . . . . . . . . . . . 51.2.3 Funzioni vettorizzate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.4 Sulle variabili . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.3 Comandi utili . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.4 Introduzione alla programmazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.5 Problemi e Algoritmi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.6 MATLAB come linguaggio di programmazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131.7 Primo script e prima function in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . 131.8 I predicati elementari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151.9 Strutture e cicli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.9.1 Ciclo if . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.9.2 Ciclo while . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241.9.3 Ciclo for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 331.9.4 Introduzione a file di scrittura dati . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 351.9.5 Modificare le function di punto fisso per creare nuove function . . . . . . . . . . . . . . 39

2 Su interpolazione e approssimazione 432.1 Interpolazione monomiale e con i polinomi di Lagrange . . . . . . . . . . . . . . . . . . . . . . . 442.2 Interpolazione con la tabella delle differenze divise di Newton . . . . . . . . . . . . . . . . . . . . 502.3 Interpolazione lineare a tratti e spline cubica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

2.3.1 Interpolazione lineare a tratti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 542.4 Curve parametriche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

2.4.1 Il ciclo switch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 602.5 Approssimazione di dati . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

2.5.1 Retta di regressione sugli scarti verticali . . . . . . . . . . . . . . . . . . . . . . . . . . . . 622.5.2 Modello potenza e modello esponenziale . . . . . . . . . . . . . . . . . . . . . . . . . . . 642.5.3 Caricare i dati di input da file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

3 Metodi diretti per sistemi lineari 693.1 Prodotti matrice-vettore, matrice-matrice, scalari . . . . . . . . . . . . . . . . . . . . . . . . . . . 713.2 Metodi di sostituzione allindietro e in avanti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 723.3 Metodi diretti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

ii

Indice

3.4 Norme di vettori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

4 Metodi iterativi per sistemi lineari 794.1 Metodo di Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.2 Modificare jacobi.m per implementare gli altri metodi . . . . . . . . . . . . . . . . . . . . . . . 83

4.2.1 Metodo di Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.2.2 Metodo di rilassamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

5 Problemi non lineari 875.1 Metodo di Newton per sistemi non lineari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875.2 Minimi quadrati non lineari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

6 Formule di integrazione numerica 956.1 Formule di quadratura numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 956.2 La formula dei trapezi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 956.3 La formula di Cavalieri-Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

7 Equazioni Differenziali Ordinarie 1037.1 Metodi per ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

7.1.1 Confronto tra i metodi sullequazione test . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.2 I metodi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1057.3 Sistemi di equazioni differenziali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

7.3.1 Function di MATLAB e Octave per ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1097.3.2 Loscillatore semplice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1107.3.3 Equazioni del moto di Keplero . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

iii

CAPITOLO 1

Introduzione

Luniverso non potr essere lettofinch non avremo imparato illinguaggio ed avremo familiarizzatocon i caratteri con cui scritto. scritto in linguaggio matematico, ele lettere sono triangoli, cerchi edaltre figure geometriche, senza lequali umanamente impossibilecomprendere una singola parola.

Galileo Galilei

1.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Avvio di MATLAB e Octave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2.1 Le variabili . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.2 Prima introduzione (breve) su matrici e vettori in MATLAB . . . . . . . . . . . . . . . . . 51.2.3 Funzioni vettorizzate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.4 Sulle variabili . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.3 Comandi utili . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.4 Introduzione alla programmazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.5 Problemi e Algoritmi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.6 MATLAB come linguaggio di programmazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131.7 Primo script e prima function in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131.8 I predicati elementari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151.9 Strutture e cicli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.9.1 Ciclo if . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.9.2 Ciclo while . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241.9.3 Ciclo for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 331.9.4 Introduzione a file di scrittura dati . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 351.9.5 Modificare le function di punto fisso per creare nuove function . . . . . . . . . . . . . . . 39

1.1 Introduzione

MATLAB un acronimo per MATrix LABoratory ed uno strumento utile per il calcolo numerico e per lavisualizzazione. MATLAB , infatti, pu essere considerato un linguaggio di programmazione ad alto livelloe, allo stesso tempo, un ambiente interattivo per lo sviluppo di algoritmi, per la visualizzazione e lanalisi didati, e per il calcolo numerico.

1

1. INTRODUZIONE

Quindi, per imparare a usare bene MATLAB , dobbiamo imparare le tecniche di programmazione (cosa un programma, quali sono le strutture di un programma, come si scrive, come si esegue). Ma dobbiamo an-che tenere in conto che MATLAB ha le sue proprie funzioni (functions), che permettono di eseguire calcolimolto sofisticati, ed capace di fare grafici e visualizzare i risultati in maniera diretta.

Molto spesso la conoscenza di MATLAB un requisito richiesto per molte posizioni lavorative in ambi-to ingegneristico proprio per la sua implementazione in numerosi problemi (dallingegneria aerospaziale aquella chimica, alla meccanica ...).

MATLAB un programma complesso (originariamente fu scritto in FORTRAN e successivamente riscrit-to in linguaggio C) che va bene per risolvere programmi delaborazione numerica, per lavorare con matricie per la grafica. Questo software un prodotto della The MathWorks : vedremo in particolare la versioneMATLAB 8.0.0.783 (R2012b) lavorando in ambiente Linux.

Un prodotto simile a MATLAB ma open source GNU Octave (si vada sul sito http://www.gnu.org/software/octave/). In queste dispense, faremo riferimento alla versione 3.6.2 di Octave.

Nel seguito vedremo come si lavora in MATLAB (e Octave) tenendo presente che quasi tutti i coman-di che daremo in MATLAB valgono alla stessa maniera anche per Octave con alcune differenze, comelinterfaccia grafico o altri comandi, che sottolineeremo.

1.2 Avvio di MATLAB e Octave

Tralasciamo la parte relativa alla procedura di installazione di MATLAB (e di Octave), che dipende dalcomputer e dal sistema operativo in uso, e diamo per scontato che il programma sia presente sul propriocomputer.

In ambiente Linux, per avviare MATLAB basta digitare il comando matlab da una finestra di termina-le (analogamente, basta digitare octave per entrare in ambiente Octave). Una volta avviato MATLAB , siaprir un ambiente che contiene molte finestre (la Command Window, la Workspace, la Command Histo-ry, la Current Folder). Il prompt dei comandi dato dal simbolo >> (> in Octave) e comparir il promptdei comandi in una finestra come quella mostrata in Figura 1.1 (si veda Figura 1.2 per lanalogo di Octave).Per uscire dallambiente si digita exit o quit dalla finestra dei comandi oppure dal menu File o ancoradallicona di chiusura (x in alto a destra della finestra).

Per imparare a usare MATLAB , conviene prendere familiarit con la finestra dei comandi eseguendocalcoli come se fosse una calcolatrice. Vediamo che il risultato viene assegnato ad una variabile detta ans eche nella Workspace si trovano informazioni su di essa (si veda Figura 1.3).

Le principali operazioni aritmetiche sono date da

elevamento a potenza moltiplicazione/ divisione+ addizione sottrazione

Nel fare le operazioni aritmetiche, viene data la precedenza agli elevamenti a potenza, poi alle moltiplica-zioni e divisioni, infine alla sottrazione e addizione. Quando due operatori hanno la stessa prorit vengonoeseguite le operazioni partendo da sinistra e andando verso destra. Tuttavia, quando le espressioni sonoabbastanza complicate e c il rischio di non capire bene quali operazioni vanno fatte prima e quali dopo,conviene mettere sempre le parentesi.

Le funzioni di base (quelle che troviamo anche su una calcolatrice) sono date da

2

http://www.gnu.org/software/octave/http://www.gnu.org/software/octave/

1.2. Avvio di MATLAB e Octave

Figura 1.1: Lambiente MATLAB consiste di una serie di finestre, alcune delle la Current Folder (la di-rectory corrente in cui si sta lavorando), la Workspace (lo spazio di lavoro) e la Command History (lafinestra con la storia dei comandi dati).

Figura 1.2: Avvio di Octave: la finestra che si apre corrisponde alla Command Window di MATLAB .

3

1. INTRODUZIONE

Figura 1.3: Primi calcoli in MATLAB : osservare cosa succede nella Command Window, nella Workspacee nella Command History.

Funzione MATLAB (Octave) Significatosin(x) sin(x)cos(x) cos(x)tan(x) tan(x)asin(x) asi n(x)acos(x) acos(x)atan(x) at an(x)exp(x) ex

log(x) ln(x)log10(x) log(x)sqrt(x)

px

abs(x) |x|

La cosa nuova, da sottolineare, che ora queste operazioni si possono fare su matrici e vettori e che inMATLAB (Octave) si lavora soprattutto su matrici e vettori.

Per lavorare meglio, nel seguito introduciamo il concetto di variabile e, successivamente, le matrici e ivettori.

1.2.1 Le variabili

Abbiamo visto che, se scriviamo nella finestra di lavoro, sin(10) il risultato ans =-0.5440. Il ri-sultato della funzione viene assegnato ad una variabile che, di default, viene chiamata ans. Se cambiamo ilcalcolo, il risultato sar del tipo ans = cui associato il valore numerico del risultato. Tuttavia, noi possiamodecidere di assegnare il risultato di una certa operazione o di una funzione ad una variabile a cui noi stes-si assegniamo un nome ben preciso. Per esempio, vogliamo che il risultato di sin(10) sia assegnato ad unavariabile che chiamiamo b. Per fare ci, basta scrivere

>> b=sin(10)

4


Figura 1.4: Variabili nella Workspace

b =

-0.5440

In questo modo il valore della funzione non pi assegnato alla variabile ans ma alla variabile b. Anche10 pu diventare il valore di una variabile, in cui poi valutare la funzione seno.

>> a=10

a =

10

>> b=sin(a)

b =

-0.5440

Osserviamo che ora nella Workspace, troviamo il nome delle variabili introdotte e il loro valore (si vedafigura ??).

1.2.2 Prima introduzione (breve) su matrici e vettori in MATLAB

Ci che diremo ora su vettori e su matrici sar qualcosa di breve, (quando parleremo di vettori e matricinel corso di Calcolo Numerico, li rivedremo sotto unaltra luce) utile per alcune semplici e prime applicazionidi MATLAB (o Octave).

I vettori si possono scrivere come vettori riga o vettori colonna.Per esempio:

>>x=[1 2]>> x=[1 2]

x =

1 2

5

1. INTRODUZIONE

>> y=[10 20]

y =

1020

Un vettore rigax viene scritto usando le parentesi quadre e scrivendo le componenti del vettore una dopolaltra. Un vettore colonna pu essere creato facendo la trasposta di un vettore riga (mediante il simbolo ),oppure mettendo un punto e virgola dopo ciascuna componente:

>> y=[10;20]

y =

1020

In MATLAB i vettori altro non sono che un caso particolare di matrice a n righe e 1 colonna (vettorecolonna) o a 1 riga e n colonne (vettore riga). Quindi per scrivere una matrice si scrivono i valori della matriceriga per riga andando a capo con il punto e virgola:

>> A=[1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16]

A =

1 2 3 45 6 7 89 10 11 12

13 14 15 16

Per indicare un valore della matrice A basta specificare lindice di riga e di colonna: per esempio

>> A(2,2)

ans =

6

Per indicare gli elementi di tutta la colonna i si usa A(:,i), mentre A(i,:) indica gli elementi della rigai.

>> A(:,2)

ans =

26

1014

>> A(2,:)

ans =

5 6 7 8

6


Loperatore due punti pu dunque essere usato per estrarre unintera riga o colonna da una matrice ounintera sottomatrice. Se vogliamo estrarre le ultime due righe e colonne della matriceA, si digita il comando

>> M=A(3:4,3:4)

M =

11 1215 16

In questo modo abbiamo creato la matriceM che ha come elementi quelli della sottomatrice diA con le ultimedue righe e colonne.

Per scrivere una matrice vuota A, invece, si usa listruzione A=[ ].

1.2.3 Funzioni vettorizzate

Le funzioni che abbiamo visto prima si possono applicare direttamente su matrici e vettori, sono ciovettorizzate. Cosa significa questo? Significa che se abbiamo un vettore, ad esempio x = [1,2,3,4] e valutiamosin(x), il risultato sar un vettore che ha come componenti il valore della funzione nelle singole componentidel vettore, in questo caso [sin(1),sin(2),sin(3),sin(4)]. Stessa cosa vale per le matrici.

Vediamo un esempio

>> x=[2 4 6 8 10]

x =

2 4 6 8 10>> log(x)

ans =

0.6931 1.3863 1.7918 2.0794 2.3026>> A=[0.1 0.2 ; 0.3 0.4]

A =

0.1000 0.20000.3000 0.4000

>> sin(A)

ans =

0.0998 0.19870.2955 0.3894

Per quanto riguarda le operazioni aritmetiche, laddizione e la sottrazione sono vettorizzate per definizio-ne. Invece, le operazioni di moltiplicazione, divisione ed elevamento a potenza vengono vettorizzate facendoprecendere il simbolo di moltiplicazione, divisione e elevamento a potenza dal simbolo del punto (.*, ./,. ) permettendo, in tal modo, che le operazioni vengano fatte componente per componente delle matrici (odei vettori) su cui sono applicate.

>> x=[2 4 6 8 10]

x =

2 4 6 8 10

7

1. INTRODUZIONE

>> y=[1 3 5 7 9]

y =

1 3 5 7 9>> x.*y

ans =

2 12 30 56 90

Se facessimo semplicementex*y avremmo un messaggio di errore perch richiederemmo unoperazionedi prodotto tra due matrici che hanno n righe e 1 colonna ciascuna, e questa operazione non si pu fare(vedremo in seguito che si pu fare il prodotto tra due matrici A e B se il numero di colonne di A corrispondeal numero di righe di B). Vediamo un altro esempio

>> A=[ 1 2; 3 4]

A =

1 23 4

>> B=[2 4; 6 8]

B =

2 46 8

>> A*B

ans =

14 2030 44

>> A.*B

ans =

2 818 32

Questa volta, abbiamo considerato due semplici matrici A B . Se facciamo A*B, abbiamo come risulta-to la matrice prodotto (la studieremo pi avanti nel corso di Calcolo Numerico), se invece facciamo A.*B,abbiamo come risultato la matrice in cui ogni componente dato dal prodotto dei corrispondenti elementisulle due matrici di partenza.

1.2.4 Sulle variabili

Negli esempi precedenti, abbiamo dati a matrici e vettori un nome, x,y,A,B... Abbiamo quindiassegnato delle variabili.

Una variabile ha, quindi, un significato simile a quello che conosciamo gi nelle espressioni algebriche:nel nostro caso, per, ad una variabile sempre associato un valore numerico, se si tratta di una variabilescalare, oppure pi valori, se si tratta di vettori o di matrici. In figura 1.5 vediamo come nella Workspace,

8


Figura 1.5: Workspace con diversi tipi di variabili.

ciascuna variabile con cui si sta lavorando nella finestra dei comandi, sia elencata con il suo nome e con ilsuo valore (che il valore che le stato assegnato se si tratta di una variabile scalare, mentre dichiarato iltipo di variabile se di tipo matriciale). Dalla figura, possiamo osservare comeG ci sono due variabili, una chiamata A e una a: questo perch MATLAB (Octave) sensibile alle lettere

minuscole e maiuscole: una variabile chiamata A diversa da a!G il valore per A : vuol dire che si tratta di una matrice 10 10 i cui valori sono

memorizzati in doppia precisione (double)G la matrice B, invece, 1010 in singola precisione (single)G la matrice C 44 ed uint8 (che significa?, ora ci arriviamo)G abbiamo poi due scalari, a e b.G ci sono due vettori, uno riga e uno colonna, in doppia precisione

Ciascuna variabile, una volta che viene generata, di un certo tipo e viene classificata di conseguenza. Didefault, ogni variabile (scalare o matriciale) viene considerata in doppia precisione (quindi di tipo double).Mediante apposite functions di MATLAB (Octave) possibile convertire il valore di una variabile in un altrotipo: se scriviamo A=single(A), la matrice A viene convertita in una matrice che ha lo stesso nome, i cuivalori per sono ora memorizzati in singola precisione. La function uint8 converte una variabile (che puessere scalare o matriciale) in una variabile il cui valore ora un intero a 8 bit senza segno (con valori chevanno da 0 a 255). Se la variabile gi di tipo intero con valori che si trovano in questo intervallo, non ci sonoproblemi, ma se la variabile assume valori al di fuori di questo intervallo, allora i valori sono arrotondati alpi vicino valore dellintervallo ammissibile. Questo tipo di dato ha senso quando si lavora con numeri interie per fare certi tipi di operazioni. Vediamo esempi semplici

>> x=[0 1 2 3]

x =

0 1 2 3

>> uint8(x)

ans =

9

1. INTRODUZIONE

0 1 2 3

>> x=[0 -1 -2 3]

x =

0 -1 -2 3

>> x=uint8(x)

x =

0 0 0 3

Allinizio il vettore x ha tutti valori positivi e uint8 non cambia i valori del vettore (viene generato unvettore ans identico a x). Nel secondo caso il vettore x ha elementi negativi che vengono trasformati in 0.Questa volta il risultato della funcion uint8 assegnato alla stessa variabile x che viene quindi sovrascrittae il vettore non pi quello di prima.

Se usiamo le function double e single non osserviamo nessun cambiamento a prima vista, ma inrealt lo spazio di memoria per la variabile che cambia.

1.3 Comandi utili

Per lavorare meglio sia nella finestra dei comandi sia, successivamente, per scrivere ed eseguire dei veri epropri programmi, risultano utili le seguenti funzioni:G who fornisce lelenco di tutte le variabili presenti nella finestra dei comandi (lo si pu vedere anche

nella Workspace;G whos fornisce lelenco di tutte le variabili insieme allo spazio in memoria da esse occupato (con que-

sto comando si pu osservare come una stessa variabile, memorizzata comedouble o come single,occupa un diverso spazio di memoria);

G help pu essere usato da solo o insieme al nome di una function propria di MATLAB (Octave)o creata dallutente, e mostra tutte le informazioni utili per capire come usare MATLAB o quellafunction;

G clear se usata da sola cancella tutte le variabili presenti nella finestra dei comandi, se invece seguita da un elenco di variabili (messe una di seguito allaltra senza virgole) cancella tutte le variabilidellelenco;

G il punto e virgola ; messo dopo unistruzione non fa vedere il risultato dellistruzione nella finestradei comandi;

G il simbolo % il simbolo per scrivere commenti: ci che viene scritto dopo % rappresenta uncommento;

G clc ha la funzione di pulire la Command Window e di portare il cursore di nuovo in alto a sinistra nellafinestra (la dove era nel momento in cui abbiamo avviato MATLAB (Octave). Le variabili non vengonocancellate.

G diary permette di salvare, su un file che viene chiamato diary, il contenuto di ci che viene scrittonella finestra dei comandi. Il file diary si trova nella directory corrente in cui si sta lavorando;

G diary off chiude il file aperto mediante listruzione diary;G diary filediary comando analogo a diary, ma il file non si chiamer diary bens

filediary, dove filediary rappresenta il nome del file definito dallutente;G save filesave salva tutte le variabili presenti nel Workspace, nel file filesave.mat con

filesave nome scelto dallutente, in modo da conservare tutte le variabili con cui si stava lavorando.Questo comando utile se si deve chiudere la sessione di lavoro ma si vuole poi ripartire a lavorareproprio da quelle variabili;

10

1.4. Introduzione alla programmazione

G load filesave ripristina lo stato del Workspace, caricando tutte le variabili che erano statesalvate in filesave.mat in modo da poterci lavorare di nuovo. Si era chiusa la sessione di lavoro eora si pu ripartire a lavorare con le stesse variabili di prima.

Ce una notevole differenza di significato tra i files generati dallistruzione diary e quelli generati me-diante save. Ci che salviamo con il comando diary paragonabile a ci che pu essere scritto alla lava-gna e di cui prendiamo appunti: scriviamo un file di testo che possiamo elaborare o da cui possiamo trarredei risultati per i nostri problemi. Al contrario, quello che viene salvato mediante il comando save costi-tuito da una o pi variabili, che possono essere riutilizzate nuovamente in MATLAB (Octave) nella finestradei comandi, in sessioni successive a quella in cui sono state create, senza dover rieseguire tutti i comandiche hanno portato alla creazione delle variabili stesse. Infatti, quando si chiude MATLAB (Octave), tutte levariabili e i dati vengono persi!

Il comando save utile anche per salvare singole matrici o una lista di matrici nella directory corrente.Ci soffermiamo sulla possibilit di salvare i dati in formato ASCII, in modo da poter utilizzare i dati che sal-viamo non solo in MATLAB (Octave) ma anche in programmi scritti con altri linguaggi di programmazione.Supponiamo di avere una matrice A e di volerla salvare nel file matriceA.dat. Scriveremo il comando

save matriceA.dat A -asciiQuando vogliamo caricare in MATLAB la matrice dal file useremo il comandoload(matriceA.dat)e in questo modo avremo la matrice nella variabile matriceA (il nome del file senza lestensione dopo il

punto).

1.4 Introduzione alla programmazione

MATLAB (e di conseguenza Octave, dora in poi quando parleremo di MATLAB ometteremo di scrivereanche Octave, intendendo che il discorso vale anche per Octave e sottolineeremo solo le differenze impor-tanti) si pu intendere come un linguaggio di programmazione, anche se non rientra in nessuna delle duecategorie di base dei linguaggi di programmazione.

Ci sono infatti due categorie di linguaggi per scrivere programmi:G linguaggi di alto livello (come FORTRAN, C, C++)G linguaggi di basso livello (come assembler).

Il linguaggio macchina, invece, il linguaggio specifico dellelaboratore e consiste di cifre binarie 0 e 1 cheidentificano istruzioni elementari. Ogni elaboratore ha il suo linguaggio macchina.

I linguaggi di basso e alto livello si possono tradurre in un qualsiasi linguaggio macchina e possono essereutilizzati su qualsiasi elaboratore (sono portabili). Con i linguaggi di basso livello, i codici binari (fatti da0 e 1) vengono sostituiti da opportune parole chiave che permettono una corrispondenza uno-uno tra illinguaggio di basso livello e il linguaggio macchina.

Con i linguaggi di alto livello, una singola istruzione pu tradurre pi istruzioni di linguaggio macchina.Inoltre un linguaggio di alto livello facile da capire (nella sua sintassi, nelle sue regole, nel suo modo diessere utilizzato) da chi programma.

Perci, quando si programma con un linguaggio di alto livello, il programmatore deve conoscere e sa-pere bene la grammatica, la sintassi, il vocabolario del linguaggio stesso. Uno volta scritto il programmain modo corretto, occorre fare capire al computer che c un programma da eseguire: bisogna quindi tra-durre in linguaggio macchina il programma scritto con il linguaggio di alto livello: questo passaggio lo si famediante la compilazione: un comando particolare chiama un programma di compilazione che viene appli-cato al programma scritto nel linguaggio di alto livello (che chiamamo programma sorgente) in un program-ma equivalente scritto in linguaggio macchina e che potr essere eseguito dal calcolatore, detto programmaeseguibile.

In MATLAB , noi scriveremo qualcosa di equivalente al programma sorgente ma non lo compileremoperch ogni volta che lo vorremo eseguire, esso verr automaticamente compilato ed eseguito.

Abbiamo due tipi di files eseguibili:G lo script, che presenta una lista di comandi da eseguire in ambiente MATLAB . Al suo interno ci

pu essere la chiamata a functions proprie di MATLAB o a functions definite dallutente;

11

1. INTRODUZIONE

G le functions (di MATLAB o definite dallutente), che possono avere uno o pi variabili di ingresso(dati di input) e una o pi variabili in uscita (dati di output).

Sia lo script sia le functions si scrivono su files chiamati M-files perch sono files che vanno salvaticon lestensione .m (ad esempio prova.m, fun.m, . . . ). I files vanno scritti con leditor di testo di MA-TLAB (nella barra dei comandi di MATLAB in alto a sinistra compare la scritta New Script o New) o conun altro editor di testo in formato ASCII (in Octave non c un editor e quindi i files vanno scritti con un altroeditor di testo come gedit, ad esempio).

Per eseguire uno script salvato in nomefile.m basta scrivere il nome del programma (senza lesten-sione .m) nella finestra dei comandi o mediante il comando run(nomefile.m) o run(nomefile)oppure posizionandosi con il mouse sul nome del programma nella lista che compare nella finestra di Cur-rent Folder e con il tasto destro del moouse andare sulla scritta run che compare sulla finestra che si apre (pidifficile a scriversi che a farsi!).

Ovviamente, sia lo script sia le function definite dallutente devono trovarsi nella directory in cui si aperto MATLAB , altrimenti bisogna spostarsi nella directory con i files che servono (da MATLAB lo si fadirettamente o tramite la Current Folder o utilizzando la riga che mostra la directory in cui ci si trova o con icomandi dellambiente Linux, vale a dire cd e il percorso da fare, invece, in Octave, ci si sposta di directorysolo con i comandi dellambiente Linux).

Per quanto riguarda le function, possono essere eseguite o tramite uno script oppure direttamentenella finestra dei comandi, specificando le variabili di ingresso e di uscita.

1.5 Problemi e Algoritmi

Gli script e functions che scriveremo ed eseguiremo in ambiente MATLAB serviranno a risolveredei problemi.

Cosa intendiamo per problema? Se vediamo alletimologia del nome, esso deriva dal verbo grevo pro-ballein, dove pro significa davanti e ballein significa mettere: mettere davanti... problema ha il significatodi ostacolo, promontorio, impedimento. E difatti, se si ha un problema, la prima cosa che si pensa comerisolvere il problema (come superare lostacolo!).

Per risolvere un problema, bisogna conoscerlo e saperlo descrivere bene. A questo scopo, noi ci serviamodel linguaggio. Ma il linguaggio naturale va usato con cautela: lo stesso evento descritto, infatti, in mododiverso da un bambino rispetto ad un adulto; e lo stesso evento viene descritto con un linguaggio diverso daun matematico, da un fisico, da un ingegnere e da un letterato.

Accanto al linguaggio naturale c il linguaggio artificiale, che si presenta sotto forma di termini tecnici enotazioni (pensiamo a termini tecnici della fisica e della chimica o ai simboli matematici).

A volte, per, le stesse parole possono essere usate per dire cose completamente diverse! Basti pensa-re a termini come temperatura, pressione, massa...: hanno significati ben diversi se usati in unofficina perriparare il motore di unautomobile o se usati in un ambulatorio medico nei confronti di un paziente...

Una volta che, cono attenzione e precisione, si capito il problema, si pu passare alla fase di soluzione.Concentrando ora la nostra attenzione su problemi da risolvere con laiuto del calcolatore, dobbiamo farecapire al calcolatore quale il problema da risolvere e come lo deve risolvere. Ma siamo noi che dobbiamodare precise istruzioni da eseguire, quindi dobbiamo sapere quale tecnica di soluzione dobbiamo applicare.Queste precise istruzioni da dare al calcolatore rappresentano quello che chiamiamo algoritmo.

Un algoritmo una sequenza di passi che permettono di risolvere tutto o parte di un problema.1

Un algoritmo si pu paragonare ad una ricetta di cucina. La ricetta fatta di due partiG la lista degli ingredienti da utilizzareG la sequenza dei passi da fare per realizzare la ricetta.

1La parola algoritmo entrata in uso negli anni 50 in sostituzione di algorismo, termine con cui si indicava il processo di calcoloutilizzando i numeri arabi. Il termine algoritmo deriva dal nome di al-Khwarizmi, importante matematico arabo del nono secolo grazieal quale lEuropa impar ad usare i numeri arabi, oltre alla notazione in base 10. Le procedure che permettevano di effettuare calcoli innotazione decimale presero il nome di algorismi o algoritmi. Anche la parola algebra viene da questo matematico e, in particolare, dallasua opera Al-mukhtasar fi hisab al-jabr wa al-muqabala.

Nel medioevo (e forse anche per qualche studente di oggi!!!), si pensava che questa parola derivasse dal greco algiros (panico, dolore)e arithmos (numero).

12

1.6. MATLAB come linguaggio di programmazione

Nel realizzare una ricetta si possono avere problemi nel cercare alcuni ingredienti. Una parte dellesecuzionerichieder poco tempo, altre parti richiederanno pi tempo. Alcuni passi non richiedono che si segua uncerto ordine, altri passi richiedono che si mantenga lordine scritto sulla ricetta. . . Le stesse analogie troviamoin un algoritmo.

1.6 MATLAB come linguaggio di programmazione

Per imparare a usare MATLAB come uno strumento di programmazione, dobbiamo imparare ad usareil linguaggio che ci permetter di scrivere script e function. E, come tutti i linguaggi, dovremo impararealcune regole ben precise.

Un linguaggio di programmazione costituito da dichiarazioni (statements):G descrizione dei dati: che genere di dati abbiamo? Se vogliamo scrivere un programma in cui per

ogni studente calcoliamo la media dei suoi voti, c differenza tra i nomi degli studenti (una stringadi caratteri) e la loro media (un numero). Si parla di tipo di dati. In MATLAB non si dichiarano levariabili, come avviene invece per altri linguaggi di alto livello. Ci facilita la programmazione ma ponepi rischi (se ci si dimentica il tipo di una variabile!). Conviene, perci, prestare sempre attenzione altipo di variabili utlizzate.

G strutture di controllo. Un programma pu essere visto come una sequenza di istruzioni per risol-vere un certo problema. Alcune istruzioni vanno ripetute, altre vanno eseguite solo in determinatecondizioni. . .

G lavorazione dei dati (data processing). In un programma occorre poter lavorare sui dati e poterlimanipolare opportunamente.

G dichiarazioni di ingresso e uscita (input e output). In genere, un programma scritto in modo che idati su cui si lavora possano esistere fuori del programma stesso.

1.7 Primo script e prima function in MATLAB

Con un editor di testo scriviamo finalmente il nostro primo script. Uno script un programma (unaserie di istruzioni) che si differenzia dai programmi scritti con altri linguaggi di programmazione perch nonsi devono dichiarare le variabili. Nel programma che scrivamo, introduciamo due vettori ma non dobbiamousare istruzioni che dichiarano il tipo di variabili (vettori, singola o doppia precisione o variabile di tipo intero,...).

% questa e una riga di commentox=[1 , 0.5 , 0 , 0 . 5 , 1]y=sin ( x )disp (ho eseguito il mio primo script) %disp e una function per la

% visualizzazione% in questo caso la str inga% ho eseguito i l mio primo s c r i p t

Gi da questo primo script impariamo cose nuove. Intanto per eseguirlo, dopo averlo salvato, scriviamoil nome con cui abbiamo salvato il file (in questo caso scegliamo il nome primoscript.m) sulla CommandWindow. Vediamo cosa succede

>> primoscript

x =

-1.0000 -0.5000 0 0.5000 1.0000

y =

-0.8415 -0.4794 0 0.4794 0.8415

13

1. INTRODUZIONE

ho eseguito il mio primo script

Nello script abbiamo assegnato a due variabili x e y il valore di due vettori riga di 5 componenti eabbiamo stampato una scritta sulla Command Window. Le due istruzioni x = ... ey= ... si chiamanoistruzioni di assegnazione. Abbiamo assegnato, in questo caso, i valori delle componenti del vettore x e delvettore y. Nel file, ci che segue il simbolo % una riga di commento. La function disp una functiondi MATLAB che serve a visualizzare il contenuto di una variabile che pu essere uno scalare, una matrice,ma anche una stringa di caratteri, come in questo caso (una stringa di caratteri pu essere considerata asua volta una variabile che, per, non numerica). Modifichiamo il nostro script mettendo un ; dopoogni istruzione. Se eseguiamo di nuovo lo script, i due vettori non vengono visualizzati e compare solo ilmessaggio di esecuzione del primo script. I due vettori, tuttavia, ci sono (basti vedere nella Workspace) e sesi digita x e si schiaccia il tasto di Invio nella Command Window, ci viene visualizzata la variabile x. Stessacosa vale pery. Si provi a vedere invece cosa succede se facciamodisp(x) edisp(y). Quando si eseguonoprogrammi pi lunghi, conviene usare il ; dopo ogni istruzione di assegnazione, in modo da non vedere valoriche scorrono e che comunque non riusciamo a leggere (si pensi a matrici di dimensioni elevate). In tal modosi velocizzano i tempi di esecuzione dello script.

E se vogliamo cambiare il vettore iniziale x, come facciamo? Usiamo listruzione input. Modifichiamolo script

% questa e una riga di commentox=input ( scrivi il vettore x ) ;y=sin ( x ) ;disp (ho eseguito il mio primo script ) ; %disp e una function per la

% visualizzazione% in questo caso la str inga% ho eseguito i l mio primo s c r i p t

Quando eseguiamo questo script, sulla Command Window compare la frase che abbiamo scritto tra apici enoi dovremo scrivere il vettore che vogliamo (ad esempio, [ 1 2 3]) e schiacciare il tasto di invio.

>> primoscriptscrivi il vettore x [1 2 3]

ho eseguito il mio primo script

Proviamo ora a voler scrivere un programma che, oltre a darci il valore di sin(x) ci dia il segno della funzio-ne (1 se positivo, 0 se uguale a zero, 1 se negativo). Per fare ci noi possiamo usare la funzionesign, propriadi MATLAB . Ma siccome vogliamo imparare a scrivere noi una function, la scriviamo noi: una function chedia come output sia il valore di sin(x) sia il suo segno.

Scriviamo dunque un file, che salveremo con il nome mysin.m, in questo modo

function [ y , b]=mysin ( x )% function [ y , b]=mysin ( x )% input : x ( scalare , vettore , matrice )% output : y=sin ( x )% b= sign ( y )y=sin ( x ) ;b=sign ( y ) ;end

La prima riga fondamentale per dire che stiamo scrivendo una function: la regola da rispettare la seguente

function [variabili output] = nomefunction(variabili input)

Se le variabili di output si riducono ad una sola, possiamo evitare di mettere le due parentesi quadre. Levariabili di input e di output possono essere una o pi di una (e a volte anche nessuna). Subito dopo questa

14

1.8. I predicati elementari

istruzione, conviene mettere delle righe di commento per dire cosa fa questa function. Questo perch, facen-do help nomefunction sulla Command Window, avremo disponibili tutte le informazioni che abbiamoscritto in quelle righe di commento e ci sar utile per ricordare cosa fa quella determinata funzione.

Una function termina con listruzione end.Una volta salvata, la funzione pu essere chiamata direttamente nella Command Window, oppure allin-

terno di uno script. Se scriviamo semplicemente y=mysin(x), la function restituisce solo il primo vet-tore che viene generato. Se invece scriviamo [y,b]=mysin(x), allora ci vengono restituiti i due vettori.Possiamo usare anche nomi diversi rispetto a y b.

La function pu essere richiamata anche allinterno di uno script. Modifichiamo lo scriptprimoscriptin questo modo:

% questa e una riga di commentox=input ( scrivi il vettore x ) ;[ y , b]=mysin ( x )disp (ho eseguito il mio secondo script)

Osserviamo che, in questo caso, non abbiamo messo il ; dopo la chiamata alla function, in modo davisualizzare subito il risultato.

Dagli esempi che abbiamo visto, possiamo chiederci quando serve scrivere uno script piuttosto di unafunction. Vediamo, intanto, le principali differenze.G Le variabili usate allinterno di uno script non vengono cancellate a meno che non venga fatto tra-

mite la funzione clear. In una function tutte le variabili usate al suo interno e che non compaiononella lista delle variabili di output vengono cancellate una volta che la function viene eseguita.

G Una function pu essere usata allinterno di uno script, per snellirne la sua struttura (se determi-nate istruzioni vanno ripetute pi e pi volte, la function ci permette di scriverle una volta sola e dirichiamarle con parametri di input diversi). In tal caso, la function agisce da sottoprogramma.

1.8 I predicati elementari

Le istruzioni che traducono algoritmi di soluzione di problemi numerici non sono sempre istruzionidi tipo sequenziale, cio istruzioni da eseguire luna dopo laltra (fai questo, poi questo, poi questo anco-ra...esempioponi a=10, b=sin(a), c=a*b...). Molte volte alcune istruzioni vanno eseguite se sonovere determinate condizioni piuttosto che altre, o fino a quando vera una determinata condizione...

Esempio

Esempio 1.8.1 Per calcolare le radici di unequazione di secondo grado, ax2 +bx + c = 0 sappiamo chedobbiamo calcolare il discriminante = b4 4ac e se = 0 le due radici coincidono, se > 0 abbiamodue radici distinte, se invece < 0 le due radici non sono reali.Abbiamo una condizione su e, a seconda della condizione, facciamo una cosa piuttosto di un altra.

I simboli =, >, < sono degli operatori logici. La proposizione < 0 pu essere vera o falsa a seconda delvalore assunto dalla variabile. Gli operatori logici, quindi, ci permettono di costruire dei predicati. Vediamoquali sono gli operatori logici che useremo

15

1. INTRODUZIONE

Operatore Significato Esempio Valore> maggiore (a >b) Vero se a > b

Falso se a b>= maggiore o uguale (a>=b) Vero se a b

Falso se a < b< minore (a a=12;>> b=5;>> (a> (a==b)

ans =

0

>> (a>=b)

ans =

1

Il predicato esso stesso una variabile (si osservi che il risultato associato alla variabile ans): vale 1 se vero, 0 se falso.

Un operatore logico pu essere applicato anche a matrici e vettori: se A e B sono due matrici, il predicato(A>B) restituisce una matrice i cui elementi sono il risultato del confronto tra ciascun elemento della matriceA con il corrispondente elemento della B.

>> A=[1 2; 8 9]

A =

1 28 9

>> B=[ 1 4; 5 10]

B =

1 45 10

16

1.8. I predicati elementari

>> A>=B

ans =

1 01 0

In molti algoritmi, tuttavia, non basta vedere se un singolo predicato sia vero o falso. A volte ci servonopi predicati da mettere insieme per formare una proposizione logica.

Esempio

Esempio 1.8.2 Supponiamo di dover implementare un algoritmo in cui vanno eseguite alcune istru-zioni fintantoch vero che (1) il numero di iterazioni, it, minore di un numero massimo prefissato,itmax, e (2) il valore assoluto della differenza di due valori che sono generati dallalgoritmo (chiamia-mo questa quantit sc) maggiore di una tolleranza prefissata, tol. Abbiamo una proposizione fattadi due predicati (it tol). I due predicati sono messi insieme tramite la con-giunzione e. Fintantoch vera tutta la proposizione (quindi sono veri tutti e due i predicati) allora cisaranno da eseguire determinate istruzioni, altrimenti non si eseguiranno. Quando la proposizione non pi vera? Diventa falsa quando il primo o il secondo predicato diventano falsi, cio quando le itera-zioni diventano maggiori del massimo prefissato oppure quando la variabile sc diventa minore dellatolleranza.

Vediamo allora gli operatori logici per unire due proposizioni:

Operatore Significato& congiunzione di pi predicati| disgiunzione di pi predicati negazione di un predicato o di una proposizione

Esempio

Esempio 1.8.3 Riprendiamo i predicati di prima it < itmax e sc > tol. La proposizione (ittol) un proposizione che vera quando entrambi i predicati sono veri,diventa falsa quando uno dei due predicati falso.La proposizione (it tol) vera quando uno dei due predicati vero (pu es-sere uno vero e uno falso, la proposizione comunque vera perch | significa o, oppure). La proposizione falsa quando entrambi i predicati sono falsi.La proposizione (it = itmax. Quindi vera se it >= itmax e sar falsa se it

1. INTRODUZIONE

1.9 Strutture e cicli

Vediamo, in generale, adesso come si possono eseguire alcuni tipi di istruzioni.

1.9.1 Ciclo if

Supponiamo di voler calcolare la radice quadrata positiva di un numero reale : +px. Sappiamo che selargomento x negativo, la funzione non d un valore reale (si passa infatti nel campo dei numeri complessi:provare la function di MATLABsqrt con un numero negativo e vedere il risultato che d). Proviamo alloraa scrivere una nostra function che calcoli la radice quadrata solo per valori positivi dellargomento.

function [ y ]= mysqrt ( x )% function [ y ]= mysqrt ( x )% input x% output y= s q r t ( x ) per x >=0% altrimenti s i ha un messaggio di e r r o r e% e y = [ ] v a r i a b i l e vuotai f ( x >=0)

y=sqrt ( x ) ;else

disp (largomento e negativo)y = [ ] ;

end

Abbiamo usato il cosiddetto ciclo if. Se (if) (x>0) allora fai una determinata cosa, altrimenti (else)fai unaltra cosa. Osserviamo che per visualizzare una stringa di caratteri in cui c lapostrofo o laccento,abbiamo usato due volte gli apici.

Il ciclo if si pu schematizzare nei seguenti modi a seconda che si abbia una o pi condizione.

Ciclo if con una condizione

i f ( espr . logica ){ istruzione 1a }{ istruzione 2a }{ . . . . }

else{ is truzione 1b }{ istruzione 2b }{ . . . . }

end

Ciclo if con pi condizioni

i f ( espr . logica1 ){ istruzione 1a }{ istruzione 2a }{ . . . . }

e l s e i f ( espr . logica2 ){ istruzione 1b }{ istruzione 2b }{ . . . . }

. . . .else

{ is truzione 1z }{ istruzione 2z }{ . . . . }

end

Osserviamo che in MATLABelseif va scritto tutto attaccato altrimenti si apre un altro ciclo if concatenatoal precedente.

Se non c nessuna istruzione in alternativa alle condizioni poste, si omette lelse.Scriviamo uno script che risolve unequazione di secondo grado ax2 +bx + c = 0. Usiamo un ciclo if per

controllare il segno del discriminante =p

b2 4ac e, di conseguenza, scrivere le radici.a=input (valore di a ) ;b=input ( valore di b ) ;c=input ( valore di c ) ;delta=sqrt (b^2 4*a* c ) ;

18

1.9. Strutture e cicli

i f ( delta ==0)x1=b/(2* a ) ;x2=x1 ;disp (le radici sono coincidenti)disp ( [ x1 x2 ] )

e l s e i f ( delta >0)x1= (b + delta ) / ( 2 * a ) ;x2= (b delta ) / ( 2 * a ) ;

disp (le radici sono)disp ( [ x1 x2 ] )

elsedisp (non ci sono radici reali)

end

Salviamo lo script con il nome radici.m. Per visualizzare su una sola riga il valore delle due radici abbia-mo creato il vettore che ha come componenti x1 e x2 e lo abbiamo visualizzato mediante la function disp.Osserviamo che, quando eseguiamo questo script, dobbiamo dare noi in input i valori dei coefficienti delle-quaazione. Tutte le variabili che vengono generate nello script rimangono nella Workspace (a meno che nondecidiamo noi di cancellarle usando il comando clean). Per capire meglio la differenza tra script e function,scriviamo ora una function che fa la stessa cosa. Chiamiamo la function funradicieq2.m.

function [ x1 , x2 ]= funradicieq2 ( a , b , c )% function [ x1 , x2 ]= radicieq2 ( a , b , c )% input : a , b , c v a r i a b i l i s c a l a r i del l equazione di secondo grado% ax^2+bx+c=0% output : x1 e x2 , l e due r a d i c i del l equazionedelta=sqrt (b^2 4*a* c ) ;i f delta ==0

x1=b/(2* a ) ;x2=x1 ;

e l s e i f delta >0x1= (b + delta ) / ( 2 * a ) ;x2= (b delta ) / ( 2 * a ) ;else

disp (non ci sono radici reali)x1 = [ ] ;x2 = [ ] ;

endend

I comandi sono quasi gli stessi dello script. Per avere sempre dei valori in uscita, abbiamo posto x1, x2uguali a due variabili vuote (luso delle parentesi quadre senza nessun valore allinterno []) insieme ad unmessaggio di spiegazione, per discriminanti negativi. Questa volta possiamo eseguire la function allinternodella Command Window o mediante uno script.

>> [x1, x2]=funradicieq2(1,-2,1)

x1 =

1

x2 =

1

Tutte le variabili che sono usate allinterno della function non sono salvate nella Workspace, ad eccezionedelle variabili di output.

Sia nello script che nella function abbiamo visto un esempio di uso del ciclo if.

19

1. INTRODUZIONE

Proviamo a complicare le cose, inserendo altri cicli if allinterno del ciclo if.Nello studiare la propagazione degli errori, abbiamo visto che, per effetto della cancellazione numerica,

le due radici di unequazione di secondo grado possono non essere arrotondate correttamente se b2 >> 4ac.In tal caso, applicando la formula x1x2 = c

apossiamo evitare il fenomeno della cancellazione. Consideriamo

due casi, assumendo x1 = b +p

b2 4ac2a

e x2 = b p

b2 4ac2a

:

G b > 0: per effetto della cancellazione numerica x1 affetta da errore. In tal caso dobbiamo correggerex1 mediante la formula x1 = c/(x2a).

G b < 0: ora x2 affetto da errore (b positivo e va sottratto ad una quantit molto vicina a b stesso).Correggiamo x2 mediante la formula x2 = c/(x1a).

Andiamo a correggere la function (stessa cosa si pu fare con lo script precedente).

function [ x1 , x2 ]= funradicieq2nocanc ( a , b , c )% function [ x1 , x2 ]= radicieq2 ( a , b , c )% input : a , b , c v a r i a b i l i s c a l a r i del l equazione di secondo grado% ax^2+bx+c=0% output : x1 e x2 , l e due r a d i c i del l equazione% s i usa la formula per e v i t a r e la cancellazione numericadelta=sqrt (b^2 4*a* c ) ;i f delta ==0

x1=b/(2* a ) ;x2=x1 ;

e l s e i f delta >0i f b> [x1, x2]=funradicieq2(1,1e+6,0.1)

x1 =

-1.000007614493370e-07

x2 =

-9.999999999999000e+05

>> x1*x2

ans =

20


0.100000761449327

>> [x1, x2]=funradiceeq2nocanc(1,1e+6,0.1)

x1 =

-1.000000000000100e-07

x2 =

-9.999999999999000e+05

>> x1*x2

ans =

0.100000000000000

Lerrore sul prodotto x1x2 dellordine di 7.6 107, ma lerrore sulla radice x1 comunque un erroremolto basso. Ci dovuto al fatto che le nostre variabili sono tutte in doppia precisione e MATLAB eseguele operazioni con precisione di macchina.

Se inseriamo i valori di a,b e c in singola precisione mediante la function single, si vede qualchedifferenza in pi. Proviamo:

>> a=single(1);>> b=single(1.e+6);>> c=single(0.1);>> [x1, x2]=funradicieq2(a,b,c)

x1 =

0

x2 =

-1000000

>> x1*x2

ans =

0

>> [x1, x2]=funradicieq2nocanc(a,b,c)

x1 =

-1.0000000e-07

x2 =

-1000000

>> x1*x2

21

1. INTRODUZIONE

ans =

0.1000000

In questo caso la prima function ci d x1 = 0 al posto di x1 =1 107. Lerrore pi alto rispetto al casoprecedente ma non cos appariscente.

Proviamo allora a creare una function che ci permetta di lavorare solo con un certo numero di cifre de-cimali, come se avessimo una macchina in virgola mobile normalizzata a n cifre decimali, con n stabilitodallutente (pi piccolo n, maggiore sar il fenomeno di cancellazione numerica). Attenzione, per: leoperazioni che eseguiremo con MATLAB saranno sempre eseguite in doppia precisione mentre il risultatodovr essere convertito nella precisione che vogliamo noi.

Per capire la function che scriveremo, partiamo con un esempio. Supponiamo di avere il numero x =12.987654. Vogliamo rappresentare questo numero in una macchina in virgola mobile normalizzata conn = 4 cifre decimali. Quindi passeremo al numero x = 1.2988 101.

Lalgoritmo che implementeremo si basa sui seguenti passiG considero il numero in valore assoluto, cio 12.987654;G valuto log10 del numero: in questo caso ottengo il valore 1.113530710205556, che, se viene arrotondato

allintero pi vicino per difetto, mi d il valore 1, cio proprio lesponente della base (101);G adesso divido il valore assoluto del numero per 101 in modo da avere la parte frazionaria del numero

stesso (eventualmente, come in questo caso, in forma normalizzata): 12.987654/101 = 1.2987654;G moltiplico questo valore per 10n , dove n il numero di cifre decimali che ho deciso di avere nella

rappresentazione finale, in questo caso 4: ottengo 1.2987654e +04 (cio 12987.654)G arrotondo il valore allintero pi vicino: ho 12988;G divido per 10n : ricavo 1.2988;G ora prendo questo valore, lo moltiplico per il segno del numero di partenza e per 101 (lesponente

che avevo trovato al secondo passo): arrivo al valore 12.988. Il numero non espresso in formanormalizzata, ma in sostanza il numero che andavo cercando (12.988 =1.2988 101).

Per completare la function ci serve qualcosa che ci permetta di arrotondare allintero pi vicino per difettoe di arrotondare allintero pi vicino. Fortunatamente, in MATLAB , esistono gi due function che fannoqueste cose:G la function floor arrotonda allintero pi vicino per difetto;G la function round arrotonda allintero pi vicino con la regola dellarrotondamento simmetrico.

Queste function (come altre simili ceil, fix per la cui spiegazione rimandiamo allhelp on line)possono essere applicate a variabili matriciali (vettori e scalari vengono visti come casi particolari di matrici).

>> floor(2.5)

ans =

2

>> round(2.5)

ans =

3

>> floor(1.1)

ans =

1

>> round(1.1)

22


ans =

1

Per avere il segno di una variabile, usiamo la function sign (che vale +1 o -1 a seconda di segno positivoo negativo). Ecco allora la function, che chiamiamo roundnorm

function y=roundnorm( x , n)% function y=roundnorm( x , n)% input x : puo e s s e r e scalare , v e t t o r e o matrice% n : numero di c i f r e decimali in un sistema normalizzato% pensiamo i l numero s c r i t t o come a . b_1b_2b_3 . . . b_n x 10^p% con a ~= 0% output y : arrotondamento di x in un sistema normalizzato a n c i f r e decimali% di default , in MATLAB i l numero non viene v i s u a l i z z a t o in virgola% mobile normalizzata .%%% Esempi% x=sin (3)=0.141120008059867% n=4 % 4 c i f r e decimali% y=roundnorm( x , n)= 0.141120000000000% ( in virgola mobila normalizzata , y= 1.4112 10^(1)%% x=exp(1)=2.718281828459046% y=roundnorm( x , n)= 2.718300000000000 (=2.7183 10^0)%% x=exp(10)=2.202646579480672 e+04% y=roundnorm( x , n)= 22026 (= 2.2026 10^4)%% s i possono avere c a s i p a r t i c o l a r i% x=55.981999999999999% in MATLAB y=roundnorm( x , n)= 55.981999999999999 (= 55.982)% in questo caso , i n f a t t i ,% f *10^e = 5.598200000000000*10=55.981999999999999% in Octave y=roundnorm( x , n)= 55.9820000000000

xx = abs ( x ) ; % s i considera i l valore assoluto di xe = f loor ( log10 ( xx ) ) ; %esponente di xf = xx . / 1 0 . ^ e ; %parte frazionaria di xf =round ( f *10^n)/10^n ;y=sign ( x ) . * f .*10^ e ;end

Proviamo a vedere cosa succede se applichiamo questa function per risolvere unequazione di secondogrado. Nelle due function precedenti sulla soluzione di unequazione di secondo grado, dobbiamo andarea convertire le variabili delta, x1 e x2 applicando la function roundnorm: aggiungeremo delle istru-zioni del tipodelta=roundnorm(delta,n); x1=roundnorm(x1,n); x2=roundnorm(x2,n);,

23

1. INTRODUZIONE

mentre lintero n viene dato in input.

Esercizio 1.9.1 Come esercizio si creino, dunque, due nuove function, una che risolve lequazione di se-condo grado senza andare a considerare la correzione del fenomeno di cancellazione numerica e laltra,invece, che corregge la cancellazione numerica, in modo che, in input, oltre ai parametri a,b,c vengachiamato anche lintero n. Allinterno delle function si aggiungano le istruzioni di chiamata alla functionroundnorm.Si trovino, quindi, le radici delle equazioniG 0.1x2 42x +0.3 = 0G x2 56x +1 = 0G x2 18x +0.03 = 0

Si confrontino i risultati tra le due function, lavorando con n = 4 cifre decimali e utilizzando anche lefunction che lavorano in doppia precisione.

1.9.2 Ciclo while

Unaltra struttura molto importante data dal ciclo while. Ci possono essere situazioni in cui determinateistruzioni sono da ripetere fintantoch vera una determinata condizione. Prendiamo in esame lo schemadi bisezione: ad ogni passo del metodo si prende il punto medio dellintervallo, si confronta il valore dellafunzione nel punto medio con il valore agli estremi e si prende, come nuovo sottointervallo, lintervallo cheha come estremi il punto medio e uno dei estremi dellintervallo precedente. Queste istruzioni vanno ripetutefino a quando lampiezza dellintervallo su cui si lavora diventa molto piccola (minore di una certa tolleranzaprefissata) oppure siamo cos fortunati che il punto medio proprio la radice della funzione, oppure nonstiamo arrivando a convergenza ma abbiamo fatto un numero elevato di iterazioni.

La struttura del cilclo while fatta nel seguente modo

Ciclo while

while ( espressione logica ){ istruzione 1 }{ istruzione 2 }{ . . . }{ is truzione n }

end

Le istruzioni 1, 2, . . . , n vengono ripetute ciclicamente e non una sola volta come accade nel ciclo if: infattiogni volta che si esegue lultima istruzione allinterno del ciclo while, si torna poi allespressione logica postasubito dopo la parola while e se vera si eseguono di nuovo le istruzioni 1, 2, . . . , n, se falsa, invece, si escedal ciclo. importante, dunque, avere unespressione logica che eviti di creare cicli infiniti, che farebberoentrare in stallo lesecuzione di uno script o di una function!

Vediamo un esempio di ciclo while, scrivendo una function che implementa lo schema delle bisezioni.Vedendo questa function, impareremo nuove potenzialit di MATLAB .

function [ c , i t e r , scarto ]= funbisez ( fun , a , b , t o l l , itmax )% function [ c , i t e r , scarto ]= funbisez ( fun , a , b , t o l l , itmax )% function che calcola una radice del la funzione fun d e f i n i t a% nell i n t e r v a l l o [ a , b ]% input : function fun : puo e s s e r e una function inline , una% function .m ( in questo caso i l nome del la function% va s c r i t t o tra apici ) oppure una function_handle% estremi del l i n t e r v a l l o a e b

24


% tolleranza t o l l% numero massimo di i t e r a z i o n i itmax% output : approssimazione del la radice data dalla v a r i a b i l e c% numero di i t e r a z i o n i e f f e t t u a t e i t e r% scarto f i n a l e%% Esempi con i n l i n e function e handle function%% f = i n l i n e ( 2 . * xcos ( x ) + 1 ) ;% a=6; b=4; t o l l =1. e10; itmax=100;% [ c , i t e r , scarto ]= funbisez ( f , a , b , t o l l , itmax )% r e s t i t u i s c e% c = 1.4552e11% i t e r =36% scarto =7.2760e11%% f = i n l i n e ( x.^2 sin (3x ) + 2 ) ;% a=2; b=1; t o l l =1. e10; itmax=100;% [ c , i t e r , scarto ]= funbisez ( f , a , b , t o l l , itmax )% r e s t i t u i s c e% estremi d e l l o s t e s s o segno% c = [ ]% i t e r = [ ]% scarto = [ ]%%% a = 0 . ; b=3; t o l l =1. e10; itmax=100;% [ c , i t e r , scarto ]= funbisez ( @cos , a , b , t o l l , itmax )% r e s t i t u i s c e% c =1.5708% i t e r =34% scarto =8.7311e11%aux= f ev a l ( fun , a ) * f ev a l ( fun , b ) ; %aux= v a r i a b i l e d appoggio per

% c o n t r o l l a r e i l segno% del la funzione a g l i estremi a e b%f e v a l e una function che permette di valutare% i l valore del la funzione fun nella v a r i a b i l e a% ( f e v a l ( fun , a ) ) oppure fun ( b ) in f e v a l ( fun , b ) )

i f aux>=0disp (estremi dello stesso segno)c =[ ] ;i t e r =[ ] ;scarto =[ ] ;return % l i s t r u z i o n e return interrompe l esecuzione del la function

% se g l i estremi sono d e l l o s t e s s o segno non s i puo % applicare lo schema d e l l e b i s e z i o n i

endi t e r =0; % i t e r e la v a r i a b i l e intera che conta l e i t e r a z i o n i che vengono

% e f f e t t u a t ec =(a+b ) * 0 . 5 ;scarto=abs (ba ) * 0 . 5 ; % v a r i a b i l e che permette di c o n t r o l l a r e di quanto c i

% s i sta avvicinando al l a radicewhile i t e r t o l l

i t e r = i t e r +1;aux= f ev a l ( fun , a ) * f ev a l ( fun , c ) ;

i f aux

1. INTRODUZIONE

b=c ;else

a=c ;endc =(a+b ) * 0 . 5 ;scarto=abs (ba ) * 0 . 5 ;

endi f ( i t e r >itmax )

disp (raggiunto il numero massimo di iterazioni)end

Oltre alluso del ciclo while (e di un ciclo if allinterno del ciclo while), la novit nella function funbisez dato dal fatto che, tra i parametri di input, compare la funzione di cui si vuole approssimare una radice. Inlinea di massima il nostro schema pu essere applicato su una generica funzione f (x) (funzione reale di unasola variabile reale), al fine di poter approssimare (se esiste) una sua radice. Possiamo scrivere la funzione delnostro problema utilizzando una function (quindi scrivere un file e salvarlo con lestensione .m). Esempio

function y=fun ( x )% function y=fun ( x )% function da usare nel metodo di b i s e z i o n i% function di esempioy =2.* xcos ( x)+1end

Questa function la salviamo come fun.m e la chiamiano tra i parametri di input della function funbiseztra apici [c, iter ,scarto]=funbisez(fun,a,b,toll,itmax) In tal modo, noi diamo il nome della function su cuiapplicheremo il metodo di bisezione. Potremmo chiamare la nostra function anche in modo diverso, nonfun.m ma, ad esempio, g.m: limportante mettere tra apici il nome della function. Daltra parte, non importante il nome dato alla variabile ma il significato di ciascuna variabile e lordine con cui compare.

corretto chiamare la function in questo modo [x, it ,sc]=funbisez(fprova,extra,extrb,toller,maxit) Cisignifica che fprova.m la function di cui vogliamo calcolare lo zero, extra e extrb sono gli estremidellintervallo, etc etc. Ogni variabile ha un nome diverso rispetto a quello che abbiamo scritto nella functionma il significato corrisponde.

Non corretto invece scrivere [c, iter ,scarto]=funbisez(a,b,fun,toll,itmax) perch ora abbiamo lavariabile a come prima variabile e non la funzione su cui applicare il metodo!

Ci sono anche altri modi per inserire una funzione tra i parametri di ingresso di una function. Se la fun-zione gi predefinita si pu scrivere il nome della function tra apici: per esempio, vogliamo calcolare gli zeridella funzione f (x) = cos(x): basta scrivere [c, iter ,scarto]=funbisez(cos,a,b,toll,itmax)

Se la funzione pi complicata si pu usare la cosiddetta inline function. Sulla Command Window (oallinterno di uno script che poi chiamer la funbisez) scriviamo il seguente comando

>> f=inline(2.*x-cos(x)+1)

f =

Inline function:f(x) = 2.*x-cos(x)+1

La function di MATLAB inline ci permette di creare una function che dipende da una o pi variabili.Esempi:

f=inline(cos(x)-2) (dipende da una sola variabile),f=inline(cos(x)+sin(y)) (dipende da due variabili).

Rispetto alle function .m, questa function caratterizzata da una sola istruzione (non possiamo costrui-re function che dipendono da cicli o da pi righe di istruzioni). Una volta che chiudiamo la session diMATLAB , perdiamo questa function (se ci serve in unaltra sessione dovremo riscriverla).

26


Un altro modo per passare una funzione tra i parametri di una function di usare una cosiddetta functionhandle: una function handle si costruisce mettendo il simbolo @ prima del nome della function gi definita(ad esempio @sin, @cos) oppure definita dallutente come M-file (ad esempio @myfun). Esempi:

[c, iter ,scarto]=funbisez(@cos,a,b,toll,itmax)[c, iter ,scarto]=funbisez(@myfun,a,b,toll,itmax) ,

dove myfun.m una function che abbiamo definito noi. In maniera equivalente, si pu assegnare ad unavariabile il nome di una function handle nella Command Window o in uno script. Esempio:

>> f=@fun

f =

@fun

>> [c,iter,scarto]=funbisez(f,a,b,toll,itmax);

Ora che abbiamo visto che una funzione pu essere chiamata tra i parametri di input di una function,vediamo come valutarla al suo interno Nellesempio che stiamo considerando, come valutare allinterno dellafunbisez la funzione fun (che pu essere un M-file, una function inline o una function handle)? Si usa laseguente istruzione (che possiamo leggere nella function funbisez) feval(fun,a) dovefun il nome che diamoalla funzione allinterno della function nellesempio, funbisez, mentre a la variabile in cui dobbiamovalutare lafun. Quindifeval una function di MATLAB che permette di valutare il valore di una function:corrisponde a fun(a). Se la funzione una funzione di pi variabili, vanno messe tutte le variabili da cuidipende: feval(fun,x1,x2 ,...)

Esercizio 1.9.2 Si applichi il metodo di bisezione per trovare gli zeri di funzione delle seguenti funzioni,imponendo una tolleranza pari a 1010 e fissando itmax=100:G f (x) = ln(x)+x2 x (ammette come radice = 1)G f (x) = e1x2 x sin(x) (ammette pi di una radice)

Se risolviamo lesercizio precedente, ci accorgiamo che non sempre facile trovare un intervallo in cui lafunzione assegnata abbia segno opposto agli estremi. utile quindi avere unidea del grafico della funzionestessa. In questo senso MATLAB ci permette di visualizzare facilmente una funzione tramite la functionezplot. La function ezplot pu avere uno o pi argomenti di input. Si pu dare solo la funzione (inline,m-file o function handle) da visualizzare, e il grafico viene fatto di default per x che varia in [6,6] (se definita per valori negativi) o in [0,6] se non definita per valori negativi. Possiamo cambiare gli estremidellintervallo, aggiungendo come variabile lintervallo stesso.

Esempiezplot(f) %f e una function inline o una function handle

ezplot(f) %f e una function Mfile definita dallutente o di MATLABezplot(f ,[ a,b]) %si visualizza il grafico della f per x in [a,b]

Per visualizzare tutte gli usi della function ezplot si controlli mediante lhelp sulla Commans Window.Del tutto simile a ezplot la function fplot in cui, tuttavia, va specificato lintervallo [a,b]. Vedremo,comunque, che ci sono altri modi per fare grafici.

Ora, seguendo le stesse linee della function funbisez, scriviamo una function che applichi lo schema delpunto fisso. Chiamiamo questa function pfisso0.m (capiremo poi perch abbiamo scritto quello zero nelnome della function).

Questa volta, per vedere se ci stiamo avvicinando a convergenza, controlliamo la differenza in valoreassoluto tra due approssimazioni successive. Quindi, la variabile scarto assume davvero il significato delnome (al contrario di quanto abbiamo fatto in funbisez). Nella function diamo anche una stima del fattoredi convergenza (o costante asintotica di convergenza) M tramite le variabili asint1 e asint2. Vediamoquesta function

27

1. INTRODUZIONE

function [ xkp1 , i t e r , scartokp1 , asint1 , asint2 ]= pfisso0 ( fun , dfun , x0 , t o l l , itmax )%function [ xkp1 , i t e r , scartokp1 , asint1 , asint2 ]= p f i s s o 0 ( fun , dfun , x0 , t o l l , itmax )% s i g n i f i c a t o d e l l e v a r i a b i l i% i t e r : i t e r a z i o n e del metodo del punto f i s s o% itmax : numero massimo di i t e r a z i o n i% t o l l : tol leranza p r e f i s s a t a per l approssimazione del punto f i s s o% x0 : punto i n i z i a l e del la successione% xk : approssimazione al passo k% xkp1 : approssimazione al passo k+1% scartok : scarto al l i t e r a t a precedente% scartokp1 : valore assoluto tra l i t e r a t a corrente e quella al% passo precedente% asint1 : scartokp1 / scartok approsimazione del f a t t o r e di convergenza M% asint2 : abs ( dfun ( xkp1 ) ) approssimazione del f a t t o r e di convergenza M% fun : funzione di cui s i vuole c a l c o l a r e un approssimazione del punto f i s s o% dfun : derivata del la funzione fun%% Esempio% x0 = 0 . 1 ; t o l l =1. e12; itmax=100;% dfun= i n l i n e ( sin ( x ) ) ;% [ xkp1 , i t e r , scartokp1 , asint1 , asint2 ]= p f i s s o 0 ( @cos , dfun , x0 , t o l l , itmax )% r e s t i t u i s c e% xkp1=0.7391% i t e r =70% scartokp1 =9.4613e13% asint1= 0.6735% asint2 =0.6736scartokp1 =2.0* t o l l ;scartok=scartokp1 ;i t e r = 0 ;xk=x0 ;while ( scartokp1 >= t o l l )&& ( i t e r itmax )

disp (raggiunto il numero massimo di iterazioni)endend

Con questa function implementiamo il metodo del punto fisso e conosciamo lapprossimazione finale delpunto fisso, il numero di iterazioni effettuate, lo scarto finale, una stima della costante asintotica. Con questafunction, tuttavia, non possiamo fare nessun grafico di convergenza perch perdiamo tutte le informazioniintermedie (i valori delle approssimazioni, degli scarti, etc, ad ogni iterazione). In effetti, per alcune variabilipossiamo farci stampare su un file il loro valore in modo da sapere cosa successo ad ogni passo di iterazione.Per altre, invece, utile metterle in un grafico.

Esaminiamo ciascuna variabile:G i valori x0,xk,xkp1 rappresentano i valori dellapprossimazione iniziale al passo 0, al passo prece-

dente iter e al nuovo passo di iterazione iter+1. Questi valori ci servono per calcolare lo scarto adogni iterazione e per approssimare il fattore di convergenza tramite la variabileasint2. Potremmo so-lo stampare questi valori ad ogni passo oppure salvarceli in un vettore (e calcolare la variabile asint2

28


come vettore applicando il valore assoluto della derivata della funzione di punto fisso al vettore delleapprossimazioni).

G gli scarti scartok e scartokp1 servono per calcolare asint1 ma soprattutto se vogliamo fare ungrafico di convergenza! Quindi ci servono i valori degli scarti ad ogni iterazione. Possiamo creareun vettore in cui la componente i -sima corrisponde allo scarto al passo i . In tal modo potremo siacalcolare asint1 in modo vettoriale, sia fare un grafico semilogaritmico di convergenza.

G la variabileiter un contatore scalare. Quando occorre, creeremo il vettore che ha come componentigli interi 1,2, . . . , i ter dove i ter rappresenter il numero finale delle iterazioni effettuate.

Modifichiamo quindi la function lavorando su vettori (la possibilit di scrivere alcuni risultati su file lalasciamo per il seguito). Scriviamo il seguente file pfisso.m.

function [ x , i t e r , scarto , asint1 , asint2 ]= pf isso ( fun , dfun , x0 , t o l l , itmax )%function [ x , i t e r , scarto , asint1 , asint2 ]= p f i s s o ( fun , dfun , x0 , t o l l , itmax )% s i g n i f i c a t o d e l l e v a r i a b i l i% i t e r : i t e r a z i o n e del metodo del punto f i s s o% itmax : numero massimo di i t e r a z i o n i% t o l l : tol leranza p r e f i s s a t a per l approssimazione del punto f i s s o% x0 : punto i n i z i a l e del la successione% x : v e t t o r e d e l l e approssimazioni del punto f i s s o% scarto : v e t t o r e d e g l i s c a r t i% asint1 : v e t t o r e che approssimazione i l f a t t o r e di convergenza M usando% i l v e t t o r e d e g l i s c a r t i% asint2 : v e t t o r e che approssima i l f a t t o r e di convergenza M usando% i l valore assoluto del la derivata prima del la funzione di punto% f i s s o% fun : funzione di cui s i vuole c a l c o l a r e un approssimazione del punto f i s s o% dfun : derivata del la funzione fun% Esempio% x0 = 0 . 1 ; t o l l =1. e12; itmax=100;% dfun= i n l i n e ( sin ( x ) ) ;% [ x , i t e r , scarto , asint1 , asint2 ]= p f i s s o ( @cos , dfun , x0 , t o l l , itmax )% r e s t i t u i s c e% i l v e t t o r e x di lunghezza 70% i t e r =70% i l v e t t o r e scarto% i l v e t t o r e asint1% i l v e t t o r e asint2% per f a r e i l g r a f i c o di convergenza semilogaritmico basta s c r i v e r e% semilogy ( [ 1 : 1 : i t e r ] , scarto )

scarto =2* t o l l ;i t e r = 1 ;x=x0 ;while ( scarto ( i t e r )>= t o l l )&& ( i t e r itmax )

disp (raggiunto il numero massimo di iterazioni)end

29

1. INTRODUZIONE

end

Osserviamo in che modo abbiamo creato i due vettori asint1 e asint2.

Per il vettore asint1 la generica componente i -sima deve essere uguale ascar toi

scar toi1. Supponiamo

che il vettore scarto abbia queste componenti: scarto=[ 2 0.90 0.45 0.31 0.20 0.14] do-ve scarto(1) corrisponde al valore che si ha alliterazione 0 (il valore che, nella function, abbiamo postouguale a 2*tol in modo da poter entrare nel ciclo while). Per approssimare il fattore di convergenza al pri-mo passo di iterazione dobbiamo fare il rapporto tra lo scarto alliterazione 1 e lo scarto alliterazione 0, cio0.90

2. Alla seconda iterazione dobbiamo fare il rapporto tra lo scarto corrente e lo scarto al passo precedente,

cio0.45

0.90e cos via. Nel vettore scarto abbiamo 6 componenti, la prima corrisponde alliterazione 0, ci vuol

dire che sono state effettuate 5 iterazioni. Dobbiamo fare dunque 5 rapporti per approssimare il fattore diconvergenza ad ogni iterazione:

0.90

2,

0.45

0.90,

0.31

0.45,

0.20

0.31,

0.14

0.20.

A tal fine costruiamo il vettore (a cui ora per comodit diamo il nome sc1) che ha come componen-ti quelle del vettore scarto dalla seconda componente in poi, cio sc1=[0.90 0.45 0.31 0.200.14]. Costruiamo poi il vettore (che chiamiamo sc2) che ha come componenti quelle del vettore scartodalla prima alla penultima componente, vale a dire sc2=[2 0.90 0.45 0.31 0.20]. Ora, se di-

vidiamo i due vettorisc1

sc2componente per componente, abbiamo proprio i rapporti che abbiamo scrit-

to prima. Questa operazione si pu fare senza dover introdurre nuovi vettori. Possiamo scrivere, infatti,asint1= scarto(2:iter )./scarto(1:iter1);

Per quanto riguarda asint2, abbiamo semplicemente valutato la funzione dfun nel vettore delleapprossimazioni x.

Poi, poich la prima componente che abbiamo dato al vettore scarto non ci interessa, la togliamoandando ad aggiornare il vettore scarto=scarto(2:iter);

Questa istruzione ci permette di prendere il vettore che ha come componenti quelle del vettore scartodalla seconda componente in poi e di chiamare questo vettore di nuovo con la stessa variabile.

Infine, aggiorniamo il numero delle iterazioni, togliendo ununit perch eravamo partiti dal contatoredelle iterazioni pari a uno in corrispondenza delliterazione 0.

A questo punto, in uscita abbiamo dei vettori.Per fare il grafico di convergenza dello schema, ci serve fare un grafico semilogaritmico dove, sullas-

se delle ascisse ci sono le iterazioni effettuate e sullasse delle ordinate i logaritmi (in base 10) degli scar-ti. La function di MATLAB che ci permette di fare grafici in scala semilogaritmica sullasse delle ordinateprende il nome di semilogy. Nel nostro caso, basta creare il vettore delle ordinate, che deve avere i valori1, 2, 3, . . . , i ter e applicare questa function. Per creare un vettore che parte dal valore 1 e arriva al valore itercon passo unitario, basta scrivere it =[1:1: iter ];

Perci possiamo usare il comando semilogy(it, scarto), oppure, in modo equivalente,semilogy([1:1:iter], scarto) (si veda figura 1.6).

Introduzione alla grafica

Abbiamo visto le function ezplot, fplot, semilogy. Vediamo ora le altre function di base perpoter fare un grafico. Abbiamo semilogx (analogo di semilogy, il grafico semilogaritmico sullassedelle ascisse), loglog (il grafico in scala logaritmica su entrambi gli assi) e infine plot.

Le function plot, semilogx, semiloy e loglog hanno le stesse caratteristiche. Vediamo quindigli elementi principali di come si usa la function plot. In maniera analoga si usano le altre function cheabbiamo detto.

plot(x,y) dove x e y sono dei vettori della stessa dimensione, produce il grafico in cui sullasse delle ascissec il vettore x e sullasse delle ordinate il vettore y.

Esempio:

30


Figura 1.6: Grafico generato dallistruzione semilogy(it, scarto). per lo schema di punto fisso applicato allafunzione di punto fisso g (x) = cos(x), con x0 = 0.1.

>> x=[0, 0.1 , 0.2, 0.3, 0.4, 0.5];>> y=sin(x);>> plot(x,y)

I comandi precedenti producono il grafico di figura 1.7Si pu anche usare semplicemente listruzione plot(y) che produce il grafico del vettore in funzione degli

indici delle componenti del vettore (quindi sullasse delle ascisse ci saranno gli indici 1,2, . . . ,n del vettore ymentre sullasse delle ordinate le componenti del vettore y di lunghezza n).

Si possono fare pi grafici sovrapposti in diversi modi. Uno di questi di avere pi coppie di vetto-ri x1,y1, x2, y2 e cos via e scrivere listruzione plot(x1,y1,x2,y2, ...) Un altro modo di sovrascrive-re i grafici luno sullaltro mediante il comando hold on Abbinato poi il comando hold off che chiude ilprocesso.

Vediamo degli esempi

>> x1=[0, 0.1 , 0.2, 0.3, 0.4, 0.5];>> y1=sin(x1);>> x2=[-0.5, -0.2 0.2 0.5];>> y2=cos(x2);>> plot(x1,y1,x2,y2)

Si faccia poi

>> plot(x1,y1) % il grafico precedente non ce piu ma si crea questo>> hold on>> plot(x2,y2)>> hold off

Per esercizio si provino i comandi precedenti. La figura generata nel primo caso mostra le due curve dicolore diverso, mentre nel secondo caso le curve hanno lo stesso colore.

31

1. INTRODUZIONE

Figura 1.7: Esempio di uso della function plot con x=[0, 0.1 , 0.2, 0.3, 0.4, 0.5] ey=sin(x).

Per cambiare colore e tipo di linea della curva, oppure fare un grafico per punti, si possono fare le mo-difiche direttamente dal menu che si presenta sulla finestra del grafico oppure tramite comandi in linea.Vediamo questi ultimi (comuni ad Octave).

Listruzione plot(x1,y1,s1,x2,y2,s2, ...) dove s1, s2, ... sono delle stringhe di caratteri (da scriveretra apici, quindi) permette di cambiare colore, linea o punto della curva.

Esempi: plot(x,y,g:) produce il grafico in cui la curva di colore verde (green) e la curva a puntini.plot(x,y,r-.) d come risultato una curva di colore rosso con linea fatta a tratteggio del tipo -.-.-..plot(x,y,yo) produce una curva dove i punti non sono legati tra loro mediante una linea ma per punti,

ciascuno dei quali ha la forma di un cerchietto (o) mentre il colore giallo (yellow)plot(x,y,c--*) combina la linea di tipo tratteggio con i punti a forma di asterisco * mentre il colore il

celeste (cyan). Per vedere tutti i colori, punti e linee a disposizione basta fare help plot dalla CommandWindow.

Altri comandi utili sonoG inserire un titolo al grafico, mediante la function title.

Esempio: title (grafico di prova)G inserire una descrizione sullasse delle x (e delle y) mediante la function xlabel (ylabel).

Esempio: xlabel(iterazioni), ylabel(scarti);G Inserire una legenda mediante la function legend.

Esempio: abbiamo fatto un grafico con listruzione plot(x,y1,x,y2) dove y1 corrisponde aduna certa funzione f mentre y2 corrisponde ad una certa funzione g . Con listruzione

legend(funz f, funz g) sul grafico compare la legenda che abbiamo scritto unita al tipo di lineao punti associati a quel vettore.

G inserire una griglia di fondo al grafico o toglierla, mediante il comando gridEsempio: grid on introduce linee di griglia; grid off le toglie.

Per ciaascuno di queste function fare delle prove e controllare lhelp in linea.Negli esempi fatti, tuttavia, abbiamo creato dei vettori di lunghezza molto ridotta. Se avevamo in mente

di fare il grafico della funzione sin(x) o cos(x), il risultato che abbiamo avuto stato molto povero perch

32


avevamo pochi punti a disposizione. Potevamo certamente usare le function ezplot o fplot per questoscopo, ma poich stiamo imparando ad usare MATLAB ci conviene prendere spunto da questa occasioneper imparare qualcosa in pi sui vettori.

Nel fare il grafico semilogaritmico, abbiamo visto un modo per creare un vettore, mediante listruzione[1:1: iter ].

Listruzione x=[x0: incr: xf] produce un vettore che ha come valore iniziale x0 e come valore finale xfe le altre componenti sono date a partire da x0 con incremento dato dal valore incr: x0, x0+incr,x0+2incr, ...., xf. Con questa istruzione possiamo creare dei vettori che hanno molte componentisenza doverle scrivere a mano. Per vedere la lunghezza del vettore, basta scrivere il comando size(x) (chepu essere applicato anche a matrici) : restituisce il numero di righe e di colonne della variabile x. Si puanche usare listruzione length(x) che restituisce, invece, la lunghezza del vettore.

>> x=[0:0.1:1]

x =

Columns 1 through 6

0 0.1000 0.2000 0.3000 0.4000 0.5000

Columns 7 through 11

0.6000 0.7000 0.8000 0.9000 1.0000

>> size(x)

ans =

1 11

>> length(x)

ans =

11

Un altro modo di creare vettori dato dalla function linspace, il cui uso il seguente: x=linspace(a,b);crea, di default, un vettore (in questo caso x) che ha 100 componenti,con valori equidistanti a partire da aper finire a b: quindi il vettore ha come prima componente a, poi a + (b-a)/99, ..., fino ad arrivare allacentesima componente che vale b. Si applica la formula xi = x1 + (i 1)h, per i = 2,3, . . . ,100 , dove x1 = a eh = b a

99(in tal modo: x2 = x1 +h; x3 = x1 +2h = x2 +h, fino ad arrivare a x100 = x1 +99h = a + (b a) = b.

Se vogliamo un vettore con un certo numero n di punti equidistanti tra a e b, si scrive listruzionex=linspace(a,b,n);

A questo punto, vediamo anche un altro modo per creare dei vettori, utilizzando un ciclo che nonabbiamo visto fino ad ora, e che prende il nome di ciclo for.

1.9.3 Ciclo for

Il ciclo for strutturato nel modo seguente

33

1. INTRODUZIONE

Ciclo for

for ind= v a l i n i z : incr : v a l f i n{ i s t r u z i o n i }end

Un esempio di ciclo for per creare un vettore di 11 elementi con componenti equidistanti tra 0 e 1 puessere il seguente

x ( 1 ) = 0 ;for i =2:11

x ( i )= x ( i 1)+0.1;end

Osserviamo che usare un ciclo for per creare vettori non conveniente n ottimale. Questo non vuol direche non lo useremo mai e che possiamo dimenticarci che esiste! Dobbiamo solo prestare molta attenzionequando lo usiamo e usarlo solo quando effettivamente serve!!!!

Vediamo un esempio, in cui misuriamo il tempo che impiega MATLAB per creare un vettore di lun-ghezza n con n grande, usando la forma compatta x=[x0:incr:xn], la function linspace e il ciclofor.

Per misurare il tempo che impiega ciascuno procedimento, applichiamo due functions di MATLABtice toc che sono da usare combinate insieme perch permettono di misurare il tempo per effettuare leistruzioni che si trovano nel mezzo.

Scriviamo uno script allo scopo di creare un vettore x di lunghezza n=100001 tra a=0 e a=1 applicando

i diversi approcci che abbiamo visto. La distanza tra due componenti del vettore data da h = 1100000

= 105.

cleart i cx = [ 0 : 1 . e5 : 1 ] ;tocdisp (lunghezza vettore)disp ( length ( x ) )clear xt i cx=linspace (0 ,1 ,100001) ;tocdisp (lunghezza vettore)disp ( length ( x ) )clear xt i cx ( 1 ) = 0 ;for i =2:100001

x ( i )= x ( i 1)+1.e5;endtocdisp (lunghezza vettore)disp ( length ( x ) )

Eseguiamo questo script (salviamolo prima dandogli un nome con lestesione .m, in questo casoconfrontotictoc.m). In MATLAB otteniamo

>> confrontotictocElapsed time is 0.000631 seconds.lunghezza vettore

100001

34


Elapsed time is 0.003618 seconds.lunghezza vettore

100001


100001

Vediamo cosa succede in ambiente Octave:

octave:3> confrontotictocElapsed time is 2.90871e-05 seconds.lunghezza vettore100001

Elapsed time is 0.001910925 seconds.lunghezza vettore100001


100001

Come possiamo osservare dai risultati, il primo approccio per scrivere i vettori quello che vadecisamente meglio.

1.9.4 Introduzione a file di scrittura dati

Riprendiamo a studiare lo schema di punto fisso e i risultati ottenuti. Immaginiamo di volerli stamparesu un file per poterli studiare ed esaminare con calma. La function pfisso ci ha dato dei vettori, noi abbia-mo fatto il grafico di convergenza. Ora vogliamo salvare, per ogni iterazione, i valori delle approssimazioni,degli scarti e delle stime asint1 e asint2 inserendoli in una tabella con una riga di intestazione in cuiscriviamo il significato della colonna corrispondente. Vediamo un modo compatto per fare tutto ci. Creia-mo una matrice in cui ciascuna riga corrisponde a ci che vogliamo mettere in tabella: la prima riga avr leiterazioni fatte, la seconda i valori approssimati, la terza gli scarti, la quarta le stime asint1 e la quinta lestime asint2. La matrice deve avere lo stesso numero di colonne, quindi dobbiamo stare attenti al fatto chealcuni vettori hanno come prima componente un valore che corrisponde alliterazione 0.

Per scrivere su un file di testo, dobbiamo aprire da MATLAB un file e dare ad esso un nome: ci vienefatto con la function fopen: questa function ha come parametri il nome del file su cui andremo a scriverei dati, scritto tra apici, e aggiungeremo poi dei permessi (ad esempio r per file di sola lettura, w per file discrittura,... si faccia help fopen). Loutput della function una variabile scalare di tipo intero, a cui diamonome fid (file identifier). Questa variabile sar usata come identificatore del file su cui dobbiamo andarea scrivere mediante la function fprintf. Una volta che si scritto, si chiude il file mediante la functionfclose, scrivendo semplicemnte fclose(fid).

La functionfprintfpermette di scrivere su un file utilizzando un certo formato. Ci che diciamo adessoserve in maniera analoga per la functionsprintf che permette di scrivere con un certo tipo di formato sullaCommand Window.

I dati di input di fprintf sono la variabile fid, alcune istruzioni che riguardano il formato, da scriveretra apici, e una o pi matrici. Se, si pone fid=1 loutput sar lo schermo della Command Window. Quindilistruzione da dare fprintf(fid, formato, A, B, ...) conA,Bmatrici o array (un array puessere un vettore o una matrice fatta di numeri o di stringhe di caratteri). Bisogna prestare attenzione al fattoche la scrittura avviene colonna per colonna, (vengono letti gli elementi della prima colonna, poi gli elementidella seconda colonna, ...).

Il formato una stringa di caratteri da scrivere tra apici che possono avere diversi significati: c il formatoassociato al tipo di rappresentazione della variabile (se in formato esponenziale, fisso, di tipo intero, di stringa

35

1. INTRODUZIONE

di caratteri...) che si scrive nella forma %e (formato esponenziale) o %12.6e (formato esponenziale cui sonoriservati 12 spazi di cui 6 per la mantissa)...; ci sono i caratteri che hanno il significato di riga bianca, andarea capo, ... che sono preceduti dal simbolo \ ( \n, nuova linea, \r andare a capo,...); ci possono essere stringhedi caratteri. In tabella vediamo i principali simboli che ci interessano.

formato Significato%s stringhe di caratteri%d formato intero%f formato decimale in virgola fissa%e formato esponenziale (del tipo 3.5e +00)%E formato esponenziale (del tipo 3.5E +00)%g formato che una via di mezzo tra %f e %e senza gli zeri che non servono\ n nuova linea (il risultato come schiacciare il tasto di invio)\ r per andare a capo\\ \

Il significato associato a \ n molto importante, perch permette di andare a capo ogni volta che stataletta una riga, altrimenti viene scritta solo una riga sul file, lunghissima a seconda del contenuto da scrivere!Al

dispensamatlab

Documents

Transcript of dispensamatlab