Archivi Categorie: math

Siamo tutti un po matematici, no?

log – qual è la base? – 2

Boltzmann’s equation — carved on his gravestone

Continuo da qui.

Roberto Zanasi che già aveva contribuito al post precedente mi ricorda in un commento “Il logaritmo in base 10 lo usi se hai a che fare coi decibel, per esempio“. Vero, come dice la Wiki. Notare che la versione inglese è diversa, più approfondita ma meno pratica.

A me sono venuti in mente due casi da esaminare.

La sensibilità (velocità) della pellicola che può essere espressa in scala lineare (ASA) o logaritmica (DIN).

La conversione da un valore in scala logaritimca () al corrispondente valore in scala lineare S è data da S = 10^((S°-1)/10) con arrotondamento al valore più vicino fra quelli indicati nella prima colonna della tabella sopra [non la riporto, è là].

Verifica per un valore tipico di quando ero giovane 100 ASA corrisponde a 21 DIN e 400 ASA a 27 DIN:

$ py3
>>> 10**((21-1)/10)
100.0
>>> 10**((27-1)/10)
398.1071705534973

OK, quel 10 è un indizio…

La conversione inversa (cioè da un valore di rapidità in scala lineare al corrispondente valore in scala logaritmica) è data da S° = 10 log(S) + 1 con arrotondamento all’intero più vicino.

>>> from math import log
>>> 10*log(100)+1
47.05170185988092
>>> 10*log(100)/log(10)+1
21.0
>>> 10*log(400)/log(10)+1
27.020599913279618

Sì! sono logaritmi base 10. E siamo 2 a 1 nella gara tutti contro me 😡

Ma non m’arrendo ancora, c’è Bolzmann, la formula dell’entropia. Il suggerimento me l’ha dato Carlo Rovelli (a sua insaputa), cap. 12, Informazione, p.212 di La realtà non è come ci appare. Presto la Wiki.

In statistical mechanics, Boltzmann’s equation is a probability equation relating the entropy S of an ideal gas to the quantity W, the number of real microstates corresponding to the gas’ macrostate: S = kB ln ⁡W, where kB is the Boltzmann constant (also written as simply k) and equal to 1.38065 × 10−23 J/K.

OK, ln è l’altro modo per indicare log quando questo potrebbe generare incertezze.

Allora a che punto siamo? Rileggendo la nota di Roberto riportata nel post precedente è chiaro che anche lui è con Bolzmann e me; e anche .mau. (non so se posso citarlo, se del caso scancello) è nel giusto. Insomma abbiamo vinto, ma nel dubbio meglio specificare sempre.


L’immagine viene dalla Wiki, una delle cose belle che giustificano il Web.

log – qual è la base?

Recentemente mi sono imbattuto in quella che è diventata una delle mie ossessioni preferite. Adesso vi conto, così la supero, forse. C’era da indovinare come inserire i dati in una formula di cui avevamo una descrizione incompleta e una tabella di (pochi) valori. Si può fare? proviamo in due, indipendentemente, e otteniamo una buona correlazione ma con parametri differenti. Dal confronto salta fuori che è colpa del logaritmo. Cioè di che base usare.

Io in questi casi mi rivolgo alla Wiki. Ottima, esauriente, come di solito.

Ho chiesto agli amici di Facebook:

Se sono Veri Matematici è il logaritmo naturale, quindi in base e. Se c’è gente relativamente normale (relativamente, perché chi parla di logaritmi dal salumiere?) è il logaritmo in base dieci (e in questo caso il logaritmo naturale si scrive ln come sulle calcolatrici. Nota che gli informatici non parlano di logaritmo in base 2, perché dicono semplicemente bit.

e

Già, si dice che chi usa log per indicare logaritmo naturale poi usi Log per indicare il logaritmo in base 10, ma io non l’ho mai visto fare (perché chi usa log al posto di ln in realtà usa solo log, le altre basi sono per gli ingegneri).

Ah! sì, ricordo, quando usavo le calcolatrici, come nell’immagine lassù. E prima ancora il Bruhns per topografia (alle superiori ho frequentato l’istituto per geometri, il Guarini a Torino).

Ma me l’ero scordato; tutti (quasi) i linguaggi di programmazione usano base e. Ma allora perché a scuola (tutte?) si usa base 10? Non credo ci sia ancora chi li usi per moltiplicazioni, divisioni, radici ed elevazione a potenza. O anche qui sono male informato?

A proposito dei linguaggi c’è chi è davvero smart, ecco Python:

math.log(x[, base])

With one argument, return the natural logarithm of x (to base e).
With two arguments, return the logarithm of x to the given base, calculated as log(x)/log(base).

$ py3
>>> import math
>>> math.log(10)
2.302585092994046
>>> math.log(100)
4.605170185988092
>>> math.log(10, 10)
1.0
>>> math.log(100, 10)
2.0
>>> math.log(math.exp(1))
1.0
>>> math.log(9, 3)
2.0
>>> math.log(27, 3)
3.0
>>> math.log(42, 4)
2.6961587113893803
>>> 4**_
42.000000000000014

Bravo Python ma non è il solo, ecco Racket:

(log z [b]) → number

  z : number
  b : number = (exp 1)

Returns the natural logarithm of z. The result is normally inexact, but it is exact 0 when z is an exact 1. When z is exact 0, exn:fail:contract:divide-by-zero exception is raised.

If b is provided, it serves as an alternative base. It is equivalent to (/ (log z) (log b)), but can potentially run faster. If b is exact 1, exn:fail:contract:divide-by-zero exception is raised.

$ rkt
Welcome to Racket v6.11.
> (log 10)
2.302585092994046
> (log 100)
4.605170185988092
> (log 10 10)
1.0
> (log 100 10)
2.0
> (log (exp 1))
1.0
> (log 9 3)
2.0
> (log 27 3)
3.0
> (log 42 4)
2.6961587113893803
> (define le (log 42 4.2))
> le
2.6044944060210176
> (expt 4.2 le)
42.0

Per tutti gli altri la regola è semplice log_b(z) = log(z) / log(b).

Ossessione superata 😁 E tutti i link che ho trovato li butto, fatto 😁

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:

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)