Iterative Methods and Preconditioning for Large and Sparse Linear Systems with Applications

È stata una fatica improba, ma finalmente è uscito!

Dalla quarta di copertina:

This book describes, in a basic way, the most useful and effective iterative solvers and appropriate preconditioning techniques for some of the most important classes of large and sparse linear systems.

The solution of large and sparse linear systems is the most time-consuming part for most of the scientific computing simulations. Indeed, mathematical models become more and more accurate by including a greater volume of data, but this requires the solution of larger and harder algebraic systems. In recent years, research has focused on the efficient solution of large sparse and/or structured systems generated by the discretization of numerical models by using iterative solvers.

Potete acquistarlo su Amazon

Matematica e modelli biomedici, un cammino appena iniziato – INCA-ABACO

Ecco il secondo seminario del ciclo:

Infrastrutture di Calcolo A Basso Costo

Negli ultimi anni le interazioni tra i matematici e le scienze biomediche sono cresciute enormemente. Le sfide poste dalla ricerca di cuincontro22marzore per malattie come i tumori o i disturbi cardiovascolari necessitano di modelli sempre più precisi per migliorarne la comprensione e ottimizzare le possibili risposte terapeutiche. Inoltre, le possibilità aperte dall’utilizzo di cellule staminali, apre la porta a nuovi scenari in cui i modelli matematici potranno aiutare i medici nella rigenerazione dei tessuti lesi.
In questa presentazioni, alcuni casi concreti di interazione tra modelli matematici e applicazioni biomediche verranno descritti.

Roberto Natalini
Direttore dell’Istituto per le Applicazioni del Calcolo “M. Picone” – C.N.R.
22 Marzo 2016
ore 16.00
Aula 12
So.Ge.Ne. – Via Della Ricerca Scientifica 1, Roma, RM

Buon 3/14 alle 9.26.53 (PM)!

Ecco un paio di modi per approssimare il valore di \pi.

Partiamo da un classico, ma non troppo, che è detta serie di Leibniz-Gregory-Madhava (a seconda delle latitudini e delle epoche). Ovvero lo sviluppo in serie di potenze della funzione \arctan(x), per cui si ha che:

(1)   \begin{equation*} \pi = 4 \sum_{k=0}^{+\infty} \frac{(-1)^k}{2k+1}. \end{equation*}

Che possiamo calcolare e stampare con il seguente codice OCTAVE.

printf('Approssimazione n° 1\n');
N = input('Sommare fino a N = ');
piapp1 = 0;
figure(1)
for k=0:N
  piapp1 = piapp1 + 4*((-1)^k)/(2*k+1);
  subplot (211);
  hold on
  plot(k,piapp1,'*');
end
plot(0:0.1:N,pi*ones(1,length(0:0.1:N)),'k--');
ylabel('Valore approssimato di \pi');
xlabel('Numero di termini della serie');
title(sprintf('Approssimazione n° 1, N = %d',N));
hold off

Che ha una convergenza molto, molto, lenta … per vederlo servono i numeri di Eulero, ma si può comunque fare qualche esperimento con il codice (un po’ goffo sotto molti punti di vista) riportato sopra.

Possiamo scegliere un’altra via, se sappiamo che

(2)   \begin{equation*} \int \sqrt(1-t^2)dt = \frac{1}{2}\left(\sqrt{1-t^2}+\arcsin(t)\right) + C, \end{equation*}

abbiamo che

(3)   \begin{equation*} \pi = 4 \int_{0}^{1} \sqrt{1-t^2}dt. \end{equation*}

Cosa ce ne facciamo? Be’ possiamo approssimare l’integrale con la regola dei trapezi … Ovvero discretizziamo l’intervallo [0,1] in N+1 punti equispaziati \{x_k\}_{k=1}^{N+1}, ovvero fissiamo una griglia su [0,1] di ampiezza h = (1-0)/N, per cui l’approssimazione dell’integrale diventa:

(4)   \begin{equation*} \pi = 4 \int_{0}^{1} \sqrt{1-t^2}dt \approx \frac{1}{2N} \sum_{k=1}^N \left(\sqrt{1-(x_{k+1})^2}-\sqrt{1-(x_k)^2}\right). \end{equation*}

Si lo so, bisognerebbe pure sostituire quelle radici con un’approssimazione, ma facciamo finta di averlo fatto, o di averlo delegato ad OCTAVE.

printf('Approssimazione n° 2\n');
M = input('Numero dei punti in cui suddividere [0,1], M = ');
for h=1:M
  x = 0:1/h:1;
  y = sqrt(1-x.^2);
  piapp2 = 4*trapz(x,y);
  subplot (212);
  hold on
  plot(h,piapp2,'*');
end 
plot(0:0.1:M,pi*ones(1,length(0:0.1:M)),'k--');
ylabel('Valore approssimato di \pi');
xlabel('Numero di intervalli di [0,1]');
title(sprintf('Approssimazione n° 2, M = 1,...,%d',M));
hold off

Il risultato dei due codici precedenti è questo:

 

Ah già, buon piday!


Codice dell’esempio: pigreco.m.tar.gz.

L’immagine di copertina viene dall’illustrazione all’inizio degli Elementi di Euclide nella traduzione di Adelardo di Bath (fonte: Wikipedia/Wikimedia).

Numeri primi – #2

Non solo la matematica è reale, ma è l’unica realtà.
Martin Gardner

Dunque dunque, con il numero 1 avevamo messo da parte alcune dimostrazioni dell’infinità dei numeri primi. A questo punto, direi, che quel materiale ha fermentato abbastanza e possiamo aggiungere altra carne al fuoco (apprezzare la doppia metafora non coerente, con slittamento banale a vanvera).

Dopo la dimostrazione di Erdös, facciamo un salto indietro nel tempo e torniamo ad una dimostrazione in qualche modo classica, prelevata dalle idee di Pierre de Fermat. Per farlo cominciamo a dare la seguente definizione:

Def 3. Per \forall n \in \NN definiamo Numero di Fermat l’intero positivo della forma:

(1)   \begin{equation*} F_n = 2^{2^n} + 1 \end{equation*}

Aprendo una piccola cornice storica, Fermat ipotizzo che tutti i numeri di questa forma fosse primo, ma un matematico a caso, Eulero, tanto per dirne uno, osservò che:

  • Per n = 0, F_0 = 3, è primo!
  • Per n = 1, F_1 = 5, è primo!
  • Per n = 2, F_2 = 17, è primo!
  • Per n = 3, F_3 = 257, è primo!
  • Per n = 4, F_4 = 65537, è primo!
  • Per n = 5, F_5 = 4294967297, non è primo neanche per niente: F_5 = 641 \cdot 6700417

Chiusa la parentesi storica, torniamo al nostro problema ed enunciamo (e dimostriamo) il seguente teorema:

Teorema [Goldbach]. \forall n,n' \in N con n \neq n' si ha \operatorname{M.C.D.}(F_n,F_{n'})=1, ovvero due numeri diversi di Fermat sono sempre coprimi.

Dimostrazione. Proviamo in primo luogo, per induzione, che vale la relazione:

(2)   \begin{equation*} F_n = F_0 \cdot \ldots \cdot F_{n-1} + 2 \;\; \forall n \geq 1 \end{equation*}

per n = 1 abbiamo che:

    \begin{equation*} F_1 = 2^{2^1} + 1 = 5 = 2^1 + 1 + 2 = 5 . \end{equation*}

questo prova la base, supponiamo ora che sia vera la relazione per n e dimostriamolo per n+1, cioè vogliamo provare che:

    \begin{equation*} F_0 \cdot F_1 \cdot F_2 \cdot \ldots \cdot F_{n-1} \cdot F_n + 2 = F_(n+1). \end{equation*}

Per farlo poniamo:

    \begin{eqnarray*} x = F_0 \cdot F_1 \cdot F_2 \cdot \ldots \cdot F_{n-1}\\ y = F_0 \cdot F_1 * \cdot \ldots \cdot F_(n-1) \cdot F_n \end{eqnarray*}

Ci siamo dunque ridotti a provare che:  y + 2 = F_{n+1}. Osserviamo immediatamente che, per l’ipotesi induttiva si ha che:  y = x (x+2). E dunque ci siamo ridotti a mostrare che:  x (x+2) + 2 = F_{n+1}. Adesso, applicando di nuovo l’ipotesi induttiva abbiamo che:  x = F_{n - 2}, allora:

    \begin{equation*} \begin{split} x  (x+2) + 2 = & [(F_{n - 2})  ((F_{n - 2} + 2)] + 2 = \\ = & [(F_{n - 2}  F_n] + 2 =  \\ = & (F_n)^2 - 2F_n + 2 = \\ = & (2^{2^n} + 1)  (2^{2^n} + 1) - 2  (2^{2^n} + 1) + 2  =  \\ = & 2^{2^n}  2^{2^n} + 2\cdot 2^{2^n} + 1 - 2 \cdot 2^{2^n} - 2 + 2 = \\ = & 2^{2^n} \cdot 2^{2^n} + 1 = 2^{2^n + 2^n} + 1 = \\ = & 2^{2\cdot (2^n)} + 1 = \\ = & 2^{2^{n+1}} + 1 = F_{n+1} \end{split} \end{equation*}

Con questo abbiamo praticamente completato la dimostrazione, ora ci resta solo da supporre,per assurdo, che per 0 \leq i < j \exists \; a > 1 tale che a | F_i e a | F_j, per l’identità che abbiamo appena provato si ha che a deve necessariamente dividere F_0\cdot \ldots \codt F_{j-1}, sia F_j, cioè deve dividere la loro differenza, ovvero a deve dividere 2, ma questo implica che a = 2, ma questo è assurdo! I numeri di Fermat sono, chiaramente, tutti dispari.

Q.E.D.

A questo punto la dimostrazione dell’infinità dei primi è presto ottenuta:

Dimostrazione 3. Per il teorema fondamentale dell’aritmetica abbiamo che \forall n \in N esiste un primo p_n | F_n e per il teorema di Goldbach, questo p_n non divide nessun F_{n'} per n \neq n', dunque il numero dei p_n e maggiore o uguale del numero degli F_n, ma questi sono infiniti, dunque abbiamo concluso.

Q.E.D.

Proseguiamo con l’ultima dimostrazione dell’infinità dei primi di questa carrellata. Questa è a base di topologia ed è stata proposta da Harry Fürstenberg nel 1955 [1], di tutte quelle viste è l’unica che probabilmente richiede strumenti di matematica avanzata (diciamo extra-liceo):

Dimostrazione 4. Sia \ZZ l’insieme dei numeri interi, positivi, negativi e lo 0. Consideriamo gli elementi a,b \in \ZZ con b>0 e definiamo gli insiemi a due indici:

    \begin{equation*} N_{a,b} = \{ a + nb \;|\; n \in \ZZ\} \end{equation*}

Osserviamo immediatamente che N_{a,b} è una progressione aritmetica che si sviluppa in ambo i versi e che, inoltre, \forall b si ha che a \in N_{a,b}, il che dovrebbe iniziare a ricordarci l’idea di un intorno in senso topologico, infatti abbiamo che ogni insieme del tipo N_{c,d} \supset N_{a,b} contiene ancora a ed è dunque ancora un suo intorno. Inoltre ogni intorno di a contiene un intorno di a che è anche intorno di ognuno dei suoi punti. Resta solo da verificare che l’intersezione di due intorni di un medesimo a è ancora un intorno di a, ma anche questo è semplice da verificarsi, infatti dati a,b,c \in \ZZ con b,c > 0 abbiamo che:

    \begin{equation*} N_{a,b} \cap N_{a,c} = N_{a,\operatorname{m.c.m.}(b,c)} \end{equation*}

Possiamo quindi costruire su \ZZ la seguente topologia:

    \begin{equation*} \mathcal{T} = \{U \subseteq \ZZ \,|\, \exists a \in U, \; b \in \ZZ_{>0} \; N_{a,b} \subseteq U \}\cup\emptyset \end{equation*}

di cui osserviamo le due seguenti proprietà:

  1. Ogni insieme aperto U \neq \emptyset è infinito.
  2. Ogni insieme aperto U è anche chiuso, infatti:

        \begin{equation*} N_{a,b} = \ZZ \setminus \bigcup_{i=1}^{b-1}N_{a+i,b} \end{equation*}

    cioè è chiuso poiché complementare di un aperto (dell’unione di aperti …).

possiamo a questo punto concludere, infatti per il teorema fondamentale dell’aritmetica ogni intero, ad eccezione di \pm 1 e 0 ha dei fattori primi. Dunque ogni intero è contenuto in uno o più N_{0,p} con p primo. Cioè abbiamo ottenuto l’identità:

(3)   \begin{equation*} \ZZ \setminus \{1,-1\} = \bigcup_{p \text{ primo}}N_{0,p} \end{equation*}

Se l’insieme dei numeri primi fosse, per assurdo, finito l’insieme a destra dell’uguaglianza sarebbe chiuso, poiché unione finita di chiusi (proprietà 2), dunque l’insieme \{-1,1\} sarebbe aperto, ma questo contraddice la (proprietà 1). L’assurdo deriva dall’aver supposto finito l’insieme dei numeri primi.

Q.E.D.

E pure questa puntata è finita. Sto meditando su un eventuale “Numeri primi – #3“, ma ancora non ho raggiunto una decisione. Prendiamo pure questo come mezzo annuncio e salutiamoci alla prossima!

Bibliografia.

  1. Harry Fürstenberg, On the infinitude of primes., Amer. Math. Monthly 62 (5): 353 (1955). DOI:10.2307/2307043.

Numeri primi – #1

I Matematici hanno provato fino ad oggi, in vano, di scoprire un qualche ordine nella sequenza dei numeri primi, e ho ragione di credere che questo è un mistero in cui la mente umana non potrà mai penetrare.

Leonhard Euler (1707-1783) 

Per oscuri motivi universitari, o forse non poi così oscuri, mi sono messo a cercare tra i risultati di Teoria dei Numeri sui primi e visto che ho messo da parte qualcosa di interessante, almeno secondo il mio standard, ho deciso di proporveli, a rate chiaramente.

Iniziamo con il ricordare la definizione di numero primo:
Def. 1 Un intero n \in \NN è detto primo se n > 1 e se gli unici divisori positivi di n sono 1 ed n stesso. Altrimenti il numero è detto composito.

La prima questione da investigare è: ma quanti sono in realtà i numeri primi? Finiti? Infiniti? 42? Prima di poter rispondere, a fondo, a questa domanda, come si fa quasi sempre in matematica, abbiamo bisogno di introdurre un paio di risultati preliminari che ci aiuteranno a risolvere.

Lemma 1. (La serie armonica) La serie \sum_{i=1}^{+\infty}\frac{1}{n} è divergente.

Dim. È sufficiente osservare la seguente relazione per le somme parziali:

    \begin{equation*} \begin{split} S_N = & \sum_{i=1}^{N}\frac{1}{n} = \sum_{i=1}^{N}\int_{n}^{n+1}\frac{1}{n}dx > \\      >& \sum_{i=1}^{N}\int_{n}^{n+1}\frac{1}{x}dx = \int_{1}^{N+1}\frac{1}{x}dx = \log(N+1) \end{split} \end{equation*}

e quindi S_N > \log(N+1) \rightarrow +\infty, ovvero la serie diverge.

Q.E.D.

Osserviamo che, inoltre, questo ci dice che la serie armonica diverge almeno in modo logaritmico. In qualche puntata successiva si potrebbe far vedere che in realtà diverge proprio in modo logaritmico.

Diamo adesso il seguente risultato sui numeri naturali:
Lemma 2. Ogni numero intero n > 1 è un numero primo oppure il prodotto di numeri primi.

Dim. Dimostriamolo per induzione su n. Il teorema è banalmente valido per n=2, supponiamolo vero per tutti gli interi <n e dimostriamolo per n. Se n è primo abbiamo finito, altrimenti n ha un divisore proprio d con d\neq 1, quindi n = dc con d<n e c<n, inoltre sia c > 1 sia d > 1, dunque possiamo applicare l’ipotesi induttiva e concludere.

Q.E.D.

Dati questi due ingredienti possiamo finalmente proporre un po’ di dimostrazioni del seguente teorema:

Teorema 1. I numeri primi sono infiniti.

Dimostrazione 0. [Euclide] Supponiamo per assurdo che i numeri primi siano finiti e siano rappresentanti dalla lista \{p_1,p_2,\ldots,p_n\}, consideriamo allora il numero M = p_1\cdot p_2 \cdot\ldots\cdot p_n, chiaramente questo numero è divisibile per ognuno dei p_i. Consideriamo invece il numero N = M + 1, questo non può essere divisibile per nessuno dei \{p_1,p_2,\ldots,p_n\}, ma questa era una lista esaustiva dei primi. Abbiamo raggiunto l’assurdo e l’assurdo deriva dall’aver supposto finiti i numeri primi.

Q.E.D.

Dimostrazione 1. [Euler]  Supponiamo per assurdo che i numeri primi siano finiti e siano rappresentanti dalla lista \{p_1,p_2,\ldots,p_n\}. A questo punto consideriamo la serie armonica:

(1)   \begin{equation*} \begin{split} \sum_{n=1}^{+\infty}\frac{1}{n} = & 1 + \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + \frac{1}{5} + \cdots = \\ = & \left(1 + \frac{1}{p_1} + \frac{1}{p_1^2} + \frac{1}{p_1^3} + \cdots + \frac{1}{p_1^n} + \cdots \right)\cdot \\ & \left(1 + \frac{1}{p_2} + \frac{1}{p_2^2} + \frac{1}{p_2^3} + \cdots + \frac{1}{p_2^n} + \cdots \right)\cdot \\ & \vdots \\ & \left(1 + \frac{1}{p_n} + \frac{1}{p_n^2} + \frac{1}{p_n^3} + \cdots + \frac{1}{p_n^n} + \cdots \right) \end{split} \end{equation*}

abbiamo scritto la serie armonica come il prodotto delle serie geometriche di ragione 1/p_i \forall i=1,\ldots ,n poiché ogni elemento della serie armonica è il reciproco di un naturale e, per il lemma 2, ogni naturale è il prodotto di numeri primi. Dunque svolgendo tutti i prodotti dal secondo lato dell’uguaglianza otteniamo tutti gli elementi della serie armonica. A questo punto è sufficiente osservare che tutte le serie geometriche a secondo membro sono convergenti, dalla nota formula abbiamo che:

(2)   \begin{equation*} \sum_{n=1}^{+\infty}\frac{1}{n} = \frac{1}{1 - \frac{1}{p_1}} \cdot \frac{1}{1 - \frac{1}{p_2}} \cdot \cdots \cdot \frac{1}{1 - \frac{1}{p_n}} \end{equation*}

ma questo è assurdo, infatti dal lemma 1 abbiamo che la serie armonica è divergente, mentre dall’equazione 2 abbiamo che è uguale ad un prodotto finito. L’assurdo deriva dall’aver supposto finiti i numeri primi.

Q.E.D.

Un’altra, elegante, dimostrazione di questo fatto è quella data da Erdös. Introduciamo prima la seguente funzione:
Def 2. Funzione di distribuzione dei numeri primi:

(3)   \begin{equation*} \pi(x) = \{\text{numero dei numeri primi } < x \} \end{equation*}

Possiamo dunque dare la prova:

Dimostrazione 2. [Erdös]  Osserviamo, da principio, che ogni numero intero n \in \NN può essere rappresentanto nella forma n = rs^2 con r \in N privo di quadrati ed s intero qualsiasi. Ad esempio basta considerare s il più grande intero tale s^2 | n e porre r = n/s^2. Chiaramente questa rappresentazione non è unica, per ottenere una sovrastima del numero di partizioni di questo tipo abbiamo, in primo luogo, bisogno di sapere quanti sono i numeri r < n privi di quadrati. Ognuno di questi corrisponde ad una serie di primi distinti il cui prodotto è < n e di questi primi una soprastima è data da \pi(n), dunque ci sono al più 2^{\pi(n)} numeri privi di quadrati \leq n. La seconda domanda a cui dobbiamo rispondere per ottenere la stima è più semplice, infatti il numero interi s il cui quadrato è minore di n è semplicemente \sqrt{n}. Riassumendo se fattorizziamo ogni numero \leq n nella forma rs^2 abbiamo al più 2^{\pi(n)} possibilità per r e \sqrt{n} possibilità per s, ovvero abbiamo ottenuto che:

(4)   \begin{equation*} 2^{\pi(n)}\sqrt{n} \geq n \end{equation*}

dividendo per \sqrt{n} ambo i membri ed estraendone i logaritmi in base 2:

(5)   \begin{equation*} \pi(n) \geq \frac{\log n}{2} \end{equation*}

dunque \pi(n) \rightarrow +\infty per n\rightarrow \infty, ovvero i numeri primi sono infiniti.

Q.E.D.

Si osservi che questa dimostrazione ci ha fornito anche una stima sulla distribuzione dei numeri primi data dall’eq. 5, a questo punto conviene farne un disegnino tanto per capire ad occhio come è la stima:

[singlepic id=454 w=430 h=326 float=center]

si può fare di meglio, molto di meglio in effetti.

A questo punto direi che per il #1 basta e avanza. Al prossimo giro con un altro paio di dimostrazioni del Teorema 1. e qualche altra cosuccia.

Bibliografia.

  1. Euclide, Elementi, Proposizione IX.20
  2. Apostol, T.M., Introduction to Analytic Number Theory, Springer, 1976.
  3. P. Erdös, “Uber die Reihe \sum \frac{1}{p}“, Mathematica, Zutphen B 7 (1938).