Manuale di Mathematica 04

Daniele Lupo Appunti di Mathematica

Calcolo numerico

Introduzione

Quello che abbiamo visto finora riguarda principalmente il calcolo simbolico. Anche se è uno degli elementi che contraddistinguono questo programma dalla massa, sia per il numero di funzioni che per la loro potenza, le caratteristiche del programma non si fermano certo in quest'aspetto, anzi....

Mathematica contiene una quantità spropositata di algoritmi dedicati per poter risolvere la totalità dei problemi numerici che possono capitare a questo mondo, ed in questo universo. Le caratteristiche fondamentali sono la scelta automatica dell'algoritmo e la precisione arbitraria. Con il primo punto, possiamo stare (quasi) sempre tranquilli e lasciare decidere al programma il metodo migliore per la risoluzione di un problema: per esempio, per un integrale numerico ci sono i metodi di Simpson, quello dei trapezi e molti altri che neanche io conosco... Mathematica riesce ad analizzare la funzione, ed è in grado di decidere da solo quale metodo utilizzare per avere la soluzione migliore al problema.

Chi ha usato qualche volta, in elettronica il programma SPICE (quello originale della Berkley, non uno dei tanti programmi presenti con grafici etc), sa che vengono usati due principali algoritmi per la risoluzione delle formule integrali: il metodo trapezoidale e quello di Gear. Sono due differenti metodi per trovare integrali numerici. Di solito le formule di Gear danno risultati più precisi, ma hanno problemi di stabilità: per sistemi stiff, cioè con poli lontani dall'origine, può capitare che il metodo diverga anche quando in realtà converge, e, in minor parte, convergere quando in realtà diverge, dato che il sistema teorico e quello numerico hanno regioni diverse di stabilità, dando quindi dei risultati errati, mentre il metodo trapezoidale è a volte meno preciso, ma rispetta rigorosamente il concetto di stabilità di un sistema, divergendo solo quando effettivamente siamo di fronte a delle instabilità.

Questo porta a dover giocare con le opzioni per essere sicuri di trovare il metodo che più ci aggrada, e quello che serve per farci avvicinare con maggior precisione al risultato. Mathematica, d'altronde, per ogni funzione numerica è in grado di trovare da sola l'algoritmo ideale, applicandolo in maniera del tutto trasparente. Quando andremo a risolvere un'equazione differenziale, dovremo 'solo' preoccuparci di porre il problema nella maniera corretta, con le giuste condizioni al contorno, e Mathematica farà il resto. Ci concentriamo, insomma, sul cosa, non sul come.

Questo non toglie il fatto che siamo perfettamente in grado di dire al programma come risolvere un problema, forzando ad esempio un determinato algoritmo, ma non sono cose di cui ci preoccupiamo, almeno per ora. Personalmente, non ho mai detto a Mathematica quale algoritmo utilizzare, e non credo che lo farò mai, a meno di curiosità.

In aggiunta a questo, possiamo anche sfruttare la precisione arbitraria di Mathematica. Non solo il programma è in grado di decidere quale algoritmo applicare, ma lo fa con una precisione che possiamo decidere da soli. Non siamo costretti ad accontentarci della risoluzione dovuta alla precisione macchina, se 32 o 64 bit. Se vogliamo una soluzione con 10000 cifre significative, con 107

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

Mathematica possiamo farlo tranquillamente; il programma basa i suoi calcoli su potenti algoritmi a precisione arbitraria, che penalizzano quasi di niente la velocità di esecuzione degli algoritmi. Questi ultimi, fra l'altro, sono altamente ottimizzati, dandoci possibilità di eseguire un numero elevato di operazioni in relativamente poco tempo.

Adesso, andiamo a vedere le funzioni numeriche di Mathematica. Sono per lo più simili a quelle che abbiamo già visto per quanto riguarda il calcolo simbolico, almeno nel loro funzionamento base, e quando non andiamo a toccare le opzioni di default. Se poi vogliamo, possiamo spulciarle per risolvere anche l'impossibile, utilizzando funzioni quali NSolve ed NDSolve che sono sicuramente fra i comandi matematici più potenti in assoluto non soltanto in Mathematica, ma in qualsiasi altro software scientifico che possiate concepire.

Sommatorie, produttorie, integrali

Abbiamo visto come si risolvono le equazioni simboliche: se ci serve un risultato numerico, oppure il programma non è in grado di trovare la soluzione simbolica, possiamo utilizzare delle funzioni che sono molto simili, in quanto a sintassi: ⁄

¶

¤

i approssimazione numerica di f

NSum[f, {i, , Infinity}]

min i

min

¶

approssimazione numerica di f

NProduct[f, {i, i

min

i Ÿ

, Infinity}]

min x max

x x approssimazione numerica di f d x

NIntegrate[f, {x, , }] Ÿ Ÿ

min max x

min

x y

max max

x approssimazione di d x d y f

NIntegrate[f, {x, ,

min x y

min min

x y y

}, {y, , }]

max min max

Proviamo a effettuare qualche calcolo:

ê 8x,

NSum@Log@xD x ^ 5, 2, Infinity<D

In[1]:= +

0.0285738 0.

Out[1]= ê H1 8x,

−

NIntegrate@1 Sqrt@Gamma@xD xLD, 0, 1<D

In[2]:= 1.61784

Out[2]=

Gli integrali possono essere anche calcolati in regioni infinite, naturalmente:

8x, −Infinity,

NIntegrate@Cos@xD Exp@−Abs@x ^ 2DD, Infinity<D

In[3]:= 1.38039

Out[3]=

Il risultato ottenuto è numerico: è anche possibile lavorare sulla precisione, ma ne parleremo poco più avanti, dato che bisogna ragionarci un attimino per poter ottenere il risultato voluto.

108

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

Soluzione di equazioni

Anche in questo caso, le funzioni sono simili, anche se si aggiunge una funzione interessante, fatta apposta per il calcolo numerico: x] risolve un'equazione polinomiale numericamente

NSolve[lhs==rhs, risolve un sistema di equazioni polinomiali numericamente

N S o l v e [ {lhs = =rhs ,

1 1

lhs … y, …

==rhs , }, {x, }]

2 2 x cerca una soluzione numerica con seme x x

FindRoot[lhs==rhs, {x, }] =

0 0

cerca soluzioni numeriche per un sistema di equazioni

F i n d R o o t [ {lhs = =rhs ,

1 1

lhs … x

==rhs , }, {{x, },

2 2 0

y …

{y, }, }]

0

Abbiamo visto che non esistono soluzioni generali per equazioni polinomiali per un grado superiore a 4. Tuttavia, possiamo ancora, con il comando NSolve, trovare il risultato numerico per tutte le radici: = + − + − − +

eq x ^ 7 5 x ^ 6 34 x ^ 5 3 x ^ 4 18 x ^ 3 64 x 2 0;

In[4]:= NSolve@eq, xD

In[5]:= 88x 8x <,

8x <, 8x

→ −8.88949<, → −0.73661 − 0.860533

Out[5]= 8x <, 8x <, 8x

→ −0.73661 + →

0.860533 0.0312415<,

→ − → + →

0.716525 0.962437 0.716525 0.962437 3.89842<<

Ovviamente, anche in questo caso noi siamo in grado, tramite questo comando, di trovare soluzioni di un sistema di equazioni algebrico:

8x

= −

sis ^ 3 4 y ^ 2 3 z ^ 5,

In[6]:= −

x y z ^ 2,

+ +

x ^ 2 y ^ 2 z ^ 1 0<;

8x, êê

NSolve@sis, y, z<D TableForm

In[7]:= → −43.596 − → −36.4276 + → −6.0454 +

x 36.3488 y 43.5925 z 6

→ −43.596 + → −36.4276 − → −6.0454 −

x 36.3488 y 43.5925 z 6

→ → −1.028 → −1.0709

x 0.118835 y z

→ −0.276836 − → − → −0.088803

x 0.943164 y 1.19596 0.727046 z

→ −0.276836 + → + → −0.088803

x 0.943164 y 1.19596 0.727046 z

→ → → −0.568996

x 0.670104 y 0.346348 z

Out[7]//TableForm= +

→ −0.607146 + → + →

x 0.662123 y 0.000898087 0.351493 z 0.193327

→ −0.607146 − → − → −

x 0.662123 y 0.000898