Robitex's Blog

Ideas in the web

Archivi Categorie: Math

Plot della funzione Gamma con matplotlib


Scarica l’articolo in formato pdf

Abstract

Utilizzando alcune librerie Python di calcolo numerico utilizzeremo il pacchetto matplotlib per disegnare il grafico della funzione Gamma di Eulero. Questa infatti è la funzione scelta per introdurre questo interessante ambiente di computazione e visualizzazione, integrato e circoscritto all’interno del linguaggio Python.

Lo scopo è quello di mostrare le capacità di questo software libero e dare un idea delle caratteristiche di un linguaggio testuale basato su Python per il calcolo e la visualizzazione tecnico-scientifica.

Numeri, numeri ed ancora numeri per Python

Il linguaggio Python è particolarmente portato per i numeri. Lo dimostra questo piccolo e classico esempio in cui si calcola il fattoriale di un numero intero elevato.
Il fattoriale non è altro che il prodotto di tutti i numeri da 1 ad n e può essere calcolato da una funzione ricorsiva. Ecco l’esempio di codice:

def fact(n):
    if n==1:
        return 1
    return n * fact(n-1)

print(fact(100))

Salvato il file contenente il codice con il nome fact.py , per eseguirlo aprire una finestra di terminale e digitare $ python fact.py. Sarà immediatamente mostrato a video un numero di 158 cifre. Nello screenshot seguente vi ho riportato l’esecuzione di fact.py posizionato nel mio caso nella cartella Scrivania:

Esecuzione dello script Python per il fattoriale di 100

Esecuzione dello script Python per il fattoriale di 100

È preferibile riscrivere il codice senza copiarlo per rendersene meglio conto, inoltre ricordatevi che se non vogliamo digitare il percorso completo del nome del file da eseguire, la directory di lavoro del terminale o della console di Windows, deve essere la stessa di quella in cui si trova il file. Per settare la directory di lavoro, e quindi spostarsi nell’albero del file system, basta usare il comando cd (che sta per change directory) seguito dal nome della directory.

Ultima cosa da ricordare, questa volta agli utenti Linux, è l’assegnazione agli script Python dei permessi di esecuzione, per esempio semplicemente spuntando il flag Esecuzione nella scheda Permessi del dialogo Proprietà del file stesso.

Preparazione del sistema

In Ubuntu installare i file necessari per svolgere l’esercizio proposto è semplice perché l’interprete Python è pre-installato di default mentre le librerie matplotlib e scipy si installano sul sistema con i due comandi (la libreria numerica numpy è una dipendenza obbligatoria del pacchetto python-matplotlib e non è necessario richiederne l’installazione esplicitamente):

$ sudo apt-get install python-matplotlib
$ sudo apt-get install python-scipy

Per gli utenti Windows serve un po’ più di tempo ma mutatis mutandis, non è difficile seguire le procedure d’installazione di Python e delle librerie numeriche indicate. La documentazione è ottima.

matplotlib plotting example

La libreria matplotlib funziona allo stesso modo di Matlab: la curva di una funzione è disegnata unendo un insieme di punti del piano le cui coordinate sono definite da un set di valori del dominio della funzione e dai corrispondenti valori calcolati nel codominio.

Quello che dobbiamo fare è scrivere un programma in Python che utilizzi motori numerici e di visualizzazione sotto forma di librerie

Consideriamo un esempio di plot di funzione per f(x)=sen^2(x). Il set di valori sull’asse x saranno i 100 punti egualmente intervallati da 0 a 2\pi: setx = np.linspace(0, 2*np.pi, 100), dove np è l’oggetto che fa riferimento alla libreria numpy definito con un precedente import numpy as np (il nome assegnato è una sorta di scorciatoia).
Il plottaggio avviene utilizzando il metodo plot dell’oggetto matplotlib.pyplot, a cui passiamo il set di valori sull’asse x ed i corrispondenti ordinate sull’asse y calcolate con la funzione seno al quadrato: plt.plot(setx, np.sin(setx)**2) dove plt è l’oggetto che fa riferimento alla libreria matplotlib.pyplot definito con import matplotlib.pyplot as plt.

Non rimane che esercitarci sul codice completo del programma che genererà un file pdf, quindi un file di qualità vettoriale che scalato non perderà risoluzione:

#!/usr/bin/python

import numpy as np
import matplotlib.pyplot as plt

setx = np.linspace(0.0, 2*np.pi, 100)
plt.plot(setx, np.sin(setx)**2)

plt.savefig('sen2plt.pdf')
Il grafico della funzione seno al quadrato (click to download the pdf)

Il grafico della funzione seno al quadrato (click to download the pdf)

Sul manuale troverete le usuali proprietà per personalizzare il grafico come etichette compreso testo in codice LaTeX per il typesetting matematico, spessori e colori, griglie, ecc. Le lavorazioni sui sorgenti assomigliano operativamente a quelle degli altri tool di visualizzazione come Matlab stesso o Gnuplot, ma con i vantaggi di un linguaggio ad oggetti intuitivo. Occorre solo l’abitudine ad orientarsi tra i moduli delle librerie.

Esiste anche la possibilità di visualizzare il grafico in una finestra grafica. Non seguiremo questa strada perché interattiva per l’utente e quindi poco adatta ad essere descritta in una guida. Dico solo che per usufruire di questa modalità basta sostituire nell’ultima istruzione il metodo savefig(‘nomefile’) con show().

The Gamma function

La funzione \Gamma fu trovata da Eulero su sollecitazione di Christian Goldbach e rappresenta il fattoriale di un numero reale (od anche complesso). Si tratta quindi di una sorta di estensione del calcolo del fattoriale di un numero intero infatti per essi sussiste la seguente relazione:

\displaystyle \Gamma(n+1)=n!

La definizione standard della funzione Gamma di Eulero è la seguente:

\displaystyle \Gamma(z)=\int_0^\infty e^{-t} t^{z-1}\, dt

Per conoscere ulteriori dettagli matematici visitate la classica pagina di Wikipedia o quella di MathWorld.

matplotlib plots Gamma

La funzione \Gamma è stata implementata nel modulo speciale della libreria scipy. Il nostro obbiettivo meno ambizioso quanto a visualizzazioni 3D nel piano complesso, è tracciare il grafico per un dominio reale ed evidenziare i punti sulla curva corrispondenti ai fattoriali dei numeri interi.
Ma la funzione \Gamma non si domina facilmente. Giustamente matplotlib non disegna punti a distanza infinita, quindi occorre gestire le discontinuità ovvero i valori infiniti che essa assume per 0 e per i numeri interi negativi, appiattendo i valori della funzione quando superano un range scelto opportunamente e molto al di fuori dei limiti del grafico (questo per la verità va fatto per qualsiasi altra funzione con asintoti verticali nell’intorno dei quali i valori diventano ovviamente elevatissimi).
Poi occorre infittire molto i punti in cui calcolare \Gamma altrimenti nel grafico compariranno dei rami parassiti nei pressi degli asintoti verticali.

Passo infine la parola al codice:

#!/usr/bin/python

# import delle librerie necessarie
import numpy as np
from scipy.special import gamma as Gamma
import matplotlib.pyplot as plt

# definizione dei limiti del grafico
# per il facile 'aggiustamento'
x0, x1 = -5.0, 5.2
y0, y1 = -6.5, 7.0

# dominio lineare sull'asse x
x = np.linspace(x0, x1, 50000)

# dominio numeri interi
xx = [2,3,4]

# creazione dei vettori di ordinata
y = Gamma(x)
yy = Gamma(xx)

# 'gestione' delle discontinuita
ulim = 100*y1
llim = 100*y0
y[y>ulim] =  np.inf
y[y<llim] = -np.inf
# plottaggio della Gamma
plt.plot(x, y)

# aggiunta dei punti notevoli del fattoriali
# l'opzione 'ro' si riferisce a singoli caratteri:
# r = red , o = pallino
plt.plot(xx, yy , 'ro')

# limiti sugli assi
plt.axis([x0, x1, y0, y1])

# tracciamento di una griglia base
plt.grid()

# creazione del file pdf con il grafico
plt.savefig('gammaplt.pdf')
Grafico della Gamma di Eulero (click to download the pdf file)

Grafico della Gamma di Eulero (click to download the pdf file)

Conclusioni

Secondo me bene si fa a tener conto delle librerie Python introdotte in questo post. Basta dare un occhiata ai numerosi talk ai meeting europei di Scipy.
Impara l’arte e mettila da parte…

Annunci

Un brindisi matematico…


Una risposta estremamente difficile

Ecco la domanda:

Se brindano n persone quanti cin faranno i loro bicchieri?

Una domanda certamente fondamentale per il proseguo la cui risposta necessita di impegno estremo.

Sia B(n) la funzione che fornisce il risultato, allora considerando il punto di vista di un singolo partecipante al brindisi, egli effettuerà un cin cin con tutti gli altri ovvero n-1 tocchi e lui potrà dirsi soddisfatto. Rimarranno allora i brindisi seguenti:

\displaystyle B(n) = n-1 + B(n-1)

Applicando lo stesso ragionamento al gruppo rimanente, di fatto un ragionamento ricorsivo, dovranno essere ancora compiuti B(n-1)=n-2+B(n-2) brindisi.

Si finirà di brindare quando rimarranno due sole persone che faranno un unico tocco, allora:

\displaystyle B(n) = \sum_{i=1}^{n-1} i

Il risultato è che il numero dei brindisi in un gruppo di n persone è la somma dei numeri da 1 a n-1. Per esprimere la somma con una semplice espressione seguiamo il ragionamento del giovanissimo Gauss.
Su n, sommiamo il primo termine, ovvero il numero 1 e l’ultimo termine quindi n ed otteniamo naturalmente 1+n. Così facciamo per il secondo ed il penultimo termine della somma 2+(n-1)=1+n. Alla fine la somma è sempre il termine n+1 moltiplicato per le coppie il cui numero è n/2 ovvero:

\displaystyle (n+1)\frac{n}{2}

Tornando al numero dei brindisi fra n persone, dovremo sostituire alla relazione precedente ad n il termine n-1 ottenendo il risultato tanto atteso:

\displaystyle B(n)=n\frac{n-1}{2}

Altro percorso di dimostrazione

La cosa si fa più interessante se immaginiamo i brindisi tra le persone come rami di relazioni. questo punto di vista non si riferisce alla singola persona come quello precedente ma all’insieme intero delle relazioni.

Tutte le relazioni tra n persone, possono essere rappresentate in una griglia quadrata n \times n dove le righe e le colonne sono riferite a tutte le singole persone.
Così per esempio la cella della riga relativa alla persona A e della colonna relativa alla persona B, rappresenta il brindisi tra la A e B.

Ogni cella rappresenta una relazione di brindisi ma dobbiamo escludere quelle che si trovano sulla diagonale principale della griglia (quelle in cui la persona sulla riga è la stessa di quella della colonna) perché è escluso che si possa contare anche il brindisi con se stessi.

Delle celle rimanenti ne possiamo considerare tuttavia solo la metà perché se esiste la relazione di brindisi A \to B esisterà anche la relazione B \to A sulla griglia ma delle due ne possiamo conteggiare solo una perché si tratta delle stesso brindisi (ramo di relazione non orientato).

Una rappresentazione grafica della griglia è riportatata di seguito per 5 persone denominate da A a E, dove col pallino rosso sono evidenziate le celle di auto-relazione mentre con i quadratini verdi e blu sono evidenziate le relazioni doppie da contare solo una volta.

La griglia dei brindisi

La griglia dei brindisi

Alla fine ritroviamo elegantemente la relazione per il conteggio dei brindisi tra n persone:

\displaystyle B(n)=\frac{n^2-n}{2}.

Con questo punto di vista il numero dei brindisi non è altro che il numero delle celle della griglia che si trovano sotto oppure sopra alla diagonale principale (il numero dei quadratini di uno stesso colore se ci riferiamo all’immagine presentata che potete scaricare nel formato pdf da qui, se invece vi interessa il codice che l’ha generata consultare questo post sul forum del GuIT).

Epilogo

Se alla prossima festa siete in 12 preparatevi dunque a 66 rintocchi di bicchieri!
Alla salute!!!

Potenze di somme


Scarica l’articolo in formato PDF per la stampa

L’elegante espressione

Stavo leggendo le prime pagine del libro di Lev Landau dedicato alla teoria dell’elasticità (il volume 7 del suo corso di Fisica pubblicato in Italia da Editori Riuniti), e mi sono imbattuto in un’espressione matematica interessante. Si tratta del quadrato di un trinomio espresso come somma di prodotti.

Sappiamo che l’espressione del quadrato della somma di tre termini è:

\displaystyle (x_1+x_2+x_3)^2= x_1^2+x_2^2+x_3^2+2x_1 x_2+2x_1 x_3+2x_2 x_3.

Nel libro di Lev Landau la potenza della somma è rapprensentata come:

\displaystyle (x_1+x_2+x_3)^2=\sum_{i=1}^3 \sum_{j=1}^3 x_i x_j.

La mia prima impressione è stata quella che qualcosa non tornasse perché quando nella sommatoria gli indici si incrociano ottengo i termini del doppio prodotto con un appropriato fattore 2 mentre i termini quadrati ottenuto con indici uguali assumevano tale fattore ingiustamente. In realtà l’espressione è perfettamente corretta perché i due indici i e j sono sempre unici nella combinazione, solo che quando sono diversi danno luogo a due termini uguali. Per esempio per la coppia 1,2 ottengo il termine x_1 x_2 che si somma con il termine x_2 x_1 quando gli indici diventano appunto 2, 1, mentre esiste una sola combinazione che genera il termine x_1^2, quando gli indici sono entrambe pari ad uno.

In generale, se n è il numero dei termini che sommati sono elevati al quadrato, il numero dei termini dello sviluppo risultante della doppia sommatoria è n^2 di cui n sono i termini quadratici e n^2-n sono i termini misti ma solo la metà di essi risultano distinti. In definitiva la doppia sommatoria produce gli (n^2+n)/2 termini dello sviluppo del quadrato. Nel nostro caso n vale 3 per cui lo sviluppo conta 6 termini.

Facile è dimostrare l’elegante espressione infatti il quadrato della somma si può scrivere come:

\displaystyle \left(\sum_i x_i\right)^2=\sum_i x_i \sum_j x_j,

dove abbiamo introdotto un secondo indice j in sostituzione del primo in uno dei due termini uguali del prodotto. A questo punto otteniamo la dimostrazione riscrivendo l’espressione precedente spostando il termine x_i all’interno della sommatoria con indice j essendo in effetti una costante nei riguardi di quest’ultima somma:

\displaystyle \sum_i x_i \sum_j x_j= \sum_i \sum_j x_i x_j.

Estensione

Il libro sull’elasticità utilizza la simbologia compatta detta notazione di Einstein, per la quale si ommettono nelle espressioni matematiche le sommatorie intendendo che termini con indici vanno sommati rispetto ad essi. L’espressione in forma compatta che abbiamo incontrato prima si scrive semplicemente come (x_i)^2=x_i x_j.

La regola si può estendere per potenze intere qualsiasi di un numero qualsiasi di termini. Per esempio il cubo del trinomio può essere scritto come una tripla sommatoria del prodotto ternario:

\displaystyle (x_1+x_2+x_3)^3=x_i x_j x_k.

Il numero di sommatorie ed il numero dei termini del prodotto è pari a quello dell’esponente della potenza, mentre il numero dei termini della somma compare come indice finale delle sommatorie (implicite) stesse.

Come ulteriore esempio, la quinta potenza del binomio sarà quindi espressa da:

\displaystyle (x_1+x_2)^5=x_i x_j x_k x_l x_m

Conclusione

A parte il fatto che l’argomento del post sembra essere interessante, vorrei concludere dicendo che le scienze applicate come la fisica o l’ingegneria sanno essere un continuo generatore di idee e di problemi per la matematica e, naturalmente, che vale anche l’opposto. Non si è quindi mai sicuri nel ritenere un’astrusa teoria matematica come perfettamente inutile o nel pensare che una nuova teoria fisica non richieda la generazione di nuovi risultati assolutamente inattesi per la scienza esatta.

Il triangolo di Tartaglia (Pascal triangle)


Il triangolo di Tartaglia

Eravamo nel lontano 1983 (ricordo l’uscita dell’ultimo album “The final cut” dei Pink Floyd.), quando frequentai il mio primo corso di informatica.
Era da poco uscito l’avanzatissimo personal computer di Olivetti chiamato M20, ed il corso verteva sulla programmazione in linguaggio Basic in ambito gestionale con sessioni di pratica davanti ad alcuni M20 nuovi fiammanti. In altre parole, un’esperienza da ricordare.

Bé direte voi, cosa centra il Triangolo di Tartaglia?
Fatemi continuare!

Al corso ricevetti una serie di dispense sul Basic (che conservo ancora), e la possibilità di effettuare esperimenti reali su macchina.

Il docente era un giovane da non molto laureato in informatica, dai modi di fare gentili e pacati, che esponeva gli argomenti con molta precisione.
Ci faceva ragionare sulla sintassi del Basic disegnando alla lavagna dei diagrammi che rappresentavano vari percorsi alternativi che alla fine producevano l’espressione da inserire nel programma.

Il primo esercizio di programmazione, fu appunto il calcolo e la stampa del triangolo di Tartaglia.

Si doveva utilizzare una struttura dati chiamata array e riempirla della serie numerica del triangolo che, diciamolo, fornisce, tra le molte altre cose, i coefficienti binomiali dell’espressione algebrica (alla riga n):

\displaystyle (a+b)^n

All’epoca del corso non era quindi neanche immaginabile lo sviluppo informatico successivo, ad un corso dove i programmi dovevano essere inseriti ed editati con i numeri di linea, il processore girava a 4 MHz, e l’Hard Disk (quando c’era) poteva superare di poco i 10 MB. Il tutto al costo di ben 16.000.000 di Lire.

Poteva non entrare in Azione Metapost?

L’idea base di Metapost fu resa nota alcuni anni più tardi, nel 1989, e la prima release pubblica risale al 1994.
Questa è un occasione davvero imperdibile: scrivere oggi, con un salto nel tempo, il programma del Triangolo di Tartaglia in Metapost.

L’idea da implementare deve essere, per forza di cose, elegante, almeno per confrontarci con come eravamo allora… (ma anche con il tema del Blog che è “idee nel web”).

   1

  1 1

 1 2 1
 \/
1 3 3 1
...

Un numero nel triangolo si costruisce sommando i due numeri ad esso superiori. Per esempio alla riga quattro, che corrisponde al cubo del binomio, il secondo numero è la somma dei due numeri segnalati dagli slash.

Scartiamo la ricorsione per mantenerci semplici, e valutiamo la seguente idea:
ogni elemento del triangolo si compone di due informazioni: la prima è il valore numerico, la seconda è la coordinata per il posizionamento grafico.

Ogni valore numerico dipende da quelli della riga precedente, ed ad ogni nuova riga troviamo un elemento in più.

Basta dunque lavorare con un array (ancora un array…) chiamato values[], a cui aggiungiamo un elemento alla fine, e lo rielaboriamo tante volte quante sono le righe volute, per trasformarne i valori in quelli corretti di riga.

Abbiamo decritto dunque un classico, due cicli for nidificati, dove nel più interno a partire dalla fine, vengono modificati i valori sommando quelli attuali nell’array (che ancora detengono quelli della riga precedente), e disegnando gli elementi nei punti corretti.

L’informazione della posizione di tracciamento è mantenuta da una variabile punto chiamata pos, che localizza l’ultimo numero sulla riga corrente, mentre la posizione sulla stessa riga è ottenuta sottraendo da pos il vettore (xstep,0) il giusto numero di volte.

Regolando i valori di xstep ed ystep, è possibile aggiustare facilmente l’aspetto del triangolo.

outputtemplate:="tartaglia.mps";

beginfig(1);
   % Tartaglia array values
   numeric values[];
   values[0]=0;

   % triangle deep
   numeric deep;
   deep := 13;

   % position of the triangle top number
   pair pos;
   pos := origin;

   % step among row
   numeric ystep;
   ystep:= 18;
   % step among column
   numeric xstep;
   xstep := 36;

   % main loop
   for k:= 1 upto deep:
      values[k] := 1;
      label("1",pos);

      for i:= k-1 downto 1:
         values[i]:=values[i]+values[i-1];
         label(decimal values[i], pos - (k-i)*(xstep,0));
      endfor;

      pos := pos + (xstep/2,-ystep);
   endfor;
endfig;
end
Il triangolo di Tartaglia

Il triangolo di Tartaglia

Se avete compreso il codice, e quindi probabilmente l’avete anche eseguito (magari cambiando il numero di righe), e ne siete rimasti colpiti dall’eleganza, allora il merito andrà ad oggi, viceversa, andrà ad un magnifico giorno primaverile di 27 anni fa, quando alla lavagna cominciò a comparire un numero alla volta, il triangolo di Tartaglia.

Ciao, alla prossima.

METAPOST ed il prodotto scalare


Se la geometria aiuta

Recentemente riflettevo come raramente il grande pubblico televisivo sente parlare di matematica. Eppure è una scienza fondamentale per l’evoluzione della nostra specie. La tecnologia, dai telefonini alla robotica, la medicina, dagli scenari epidemici alla genetica, l’economia, dalla teoria dei giochi ai mercati finanziari, sono mosse dai numeri e dalle loro relazioni.

Prodotto scalare

Dunque, appena oltre la geometria elementare fatta di triangoli e circonferenze, si trova l’algebra lineare, ovvero il mondo dei vettori.
Una delle prime operazioni che si imparano dei vettori è quella del prodotto scalare, detto così perché riceve due vettori per fornire un numero, che è una grandezza scalare in distinzione da quelle vettoriali.

Si tratta di moltiplicare le lunghezze dei due vettori per il coseno dell’angolo da essi formato. Esiste anche una forma algebrica che consiste nella somma dei prodotti delle corrispondenti componenti dei vettori.

In formule, l’espressione geometrica è:

\displaystyle u \cdot v = uv\cos{\theta}

mentre quella algebrica è:

\displaystyle u \cdot v = u_1v_1+u_2v_2+u_3v_3

L’equazione della retta ed il prodotto scalare

Se l’angolo formato dai vettori è l’angolo retto, il coseno sarà zero e quindi anche il loro prodotto scalare.

Noto quindi un vettore n, i punti Q della retta ortogonale ad n e passante per il punto P soddisfano l’equazione:

\displaystyle n \cdot PQ=0

Se nel piano le componenti di n sono (a,b), quelle di PQ sono (x-x_p,y-y_p), l’equazione della retta si può scrivere nella forma lineare standard (assegnando a c la parte costante dell’espressione):

\displaystyle n \cdot PQ= a(x-x_p)+b(y-y_p)= ax+by+c=0

Applicazione in Metapost

La domanda è già nell’aria: perché non utilizzare le capacità di Metapost di risolvere sistemi lineari per rappresentare una retta a partire da un vettore?

(xpart r - x.P) * xpart n + (ypart r - y.P) * ypart n = 0;
 xpart r = 150;

Il codice completo è il seguente, in cui si è definito una funzione chiamata drawline, che prende un vettore (ed un colore), risolve il sistema come prodotto scalare e disegna il risultato con tanto di etichette.

outputtemplate:="scalar.mps";
beginfig(1);

% solve function to draw the line
def drawline(expr n, col)=
  pair r;

  % linear system
  (xpart r - x.P) * xpart n + (ypart r - y.P) * ypart n = 0;
   xpart r = 150;
  
  drawoptions(withcolor col withpen pencircle scaled 1.4);
  draw (-1.5)[z.P, r] -- 1.28[z.P, r];
  drawarrow origin -- n;
  drawoptions(withcolor black);
  
  % labels
  dotlabel.bot(btex $R$ etex, r);
  label.lft(btex $n$ etex,0.5[origin,n]);
enddef;

% punto di aggancio rette
z.P = (50,100);

pair a;
a := (0,75);

drawline(a rotated  20, blue);
drawline(a rotated -20, 0.7red);
drawline(a rotated -45, green);

dotlabel.bot(btex $P$ etex, z.P);

% ascissa limite
draw (150,-30)--(150,160) dashed evenly;

% draw axis
dim := 200;
dim1 := 12;

drawarrow origin -- (dim,0); % x
drawarrow origin -- (0,dim); % y

label.bot(btex $x$ etex,(dim-dim1,0));
label.lft(btex $y$ etex,(0,dim-dim1));
label.llft(btex $O$ etex,origin);

endfig;
end;
Scalar product: a Metapost example

Scalar product: a Metapost example

Punti e semipiani

Se pensiamo all’equazione della retta nella forma vista del prodotto scalare, è facile rendersi conto che è possibile determinare a quale dei due semipiani appartenga un punto qualsiasi.

Un punto al di sopra della retta (rispetto al vettore (a,b) chiamiamolo “generatore”), formerà un angolo con il vettore generatore inferiore all’angolo retto, per cui il prodotto scalare risulterà positivo. Viceversa se il punto è sotto la retta l’angolo è superiore all’angolo retto pertanto il prodotto scalare sarà negativo.

In definitiva, un punto di coordinate (x,y) sarà sopra o sotto la retta a seconda del segno del valore dell’espressione lineare ax+by+c.

Nello spazio

In “3D” l’espressione lineare ax+by+cz+d=0 è l’equazione del piano ortogonale al vettore (a,b,c), e sottoforma di prodotto scalare permette di capire come assegnare ad un punto qualsiasi il semispazio di appartenenza.
Non male per una banale moltiplicazione.

Ciao, alla prossima.

Metapost: proiezione punto su retta


Potenzialità lineari di Metapost

La sintassi di Metapost è da considerare un po’ datata ed a volte puntigliosa?
Meglio usare programmi di disegno preferibilmente vettoriale?

Oppure possiamo parlare la lingua di Metapost ragionando tranquillamente di relazioni geometriche, semipiani e vettori…

Allora, Metapost parlando, ci proponiamo in questo post di disegnare il punto proiezione di uno dato su un segmento, alla maniera di John Hobby, il creatore originale di questo potente tool basato su Postscript, entrando nei dettagli tecnici.

Un punto qualsiasi

In Metapost per dire che ci stiamo riferendo ad un punto qualsiasi (od anche ad un valore qualsiasi, dipende dal contesto), si usa il termine whatever. Con le informazioni disponibili, Metapost tenterà comunque di individuare un punto esatto (a meno delle approssimazioni numeriche), permettendoci di esprimere facilmente vincoli geometrici.

È anche possibile utilizzare una mediation expression per trovare punti su un segmento. Si premette un numero ai punti tra parentesi quadre, così il punto di mezzo del segmento AB è 0.5[A,B].
Se controllate, abbiamo già fatto uso della mediation in altri post, come per esempio quello riguardante il fiocco di Von Koch, per suddividere il lato in tre parti uguali.

Un punto qualsiasi sulla retta passante per i punti A e B è un’affermazione che si esprime con whatever[A,B].

Rette ortogonali

Che significa in Metapost sottrarre due punti? L’operazione sottrae le coordinate così l’espressione B-A possiamo pensarla come il vettore AB esattamente come se utilizzassimo la notazione comune in Meccanica Razionale, ma con il primo estremo nell’origine.

L’operazione di rotazione di 90° della differenza tra due punti, fornirà quindi il vettore uscente dall’origine e di direzione ortogonale alla retta passante per i due punti. Un breve esempio chiarificatore:

outputtemplate:="ortoa.mps";

beginfig(1);
% draw axis
numeric dim;
dim := 200;
drawarrow origin -- (dim,0); % x
drawarrow origin -- (0,dim); % y

pair a,b;
a := ( 10,32);  % se non sono indicate unità viene
b := (180,80);  % assunto il bp (big point) di Postscript

drawoptions(withcolor red withpen pencircle scaled 2);
drawarrow a--b;
drawarrow origin -- (b-a) dashed evenly;
drawarrow origin -- (b-a) rotated 90 dashed evenly;

% labels
drawoptions(withcolor blue);
dotlabel.bot("A",a);
dotlabel.bot("B",b);

label.bot(btex $x$ etex,(dim,0));
label.lft(btex $y$ etex,(0,dim));
label.llft(btex $O$ etex,origin);
endfig;
end;
A Metapost example on vector subtraction and rotation

A Metapost example on vector subtraction and rotation

Proiezione ortogonale (potenza di Metapost)

Chiamiamo C il punto da proiettare e P il punto di proiezione sulla retta AB.
Allora, se si vuole che il vettore di proiezione CP sia proporzionale al vettore della retta di proiezione scriveremo:
p-c = whatever*(b-a) rotated 90;
Nell’espressione whatever assume il ruolo di un fattore per condizione geometrica di parallelismo dei vettori.

Aggiungiamo il secondo vincolo, ovvero che il punto P appartiene alla retta AB:
p = whatever[a,b];

Adesso Metapost è in grado di risolvere il sistema di equazioni per trovare il punto P.

Codice completo e compilabile

Ricordo che l’output Metapost è direttamente utilizzabile da pdfLaTeX, quindi è possibile includerlo in documento per ottenere il pdf. Ecco quindi il codice immediatamente compilabile, dove la distanza di P dalla retta viene calcolata con la distanza pitagorica implementata in Metapost con l’operatore ++, e stampata sul terminale con il comando show:

outputtemplate:="ortob.mps";

beginfig(1);
% draw axis
numeric dim;
dim := 200;
drawarrow origin -- (dim,0); % x
drawarrow origin -- (0,dim); % y

pair a,b,c,p;
a := ( 25,32);  % se non sono indicate unità viene
b := (180,80);  % assunto il bp (big point) di Postscript
c := (80,180);

p-c = whatever * (b-a) rotated 90;
p = whatever[a,b];

drawoptions(withcolor red withpen pencircle scaled 2);
draw a--b;
draw p--c dashed evenly;

% labels
labeloffset := 6;
drawoptions(withcolor blue);
dotlabel.bot("A",a);
dotlabel.bot("B",b);
dotlabel.top("C",c);
dotlabel.bot("P",p);

label.llft(btex $x$ etex,(dim,0));
label.llft(btex $y$ etex,(0,dim));
label.llft(btex $O$ etex,origin);

numeric dist;
dist := xpart (c - p) ++ ypart (c - p);
show dist;%   -> 125.10623 bp
endfig;
end;
Metapost figure on point projection

Metapost figure on point projection

Saluti e buona estate!

La curva di Hilbert in METAPOST


Figura misteriosa…

The Hilbert Curve artwork

Si tratta di un quadro in stile geometrico di qualche artista contemporaneo da appendere a casa o al lavoro…? No, è solamente una rappresentazione delle curve di Hilbert dei primi tre livelli eseguita in METAPOST… ah ah ah.

Ancora una figura frattale…

Dopo il triangolo di Sierpinski ed il fiocco di Von Kock, consultabili qui, mi sono riproposto di disegnare la curva di Hilbert, ovvero una curva frattale in grado di riempire un intero quadrato.
Questo tipo di curve è stato studiato dal matematico italiano Giuseppe Peano ed esistono come curva limite di una successione infinita. Si può dimostrare infatti che la curva limite esiste come funzione e ricopre interamente un quadrato nel piano.

Le iterazioni della curva

Il primo passo per costruire la curva di Hilbert è disegnare tre lati consecutivi di un quadrato passanti per i quattro punti base. Su questa prima figura si disegnano gli stessi tre lati di lunghezza pari a metà della precedente, centrandoli sui quattro punti base, ruotandoli e collegandoli con tre segmenti come nella seguente figura.

The recursive drawing of the Hilbert Curve

Nell’immagine sono disegnati sovrapposti i primi due livelli, il primo in colore blu chiaro con i quattro punti base, il secondo in colore blu scuro, con i tre segmenti di collegamento in colore rosso scuro.

Per il terzo livello si procede allo stesso modo, considerando i nuovi punti basi del secondo livello e disegnando i relativi collegamenti. Nell’immagine di apertura, il primo livello è in colore verde, il secondo è in un rosso scuro ed il terzo è in blu.

Contiamo i lati ed i punti della curva

Ad ogni passo della successione i segmenti che compongono la curva si moltiplicano per quattro aggiungendosene tre per i collegamenti. In formule ciò si scrive (dove n è il numero del livello a partire da 1):

\displaystyle L(n)=3\sum_{i=0}^{n-1} 4^i = 4^n - 1.

Il numero dei punti della curva si ricava aggiungendo 1 al numero dei segmenti L(n)+1 = 4^n, dunque la crescita è esponenziale! All’ottava iterazione il numero dei punti è 65536!

Entra in scena METAPOST

Per disegnare i tre lati consecutivi dell’iterazione consideriamo il centro c ed il punto di mezzo del lato centrale x. Analizziamo il seguente codice:

beginfig(1)
  def trelati(expr c, x)=
    pair ax, bx, a, b;
    ax := c  rotatedabout(x, 90);
    bx := c  rotatedabout(x,-90);
    a  := ax shifted 2(c - x);
    b  := bx shifted 2(c - x);

    % disegno linee
    pickup pencircle scaled 5;
    draw a--ax--bx--b;
    drawarrow x--c withcolor .75blue;
  enddef;

   base = 100;
   pair cen, cx;
   cen := (base,base);  % centro
   cx  := (base,0);     % base

   trelati(cen, cx);
   trelati(cen shifted 3cx, cen shifted 4cx);
   trelati((base,-1.5base),(base,-.5base));
   trelati((4base,-1.5base),(3base,-1.5base));
endfig;
end
METAPOST Drawing test

METAPOST Drawing test

I quattro punti base ax, bx, a, b vengono definiti con trasformazioni geometriche dei punti centro e base c, x. Nel codice si mostra come la scelta del vettore che viene evidenziato dalla freccia blu, quindi l’ordine con il quale si passano gli estremi del vettore è essa stessa un informazione essenziale.

La ricorsione risolve la curva

La soluzione ricorsiva che vi presento la trova particolarmente elegante (ed adatta a chiudere l’anno), se non altro perché risolve un primo tentativo che oltre a non funzionare prevedeva altri due punti come argomento della funzione.

%
% Recursive Hilbert function
%

vardef hilbertCurve(expr c, x, level)=
    % the four basepoint
    save hilpath;
    path hilpath;
    save ax, bx, a, b;
    pair ax, bx, a, b;

    ax := c  rotatedabout(x, 90);
    bx := c  rotatedabout(x,-90);
    a  := ax rotatedabout(c ,-90);
    b  := bx rotatedabout(c , 90);

    if level = 1:
        hilpath := a--ax--bx--b;
    else:
        save atemp, btemp, axtemp, bxtemp;
        pair atemp, btemp, axtemp, bxtemp;

	atemp  := 0.25[a,b];
	btemp  := 0.25[b,a];
	axtemp := 1.25[a,ax];
	bxtemp := 1.25[b,bx];

        hilpath :=
        reverse hilbertCurve( a,  atemp, level-1) --
	        hilbertCurve(ax, axtemp, level-1) --
	        hilbertCurve(bx, bxtemp, level-1) --
        reverse hilbertCurve( b,  btemp, level-1);
    fi
    hilpath
enddef;

beginfig(1)
    pickup pencircle scaled 5;
    draw hilbertCurve(origin,(0,-100),4) withcolor red;
endfig;
end

La parola chiave reverse inverte il path e con questo si riesce a disegnare la curva di Hilbert posizionata comunque nel piano e del livello voluto (che dipende dalla memoria disponibile per METAPOST).

Scaricate pure questo file pdf della curva all’iterazione 5.

Buon anno a tutti!!!

Il fiocco di Von Koch in METAPOST


La curva di Von Koch

Il bellissimo fiocco di neve di Niels Fabian Helge von Koch, matematico svedese vissuto a cavallo del 1900, è una delle prime strutture frattali mai descritte.

Naturalmente è utile consultare la pagina web sull’argomento su wikipedia e quella dedicata da Wolfram MathWorld, ma certo non sarà difficile reperire in rete altre notizie se ancora non bastasse.

Si tratta di una curva chiusa di perimetro infinito ma di area finita, molto nota. Potevamo perdere l’occasione di generarla in METAPOST?

Effettivamente la curva si può costruire con il pacchetto METAOBJ, ma come problema ricorsivo non è male da implementare con quello che sappiamo già dai precedenti post.

La costruzione

Per ciascun lato di un triangolo equilatero si trovano i due punti che lo dividono in tre parti uguali, poi si costruisce il triangolo equilatero il cui lato è il segmento centrale che poi si elimina, e si reitera all’infinito su tutti i nuovi segmenti.

Building the Von Koch's curve

Building the Von Koch's curve

Come si vede dalla figura qui sopra il procedimento è semplice, ed avremo potuto introdurlo su un solo segmento. Si vedono le curve di Von Koch di livello 0, 1 e 2, con rispettivamente 3, 12, e 48 segmenti distinti.

METAPOST in action

Diciamo che nella variabile numerica lato memorizziamo la lunghezza del lato del triangolo di base che darà origine alla curva di Von Koch snowflake. Il codice per disegnare l’iterazione successiva deve trovare i due punti intermedi sul lato e costruire il triangolo equilatero sul terzo medio. Ecco il codice METAPOST:

beginfig(1);
numeric lato;
lato := 100;

pair a,b,p,q,r;

a := origin;
b := (lato,0);
p := 1/3[a,b];
q := 1/3[b,a];
r := q rotatedabout(p,-60);

draw a--p--r--q--b
withcolor .625red;
endfig;
end

Abbiamo fatto uso della comoda notazione per ricavare i punti P e Q che tripartiscono il segmento AB disteso sull’asse x e di lunghezza lato, e l’altrettanto comoda funzione rotatedabout() non nativa ma costruita con trasformazioni elementari, con la quale si crea il punto R, vertice del triangolo di prima iterazione.

Non rimane che costruire la funzione ricorsiva.
Faremo così: accetteremo in input i punti estremi del segmento ed il livello fino a cui la funzione disegnarà la curva di Koch, e su questo baseremo il controllo di terminazione della ricorsione:

se il livello richiesto è zero, allora disegna il segmento AB, altrimenti trova i punti intermedi P, Q ed R e lancia la funzione stessa sui quattro sotto lati AP, PR, RQ e QB.

Invece di elaborare oggetti path lavoriamo così più semplicemente sui punti vertice della curva. Lo stack delle chiamate ricorsive ci darà gratis l’ordine di disegno tra i punti fino al livello desiderato.

Il numero dei segmenti disegnati n è una funzione esponziale del livello v, si ha infatti (per il numero dei segmenti della la curva di Kock completa moltiplicare ulteriormente per 3):

\displaystyle n = 4^v

Ecco il codice completo:

beginfig(1);
pickup pencircle scaled .75pt;

vardef koch(expr aa,zz,level)=
save p,q,r;
pair p,q,r;

if level=0:
draw aa -- zz
withcolor .625red;
else:
p := 1/3[aa,zz];
q := 1/3[zz,aa];
r := q rotatedabout(p,-60);

koch(aa, p,level-1);
koch( p, r,level-1);
koch( r, q,level-1);
koch( q,zz,level-1);
fi
enddef;

numeric lato ;
numeric n;

lato := 420pt;
n := 5;
koch( origin   , (lato,0)     , n);
koch( (lato,0) , lato*dir(60) , n);
koch( lato*dir(60) , origin  , n);

endfig;
end

dove abbiamo utilizzato una macro vardef dando istruzione di utilizzo di p, q, ed r come variabili locali (comando save).

Pronti a scaricare il risultato nel formato PDF? Fate click allora sull’immagine sottostante.

A Von Koch flake with level 5

A Von Koch flake with level 5

Alla prossima.

Il triangolo di Sierpiński con METAPOST


La ricorsione è sempre la via migliore?

A not recursive algorithm for the Sierpinski fractal triangle

A not recursive algorithm for the Sierpinski fractal triangle

Sembra che il miglior modo di affrontare una struttura ricorsiva come il triangolo del matematico polacco sia appunto un algoritmo ricorsivo.

Tuttavia leggendo il codice relativo nel post predecente ho cominciato a ragionare sulla possibilità di chiamare ricorsivamente la funzione di disegno nominata sierpinski() una sola volta anziché tre volte, una per ciascun triangolo del livello inferiore.

L’idea che mi è venuta è quella di disegnare il triangolo dal basso, triplicando per ciascun livello successivo il disegno. Tuttavia questo si può elegantemente fare con un ciclo iterativo semplice!

Dunque vi propongo una seconda soluzione (scritta in METAPOST) in cui il disegno è ottenuto con un normalissimo ciclo for.

Il codice spiegato passo passo

Come d’abitudine propongo l’implementazione dell’idea sviluppandola di pari passo con il codice, seguito dal disegno che ne deriva.

Per prima cosa disegniamo il primo triangolo equilatero (per una breve e mirata introduzione a METAPOST potete far riferimento utile al post precedente su questo stesso blog).

beginfig(1);

u = 40 mm;
fill origin -- (u,0) -- u * dir(60)-- cycle
withcolor .625red;
drawarrow origin -- u * dir(60)
withpen pencircle scaled 5 pt;
endfig;
end

An example with the dir() METAPOST function

An example with the dir() METAPOST function

Impostare una variabile unità qui settata a 40 mm è molto comodo, ma quello che dovreste notare di questo codice è l’uso della funzione dir() per fornire al path le coordinate del vertice alto del triangolo equilatero. La funzione restituisce un vettore unitario che forma l’angolo specificato con l’asse x di senso antiorario con centro nell’origine.
Questo vettore viene visualizzato con il comando drawarrow seguente.

METAPOST disegna il vettore restituito da dir() come un punto corrispondente al secondo estremo dello stesso. Del resto sarebbe scomodo se nel path andasse il vettore piuttosto che il punto.

Nel prossimo esempio possiamo vedere l’uso del costrutto for impiegando la funzione dir() a riprova che dir(angolo) produce il disegno di un punto:

beginfig(2);
numeric u;
u = 40 mm;

for ang = 0 step 15 until 345:
draw u * dir(ang)
withpen pencircle scaled 12 pt;
endfor;
endfig;

A polar example with dir() METAPOST function

A polar example with dir() METAPOST function

Per andare avanti sul nostro obbiettivo Sierpinski, occorre spiegare il concetto di picture. I comandi grafici non fanno altro che plottare disegni in una picture. Ci sono sempre disponibili due oggetti picture, quella corrente detta currentpicture e quella vuota detta nullpicture. Un po’ di codice chiarirà di cosa si tratta:

beginfig(3);
% set the thin of the pen
pickup pencircle scaled 3pt;

draw fullcircle scaled 80pt
withcolor 0.5 green;

picture pic;
pic := currentpicture;
currentpicture := nullpicture;

draw pic;
draw pic shifted (80pt,0);
endfig;

The METAPOST picture concept at work

The METAPOST picture concept at work

All’oggetto pic di tipo picture viene assegnato il disegno corrente, che in quel momento contiene un cerchio verde, poi il disegno viene cancellato. In questo momento il cerchio verde si trova memorizzato in un solo luogo: la picture pic.
Allora possiamo tornare a riempire currentpicture disegnando prima pic nella sua posizione e poi in quella traslata a destra.

L’idea iterativa Sierpinski

% set the unit size
u = 10 mm;

% a picture object
picture sier;

sier:= currentpicture;
draw sier shifted (u,0);
draw sier shifted (dir(60) scaled u);

The basic idea for a iterative Sierpinski drawing method

The basic idea for a iterative Sierpinski drawing method

Questo frammento di codice è il cuore del programma che disegna il triangolo frattale in maniera iterativa e non ricorsiva.
Si presuppone che la currentpicture contenga all’inizio il triangolino di base molto piccolo. Dapprima si memorizza il disegno corrente nella variabile sier, di tipo picture naturalmente, e semplicemente l’aggiungiamo al disegno corrente stesso prima a destra e poi in alto a 60° al posto giusto.

In altre parole come si vede in figura, il triangolo rosso viene memorizzato in una picture e poi viene riaggiunto due volte nella posizione dei triangoli blu. Basterà iterare all’infinito questo procedimento!

Il ciclo for Sierpinski

beginfig(4);
% set the unit size
u = 0.65 mm;

% draw the basic smallest triangle
fill origin--(u,0)-- dir(60) scaled u --cycle
withcolor .625red;

% a picture object
picture sier;

% recursive level
level = 8;

% main drawing cycle
for i=0 upto level-1:
tmp := u * (2**i);
sier:= currentpicture;
draw sier shifted (tmp,0);
draw sier shifted (dir(60) scaled tmp);
endfor;

endfig;

Adesso il codice è completo. L’operazione precedente di copia a destra ed a 60° è inserita nel ciclo iterativo for che produce il triangolo del livello desiderato (ho scelto il livello 8).
Aggiungo che ad un dato livello i la copia del disegno deve esser fatta ad una distanza di u \cdot 2^i, poichè ogni volta che si sale di un livello, il lato del triangolo raddoppia (il doppio asterisco significa quindi elevamento a potenza).

Alla prossima!

Il triangolo di Sierpiński


Il disegno della ricorsione

Visualizzare la ricorsione è l’idea di questo post. E la figura ricorsiva più semplice da disegnare mi sembra proprio il triangolo di Sierpiński. Dopo qualche notizia essenziale sul software che userò, passeremo al codice…

Sierpinski Triangle

Sierpinski Triangle of 6 level (click to download the pdf file)

Il linguaggio grafico di Asymptote

Per le nostre sperimentazioni useremo un software chiamato Asymptote, che consiste in un linguaggio di alto livello derivato dal METAPOST, piuttosto espressivo, simile a C++ e Java, quindi fortemente tipizzato.

Si possono instanziare i classici tipi numerici come int o real, e produrne di nuovi per mezzo del costrutto struct.

Infine, occorre terminare le istruzioni con il punto e virgola, utilizzare le parentesi graffe per definire i blocchi di codice, e produrre i commenti nel formato C++, ma vediamo subito un esempio (supporremo che una versione di Asymptote sia installata sul vostro sistema Windows, Linux o Mac OS che sia, magari tramite la TeX Live 2009).

Per disegnare una spezzata scrivete in un file di testo con estensione .asy il seguente codice, dove le coordinate sono separate da un doppio trattino che significa come in molti altri linguaggi del mondo LaTeX, una linea retta:

draw((0,0)--(2.5,2.5)--(5,0.5)--(7.5,5));

Poi in una finestra di console, date il comando di compilazione: asy -f pdf nomefile.asy. Al termine dell’elaborazione verrà rilasciato un file pdf contenente il risultato grafico desiderato.

Qualcosa di più complesso

Cominciamo a disegnare quello che sarà il triangolo di livello 0, utilizzando il tipo pair ovvero il tipo che in Asymptote corrisponde ad una coppia di numeri reali:

// triangle line length
real ell = 100;

// level 0 triangle:
pair p1 = (0,0);
pair p2 = (ell,0);
pair p3 = (ell/2,sqrt(3)*ell/2);
draw(p1--p2--p3--cycle);

Il termine cycle specifica che il percorso da disegnare sarà chiuso e corrisponde di fatto al primo punto p1.

Sierpiński

Il triangolo del matematico polacco venne da lui elaborato nel 1915, quando i sovietici tentavano di impossessarsi del cuore di quel popolo.

Disegnato un triangolo (che supporremo equilatero) detto di livello 0, si disegna il triangolo i cui vertici sono i punti di mezzo dei lati del triangolo principale, poi si applica ricorsivamente lo stesso procedimento ai tre sotto triangoli nord, sud-est e sud-ovest.

La struttura è un frattale in quanto è autosomigliante. Una semplice osservazione ci dice che man mano che si scende di livello le lunghezze dei lati dei triangoli si dimezzano, per cui se l_0 è la misura del lato del triangolo del livello 0, quella del lato a livello n sarà:

\displaystyle l_n = \frac{l_0}{2^n}.

Il numero dei triangoli t in funzione del livello v si ottiene considerando che scendendo di un livello il loro numero cresce di 3^v-1. In formule:

\displaystyle t = 1 + \sum_{i=1}^v 3^{i-1}

Per esempio al livello 5 i triangoli disegnati sono 122, ma al 10° salgono a 29.525 ed al 14° sono ben 2.391.485.

Il codice (non ottimizzato)

Per terminare i cicli ricorsivi altrimenti infiniti, è possibile utilizzare un criterio grafico: se lo spessore di tracciamento delle linee è tale che i triangoli dei livelli ulteriori risultano coperti, è inutile proseguire con il disegno.

La funzione ricorsiva chiamata sierpinski() prende i tre punti vertici del triangolo ed effettua immediatamente il test seguente (s è lo spessore del tratto ed l la lunghezza del lato del triangolo):

\displaystyle \frac{s}{2}\geq \frac{l}{2\sqrt{3}}

Se il test è falso significa che possiamo procedere al disegno del sotto-triangolo con i vertici nei punti di mezzo dei lati del triangolo principale, ed incrementare l’intero che conterrà alla fine il numero totale dei triangoli disegnati.

Il passo finale è lanciare ricorsivamente tre volte la funzione stessa sui sotto-triangoli del livello successivo.

Il codice non è ottimizzato perché in realtà non serve calcolare la distanza tra due punti per ottenere la lunghezza del lato del triangolo del livello attuale perché è ottenibile semplicemente dimezzando quella corrispondente al livello precedente.

Ecco il listato:

//
// Sierpinsky triangle: draw it with a recursive function
// Copyright (c) 2009 Roberto Giacomelli
// Term of licence see: https://robitex.wordpress.com/legalese/
//

// var setup
int n = 1;
real l0 = 100;
real s = 1.0;

// set the linewidth
defaultpen(s);

// main triangle
pair p1 = (0,0);
pair p2 = (l0,0);
pair p3 = (l0/2,sqrt(3)*l0/2);
draw(p1--p2--p3--cycle);

// function med() return the midpoint
// of the line ab
pair med(pair a, pair b){
return ((a.x+b.x)/2 , (a.y+b.y)/2);
}

// function dist() computes the length
// of line ab
real dist(pair a, pair b){
return sqrt((a.x-b.x)^2+(a.y-b.y)^2);
}

// main recursive function
void sierpinski(pair a, pair b, pair c){
// check graphic stop condition
if ( s<= dist(a,b)/sqrt(3) ){return;}

// draw subtriangle of the level
pair am = med(b,c);
pair bm = med(a,c);
pair cm = med(a,b);
draw(am--bm--cm--cycle);
++n;

//well, and now recursive running
// sublevel south east triangle
sierpinski(a,cm,bm);
// sublevel south west triangle
sierpinski(cm,b,am);
// sublevel north triangle
sierpinski(c,bm,am);
}
// draw the Sierpinski figure
sierpinski(p1,p2,p3);

// write on console the number
// of drawed triangles
write(n);

Scarica l’esempio in pdf : Sierpinski Triangle Figure
Alla prossima!

%d blogger hanno fatto clic su Mi Piace per questo: