Category Archives: math

Siamo tutti un po matematici, no?

maxima, per la mate nella scuola

MaximalogoRecentemente su Facebook teh notorious proooof Zar ha chiesto su ‘na roba che poteva risultare ambigua come interpretazione.
Ho scoperto il post tardi, su sollecitazione dell’Autore, ho detto la mia:

  • scrivete chiaramente, se del caso usate le parentesi;
  • verificate con tools disponibili online, p.es. MathStudio e Wolfram|Alpha;
  • ma ho anche un’altra dritta, voglio però prima sentire altri, specie matematti e insegnanti.

Solo che non ci sono più stati contributi. Questo è il brutto del Web ultimamente (secondo me): poca interazione, colpa del mobile (secondo me).
Ora che sono passati tre giorni e sperando che non sia stato io a spoilare il post vi dico quel che volevo dire ancora 😉

Niente di nuovo qui su Ok, panico se ne è parlato ampiamente in passato:

Ah! lo so, ma OK: c’è l’installer pe Windows ed è –ovviamente– free. Non avete scuse. Trovate tutto qui e se volete date una scorsa alla Wiki.

m0

Uh! ‘n po’ di cose:

  • segalazione di espressione ambigua;
  • il due-punti non divide;
  • le parentesi sono tue amiche ma non sempre servono;
  • sì gli spazi si possono usare;
  • non dimenticate il punto-virgola, se del caso potete aggiungerlo sulla riga successiva.

Dopo le prime volte (ma vedi dopo) si può lanciare così:

m1

OK, ho semitrollato (sono vecchio), c’è l’usuale versione finestrosa; Maiusc-Invio per valutare l’espressione 🙂

m2

Poi –c’è chi ricorda tutto questo nel post di Zar– c’è SymPy (anche su OKp se ne è parlato) ma richiede la conoscenza di Python, qualcosa in più della semplice matematica che si fa a scuola (secondo me, vorrei essere smentito).

Uh! devo dire che l’immagine lassù è GPL, altrimenti mi katsiano 😳

:mrgreen:

Annunci

Maxima

logoMaxima is a computer algebra system dice il manuale che continua con
Maxima is derived from the Macsyma system, developed at MIT in the years 1968 through 1982 as part of Project MAC. MIT turned over a copy of the Macsyma source code to the Department of Energy in 1982; that version is now known as DOE Macsyma. A copy of DOE Macsyma was maintained by Professor William F. Schelter of the University of Texas from 1982 until his death in 2001. In 1998, Schelter obtained permission from the Department of Energy to release the DOE Macsyma source code under the GNU Public License, and in 2000 he initiated the Maxima project at SourceForge to maintain and develop DOE Macsyma, now called Maxima.

Nel Sotware center di Ubuntu ce ne sono 3 versioni, ognuna delle quali ha i suoi pregi.

umax
Un caso semplice:

Maxima REPL nel terminale
m0
Xmaxima
m1

in questo caso quit() non funziona, nemmeno l’usuale Ctrl-D, usare File|Exit.

wxMaxima
m2Per valutare l’espressione usare Shift-Enter. E File|Esci per uscire.

Quale versione usare? Dipende. Se voglio un output testuale una delle prime due, la terza se voglio immagini belle:

m3

I grafici vengono visualizzati attraverso Gnuplot, con le tre versioni ottengo lo stesso risultato.

m4
m5
m6E naturalmente…

plot3d (sin(sqrt(x^2 + y^2))/sqrt(x^2 + y^2),
    [x, -12, 12], [y, -12, 12])$

m7

e altrettanto ovviamente il grafico può essere esaminato, sapete Gnuplot 😀

 

m8

Insomma una volta scoperto Maxima non ne puoi più fare a meno. Grazie a shintakezou che me l’ha presentata.

Due cose ancora: sì, c’è anche per Windows e il sito di Maxima (è sempre lo stesso cambia solo l’interfaccia) dove trovate il download e tutta la documentazione, indispensabile.

:mrgreen:

Un’altra rapida e semplice incursione nel mondo di Maxima

Continuo la panoramica di Maxima, iniziata proprio con Maxima: panoramica iniziale.

Limiti

Naturalmente Maxima può anche calcolare limiti. Per esempio:

lim1

Questo è un limite notevole. In certi casi il limite “venendo da destra” è diverso dal limite “venendo da sinistra”: possiamo aggiungere un altro argomento a limit che specifica come deve essere fatto il limite, se da destra (plus) o da sinistra (minus).

lim2

Quando c’è questo tipo di ambiguità e non specifichiamo come vada fatto il limite, otteniamo und (undefined).

Naturalmente possiamo far tendere x all’infinito, inf, o meno infinito, minf.

lim3

Come risposta potete ottenere anche infinity (infinito complesso) o ind, indefinito ma limitato (come per esempio è il caso del seno e del coseno). Per qualche altra cosetta rimando al manuale.

Serie di Taylor/Laurent

A Maxima possiamo anche chiedere l’espansione (troncata) in serie di Taylor. Potete sperimentare con le espansioni più note, che ovviamente Maxima gestisce, o con altre cosucce come per esempio la seguente:

taylor1

La lettera γ è la costante di Eulero-Mascheroni. Esce fuori perché il fattoriale è stato espresso (diciamo per forza di cose… ma non dettaglio) in termini della funzione Γ, e infatti otteniamo lo stesso risultato se al posto di x! scriviamo gamma(x+1). Se calcoliamo la derivata dobbiamo ottenere, come potete leggere su Wikipedìa, ψ0(x + 1)Γ(x + 1). Se la calcoliamo a x=0, otteniamo -γ; Maxima “sa“ anche tutto questo:

taylor2

(Ho preso l’ultimo risultato tramite % e poi gli ho detto di considerare x=0.) Con ciò abbiamo sfiorato un paio di funzioni speciali, che naturalmente Maxima conosce.

Una volta che abbiamo l’espansione in serie (troncata) in un intorno di un punto, possiamo usarla per calcolare la funzione in quell’intorno; naturalmente otteniamo un’approssimazione. Ne approfittiamo per vedere che Maxima è anche un linguaggio di programmazione e può fare pure grafici (tipicamente è gnuplot in realtà a farli)…

Vogliamo avere un’idea di come migliora l’approssimazione della funzione gamma man mano che tronchiamo la sua espansione in serie di Taylor sempre più tardi (ma mi fermo presto perché già intorno al termine alla 11 comincia a metterci un po’ troppo… non aspettiamoci che Maxima abbia performance stratosferiche.)

Grafici (e Turing…)

Potremmo anche abbandonare l’ambiente interattivo, aprire un editor, e scrivere il codice come siamo abituati; questo magari sarà effettivamente più comodo in qualche puntata successiva (se ci arrivo!:-D). Per ora continuo ad usare wxmaxima.

L’obiettivo che mi sono posto alla fine del paragrafo precedente è in pratica questo: graficare un po’ di espansionsioni in serie (troncate) di Taylor della Γ, e la Γ medesima, per vedere cosa cambia man mano che aggiungiamo termini.

La Γ è definita per valori positivi della variabile; la funzione che stiamo sviluppando è in realtà Γ(x+1) (siamo partiti dal fattoriale), e la stiamo sviluppando intorno a x=0, per cui non ha senso andare più lontano di 1. Insomma il nostro intervallo sarà [0, 1].

Il codicillo che andiamo ad eseguire ha il seguente aspetto:

code0

(Ricorda: se abbiamo pasticciato un po’ e vogliamo pulire un po’ il nostro ambiente, eseguiamo reset()$ e kill(all)$.)

Abbiamo prima inizializzato una lista vuota; poi abbiamo usato un ciclo for per popolarla, tramite append. Il for parte da 2 e non supera 10 a passi di 3, cioè insomma calcoliamo taylor per 2, 5, 8. Alla fine appendiamo anche la gamma(x+1) che è la funzione “giusta e vera” e gli diciamo di fare un grafico 2D usando come intervallo per la variabile x [0, 1]. Dopo un po’ compare il nostro grafico.

Se vogliamo salvarlo invece che vederlo, possiamo scrivere qualcosa come:

plot2d(t, [x, 0, 1], [gnuplot_term, 'png], [gnuplot_out_file,"graph0.png"])$

Controllando un po’ meglio gnuplot è possibile ottenere grafici “migliori”; senza fronzoli il file generato è questo:

graph0

Se rimaniamo nell’intorno (positivo) del punto in cui abbiamo sviluppato la serie riusciamo a ottenere approssimazioni migliori man mano che aumentiamo i termini della serie; allontanandoci notiamo invece che le cose vanno diveramente e capita che il verde (che ha più termini) va persino peggio del blu (n=2)… Se giocate un po’ vi accorgete che lontano dallo zero avvengono oscillazioni piuttosto consistenti man mano che aggiungiamo termini. Del resto nessuno ci dice che abbia senso allontanarsi dallo zero…

Le possibilità grafiche di Maxima non si esauriscono qui, naturalmente, e nemmeno la capacità di esprimere algoritmi.

Vettori e matrici

Potevano mancare? Possiamo inserire matrici e calcolarne il determinante; come al solito possiamo mettere variabili (simboli) invece di numeri:

dett2x2

Se definiamo matrici più grandi possiamo dedurre la regola usata, che è proprio una di quelle che (forse) abbiamo imparato.

dett3x3

Provate da voi invert(m2), detout; per invertire la matrix 3×3. Noterete che al denominatore c’è proprio il determinante (detout serve per metterlo a fattore comune). Non a caso. Provate

invert(m2), a=0, b=0, c=0;

Quelle imposizioni rendono il determinante nullo, per cui Maxima vi risponde con un errore che può essere un po’ criptico (0 to a negative exponent), ma che dovrebbe farvi capire che c’è una divisione per 0 di mezzo o almeno questo: la matrice non è invertibile.

Possiamo calcolare il prodotto tra due matrici, e “scopriamo” che il prodotto matriciale non commuta:

AB-BA

Eccetera: anche qui ciò che possiamo fare è veramente tanto e potete esplorare e sperimentare con il capitolo del manuale Matrices and linear algebra. Magari proverò ad approfondire tramite qualche applicazione concreta (esercizi?)

Alla prossima

Nella prossima parte proverò ad utilizzare Maxima in modo un po’ più creativo (o forse mi ributto su qualche esercizietto, dipende). Nel frattempo spero che già abbiate iniziato a divorare il manuale e fare i vostri esperimenti.

Tre numeri a caso

ResearchBlogging.org

Nel post precedente ho enunciato quello che sembra un problemino abbastanza semplice. Talmente semplice che non è che uno ci stia a pensare più di tanto, perché la soluzione sembra ultrabanale. Ecco il problema:

Estrarre tre numeri reali positivi a caso in modo che la loro somma sia pari a 1. Fare in modo che le triplette di numeri siano distribuite uniformemente, ovvero ognuna delle possibili triplette deve avere la stessa densità di probabilità di essere estratta a caso.

Sembra talmente semplice che di solito uno applica un metodo molto intuitivo e semplice senza stare a pensarci più di tanto:

Si estraggono a caso tre numeri reali positivi tra 0 e 1, siano x1, x2 e x3, con densità di probabilità uniforme (ovvero ogni possibile reale tra 0 e 1 ha la stessa probabilità di uscire). Quindi si sommano, sia s=x1+x2+x3. I tre numeri cercati sono U1=x1/s, U2=x2/s, U3=x3/s.

Semplice ed intuitivo. Peccato che sia sbagliato: facendo così, certe triplette vengono fuori più frequentemente di altre. Come è possibile? E’ quello che si sono chiesti i miei colleghi e co-autori Enrico e Giorgio nel paper citato in fondo a questo post.

Il problema si presenta spesso nelle nostre ricerche: ci occupiamo di algoritmi per la schedulazione dei processi in sistemi in tempo reale, e spesso dobbiamo misurare e comparare le loro performance. Per fare queste comparazioni si generano a caso dei processi e si cerca di applicare un test di schedulabilità per ogni algoritmo. Spesso, il parametro più importante di un processo per noi è il carico di lavoro che impone sulla CPU. Il carico si esprime come una percentuale, quindi un numero reale tra 0 e 1. Inoltre, una CPU non puà essere caricata per più del 100%, per cui ecco che salta fuori il problema: dobbiamo generare a caso 3 (o più) numeri, uno per ogni processore, in modo che la loro somma sia un valore predeterminato e fisso non superiore a uno. Per semplicità, assumiamo qui che tale valore deve essere esattamente pari a 1.

Enrico stava appunto cercando di rifare degli esperimenti comparativi fra i due algoritmi più famosi: Fixed Priority e Earliest Deadine First, e guardando i lavori in letteratura si accorse che tutti utilizzavano il metodo di cui sopra. Al mio amico venne il sospetto che tale metodo non fosse “fair”, cioè tendesse a favorire uno dei due algoritmi contendendi. Per cui si mise a studiare il problema. Per prima cosa, generò un buon numero di triplette con il metodo della normalizzazione. Tali numeri si possono graficare in 3D, e vanno tutti a finire su un triangolo che ha i sui vertici nei punti (1,0,0), (0,1,0), (0,0,1). Ecco il triangolo che venne fuori:

Come vedete con questo metodo i punti tendono a finire più spesso verso il centro del triangolo. Come mai? C’è una spiegazione geometrica semplice: Guardate la seguente sezione in 2D (cioè, cosa succederebbe se i numeri da estrarre fossero solo 2):

I punti x1 e x2 generati a caso tra 0 e 1 stanno nel quadrato. Dopo la normalizzazione U1 e U2 vengono riportati sulla retta retta obliqua che va da (0,1) a (1,0). Quindi, tutti i punti gialli andranno sul segmento in alto; tutti i punti nell’area rossa andranno sul segmento nel mezzo. Poiché l’area rossa è molto maggiore dell’area gialla, ne consegue che molti più punti finiranno sul segmento centrale rispetto a quelli che finiscono sui segmenti laterali. Torna? Beh, ora immaginatelo in 3D.

Qualcun’altro ha proposto questo metodo:

genero solo due numeri: U1 come distribuzione uniforme tra (0,1); U2 come distribuzione uniforme tra (0, 1-x1). Infine, U3 = 1-U1-U2.

Purtroppo non va meglio, ecco cosa esce fuori:

Stanno più frequentemente da una parte! Questo perché, dopo aver generato il primo, il secondo non ha distribuzione propriamente uniforme tra 0 e 1…

Un metodo corretto è quello del rejection sampling, proposto da Alb su friendfeed, e da hronir su okpanico. Si generano punti su un rettangolo nel piano 2D che ha base e altezza pari alla base del triangolo; se il punto finisce fuori dal triangolo, si butta via e se ne genera un altro. Se invece finsce dentro il triangolo, si proietta sui tre assi per ottenere i tre numeri originari. In questo modo siamo sicur di generare uniformemente. Come riprova, ecco la figura ottenuta:

Bella vero?

C’è un solo problema: se invece di generare solo 3 numeri ne vogliamo generare N, allora il numero di punti che cascano fuori dall’iper-triangolo in N-1 dimensioni sono molti di pù di quelli che cascano dentro. Talmente tanti di più che il metodo diventa improponibile. Quindi, al crescere di N tale metodo non è utilizzabile.

E questo era proprio il problema di Enrico e Giorgio: come generare N numeri a caso tali che la somma faccia un numero fisso minore di 1? Ci sono altri due metodi possibili. Il primo è il seguente:

Generiamo N-1 numeri in maniera uniforme compresi tra (0,1), chiamiamoli x(1), x(2), … x(N-1). Poi li ordiniamo dal più piccolo al più grande, e aggiungiamo x(0) = 0 in testa, e x(N)=1 in fondo. A questo punto, i numeri richiesti saranno:
U(1)=x(1)-x(0)=x(1),
U(2)=x(2)-x(1),



U(i) = x(i) – x(i-1),


U(N)=x(N)-x(N-1)=1-x(N-1).

Per capire come funziona, pensate di avere un segmento lungo 1, e di “gettare” N-1 punti a caso su di esso. La lunghezza dei segmenti delimitati tra due punti consecutivi sono i numeri che cerchiamo.

Per N=3 si ottiene un grafico 3D proprio uguale a quello mostrato qui su: provare per credere! Inoltre il metodo ha complessità O(N log(N)), quindi va abbastanza bene. Perché funziona? Non è facilissimo fare la dimostrazione. Bisogna dimostrare che la probabilità che un segmento abbia lunghezza minore di un certo L è proporzionale a L. Facendo i conti con le probabilità nel caso N=3 se ne dovrebbe uscire facilmente. L’altra sera l’abbiamo fatto con un mio studente e ce ne siamo venuti via con un paio di paginette… 🙂 (se qualcuno ha tirato fuori un ragionamento migliore ce lo faccia sapere!)

L’ultimo metodo lo prendo dall’articolo, Enrico e Giorgio l’hanno chiamato UUniFast.

La distribuzione di probabilità della somma di N variabili casuali è la convoluzione delle distribuzioni delle variabili sommate. La pdf uniforme è una retta di inclinazione 1: P\{x \leq a\} = a. La convoluzione di N variabili uniformi è un poliomio di grado N: P\{s \leq a\} = a^N. Si può quindi sfruttare questo fatto come segue:

  • Per prima cosa si genera un numero tra 0 e 1 uniformemente, e se ne fa la radice (N-1)-esima. Chiamiamo questo numero s(N-1).
  • Quindi si genera un secondo mumero tra 0 e s(N-1), e se ne fa la radice (N-2)-esima. Chiamiamo questo numero s(N-2).
  • E così via fino a s(1).

I numeri richiesti sono quindi:

U(1)=S(1),

U(2)=S(2)-S(1),

U(i)=S(i)-S(i-1),


U(N)=1-S(N-1).

Ok, panico, se qui non mi avete seguito è perché magari non vi ricordate bene la teoria delle probabilità! Comunque, tutti i dettagli sono nel paper qui sotto.

A proposito: Enrico e Giorgio dimostrarono che il metodo “normalizzato” era a favore dell’algoritmo EDF, e quindi rifacendo le simulazioni con il nuovo metodo di generazione ristabilirono la giustizia tra gli algoritmi! (vince sempre EDF, ma con un distacco minore…)


Bini, E., & Buttazzo, G. (2005). Measuring the Performance of Schedulability Tests Real-Time Systems, 30 (1-2), 129-154 DOI: 10.1007/s11241-005-0507-9

Problema: numeri a caso

Tempo di quiz! Ecco un problemino per chi vuole cimentarsi nel week-end (la soluzione la darò lunedì sera).

Voglio un metodo per generare 3 numeri (non negativi) a caso in maniera che la loro somma sia 1. Attenzione: voglio una distribuzione uniforme delle triplette: non ci deve essere una tripletta più probabile di un altra. Sapete fornirmi un metodo? E sapreste generalizzarlo a N numeri?

Mettete pure le vostre soluzioni nei commenti!

Errare (non) è (solo) umano

La calcolatrice elettronica pisana al museo.

I precedenti post del collega Juhan e di Zar mi hanno ricordato i bei vecchi tempi di quando ero uno studente pischello a Ingegneria a Pisa (nei lontanissimi anni ’90) e tra le visite a piazza dei Miracoli in cerca di belle turiste straniere (sempre andate male) e rovinose partite di calcetto (la nostra squadra si chiamava “Inseminators”, ma non finì il campionato tanto bene), ero costretto a studiare il corso di “Calcolo Numerico”. Lo teneva il prof. Lombardi, uno di quelli che  partecipò alla costruzione della mitica CEP (Calcolatrice Elettronica Pisana, ma guardate anche qui), che potete oggi ammirare in un museo sempre in quel di Pisa. I tempi del prof. Lombardi dovevano essere migliori dei miei, perché l’Italia era all’avanguardia nella ricerca sull’elettronica e l’informatica, l’Olivetti di Adriano Olivetti era una gran bella azienda, e gli italiani erano pieni di belle speranze. Come avrete capito il prof. Lombardi non era giovanissimo, e si perdeva spesso (almeno una volta a lezione) nei ricordi del bel tempo che fu. Ma non divaghiamo troppo.

Il corso di Calcolo Numerico, tenuto ancora oggi ed obbligatorio in un certo numero di corsi di laurea a Ingegneria, insegna metodi numerici per il calcolo: per esempio a come risolvere numericamente equazioni, sistemi di equazioni, fare integrali, ecc. Uno dei primi argomenti è il trattamento dell’errore di rappresentazione nei numeri.  Lungi da me rifare qui il corso: vorrei solo farvi notare alcune cosine, che magari vi faranno capire come non sia poi così facile far fare i calcoli a un computer. Dato che non mi ricordo più tanto bene, magari dirò qualche cappellata e allora mi correggerete (spero).

Come tutti saprete, non tutti i numeri reali sono rappresentabili con esattezza. In particolare, ogni numero all’interno del nostro amato PC è rappresentato con un numero finito di bit. La rappresentazione standard per numeri con la virgola è la rappresentazione in virgola mobile. Praticamente ogni numero viene rappresentato come una mantissa e un esponente. Ad esempio, se devo rappresentare il numero 0.00002 il ‘puter lo rappresenta come 0.2 \cdot 10^{-4}. Un certo numero di cifre sono riservate alla mantissa, il resto sono riservate all’esponente. Per fare gli esempi semplici, supporrò che 5 cifre decimali siano dedicate alla mantissa, e 2 cifre decimali all’esponente. Si può quindi arrivare in valore assoluto da 0.00001 \cdot 10^{-99} a 0.99999 \cdot 10^{99}.

Quindi, per seguire l’esempio classico, se devo rappresentare \frac{1}{3} posso farlo solo con un certo errore, nel caso in questione l’errore sarà 0.00000333 \ldots. Qual’è l’errore massimo di rappresentazione che posso avere su un qualunque numero? E’ chiaro che se il numero che devo rappresentare non è troppo grande (cioè maggiore di 10^{99}) né troppo piccolo (cioè minore di 0.00001 \cdot 10^{-99}), l’errore si ha solo sulla mantissa. Cioè, se indico con r la rappresentazione del numero x, si ha l’errore assoluto

E = |x - r|.

Adesso scrivo x e r come mantissa ed esponente, utilizzando lo stesso esponente p:

x = m \cdot 10^{p} e r = m^\prime \cdot 10^{p},

dove ho indicato con m la mantissa esatta e con m^{\prime} la sua rappresentazione. Poi, invece di calcolare l’errore assoluto, calcolo l’errore relativo, cioè in percentuale sul numero:

\delta = \frac{|x - r|}{x} = \frac{|m - m^{\prime}| 10^{p}}{m 10^{p}} = \frac{|m - m^{\prime}|}{m}.

Supponendo che facciate i calcoli facendo un’arrotondamento (e non un troncamento), la differenza |m - m^{\prime}| dipenderà da quante cifre utilizzate per rappresentare la mantissa; nel caso ad esempio di 5 cifre, tale errore al massimo sarà pari a 0.5 \cdot 10^{-5}, e nel caso peggiore (m più grande possibile) avremo appunto \delta = |0.5 \cdot 10^{-5}|. Quindi, possiamo scrivere che:

|x (1 -\delta)| \leq |r| \leq |x (1 +\delta)|,

ovvero l’errore massimo relativo che posso fare nel rappresentare il numero x è appunto pari a \delta.

Fin qui è stata matematica noiosa e con un sacco di simboli (e probabilmente a questo punto il 99% dei lettori avrà smesso di leggere, e saremo rimasti in pochi, al limite sarò rimasto da solo. Va beh, meglio pochi ma buoni!). Era roba preparatoria, adesso però viene il bello: ci resta da vedere cosa succede quando facciamo dei calcoli. Analizzerò due operazioni: somma e moltiplicazione.

Supponiamo di dover sommare x e y. Le loro rappresentazioni sono \alpha e \beta. Naturalmente, quando facciamo i calcoli sul PC, invece di sommare x+y, sommiamo \alpha + \beta. Qual’è il massimo errore relativo sulla somma? Vediamo di fare i calcoli per bene, e consideriamo prima il caso in cui x e y abbiano lo stesso segno (per esempio positivo):

\frac{|\alpha + \beta - (x+y)|}{(x+y)} \leq \frac{(x+y)(1+\delta) - (x+y)}{(x+y)} = \delta

Ovvero, l’errore relativo sulla somma è lo stesso che l’errore sui termini originali. Bene, siamo confortati, possiamo sommare numeri senza paura di perdere precisione! O no? Aspettate, ci resta da vedere l’altro caso. Supponiamo che due termini x e y hanno segno diverso. Per semplificare i calcoli  supponiamo il primo sia positivo e il secondo negativo, e inoltre che il primo sia maggiore del secondo. Ecco come calcolare l’errore:

\frac{|\alpha + \beta - (x+y)|}{x+y} = \frac{|(\alpha - x) + (\beta - y)|}{x+y} \leq \frac{|\delta x - \delta y|}{x+y} = \frac{|\delta|(x-y)}{x+y}

Purtroppo è successa una cosa non molto simpatica. Se x e y hanno segno diverso e sono simili fra di loro, l’errore relativo può aumentare e di molto. Infatti, al denominatore avremo una quantità piccola mentre al numeratore una quantità grande (ricordatevi che x è positivo e y è negativo!). Ad esempio: se prendiamo i numeri 5.02 e -5.01, ognuno dei quali con un errore assoluto massimo di 0.00001, la somma sarà 0.1, con un errore massimo di 0.00002. Però 0.1 è adesso circa 500 volte più grande di 5.01, e quindi l’errore è relativamente più grande sul risultato di quanto non fosse sui numeri originali!

Insomma, non si può neanche sommare in pace! Ma se ci pensate è giusto che sia così. Se sottraete due rappresentazioni uguali, come si fa a sapere se la differenza originale fa zero oppure no? In percentuale potreste anche fare il 100% o il 1000% di errore! Va bene, basta saperlo, e tenerne conto opportunamente quando andate a fare i calcoli.

Vediamo adesso la moltiplicazione e speriamo ci vada meglio. L’errore si esprime come:

\frac{\alpha\beta - xy}{xy}

Adesso proviamo a fare un po’ di magie algebriche:

\frac{\alpha\beta - xy}{xy} = \frac{\alpha\beta - \alpha y + \alpha y - xy}{xy}= \frac{\alpha (\beta - y) + y (\alpha - x)}{xy} =

= \frac{\alpha(\beta-y) + x(\beta-y) - x(\beta-y)+y(\alpha - x)}{xy} = \frac{(\alpha-x)(\beta-y) + x(\beta-y) + y(\alpha-x)}{xy}

Osserviamo il primo termine della somma a numeratore: è il prodotto degli errori assoluti fra i due numeri. Se ad esempio gli errori sono entrambi 10^{-5}, il loro prodotto è molto più piccolo: 10^{-10}. Insomma, fidatevi, lo possiamo trascurare. Quindi:

\frac{\alpha\beta - xy}{xy} \approx \frac{\beta-y}{y} + \frac{\alpha-x}{x}.

Insomma, l’errore relativo risultante è la somma dei due errori relativi. Qundi, ogni volta che moltiplichiamo due numeri con lo stesso errore relativo, l’errore risultante raddoppia!

Immaginate di avere un numero con un errore relativo di 2^{-23} (ovvero circa 10^{-7}, e immaginate di elevarlo alla ventesima potenza (come nell’equazione di Zar). L’errore relativo sarà circa 20 volte superiore ovvero 20 \cdot 2^{-23}, ovvero, ovvero circa 2 \cdot 10^{-6}. Ci siamo mangiati un po’ di precisione, ma neanche poi tanta.

Voglio mettere un attimo a fuoco queste cose: non si scappa da questi problemi. Essi hanno a che fare con due fatti non facilmente eliminabili:

1) la precisione nella rappresentazione dei numeri reali all’interno dei calcolatori

2) proprietà matematiche indipendenti dall’esistenza o meno dei calcolatori.

In parole povere: gli ingegneri devono rassegnarsi: se hanno da fare un calcolo complicato devono analizzare l’algoritmo cercando di eliminare il più possibile le operazioni inutili, operando semplificazioni a priori, usando formule invece di calcoli iterativi ogni volta che è possibile. Insomma, non possiamo dare niente per scontato, neanche avendo a disposizione una bella e potente macchina di calcolo.

(siete arrivati fino a qui? complimenti per la pazienza!)

L’alternativa ortodossa a calc — bc

no! non sto parlando di questo

Credevate che il post su calc della scorsa settimana fosse l’ultimo della serie? Allora credevate male perché c’è bc che, a differenza di calc, viene installato automaticamente con Linux (e con Unix, naturalmente).
L’autore di bc è Mark Hopkins (almeno secondo questa pagina), il man di Ubuntu da invece Philip A. Nelson.

bc è un vero e proprio linguaggio di programmazione. Il manuale completo al solito è disponibile con

man bc

o online, All’indirizzo http://man.cx/bc dove potete trovarne un’altra versione in inglese, spagnolo, finlandese, ungherese, giapponese, coreano, polacco e turco –no l’italiano non c’è– diversa da quella installata su Ubuntu. Inoltre le versioni nelle varie lingue sono diverse (OK, non ne conosco nessuna di quelle elencate, ho semplicemente contato le opzioni in synopsis).

L’uso più elementare di bc è un pochino quirkoso

cioè bisogna passargli l’operazione da eseguire attraverso echo, usando la pipe. E ancora


bc liscio determina la precisione (le cifre dopo il separatore decimale) in base agli operandi che riceve; è possibile modificare questo comportamento settando la variabile scale. Notare che in questo caso dobbiamo racchiudere la stringa da echoare tra doppi apici.

È naturalmente possibile salvare il risultato in una variabile, così

o, seguendo la sintassi antica


che pare quella ancora preferita da qualcuno. A proposito: state seguendo il corso di BASH degli amici Bit3Lux e Lightuono, qui o qui vero?

Naturalmente si vede subito che questo modo d’uso è scomodo se si vogliono eseguire più calcoli in sequenza, eventualmente salvando i risultati in variabili. Ecco un esempio

Ma si può fare di meglio: con l’opzione -l (o l’equivalente lunga –mathlib) scale viene automaticamente impostato a 20 e sono disponibili le funzioni standard, come a – arcotangente dell’angolo, espresso in radianti.

Copiando da man.cx ecco la definizione di una funzione, mye, valore approssimato dell’esponenziale e suo esempio d’uso. È fatta per scopi puramente didattici in quanto la funzione e – esponenziale è predefinita quando si specifica l’opzione -l. Ecco il file test_e.bc

scale = 20
define mye(x) {
	auto a, b, c, i, s
	a = 1
	b = 1
	s = 1
	for (i = 1; 1 == 1; i++) {
		a = a*x
		b = b*i
		c = a/b
		if (c == 0) {
			return(s)
		}
		s = s+c
	}
}

/* calcolo l'esponenziale dei
   primi dieci interi
*/

for (i = 0; i <= 10; i++) {
   mye(i)
}

quit

eseguendola otteniamo


E con questo credo di aver illustrato cosa si può fare con bc.
Storicamente –ma non in Linux– bc si appoggia su dc.
Già resta da parlare di dc, prossimamente, forse 8)

Grafici con Scilab

Scilab è ottimo. Per tante cose. Io l’ho scoperto per i grafici, lo sto utilizzando e in questo post vorrei riassumere cosa ho scoperto.
Intanto è ottimo, lo uso e lo consiglio a tutti.

Pero non è facilissimo, è un prodotto grosso (proteiforme), si evolve (giustamente) e qualche istruzione diventa obsoleta (com’è giusto che sia).
Non sempre i manuali sono aggiornati e per quel che ho trovato manca un “Getting Started” o un “How To” per le singole branche in cui può essere suddiviso. In particolare per il plotting ho dovuto scorrere bazillioni di files .html e .pdf e in nessuno ho trovato tutto quello che cercavo.

Tra i tanti considero fondamentali:

  1. Scilab: Graphics di Anuradha Amrutkar dell’Indian Institute Of Technology, Bombay, 15th April, 2010, disponibile qui;
  2. Graphics: Getting started, qui e naturalmente qui cercando Graphics Library (e c’è tutto il resto del manuale online, da perdersi).

Allora ecco cosa ho ottenuto. Il mio sogno sarebbe stato quello di costruire uno script che, lanciato da terminale, generasse automaticamente il grafico richiesto; non è esattamente quello che sono riuscito a fare, difatti lanciando l’interprete di Scilab ottengo:
e, ovviamente, si apre il terminale di Scilab. Trascinando uno script, qualunque script, su quest’ultimo mi compare questo messaggio d’errore:
Ci sono due cose: 1) un messaggio solo parzialmente visibile, che fa riferimento a ecanvas"* che non so cosa sia, e 2) la variabile mostra non definita. Strano che c’entri una variabile con il nome italiano, non citata da nessuna parte.
Ho provato a inserire istruzioni relative a driver(), xinit() e usecanvas() senza risultato. Per mostra basta settarla a %t o %f (non so qualse sia quella giusta).

In genere la prima esecuzione dello script abortisce, completamente o produce solo parte del grafico, basta ridroppare lo script una seconda volta e tutto funziona.

Questi sono alcuni esempi, presi dalla documentazione citata.

mostra = %T;
t=linspace(0,%pi,20);
a=gca();a.data_bounds=[t(1) -1.8;t($) 1.8];

plot2d(t,[cos(t'),cos(2*t'),cos(3*t')],[-5,2 3]);
e=gce();
e1=e.children(1);e1.thickness=2;e1.polyline_style=4;e1.arrow_size_factor = 1/2;
e2=e.children(2);e2.line_style=4;e2.thickness=3;
e3=e.children(3);e3.line_mode='on';e3.mark_background=5;

hl=legend(['cos(t)';'cos(2*t)';'cos(3*t)']);

xs2png(0, "e4.png")

mostra = %F;
t = linspace(-%pi, %pi,30);
function z = mysurface(x,y);
    z = x * sin(x) ^ 2 * cos(y);
endfunction
contour(t,t,mysurface,10)
xs2png(0, "aa2.png")

mostra = %F;
x = [0 : 0.1 : 2 * %pi]';
plot2d(x, [sin(x) sin(2*x) sin(3*x)], ..
       [1,2,3], leg = "L1@L2@L3", nax = [2,10,2,10],rect = [0,-2,2*%pi,2]);
xs2png(0, "scp4.png")

mostra = %f;
t=-%pi:0.1:%pi; m=sin(t)'*cos(t);
clf()
xset("colormap",jetcolormap(64))
colorbar(-1,1)
Sgrayplot(t,t,m, strf="041")
xs2png(0, "scp6.png")

mostra = %f;
n = 30;
nt = 100;
x = linspace(0,2*%pi,n);
y = linspace(0,%pi,n/2);
z = sin(x')*sin(y);
t = linspace(0,4*%pi,nt);
show_window(); clf()
f=gcf();
f.color_map=jetcolormap(64);
f.pixmap='on';
colorbar(-1,1)
Sgrayplot(x,y,cos(t(1))*z, strf="042", zminmax=[-1,1])
c=gce(),e=c.children
xtitle("Kaa''s eyes")
for i = 1:nt
  e.data(:,3)=matrix(cos(t(i))*z,-1,1);
  show_pixmap()
end
f.pixmap='off';
xs2png(0, "scp7.png")

(quest’ultima è dinamica, provate…).

Bon. Per adesso Scilab lo uso così.
Ho trascurato altre cose, per esempio lo sapete che è uscita la nuova versione di newLISP, la 10.3? Lo sapevo che non interessa nessuno 😦
Ho trascurato anche Python, ma restate sintonizzati, prossimamente…
E poi tante altre cose nuove 8) e nuovi collaboratori 8)

Scilab, le funzioni e tutte le altre meraviglie

In questo post, continuazione del precedente,si coninua l’esame del pdf progscilab per un illustrazione veloce delle notevoli capacità del package. È utile per arrivare a una mia conclusione, giù in fondo in “Conclusione” 😉

In Scilab le funzioni sono estremamente flessibili. Intanto si distinguono in macro, scritte nel linguaggio di Scilab e primitive, scritte in qualche altro linguaggio (per esempio C, C++, Fortran) e compilate.
Nella maggior parte dei casi non c’è differenza tra i due tipi e per questo si usa il termine generico function. Ovviamente la cosa cambia nel debug.
Le funzioni sono variabili: non male per un progetto partito avendo come riferimento il Fortran 😀
D’altronde oggi è una cosa abbastanza comune, a me capita di farla anche senza volerlo: in Python basta dimenticare le parentesi. Ma è di una comodità unica, conosciuta come callback.
A p.57 e ssg viene spiegato molto bene, con esempi, come creare funzioni con numero di argomenti variabili, sia in input che in output, cose tecniche: se servono leggetevele e sapetele 😉
A p.81 si parla dello scope delle variabili. Diversi casi sono possibili, attenzione! La soluzione più semplice (il livello dello stack) di p.82 può creare problemi, raccontati nelle pagine successive. A me casi del genere sono capitati più volte, provo una discreta Schadenfreude venire a conoscenza che succede anche a altri, tanto spesso che vale la pena di riportarli in un manuale 8)
A p.93 e ssg viene descritto come sia possibile generare dinamicamente le funzioni come stringhe cosa che conferisce a Scilab un livello di flessibilità altrimenti irraggiungibile (ma a me comincia a frullarmi in t’a capa il Lisp).
Le tecniche per rendere l’esecuzione veloce sono descritte a partire da p.99, Performances.
Da copiare l’idea della coppia di funzioni gemelle tic() e toc(): la prima fa partire il cronometro, la seconda lo blocca. In uno dei prossimi post le faccio in Python, newLISP e Fortran, forse, se nel frattempo non me ne sono dimenticato o sono stato attratto da qualcos’altro.
Bello l’esempio d’uso del profiling, da p.102 (profilaggio? la banda del nano usa “briffare” addirittura!).
A p.110 inizia la vettorizzazione, da leggere. Poi vengono i trucchi per ottimizzare il codice, p.115. Regola aurea: usare le funzioni di Scilab e evitare come la peste i cicli for.
Poi basta. Il discorso diventa troppo tecnico e mi sono addormentato 😦

Il pdf 419922 An Introduction to Scilab (un’altra) di Satish Annigeri
(prof. per ingegneri civili) è conciso (55 pagine) e pratico: a p.35 ci dice come gestire file su disco, nella pagina successiva già si occupa di plottaggio di curve polinomiali e poi passa a leggere file di Excel.
A p.47 passa a un caso pratico da ingegneri civili, calcolo di una struttura reticolare piana arrivando a determinare le sollecitazioni nelle aste. Notevole.

Il pdf intro dello Scilab group a p.87 e ssg ci illustra come gestire i grafici.
L’esempio di p.81 mi da errore e allora ho provato a dare le istruzioni direttamente dentro al Terminale si Scilab. Comodissimo il comando xset() che apre la finestra grafica e permette di modificare il grafico. Ecco il mio risultato:

Proseguendo con gli esempi non tutti vengono, ci sono errori nei comandi, alcuni di difficile individuazione, i font sono diversi da quelli del manuale (Windoze strike again!), …
Ci sono però numerosi file di esempi, richiamabili direttamente anche all’interno dell’editor (alle volte sono davvero gnugnu).

Il capitolo 6 a p.101 ci introduce a una delle caratteristiche più avanzate di Scilab: Interfacing C or Fortran programs with Scilab. Chiaramente questo sarà usato solo da utenti skillati per task non banali ma sembra fattibile, semplice e ci sono esempi. No non ci provo, sono troppo vecchio per quest cose. Però: avete nai provato a farlo con M$ DOS (sto parlando quando ancora non si usava Windoze): frustrante. Molto meglio la filosofia *NIX/Linux.

Per la parte matematica sono naturalmente previsti collegamenti con Matlab e Maple.

Conclusione

Il manuale Programming in Scilab di Michaël Baudin, contenuto nel file progscilab-so-v0.6.pdf è fondamentale.
Riassumerlo è impossibile 8)
Ciò porta a una conclusione: questi post su Scilab non possono sostituire il manuale. Per diversi motivi:

  • il manuale è ottimo;
  • Scilab è un tool molto grande e completo, il suo linguaggio di programmazione anche;
  • compito del blog non è fare un corso bensì segnalare le chicche che trovo, anzi che mi piacciono;
  • quando possibile adoro proporre script che per questioni oggettive debbono essere piccoli, se piacciono sarà il lettore a completarli.

Prossimamente

Il modulo Plot2d mi attira. Ma ci sono anche altre cose nel mondo della programmazione. È da parecchio che non nomino il newLISP, Pyhton idem. Poi ci sarebbero altri linguaggi, non quelli predominanti C++ e Java, su quelli c’è già troppa roba, altri meno noti e pubblicizati ma intriganti. E poi ci sarebbe da considerare anche gli utenti di Windows: alcune cosa che con Linux si danno per scontate lì mancano, o co$tano.
Prossimamente, se non mi sto beccando l’influenza, sempre su questo canale 😀

Scilab – polinomi e matrici

Dopo l’intro e aver visto come usare Scilab per fare un grafico oggi racconterò di come usarlo per fare calcoli. Scilab adora contare e allora via con i calcoli.

Intanto sto ancora imparando ma è disponibile tutto quello che serve, anzi molto di più. Parecchi documenti si trovano a Scilab Wiki. In particolare Tutorials and documentations dove c’è, tra un mucchio di altre cose, Programming in Scilab di Michaël Baudin.

A p.22 introduce i polinomi. Seguendolo definiamo un polinomio in funzione delle sue radici 1 e 2: p(x) = (x-1)(x-2) con

p = poly([1,2], "x")

Si può definire il polinomio in funzione dei suoi coefficienti (in ordine crescente): q(x) = 1 + 2x sarà

q = poly([1,2], "x", "coeff")

A questo punto i polinomi si possono sommare, sottrarre, moltiplicare:

-->m = p * q
 m  =

              2    3
    2 + x - 5x + 2x   

-->

Per valutare il valore del polinomio in un dato punto si usa horner. Ecco, niente panico!, un po’ di sano terrorismo matematico, qui su wikipedia, ma con Scilab niente panico!

:-->p = poly([1,2], "x")
p  =

          2
2 - 3x + x

-->horner(p, [0 1 2 3])
ans  =

2.    0.    0.    2.

che sono, ovviamente, i valori per x = 0, 1, 2, 3. Quindi dimenticate la pagina di wiki: con Scilab i calcoli (noiosi e (dai) astrusi) li fa lui.

roots ci da le radici, verifica:

-->p = poly([1,2,3,4,5], "x")
 p  =

                     2     3     4   5
  - 120 + 274x - 225x + 85x - 15x + x   

-->roots(p)
 ans  =

    1.
    2.
    3.
    4.
    5.    

OK, viso che funziona

-->p = poly([1,2,3,4,5], "x", "coeff")
 p  =

               2    3    4
    1 + 2x + 3x + 4x + 5x   

-->roots(p)
 ans  =

    0.1378323 + 0.6781544i
    0.1378323 - 0.6781544i
  - 0.5378323 + 0.3582847i
  - 0.5378323 - 0.3582847i  

numeri brutticomplessi, ma ecco:

-->t = poly([24, -50, 35, -10, 1], "x", "coeff")
 t  =

                  2     3   4
    24 - 50x + 35x - 10x + x   

-->roots(t)
 ans  =

    1.
    2.
    3.
    4.  

Scilab è sempre più simpatico 😀
Uno (cioè a dire il vero io, fino a un istante fa) potrebbe pensare che Scilab tratti solo matrici bidimensionali: falso, c’è l’ipermatrice (hypermatrix)

This feature is familiar to Fortran developers, where the multi-dimensional array is one of the basic data structures.

Quindi la volta scorsa ho detto una cosa non vera: chiedo venia ma mi verrebbe voglia di dare la colpa al Fortran 😦

A p.51 di progscilab si inizia a parlare delle funzioni. Potrebbe essere l’argomento di un prossimo post, forse…