Manuale di Mathematica 11

Daniele Lupo Appunti di Mathematica

Eq differenziali

Prima di cominciare inizializziamo caricando il package per visualizzare i grafici:

Needs@"Graphics`Graphics`"D;

Needs@"Graphics`Colors`"D;

Needs@"Graphics`Animation`"D;

Needs@"DifferentialEquations`NDSolveProblems`"D;

Needs@"DifferentialEquations`NDSolveUtilities`"D;

ü Introduzione

In quest'appendice vedremo meglio uno degli aspetti migliori (a mio avviso) di Mathematica, ossia i

suoi comandi e le sue opzioni per la risoluzione di equazioni differenziali, sia in forma simbolica,

che in forma numerica. Per esempio, possiamo trattare con estrema semplicità sia le equazioni, che i

sistemi di equazioni, anche di tipo misto, dove alcune equazioni del sistema sono differenziali,

mentre altre sono algebriche. Inoltre, possiamo trattare anche le equazioni differenziali alle derivate

parziali. Insomma, possiamo fare veramente di tutto. Tuttavia, se le equazioni differenziali che

trattiamo sono alquanto complicate, dobbiamo cominciare a cercare di capire meglio come

funzionano questi comandi, vedendoli in maniera più dettagliata. Il fatto che risolvano da soli la

maggior parte dei problemi è sicuramente una grande cosa, ma dobbiamo comunque tener conto del

fatto che mettere mano alle opzioni permette di avere un controllo maggiore su quello che si fa, e

permette di capire anche il funzionamento di certi algoritmi: inoltre potremmo anche essere in grado

di selezionare manualmente le opzioni e gli algoritmi migliori per il problema che ci serve.

ü Tipi di equazioni

Il comanid DSolve è in grado di risolvere sia singole equazioni differenziali, che sistemi di equazioni.

Ł HxL

Equazioni differenziali ordinarie ODE: sono i tipi di equazioni in cui la variabile indipendente è

una sola, mentre quelle dipendenti y possono essere una soltanto, nel caso di una sola equazione

i

differenziale, oppure più di una, nel caso in cui abbiamo a che fare con dei sistemi di equazioni

differenziali.

Ł Equazioni differenziali alle derivate parziali PDE: questo tipo di equazioni differenziali si

distinguono dalle prime per il fatto che la loro funzione (o variabile dipendente) dipende da più di

una variabile indipendente yHt , t , ...L. Questi casi sono più difficili da risolvere simbolicamente

1 2

rispetto alle ODE: Mathematica non è oggettivamente in grado di poter trovare soluzioni a qualsiasi

equazione PDE, per il semplice fatto che non esiste una teoria che permetta di calcolarle in maniera

esatta. Il problema è talmente difficile che ci sono, per fare un esempio, dei casi in cui ci sono premi

per chi riesce a trovare soluzioni esatte per questo tipo di equazioni differenziali; un esempio è dato

dal milione di dollari messo in palio (da chi adesso non ricordo, sinceramente), per chi riesca a

476

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

trovare soluzioni in forma simbolica delle equazioni di Navier-Stokes, per la particolare importanza

che rivestono in moltissimi campi dell'ingegneria.

Mathematica è in grado di risolvere PDE in maniera simbolica nella maggior parte dei casi in cui

l'equazione è del primo ordine, ed in un minor numero di casi quando la PDE è di ordine superiore.

Ł Equazioni algebrico-differenziali DAE: si tratta in particolar modo di sistemi, in cui compaiono sia

equazioni differenziali, che equazioni algebriche, in cui non compaiono le derivate. Per questo tipo

di sistemi valgono gli stessi problemi delle PDE, cioè che non esiste una teoria tanto robusta da poter

trovare in tutti i casi le soluzioni in maniera simbolica; anche in questo caso Mathematica non è in

grado di trovare le soluzioni per casi difficili, ma molti problemi sono comunque alla portata di

questo programma. Risoluzione Simbolica

Adesso vedremo meglio, più in dettaglio, e con un maggior numero di esempi, il funzionamento del

comando DSolve. La soluzione simbolica, è sempre da preferire, quando è possibile averla, perchè

naturalmente non abbiamo in questo caso problemi tipici del calcolo numerico, come errori di

approssimazione, necessità di specificare sempre le condizioni iniziali, la scelta dell'intervallo di

risoluzione e così via. In casi avanzati le soluzioni possono occupare anche diverse pagine, ma fin

quando lasciamo tutto il lavoro di elaborazione a Mathematica, non dobbiamo preoccuparcene più di

tanto, lasciando al programma il compito di elaborare le funzioni ottenute.

ü DSolve

Come ben sappiamo, è questo il comando che dobbiamo utilizzare quando vogliamo trovare la

soluzione di un'equazione differenziale in forma simbolica. Data la natura, appunto, simbolica della

soluzione, se si omettono le condizioni iniziali, Mathematica darà il risultato in forma generale, con

le opportune costanti, che di default sono date nella forma C[n]. Per poter scrivere la derivata, basta

utilizzare il carattere ' (quello che vi trovate accanto lo zero nella parte superiore della tastiera, per

capirci), come se si trattasse del simbolo che si usa nei libri; per poter creare derivate di ordine

superiore basta mettere tanti ' per quanto è necessario. Vediamo questo esempio:

+

DSolve@y ''@xD x y@xD 0, y@xD, xD

88y@xD 1ê3 1ê3

+

→ xD C@1D AiryBi@H−1L xD C@2D<<

AiryAi@H−1L

Possiamo vedere due cose, che in fondo già sapevamo (ma ricordare non fa mai male, vero, cari miei

folletti che volano sulle favole dell'ingegneria? Ehm....).

Prima di tutto, occorre notare come viene data la soluzione, cioè sotto forma di regola, invece che

direttamente come funzione. Questo permette di andare a sostituire la funzione dove serve. Inoltre,

come caso più generale, permette anche di andare a sostituire più funzioni in una volta sola, nel caso

477

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

avessimo più soluzioni. per andare a sostituirla basta usare ReplaceAll, che come sapete si può

scrivere anche come /.

ê.

y@xD %

8AiryAi@H−1L 1ê3 1ê3

+

xD C@1D AiryBi@H−1L xD C@2D<

Viene restituita come lista, anche se è formata da un solo elemento, per manterene la generalità della

soluzione: infatti, dato che le soluzioni possono essere più di una, si racchiudono in una lista, e

questa è semplicemente il caso particolare in cui abbiamo un solo elemento.

La seconda cosa da notare è come abbiamo introdotto le costanti, che Mathematica ha aggiunto

automaticamente, in mancanza delle opportune condizioni iniziali. Se le avessimo imposte, avremmo

anche ottenuto la soluzione particolare: êê

+

DSolve@8y ''@xD x y@xD 0, y@0D 3, y@3D 0<, y@xD, xD FullSimplify

99y@xD → H−AiryAi@−xD EN í

J3 2

+

2ê3 è!!!

AiryBi@−3D AiryAi@−3D AiryBi@−xDL GammaA

3

J2 EN== 3

1

BesselJA , 2 3

3

In questo caso abbiamo posto le condizioni iniziali, ed abbiamo ottenuto la soluzione particolare:

ê. 8x, 8Orange<,

→

8

Plot@y@xD %, 0, 12<, PlotStyle

→

8Blue, êê. 8%@@1,

Epilog →

8Red, êê. 8%@@1,

PointSize@0.04D, Point@80, y@xD 1DD, x 0<<D<,

→

< PointSize@0.04D, Point@83, y@xD 1DD, x 3<<D<

D

3

2

1 2 4 6 8 10 12

-1

-2

-3

Graphics 478

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

Come possiamo vedere, abbiamo aggiunto al grafico dei punti che rappresentano la funzione

calcolata nelle condizioni iniziali, e possiamo vedere come coindidano con quelle che avevamo

83, 80,

imposto in partenza, cosa che effettivamente era naturale. Guai ad non essere così. Potevo anche

semplicemente utilizzare le coordinate 0< e 3< per rappresentare direttamente i punti, ma in

questa maniera vi ho fatto vedere come ottenere lo stesso risultato lavorando con le funzioni e con le

regole di sostituzione.

Potete anche vedere come abbiamo imposto le condizioni iniziali nel comando DSolve, cioè creando

sintatticamente un sistema di equazioni simile al DAE, dove però c'è l'importante differenza che,

invece di equaioni delle variabili dipendenti generiche, compaiono equazioni per valori particolari

della funzione. Ovviamente, le condizioni iniziali devono essere compatibili con il problema.

Sempre nell'esempio fatto, infatti, abbiamo bisogno di due condizioni iniziali per avere la soluzione

esatta; un numero maggiore porta a generare un errore:

+

DSolve@8y ''@xD x y@xD 0, y@0D 3, y@3D 0, y '@0D 3<, y@xD, xD

— DSolve::bvnul : For some branches of the general solution, More…

the given boundary conditions lead to an empty solution.

8<

Questo perchè le condizioni iniziali sono inconsistenti con li problema. Ovviamente, in questo caso,

l'errore non è del programma, ma è vostro, ed è un modo per dirvi di andare a studiare...

Un'importante aspetto del comando DSolve, specialmente dalla versione 5 di Mathematica, è la sua

capacità di lavorare con funzioni definite a tratti, che quindi possono presentare delle discontinuità,

come la funzione UnitStep:

8y@tD

DSolve@ − + −

UnitStep@t 1D t ^ 2 UnitStep@t 2D y '@tD, y@0D 3<, y@tD, tD

99y@tD − − + +

2 t 1+t 2 2 2

2 10 3 2 t t

→ − +

2

J3 N

− − + +

2 t 1+t 2 2 2

10 3 2 t t

2

−1+t + − +

UnitStep@2 tD

H3 2

−1+t

− − −

UnitStep@1 tD 3 UnitStep@2 tDL==

479

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

ê. 8t, 8Blue<D

→

Plot@y@tD %@@1DD, 0, 3<, PlotStyle

30

25

20

15

10

5 0.5 1 1.5 2 2.5 3

Graphics = =

Potete notare i punti di discontinuità che si hanno per t 1, t 2, che sono i punti dove iniziano i

gradini nella nostra equazione differenziale.

Le funzioni definite a tratti possono essere anche utilizzate nei sistemi, come in questo caso:

8

=

sistema − − + +

UnitStep@t 3D y '@tD x@tD t UnitStep@−t 3D x '@tD t ^ 2 x@tD,

y '@tD t x '@tD,

x '@1D 4,

<;

y@0D 4

Per comodità, oltre al sistema, ho specificato anche le condizioni iniziali nel sistema stesso: adesso

posso andare a risolverlo normalmente, utilizzando sempre DSolve:

8x@tD, êê

DSolve@sistema, y@tD<, tD FullSimplify

Ø

99x@tD ≤ 1

− +Log@tD t≤3

∞

≤ t

1+

± 10

− +t+Log@tD

→ True ,

2 3 H1 H−1 D ==

Ø

≤ 7

≤ − +t

− + + + + >

1

2ê3 H1 D

4 2 2 tL tL ExpIntegralEi@− t 3

3

∞ 3

≤

→

y@tD ≤

± 1

+ + + 1

1−

4 t tL ExpIntegralEi@− True

t t

Notate la potenza del comando, che risulta in grado di definire funzioni a tratti anche nell'argomento

dell'esponenziale, e come la soluzione includa funzioni avanzate come ExpIntegralEi

480

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

ê. ê. 8t, −7,

88Orange<, 8Dashing@80.05,

Plot@8x@tD %, y@tD %<, 7<,

→

PlotStyle 0.015<D, Green<<D

200

150

100

50

-6 -4 -2 2 4 6

-50

-100

-150

Graphics

In maniera analoga, possiamo risolvere anche equazioni differenziali alle derivate parziali, andando

a specificare sempre quali sono le variabili indipendenti. Considerando la seguente equazione

differenziale

= + −

eqn x y ^ 2 D@u@x, yD, xD y ^ 2 x ^ 2 D@u@x, yD, yD x y u@x, yD

@x, @x,

H0,1L H1,0L

−x + +

2 2 2

y u@x, yD x y u yD x y u yD

8x,

=

sol DSolve@eqn 0, u@x, yD, y<D

99u@x, H−x

E

è!!!!!!!!

!!!!!!!!!!!!!

è!!!!!!!!

!!!!!!!!!!!!

x

2 ArcTanA −x2 +2 1

y

→ +

2

yD C@1DA 2 yLE==

−x +2

2 y 2

Da questa scrittura possiamo notare subito un paio di cose importanti. Prima di tutto, possiamo

notare il diverso tipo di notazione, quando dobbiamo scrivere delle PDE. infatti, scrivendo qualcosa

come u'@x, tD, Mathematica (e neanche noi, a dire il vero), non è in grado di distinguere la variabile

rispetto alla quale stiamo derivando; per questo motivo, dobbiamo per forza di cose andare ad

esplicitare tramite il comando D il tipo di derivata parziale che andiamo ad utilizzare. Comparirà in

maniera corretta nell'output, come possiamo andare a vedere.

Inoltre, i coefficienti C[i], in questo caso, non rappresentano più delle costanti, ma sono a tutti gli

H-x

effetti delle funzioni. Si può notare considerando che C[1], nel nostro caso, è una funzione, avente

1

ÅÅÅÅ +

2

come argomento 2 yL, come possiamo vedere dalla soluzione che abbiamo ottenuto. Dalla

2

soluzione generale, quindi, se vogliamo arrivare ad una particolare, dobbiamo andare a sostituire a

C[i] delle funzioni (che naturalmente, a seconda delle condizioni al contorno, possono essere anche

481

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

delle costanti...). Possiamo trovare una soluzione particolare del nostro problema, ipotizzando per

esempio, di mettere come caso particolare il seno:

ê. ê. 8C@1D@t_D

= →

solpar u@x, yD %@@1DD Sin@tD<

H−x

E

è!!!!!!!!

!!!!!!!!!!!!!

è!!!!!!!!

!!!!!!!!!!!!

x

2 ArcTanA −x2 +2 1

y +

2

SinA 2 yLE

−x +2

2 y 2

Notiamo come sia stato necessario preservare la struttura dell'argomento della costante C[i], in

quanto è quello necessario per la risoluzione del problema. Adesso possiamo andare a visualizzare la

-x + ¥

2

funzione, notando che ottengo una soluzione reale soltanto quando 2 y 0, cioè quando è

reale il radicale, considerato che non ho limiti per quanto riguarda l'argomente del seno e dell'Arc-

Tan, sempre che siano ovviamente numeri reali :-). Per cui otterremmo degli errori nel caso in cui

specifichiamo durante il plottaggio, dei punti al di fuori questo dominio. Per risolvere il problema

allora andiamo a plottare una funzione Piecewise, che facciamo corrispondere alla nostra funzione se

rientra nel dominio, e la poniamo uguale a 0 se invece cade fuori i nostri limiti permessi:

8x,

= + >

dominio Reduce@−x ^ 2 2 y 0, y<D

2

x

∈ >

x Reals && y 2 482

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

88

DisplayTogetherArray@ 80,

8x, 8y,

DensityPlot@Piecewise@88solpar, dominio<, True<<D,

−4, −.1,

4<, 10<,

→

PlotPoints 200,

→

Mesh False, →

HHue@Sqrt@#DD

ColorFunctionScaling False,

→ 80,

ColorFunction &LD,

8x, 8y,

Plot3D@Piecewise@88solpar, dominio<, True<<D,

−4, −.1,

4<, 10<,

→

PlotPoints 130,

→ 8Automatic, 80,

Mesh False,

→ 8−2.004,

PlotRange Automatic, 5<<,

−> −1.853,

<< ViewPoint 2.000<D

D 10

8 5

6 4

3

2

1 4

4 10

0 2

7.5 0

5

2 -2

2.5 -4

0

0

-4 -2 0 2 4

GraphicsArray

Abbiamo visto come possiamo risolvere sia le ODE che le PDE. Naturalmente, possiamo risolvere

anche le DAE, sempre nei casi in cui è permesso. Possiamo considerare il seguente sistema:

8

=

dae

f '@xD 5 g ''@xD,

+

f@xD g@xD 3 Cos@xD,

<;

f@PiD 0, f '@0D 0 483

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

8f,

=

soluzione DSolve@dae, g<, xD

99f → FunctionA8x<,

H−5

15 π π π π

x x x x

− − πê5 + + +

+ + + −

xê5

5 5 5 Cos@xD Sin@xDLE,

5 5 5 5 5 5 5 5

26 3 π x

− −

→

g FunctionA8x<,

H25 5 5

26

π π π

x x x

πê5 + + +

− − + +

xê5

25 25 Cos@xD 5 Sin@xDLE==

5 5 5 5 5 5

ê. 8x,

+ −9,

88Blue<, 8Green<, 8Orange<<,

Plot@Evaluate@8f@xD, g@xD, f@xD g@xD< soluzioneD, 9<,

→

PlotStyle →

D

Background GrayLevel@0.9D

10

5

-7.5 -5 -2.5 2.5 5 7.5

-5

-10

-15

Graphics

Possiamo notare come la loro somma dia effettivamente 3 cosHxL, come avevamo scritto nel sistema

misto.

Il modo di lavorare di DSolve è abbastanza complicato ed oscuro per la maggior parte degli esseri

umani (diciamo tutti apparte la ventina di persone che l'hanno programmato). Tuttavia imparare ad

usarlo bene, permette di inserire l'input in una maniera che può anche agevolare il comando a

lavorare.

Quello che faremo adesso è quello di andare ad analizzare un tantino più in dettaglio quello che

possiamo ottenere con DSolve, e capire un pochino meglio le equazioni differenziali.

ODE

Sono le equazioni che si incontrano per la prima volta, perchè corrispondono a quelle più semplici

da risolvere. Il caso più diretto è un'equazione differenziale di questo tipo:

484

Printed by Mathematica for Students

Daniele Lupo Appunti di Mathematica

HtL HtL

¢ n

y f

Questo tipo di equazione differenziale si risolve semplicemente integrando la parte destra

dell'equazione stessa:

DSolve@y '@tD Sin@tD ^ 2, y@tD, tD

99y@tD 1

t + −

→ C@1D Sin@2 tD==

2 4

La soluzione è semplice e diretta. Dato la natura di questo tipo di equazione, potev