Appunti di Mathematica

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. Ł Equazioni differenziali ordinarie ODE: sono i tipi di equazioni in cui la variabile indipendente è una sola, mentre quelle dipendenti yi HxL possono essere una soltanto, nel caso di una sola equazione 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 yHt1 , t2 , ...L. Questi casi sono più difficili da risolvere simbolicamente 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 Printed by Mathematica for Students 476 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 → AiryAi@H−1L1ê3 xD C@1D + AiryBi@H−1L1ê3 xD C@2D<< 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 Printed by Mathematica for Students 477 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−1L1ê3 xD C@1D + AiryBi@H−1L1ê3 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 99y@xD → 0, y@0D 3, y@3D 0<, y@xD, xD êê FullSimplify J3 32ê3 H−AiryAi@−xD AiryBi@−3D + AiryAi@−3D AiryBi@−xDL GammaA J2 BesselJA 1 è!!! , 2 3 EN== 3 2 EN í 3 In questo caso abbiamo posto le condizioni iniziali, ed abbiamo ottenuto la soluzione particolare: Plot@y@xD ê. %, 8x, 0, 12<, PlotStyle → 8Orange<, Epilog → 8 8Blue, PointSize@0.04D, Point@80, y@xD êê. 8%@@1, 1DD, x → 0<<D<, 8Red, PointSize@0.04D, Point@83, y@xD êê. 8%@@1, 1DD, x → 3<<D< < D 3 2 1 2 4 6 8 -1 -2 -3 Graphics Printed by Mathematica for Students 478 10 12 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 semplicemente utilizzare le coordinate 83, 0< e 80, 3< per rappresentare direttamente i punti, ma in imposto in partenza, cosa che effettivamente era naturale. Guai ad non essere così. Potevo anche 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 DSolve::bvnul : For some branches of the general solution, the given boundary conditions lead to an empty solution. 3<, y@xD, xD More… 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: DSolve@ 8y@tD UnitStep@t − 1D + t ^ 2 UnitStep@t − 2D 99y@tD → − J3 2 −1+t + 2 − 10 2 2 t −3 − 10 1+t t −3 UnitStep@1 − tD H3 − 3 Printed by Mathematica for Students +2 2 1+t +2 2 2 −1+t t+ 2 2 t2 t+ 2 + t2 y '@tD, y@0D N UnitStep@2 − tD + UnitStep@2 − tDL== 479 3<, y@tD, tD Daniele Lupo Appunti di Mathematica Plot@y@tD ê. %@@1DD, 8t, 0, 3<, PlotStyle → 8Blue<D 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: sistema = 8 UnitStep@t − 3D y '@tD − x@tD t + UnitStep@−t + 3D x '@tD t ^ 2 y '@tD t x '@tD, x '@1D 4, y@0D 4 <; x@tD, Per comodità, oltre al sistema, ho specificato anche le condizioni iniziali nel sistema stesso: adesso posso andare a risolverlo normalmente, utilizzando sempre DSolve: DSolve@sistema, 8x@tD, y@tD<, tD êê FullSimplify 99x@tD → 2 Ø ≤− 1+∞ ≤ ±− 1 t +Log@tD 10 3 +t+Log@tD t≤3 True , 2ê3 + 2 − 73 +t H1 + H−1 + tL tL + Ø ExpIntegralEi@− ≤ ≤ 4−2 y@tD → ∞ 1 ≤ ≤ 4 + 1− t t H1 + tL + ExpIntegralEi@− 1 D ± t 1 3 D t>3 True == 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 Printed by Mathematica for Students 480 Daniele Lupo Appunti di Mathematica Plot@8x@tD ê. %, y@tD ê. %<, 8t, −7, 7<, PlotStyle → 88Orange<, 8Dashing@80.05, 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 y u@x, yD + x2 y2 uH0,1L @x, yD + x y2 uH1,0L @x, yD sol = DSolve@eqn 99u@x, yD → 0, u@x, yD, 8x, y<D x 2 ArcTanA è!!!!!!!! !!!!!!!!!!!!! E −x2 +2 y è!!!!!!!! !!!!!!!!!!!! −x2 +2 y C@1DA 1 H−x2 + 2 yLE== 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 come argomento ÅÅÅÅ12 H-x2 + 2 yL, come possiamo vedere dalla soluzione che abbiamo ottenuto. Dalla effetti delle funzioni. Si può notare considerando che C[1], nel nostro caso, è una funzione, avente 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 Printed by Mathematica for Students 481 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: solpar = u@x, yD ê. %@@1DD ê. 8C@1D@t_D → Sin@tD< x 2 ArcTanA è!!!!!!!! !!!!!!!!!!!!! E −x2 +2 y è!!!!!!!! !!!!!!!!!!!! −x2 +2 y SinA 1 H−x2 + 2 yLE 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 funzione, notando che ottengo una soluzione reale soltanto quando -x2 + 2 y ¥ 0, cioè quando è reale il radicale, considerato che non ho limiti per quanto riguarda l'argomente del seno e dell'ArcTan, 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: dominio = Reduce@−x ^ 2 + 2 y > 0, 8x, y<D x ∈ Reals && y > Printed by Mathematica for Students x2 2 482 Daniele Lupo Appunti di Mathematica DisplayTogetherArray@ 88 DensityPlot@Piecewise@88solpar, dominio<, 80, True<<D, 8x, −4, 4<, 8y, −.1, 10<, PlotPoints → 200, Mesh → False, ColorFunctionScaling → False, ColorFunction → HHue@Sqrt@#DD &LD, Plot3D@Piecewise@88solpar, dominio<, 80, True<<D, 8x, −4, 4<, 8y, −.1, 10<, PlotPoints → 130, Mesh → False, PlotRange → 8Automatic, Automatic, 80, 5<<, ViewPoint −> 8−2.004, −1.853, 2.000<D << D 10 8 5 4 3 2 1 10 0 6 4 4 2 7.5 0 5 2 -2 2.5 0 -4 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: dae = 8 f '@xD 5 g ''@xD, f@xD + g@xD 3 Cos@xD, f@PiD 0, f '@0D 0 <; Printed by Mathematica for Students 483 Daniele Lupo Appunti di Mathematica soluzione = DSolve@dae, 8f, g<, xD 99f → FunctionA8x<, π x 15 − π5 − x5 H−5 πê5 + 5 5 + 5 + 5 26 3 − π5 − x5 g → FunctionA8x<, 26 π x H25 πê5 − 25 5 + 5 − 25 xê5 + π 5 Cos@xD − xê5 +5 + Cos@xD + 5 x 5 π 5 + x5 π 5 + x 5 π 5 + x 5 Sin@xDLE, Sin@xDLE== Plot@Evaluate@8f@xD, g@xD, f@xD + g@xD< ê. soluzioneD, 8x, −9, 9<, PlotStyle → 88Blue<, 8Green<, 8Orange<<, Background → GrayLevel@0.9D D 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: Printed by Mathematica for Students 484 Daniele Lupo Appunti di Mathematica y¢ HtL ‹ f HtL Questo tipo di equazione differenziale si risolve semplicemente integrando la parte destra dell'equazione stessa: DSolve@y '@tD 99y@tD → Sin@tD ^ 2, y@tD, tD t 1 + C@1D − Sin@2 tD== 2 4 La soluzione è semplice e diretta. Dato la natura di questo tipo di equazione, potevamo risolverla anche in questa maniera: Integrate@Sin@tD ^ 2, tD t 1 − Sin@2 tD 2 4 Ottenendo esattamente lo stesso risultato; l'unica differenza risiede nel fatto che Integrate non mette direttamente la costante C (tutti gli integrali indefiniti dovrebbero essere in teoria come C + FHtL), mentre l'equazione differenziale l'aggiunge, perchè varia al variare delle condizioni al contorno della soluzione generica DSolveAy @xD x Sin@xD y@xD ^ 2 , y, xE è!!!!!!!!!! 3−x è!!!!!!!!!! è!!!!!!!!!! 3 − x − 2 3 +2 x 3 − x − è!!!!!!!!!! x è!!! 4 3 + x C@1D − H1 + 6 L H−1L3ê4 π Erfi@H−1L1ê4 3 − x D − è!!! è!!!!!!!!!! H1 − 6 L H−1L1ê4 6 + x π Erfi@H−1L3ê4 3 − x DLD<< 88y → Function@8x<, H4 3 + x L ë H−2 3 Adesso abbiamo invece un esempio di equazione omogenea: omo = y '@xD −Hx ^ 4 − 3 y@xD ^ 3L ∗ Hx ∗ y@xDL; sol = DSolve@omo, y, xD êê FullSimplify 99y → FunctionA8x<, − 9y → FunctionA8x<, 9y → FunctionA8x<, Printed by Mathematica for Students I2 I2 I2 x6 2 x6 2 x6 2 H−2L1ê3 Hx6 L Hx6 L1ê3 C@1D + 3 21ê3 21ê3 Hx6 L Hx6 L1ê3 C@1D + 3 21ê3 x6 2 1ê9 x2 Gamma@ 485 , x6 2 1ê9 x6 2 x2 Gamma@ H−1L2ê3 21ê3 Hx6 L Hx6 L1ê3 C@1D + 3 21ê3 1 3 x6 2 1ê3 1 3 , x6 2 DM 1 3 , x6 2 DM 1ê9 x2 Gamma@ DM 1ê3 1ê3 E=, E=, E== Daniele Lupo Appunti di Mathematica Come vediamo, abbiamo tre funzioni che risolvono questa equazione differenziale. Effettivamente, però, le tre funzioni coincidono: f1 = Table@Abs@y@xDD ê. sol@@1DD ê. C@1D → p, 8p, −5, 5, .5<D; f2 = Table@Abs@y@xDD ê. sol@@2DD ê. C@1D → p, 8p, −5, 5, .5<D; f3 = Table@Abs@y@xDD ê. sol@@3DD ê. C@1D → p, 8p, −5, 5, .5<D; f1 f2 f3 True Abbiamo visto che sono uguali, indipendentemente dal parametro, e basta quindi plottare soltanto una famiglia di equazioni, diciamo f1: colori = PlotStyle → Table@8Hue@Random@DD<, 8p, 0, 20<D; Plot@Evaluate@f1, 8x, 0, 2<, coloriDD; 3 2.5 2 1.5 1 0.5 0.5 1 1.5 2 Piccola particolarità: non è necessatio che PlotStyle definisca un numero di stili uguale al numeri di funzioni, perchè viene ripetuto ciclicamente. Detto questo, abbiamo visto come ottenere le soluzioni generali. Questo può essere più utile di quanto sembri. Infatti, possiamo avere vari rami della soluzione generale. Per qualcuno di essi, può non valere la condizione iniziale che vogliamo. Ipotizziamo di avere questa equazione differenziale: rami = y '@xD Printed by Mathematica for Students −Hx − 3 y@xD ^ 2L ê Hx ^ 2 ∗ y@xDL; 486 Daniele Lupo Appunti di Mathematica sol = DSolve@rami, y, xD 6 −6êx JC@1D + 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% 99y → FunctionA8x<, − $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ExpIntegralEiA EN% E=, x 6 −6êx JC@1D + 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% ExpIntegralEiA EN% E== 9y → FunctionA8x<, $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x r1 = Table@y@xD ê. sol@@1DD ê. C@1D → p, 8p, −5, 5, .5<D; r2 = Table@y@xD ê. sol@@2DD ê. C@1D → p, 8p, −5, 5, .5<D; DisplayTogetherArray@8 Plot@Evaluate@r1, 8x, 0, 4.5<, coloriDD, Plot@Evaluate@r2, 8x, 0, 4.5<, coloriDD <D -0.25 -0.5 -0.75 -1 -1.25 -1.5 1 2 3 4 1.5 1.25 1 0.75 0.5 0.25 1 2 3 4 GraphicsArray In questo esempio abbiamo evidenziato come, ad un'equazione differenziale, possano corrispondere più rami della soluzione, definite da due (o più, a seconda dei casi) funzioni. Vediamo quello che succede adesso: DSolve@8rami, y@2D — 2<, y, xD DSolve::bvnul : For some branches of the general solution, the given boundary conditions lead to an empty solution. More… 99y → FunctionA8x<, 6 è!!! $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% −6êx J2 3 − ExpIntegralEi@3D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% + ExpIntegralEiA EN E== 2 x In questo caso abbiamo avuto un warning... La soluzione generale, infatti. è formata da due rami ma, dal comportamento che abbiamo visto nei grafici, vediamo come solamente uno dei due rami sia in grado di soddisfare la condizione al contorno che abbiamo imposto. Il warning sta appunto a dirci che, nonostante Mathematica abbia trovato due soluzioni indipendenti per l'equazione differenziale, Printed by Mathematica for Students 487 Daniele Lupo Appunti di Mathematica soltanto una viene utilizzata per trovare quella particolare. In altre occasioni, invece, Mathematica riesce a trovare la soluzione, per così dire, in maniera parziale; questo può capitare perchè, ad esempio, durante i calcoli si ritrova un integrale che non riesce a risolvere in maniera analitica, come nell'esempio seguente: integ = x y '@xD Sin@xD ^ 4 y@xD Cos@xD; DSolve@integ, y, xD 99y → FunctionA8x<, Ÿ1 − x 8 Cos@K$101533D K$101533 H−3+4 Cos@2 K$101533D−Cos@4 K$101533DL K$101533 C@1DE== Come possiamo vedere, nella soluzione compare un integrale. Questo succede perchè Mathematica se lo è ritrovato davanti, senza essere in grado di risolverlo. La variabile utilizzata per l'integrale è una variabile locale, di quelle che il programma si crea da solo. Se vogliamo, possiamo andarlo a sostituire con un'altra variabile, per rendere la soluzione un tantinello più leggibile: % ê. K$101533 → t 99y → FunctionA8x<, Ÿ1 − x 8 Cos@tD t H−3+4 Cos@2 tD−Cos@4 tDL t C@1DE== Naturalmente, la soluzione non è stata calcolata completamente. Per avere il valore della funzione in un punto, occorre sempre valutare numericamente l'integrale, quindi l'intera espressione y@3D ê. % ê. C@1D → 1 êê N 83.39381 × 10−18 < Oltre agli integrali, nel comando DSolve occorre anche risolvere delle equazioni, in maniera simbolica, mediante il comando Solve. Quando quest'ultimo non è in grado di lavorare al meglio, Mathematica ce lo segnala: errSolve = Sin@y '@xDD x Printed by Mathematica for Students y@xD; 488 Daniele Lupo Appunti di Mathematica DSolve@errSolve, y, xD — Solve::ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. More… — Solve::tdep : The equations appear to involve the variables to be solved for in an essentially non−algebraic way. More… — Solve::tdep : The equations appear to involve the variables to be solved for in an essentially non−algebraic way. More… — Solve::tdep : The equations appear to involve the variables to be solved for in an essentially non−algebraic way. More… — General::stop : Further output of Solve::tdep will be suppressed during this calculation. SolveA‡ 1 y@xD x 1 K$112646 − ArcSin@K$112646D K$112646 More… C@1D − Log@xD, y@xDE In questo caso, non viene dato il risultato completo, ma neanche quello che avevamo in partenza. Quello che abbiamo ottenuto è un'espressione che coinvolge il comando Solve, perchè non è stato in grado di dare la soluzione cercata, restituendocelo così com'è. Anche in questo caso, come nel precedente, per avere soluzioni occorrono altre manipolazioni, e risoluzioni per via numerica, dato che non si possono trovare soluzioni in forma chiusa (dubito che possa risultare utile un risultato dove comare Solve ed un integrale non calcolabile simbolicamente). A questo punto conviene usare NDSolve, ma ne parleremo nell'opportuna sezione. Un altro tipo di equazioni differenziali che Mathematica è in grado di risolvere sono quelle di Bernoulli: bern = y '@xD + 3 x y@xD x ^ 4 y@xD ^ 4; solbern = DSolve@bern, y, xD 99y → FunctionA8x<, − 9y → FunctionA8x<, 9y → FunctionA8x<, I6 x + 18 x3 + 54 I6 x + 18 x3 + 54 I6 x + 18 x3 + 54 3 H−2L1ê3 9 x2 2 C@1D − 9 x2 2 3 21ê3 9 x2 2 C@1D − 9 x2 2 9 x2 2 C@1D − 9 x2 2 3 H−1L2ê3 21ê3 è!!!!!!! 2 π ErfA 3x è!!!! 2 EM 1ê3 E=, è!!!!!!! 2 π ErfA 3x è!!!! 2 EM 1ê3 E=, è!!!!!!! 2 π ErfA 3x è!!!! 2 EM 1ê3 E== Naturalmente, avendole trovate, sono le soluzioni dell'equazione differenziale, e possiamo verificarlo facilmente: Printed by Mathematica for Students 489 Daniele Lupo Appunti di Mathematica bern ê. solbern êê Simplify 8True, True, True< Riccati, che sono nella forma y£ HxL f HxL + gHxL yHxL + hHxL yHxL2 . Qua sotto potete vederne un Un altro tipo di equazoni differenziali, particolarmente difficili da risolvere, sono le equazioni di esempio: riccati = y '@xD + 1 ê x ^ 3 − 3 x ^ 3 y@xD ^ 2 0; DSolve@riccati, y, xD 99y → FunctionA8x<, 1 è!!! 2 è!!! è!!! 3 xD − 3 x HBesselY@1, − 3 xD − 2 1 è!!! è!!! è!!! BesselY@3, − 3 xDL + J2 x BesselJ@2, − 3 xD − 3 2 è!!! è!!! x2 HBesselJ@1, − 3 xD − BesselJ@3, − 3 xDLN C@1DN í è!!! è!!! H3 x3 Hx2 BesselY@2, − 3 xD + x2 BesselJ@2, − 3 xD C@1DLLE== −J2 x BesselY@2, − nella seguente forma: y£ HxL f HxL + gHxL yHxL + hHxL yHxL2 + kHxL yHxL3 Un altro tipo di equazioni differenziali difficili da risolvere sono quelle di Abel, che vengono poste abel = y @xD y@xD3 − x y@xD2 x−1 ; DSolve@abel, y, xD — — — InverseFunction::ifun : Inverse functions are being used. Values may be lost for multivalued inverses. More… InverseFunction::ifun : Inverse functions are being used. Values may be lost for multivalued inverses. More… Solve::tdep : The equations appear to involve the variables to be solved for in an essentially non−algebraic way. More… 1−x+ SolveA 1 y@xD −1 + x + C@1D + ExpIntegralEiA1 − x + 1 E y@xD 0, y@xDE In questo caso abbiamo diversi warning: questi sono dovuti al fatto che, durante la ricerca della soluzione esatta di questo tipo di equazioni, vengono utilizzati inverse di ODE, per cui compaiono delle funzioni inverse, di cui Mathematica non conosce direttamente la struttura, e se in questo modo si perdono alcuni valori durante l'inversione. Nel caso seguente, invece, tutto quanto va a buon fine: abel2 = y @xD Printed by Mathematica for Students y@xD x − 3 y@xD2 + x y@xD3 ; 490 Daniele Lupo Appunti di Mathematica sol = DSolve@abel2, y, xD 99y → FunctionA8x<, 9y → FunctionA8x<, 1 1 − E=, x 1 "################ ###### x2 x2 + C@1D 1 1 + E== x 1 "################ ###### x2 x2 + C@1D a1 = Table@y@xD ê. sol@@1DD ê. C@1D → p, 8p, 0, 1000, 20<D; a2 = Table@y@xD ê. sol@@2DD ê. C@1D → p, 8p, 0, 1000, 20<D; DisplayTogetherArray@8 Plot@Evaluate@a1, 8x, −2, 2<, coloriD, PlotRange → 8Automatic, 8−10, 10<<D, Plot@Evaluate@a2, 8x, −2, 2<, coloriD, PlotRange → 8Automatic, 8−10, 10<<D <D 10 7.5 5 2.5 -2 -1 10 7.5 5 2.5 1 -2.5 -5 -7.5 -10 2 -2 -1 -2.5 -5 -7.5 -10 1 2 GraphicsArray Come possiamo vedere, anche in questo caso abbiamo due rami distinti per la soluzione dell'equazione differenziale. Finora abbiamo considerato solamente alcuni esempi di equazioni ODE di primo grado. Naturalmente, però, possiamo considerare anche ordini superiori: possiamo prendere un esempio di ODE a coefficienti costanti, che sono il caso sempre più semplice: ode2 = y ''@xD − y@xD 0; DSolve@ode2, y, xD 88y → Function@8x<, x C@1D + −x C@2DD<< Naturalmente, in questo casi abbiamo due costanti di indeterminazione C[1] e C[2]. Vediamo come varia la soluzione al variare di quesi valori: Printed by Mathematica for Students 491 Daniele Lupo Appunti di Mathematica funz = FlattenATableA9y@xD ê. %@@1DD ê. 8C@1D → a, C@2D → b<, a+3 b+3 , .1E=, 8a, −3, 3, .4<, 8b, −3, 3, .4<E, 1E; RGBColorA , 6 6 Ho creato in questa maniera una lista di elementi, in cui il primo rappresenta la funzione, mentre il secondo il colore della funzione, che varia al variare dei parametri: Plot@Evaluate@funz@@All, 1DDD, 8x, −5, 5<, PlotStyle → funz@@All, 2DDD 75 50 25 -4 -2 2 4 -25 -50 -75 Graphics Per questo tipo di equazioni, la soluzione è data come combinazione lineare di esponenziali. Dato che abbiamo in questo esempio il secondo ordine, il risultato è combinazione lineare di due esponenziali. Per questo tipo di equazioni, come ben saprete, si può ottenere il risultato, a partire dall'equazione caratteristica, dalla quale si ricavano poi gli esponenti degli esponenziali che formano la base per lo spazio delle soluzioni. Prendiamo quest'altra equazione: car = y '''@xD + 4 y ''@xD − 11 y '@xD − 30 y@xD 0; A questo punto, possiamo ricavarci l'equazione caratteristica: car ê. 8Derivative@n_D@yD@xD → λ ^ n, y@xD → 1< −30 − 11 λ + 4 λ2 + λ3 0 Dato che tutto è rappresentato come espressioni, in Mathematica, lo è anche la forma y'' HxL, che viene descritta da Derivative@nD@yD@xD, dove n rappresenta il grado di derivazione, y rappresenta la funzione ed x la variabile indipendente. Abbiamo utilizzato quindi Derivative nella regola di sostituzione, per andare a sostituire tutte le derivate con la variabile l, elevata al corrispettivo esponente. Questa regola permette di ricavarci direttamente l'equazione caratteristica, senza bisogno di doverla ricopiare. Una volta ottenuta, possiamo risolverla: Printed by Mathematica for Students 492 Daniele Lupo Appunti di Mathematica Solve@%, λD 88λ → −5<, 8λ → −2<, 8λ → 3<< Ed, arrivati a questo punto, possiamo creare la funzione che è soluzione dell'equazione differenziale, andando a mapparci le soluzioni in una combinazione lineare. Il comando seguente prende le soluzioni che abbiamo trovato, e le trasforma da regole a valori. A questo punto, applichiamo MapIndexed alla lista delle soluzioni, in maniera tale da creare una lista dove ogni elemento è formato dall'esponenziale dove compare la soluzione, e dal coefficiente di ogni elemento, con opportuno indice. A questo punto, applico Plus alla lista trovata, in modo da sommare tutti gli elementi presenti: Plus @@ HMapIndexed@a#2@@1DD Exp@x #1D &, λ ê. %DL −5 x a1 + −2 x a2 + 3x a3 Notate come i comandi che ho utilizzato siano generici, il che mi permette di utilizzarli senza nessuna modifica per ODE a coefficienti costanti di qualsiasi ordine. L'unica accortezza è che, per equazioni di livello elevato, potrebbe essere necessario sostituire Solve con NSolve. Naturalmente, potevamo fare tutto direttamente tramite DSolve: DSolve@car, y, xD 88y → Function@8x<, −5 x C@1D + −2 x C@2D + 3x C@3DD<< D'altronde, spiegare la risoluzione classica credo che abbia giovato a molti di voi, impratichendovi un pochino con le manipolazioni. D'altronde, qua ormai si parla seriamente... Possiamo eseguire lo stesso procedimento mediante un'altra equazione: car2 = y ''''@xD − 10 ∗ y '''@xD + 54 ∗ y ''@xD − 130 ∗ y '@xD + 125 ∗ y@xD car2 ê. 8Derivative@n_D@yD@xD → λ ^ n, y@xD → 1< 125 − 130 λ + 54 λ2 − 10 λ3 + λ4 0 Solve@%, λD 88λ → 2 − <, 8λ → 2 + <, 8λ → 3 − 4 <, 8λ → 3 + 4 << Plus @@ HMapIndexed@a#2@@1DD Exp@x #1D &, λ ê. %DL H2− L x a1 + H2+ L x Printed by Mathematica for Students a2 + H3−4 L x a3 + 493 H3+4 L x a4 0; Daniele Lupo Appunti di Mathematica In questo caso abbiamo coefficienti immaginari, per cui possiamo trasformare l'espressione ottenuta in maniera da far comparire seni e tett... ehm, coseni: Con il comando DSolve avremmo ottenuto direttamente: DSolve@car2, y@xD, xD 88y@xD → 2x C@2D Cos@xD + 3x C@4D Cos@4 xD + 2x C@1D Sin@xD + 3x C@3D Sin@4 xD<< Lascio a voi il compito di mostrare l'equivalenza di queste due espressioni, tenendo conto del fatto che non è così immediato come possa sembrare a prima vista... Dovrete lavorarci un po', cari miei discepoli. x y HxL + a x y HxL + b yHxL Vediamo 2 ££ adesso di £ trattare l'equazione di Eulero, che ha la seguente forma: 0. Possiamo vederne un esempio: eulero = x ^ 2 ∗ y ''@xD + 5 ∗ x ∗ y '@xD + 9 ∗ y@xD 0; DSolve@eulero, y@xD, xD 99y@xD → è!!! è!!! C@2D Cos@ 5 Log@xDD C@1D Sin@ 5 Log@xDD + == x2 x2 In questo caso non abbiamo nente di particolarmente complicato, come soluzione. Possiamo invece forma Hc x + dL2 y££ HxL + a Hc x + dL y£ HxL + b yHxL considerare una generalizzazione, data dall'equazione lineare di Legendre, che sono equazioni nella 0: legendre = H9 x + 1L ^ 2 ∗ y ''@xD + 5 ∗ H3 x − 1L ∗ y '@xD + 6 ∗ y@xD 0; DSolve@legendre, y@xD, xD 99y@xD → è!!!!!! 54 I−11− 67 M 1 H−5L C@1D 2 3 J N 2 1 + 18 x + 81 x è!!!!!! è!!!!!! 1 %%%%%%%%%%% 11 67 2 67 20 $%%%%%%%%%%%%%%%%%%%%%%%% E+ Hypergeometric1F1A− − , 1− , 1 + 18 x + 81 x2 27 27 27 27 1 27 è!!!!!! I−11− 67 M 2 27 è!!!!!! I−11− 67 M 1 9 è!!!!!! I11+ 67 M 1 è!!!!!! 54 I−11+ 67 M 1 N C@2D 1 + 18 x + 81 x2 è!!!!!! è!!!!!! 1 %%%%%%%%%%% 11 67 2 67 20 $%%%%%%%%%%%%%%%%%%%%%%%% Hypergeometric1F1A− + E== , 1+ , 1 + 18 x + 81 x2 27 27 27 27 H−5L 27 1 è!!!!!! I−11+ 67 M 2 2 27 è!!!!!! I−11+ 67 M 1 39 è!!!!!! I11− 67 M J Sebbene sia lunga, il ramo della soluzione è unico anche in questo caso. Printed by Mathematica for Students 494 1 Daniele Lupo Appunti di Mathematica Come avete visto, Mathematica è in grado di risolvere equazioni differenziali che coinvolgono funzioni speciali, come quelle di Bessel: bes = x ^ 2 ∗ y ''@xD + x ∗ y '@xD + Hx ^ 2 − 5L ∗ y@xD DSolve@8bes, y@1D 0, y '@1D 0; 3<, y@xD, xD êê FullSimplify 3 è!!! è!!! π H−BesselJ@ 5 , xD BesselY@ 5 , 1D + 2 è!!! è!!! BesselJ@ 5 , 1D BesselY@ 5 , xDL== 99y@xD → Plot@y@xD ê. %, 8x, 0, 15<, PlotStyle → 88Red, Thickness@0.01D<<D 5 2.5 2 4 6 8 10 12 14 -2.5 -5 -7.5 -10 -12.5 Graphics Una leggera modifica all'equazione di sopra comporta soluzioni con più funzioni speciali: bes2 = x ^ 2 ∗ y ''@xD + x y '@xD + Hx − 5L ^ 2 ∗ y@xD 0; DSolve@bes2, y@xD, xD êê FullSimplify 99y@xD → − x x5 JC@1D HypergeometricUA C@2D LaguerreLA− 1 , 10 , 2 2 1 , 1 + 10 , 2 2 xE + xEN== Possiamo utilizzare anche altri tipi di coefficienti, che non siano semplicemente polunomi: DSolve@ y ''@xD − Hk ^ 2 + 2 ∗ Sech@xD ^ 2L y@xD == 0, y@xD, xD 99y@xD → C@1D LegendrePA C@2D LegendreQA Printed by Mathematica for Students 1 2 1 2 H + H + è!!! 7 L, k, Tanh@xDE + è!!! 7 L, k, Tanh@xDE== 495 Daniele Lupo Appunti di Mathematica Le equazioni che DSolve è in grado di calcolare possono essere sia omogenee, come abbiamo visto negli esempi precedenti, che inomogenee, come nel caso che segue: trig = DSolve@x ^ 2 y ''@xD + y@xD x ^ 4 + Log@xD, y@xD, xD 1 è!!! 1 è!!! è!!! è!!! x C@1D CosA 3 Log@xDE + x C@2D SinA 3 Log@xDE + 2 2 2 2 1 è!!! 1 i 1 è!!! j 3 Log@xDE + x4 CosA 3 Log@xDE + j13 CosA 2 2 13 k 99y@xD → 2 2 1 è!!! 1 è!!! 3 Log@xDE Log@xD + 13 SinA 3 Log@xDE + 2 2 2 2 1 1 y è!!! è!!! x4 SinA 3 Log@xDE + 13 Log@xD SinA 3 Log@xDE z z== 2 2 { 13 CosA Analogamente a quanto abbiamo fatto precedentemente, possiamo particolareggiare le soluzioni, disegnandone alcune con costanti definite mediante il comando Table: sol = Flatten@Table@8y@xD ê. trig ê. 8C@1D → i, C@2D → j<, 8Hue@i + jD<<, 8i, 0, 10, 2<, 8j, 1, 6, .3<D, 1D; Plot@Evaluate@sol@@All, 1DDD, 8x, 0, 2<, PlotStyle → sol@@All, 2DD , PlotRange → 8Automatic, 8−7, 20<<D 20 15 10 5 0.5 1 1.5 2 -5 Graphics Mathematica è in grado di risolvere equazioni non lineari, anche se ci sono casi in cui non è possibile neanche al programma trovare la soluzione in forma chiusa in nessuna maniera conosciuta fino ad oggi. Un esempio può essere il seguente: eqn = y ''@xD Printed by Mathematica for Students 4 x ∗ y '@xD ^ 2 + y '@xD ^ 2; 496 Daniele Lupo Appunti di Mathematica DSolve@eqn, y, xD 99y → FunctionA8x<, − 1+4 x 2 ArcTanA è!!!!!!!!!!!!!!!! !!!!!!!!! E −1+8 C@1D + C@2DE== è!!!!!!!!!!!!!!!!!!!!!!!!! −1 + 8 C@1D Plot@Evaluate@y@xD ê. % ê. 8C@1D → −2, C@2D → 1<D, 8x, −1, .7<, PlotStyle → 8Blue<D; 1.8 1.6 1.4 1.2 -1 -0.75 -0.5 -0.25 0.25 0.5 0.8 0.6 esplicitamente da x, y' HxL, come in questo caso: Possiamo avere a che fare, in alcuni casi, con equazioni differenziali che non dipendono DSolve@y ''@xD == 2 Exp@3 ∗ y@xDD, y@xD, xD — Solve::ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. More… 99y@xD → LogA− 9y@xD → LogA 9y@xD → LogA H−3L1ê3 J−C@1D + C@1D TanhA 3 2 22ê3 2 1ê3 è!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C@1D Hx + C@2DL2 E N !!!!!!!!!!!!!!!!!!!!! 3 è!!!!!!!!!!!!!!!! C@1D Hx + C@2DL2 E 2 22ê3 31ê3 J−C@1D + C@1D TanhA H−1L2ê3 31ê3 J−C@1D + C@1D TanhA 22ê3 3 2 N 2 1ê3 E=, E=, 2 1ê3 è!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C@1D Hx + C@2DL2 E N E== Data la difficoltà intrinseca nel risolvere equazioi differenziali di ordine superiore, e considerando che DSolve utilzza Solve nella sua pancia, può capitare che le soluzioni coinvolgano la funzione Root, dato che è il modo di dare soluzioni per un'equazione di grado elevato, quando non è possibile trovarla in forma chiusa: Printed by Mathematica for Students 497 Daniele Lupo Appunti di Mathematica DSolve@x ^ 5 y '''''@xD + 6 x ^ 4 ∗ y ''''@xD − 2 ∗ x ^ 3 ∗ y '''@xD − x ^ 2 ∗ y ''@xD + 5 ∗ x ∗ y '@xD + y@xD 99y@xD → xRoot@1−7 #1−#1 2 +#13 &,1D Root@1−7 #1−#12 +#13 &,3D x C@1D + xRoot@1−7 #1−#1 2 C@3D + x 3 2 − è!!!! 5 2 C@4D + x 3 2 + +#13 &,2D è!!!! 5 2 0, y@xD, xD C@2D + C@5D== In questo caso non possiamo trovarci una soluzione simbolica, ma comunque possiamo averla in maniera numerica. Inoltre, consideranto che abbiamo ottenuto una formula sì in maniera non chiusa, ma sempre senza coefficienti numerici, possiamo avere facilmente risultati numerici con precisione arbitraria: f@x_D = y@xD ê. %@@1DD; Notate come abbia utilizzato = invece di := perchè altrimenti avrei sempre la valutazione di % ad ogni nuovo comando, che non rappresenterebbe più la soluzione di DSolve. Se andiamo a mettere un valore alla soluzione in maniera esatta, otteniamo: f@3D 2 3Root@1−7 #1−#1 3 +#13 &,1D C@1D + 3Root@1−7 #1−#1 +#13 &,2D 3 2 3 2 Root@1−7 #1−#12 +#13 &,3D 2 C@3D + 3 − è!!!! 5 2 C@4D + 3 + C@2D + è!!!! 5 2 C@5D Il risultato è restituito sempre in forma simbolica, come nel caso generale. Andiamo invece a sostituire un valore numerico: f@3.D 0.082313 C@1D + 1.16682 C@2D + 31.2355 C@3D + 1.5214 C@4D + 17.7468 C@5D Questa volta abbiamo visto come vengano restituiti i numeri, anche se rimangono i coefficienti C che ovviamente non erano calcolati. Possiamo utilizzare una precisione maggiore, se ci serve: N@f@3D, 50D 0.082313036603068456357530701934319583558912006682794 C@1D + 1.1668222437039807891665184590870851804424075884110 C@2D + 31.235461826749516230168482745226269524617450019969 C@3D + 1.5214024193806499467234230586261890237812609816579 C@4D + 17.746783925183628648066945795420638758100659851792 C@5D Ah, la potenza di Mathematica!!!! In altri casi possiamo avere soluzioni esatte date in maniera simbolica, come questa sottostante che restituisce la soluzione con le funzioni di Airy che compaiono all'interno di integrali Printed by Mathematica for Students 498 Daniele Lupo Appunti di Mathematica ai = y '''@xD − y ''@xD + 5 x y '@xD + 5 y@xD 0; DSolve@ai, y@xD, xD êê FullSimplify i H−1L1ê3 H−1 + 20 xL j j E C@2D + j jAiryAiA 4 52ê3 k H−1L1ê3 H−1 + 20 xL H−1L1ê3 H−1 + 20 xL AiryBiA E C@3D + AiryBiA E 2ê3 45 4 52ê3 x 1 i H−1L1ê3 H−1 + 20 K$4755L y E C@1Dz jH−1L2ê3 −K$4755ê2 π AiryAiA z ‡ − 51ê3 j 2ê3 4 5 k { 1 99y@xD → xê2 K$4755 + AiryAiA ‡ x H−1L2ê3 1 −K$4051ê2 H−1L1ê3 H−1 + 20 xL E 4 52ê3 π AiryBiA H−1L 51ê3 1ê3 H−1+20 K$4051L 4 52ê3 E C@1D y z z== K$4051z z { Sistemi ODE Abbiamo visto finora come lavorare con equazioni differenziali, ma possiamo trattare alla stessa maniera i sistemi di equazioni differenziali, che sono sistemi del tipo: FHY HnL HxL, Y Hn-1L HxL, ..., Y ' HxL, Y HxLL ã 0 dove Y HxL rappresenta un vettore di funzioni incognite: il caso più semplice di sistema, naturalmente, è il caso di sistema lineare, che può essere espresso, per il caso di sistema lineare del primo ordine, nella seguente forma: Y £ HxL AHxL Y HxL +BHxL Dove è presente il vettore delle derivate prime Y ' HxL, quello delle funzioni non derivate Y HxL, ed inoltre è presente la matrice dei coefficienti AHxL, ed il vettore dei termini noti BHxL. Naturalmente il caso si semplifica ulteriormente se la matrice A è a termini costanti, invece che dipendenti dalla variabile x. Il comando DSolve, tuttavia, si aspetta come argomento una lista di equazioni, quindi, se abbiamo la matrice A ed il vettore B, dobbiamo comunque trasformarli nelle rispettive equazioni, prima di dare il sistema in pancia a DSolve. Supponiamo che la nostra matrica A sia la seguente, assieme al vettore B: A=8 8−3, 3, 0<, 85, 0, −2<, 84, −3, −1< <; Printed by Mathematica for Students 499 Daniele Lupo Appunti di Mathematica Adesso occorre creare il vettore delle incognite, che nel nostro caso sono delle funzioni: Y@x_D := 8y1@xD, y2@xD, y3@xD<; Arrivati a questo punto, possiamo crearci il sistema; sisdiff = H#1 8y1 @xD y2 @xD #2L & ∼ MapThread ∼ 8Y '@xD, A.Y@xD< −3 y1@xD + 3 y2@xD, 5 y1@xD − 2 y3@xD, y3 @xD 4 y1@xD − 3 y2@xD − y3@xD< Una volta creato il sistema, possiamo risolverlo: solsis = DSolve@sisdiff, Y@xD, xD êê FullSimplify è!!!!!! 1 − 12 I7+ 37 M x 814 è!!!!!! è!!!!!! 1 è!!!!!! è!!!!!! II333 + 15 37 + H333 − 15 37 L 37 x + 148 2 I13+ 37 M x M C@1D − è!!!!!! 1 è!!!!!! è!!!!!! 2 H74 + 7 37 L C@2D + 148 2 I13+ 37 M x H2 C@2D − C@3DL + H74 − 26 37 L è!!!!!! è!!!!!! è!!!!!! C@3D + 2 37 x HH−74 + 7 37 L C@2D + H37 + 13 37 L C@3DLM, è!!!!!! è!!!!!! 1 è!!!!!! è!!!!!! − 12 I7+ 37 M x y2@xD → I2 I−74 − 29 37 + H−74 + 29 37 L 37 x + 814 è!!!!!! 1 è!!!!!! 148 2 I13+ 37 M x M C@1D + 3 H37 + 9 37 L C@2D + è!!!!!! 1 è!!!!!! 296 2 I13+ 37 M x H2 C@2D − C@3DL + 4 H37 − 2 37 L C@3D + è!!!!!! 1 è!!!!!! è!!!!!! 37 x H3 H37 − 9 37 L C@2D + 4 H37 + 2 37 L C@3DLM, y3@xD → 814 è!!!!!! è!!!!!! è!!!!!! 1 è!!!!!! è!!!!!! − 12 I7+ 37 M x II37 − 101 37 + H37 + 101 37 L 37 x − 74 2 I13+ 37 M x M C@1D + è!!!!!! 1 è!!!!!! è!!!!!! H74 + 40 37 L C@2D − 74 2 I13+ 37 M x H2 C@2D − C@3DL + H370 − 42 37 L è!!!!!! è!!!!!! è!!!!!! C@3D + 37 x HH74 − 40 37 L C@2D + H370 + 42 37 L C@3DLM== 99y1@xD → Vediamo di disegnare queste soluzioni per alcuni valori delle costanti; visualizziamo in tre grafici distinti y1, y2, y3: funz = Flatten@Table@# ê. solsis ê. 8C@1D → a, C@2D → b, C@3D → c<, 8a, 0, 1, .3<, 8b, 0, 1, .3<, 8c, 0, 1, .3<DD & ê@ 8y1@xD, y2@xD, y3@xD<; colori = PlotStyle → Flatten@Table@RGBColor@a, b, cD, 8a, 0, 1, .3<, 8b, 0, 1, .3<, 8c, 0, 1, .3<DD; Printed by Mathematica for Students 500 Daniele Lupo Appunti di Mathematica DisplayTogetherArray@ Plot@Evaluate@#, 8x, −1, 1.3<, colori, AspectRatio → 1, PlotRange → 8−3, 6<DD & ê@ funz D 6 6 6 4 4 4 2 2 2 -1 -0.5 0.5 1 -1 -0.5 -2 0.5 1 -2 -1 -0.5 0.5 1 -2 GraphicsArray Mamma mia, quante soluzioni... Proviamo adesso a fare qualcosa di più semplice, invocando soltanto due equazioni, e tracciando la traiettoria che si ottiene utilizzando come coordinate parametriche le due funzioni ottenute: A = 887, −80<, 85, −5<<; Y@t_D := 8x@tD, y@tD< sisT = H#1 8x @tD #2L & ∼ MapThread ∼ 8Y '@tD, A.Y@tD< 7 x@tD − 80 y@tD, y @tD 5 x@tD − 5 y@tD< solT = DSolve@sisT, Y@tD, tD êê FullSimplify 99x@tD → y@tD → t C@1D Cos@2 t C@2D Cos@2 è!!!!!! 91 tD + è!!!!!! 91 tD + t t è!!!!!! H3 C@1D − 40 C@2DL Sin@2 91 tD , è!!!!!! 91 è!!!!!! H5 C@1D − 6 C@2DL Sin@2 91 tD == è!!!!!! 2 91 Arrivati a questo punto, possiamo utilizzare ParametricPlot per disegnare la traiettoria; tuttavia creerò qualcosa di più attraente, colorando la curva a con un gradiente che dipende dal parametro: traiettoria@t_D = Evaluate@8x@tD, y@tD< ê. solT ê. 8C@1D → 2, C@2D → 4<D@@1DD; Printed by Mathematica for Students 501 Daniele Lupo Appunti di Mathematica Module@ H∗ Parametri del disegno ∗L 8inizio = 0, fine = 3, passo = 1 ê 300<, H∗ Creazione della curva ∗L Show@Graphics@ H∗ Creazione delle linee colorate che compongono il grafico ∗L 8Table@8Thickness@0.006D, Hue@t ê 6D, Line@8traiettoria@t − passoD, traiettoria@tD<D<, 8t, inizio, fine, passo<D<, H∗ Opzioni del disegno ∗L D D AspectRatio → 1, Axes → True, Background → RGBColor@1, 1, .8DD 80 60 40 20 -200 -100 100 200 300 -20 -40 -60 Graphics Possiamo anche vedere come varia la spirale al variare di un singolo parametro: vediamo i due grafici distinti, dove prima variamo C[1], e dopo C[2]: Printed by Mathematica for Students 502 Daniele Lupo Appunti di Mathematica traiettoriamultiplaC1 = Table@8x@tD, y@tD< ê. solT@@1DD ê. 8C@1D → a, C@2D → 4<, 8a, 0, 20, 4<D; traiettoriamultiplaC2 = Table@ 8x@tD, y@tD< ê. solT@@1DD ê. 8C@1D → 10, C@2D → a<, 8a, 0, 20, 4<D; DisplayTogetherArray@8 ParametricPlot@Evaluate@#D, 8t, 0, 3<, PlotRange → All, PlotStyle → 8Red, Green, Orange, Blue, Cyan<D & ê@ 8traiettoriamultiplaC1, traiettoriamultiplaC2< <D 400 300 200 100 100 50 -200 200 400 -50 -100 -1000 -500 -100 -200 -300 500 1000 1500 GraphicsArray I coefficienti del sistema lineare, comunque, possono anche essere non costanti, cioè dipendere dalla variabile indipendente anch'essi, come in questo caso, in cui le due funzioni i vettori della matrice sono ortogonali fra di loro: A = 88t, Cos@tD<, 8− Cos@tD, t<<; Y@t_D := 8x@tD, y@tD< noncostante = MapThread@#1 8x @tD #2 &, 8Y '@tD, A.Y@tD<D t x@tD + Cos@tD y@tD, y @tD −Cos@tD x@tD + t y@tD< DSolve@noncostante, Y@tD, tD 99x@tD → y@tD → t2 2 t2 2 C@1D Cos@Sin@tDD + C@2D Cos@Sin@tDD − t2 2 C@2D Sin@Sin@tDD, t2 2 C@1D Sin@Sin@tDD== Anche nel caso dei sistemi di equazioni differenziali possiamo avere casi che non sono omogenei, come abbiamo trattato finora, cioè casi in cui il vettore B non è nullo: Printed by Mathematica for Students 503 Daniele Lupo Appunti di Mathematica A = 887, −8<, 81, −1<<; B = 8t, Gamma@t ê 5D<; X@t_D = 8x@tD, y@tD<; inomogeneo = MapThread@#1 9x @tD #2 &, 8X '@tD, A.X@tD + B<D t + 7 x@tD − 8 y@tD, y @tD GammaA t E + x@tD − y@tD= 5 DSolve@inomogeneo, 8x, y<, tD êê FullSimplify 99x → FunctionA8t<, è!!!! è!!!! è!!!! è!!!! 1 è!!! è!!! I I3−2 2 M t − 2 I3−2 2 M t + I3+2 2 M t + 2 I3+2 2 M t M C@1D + 2 è!!!! è!!!! è!!! 2 I I3−2 2 M t − I3+2 2 M t M C@2D + è!!!! è!!!! è!!!! è!!!! 1 è!!! è!!! I I3−2 2 M t − 2 I3−2 2 M t + I3+2 2 M t + 2 I3+2 2 M t M 2 t è!!!! è!!!! 1 è!!! −I3−2 2 M K$1025371 − 2 −I3−2 2 M K$1025371 + ‡ J2 I 1 è!!!! è!!!! è!!! −I3+2 2 M K$1025371 + 2 −I3+2 2 M K$1025371 M K$1025371 + è!!!! è!!!! K$1025371 è!!! EN 2 I −I3−2 2 M K$1025371 − −I3+2 2 M K$1025371 M GammaA 5 è!!!! è!!!! è!!! K$1025371 + 2 I I3−2 2 M t − I3+2 2 M t M ti è!!!! è!!!! I −I3−2 2 M K$1029249 − −I3+2 2 M K$1029249 M K$1029249 j j j + − ‡ j è!!! j 4 2 1 k è!!!! è!!!! è!!!! 1 è!!! I −I3−2 2 M K$1029249 + 2 −I3−2 2 M K$1029249 + −I3+2 2 M K$1029249 − 2 z K$1029249 y è!!! −I3+2 è!!!! 2 M K$1029249 z 2 M GammaA Ez z z K$1029249E, 5 { y → FunctionA8t<, − I I3−2 è!!!! 2Mt − è!!!! è!!!! 1 è!!! I I3−2 2 M t + 2 I3−2 2 M t + 2 è!!!! i 1 2Mt jI I3−2 è!!!! − I3+2 2 M t M j è!!! j 4 2 k 1 ‡ J2 I 1 t −I3−2 è!!! 2 4 è!!!! 2 M K$1025371 −I3+2 − I3+2 è!!!! 2Mt I3+2 è!!!! 2Mt è!!! 2 è!!! 2 è!!!! 2 M K$1025371 −I3−2 M C@1D − è!!! 2 + I3+2 è!!!! 2Mt è!!!! 2 M K$1025371 M K$1025371 + è!!! 2 I + M C@2D − −I3+2 −I3−2 è!!!! 2 M K$1025371 è!!!! 2 M K$1025371 y K$1025371 z EN K$1025371z z+ 5 { è!!!! è!!!! è!!!! è!!!! 1 è!!! è!!! I I3−2 2 M t + 2 I3−2 2 M t + I3+2 2 M t − 2 I3+2 2 M t M 2 −I3+2 ti è!!!! 2 M K$1025371 M GammaA è!!!! è!!!! − −I3−2 2 M K$1029249 − −I3+2 2 M K$1029249 M K$1029249 j j− I + j ‡ j è!!! j 4 2 1 k è!!!! è!!!! è!!!! 1 è!!! I −I3−2 2 M K$1029249 + 2 −I3−2 2 M K$1029249 + −I3+2 2 M K$1029249 − 2 y K$1029249 z è!!! −I3+2 è!!!! 2 M K$1029249 z M GammaA Ez 2 z z K$1029249E== 5 { Printed by Mathematica for Students 504 + Daniele Lupo Appunti di Mathematica Anche in questo caso possiamo visualizzare soluzioni partcolari nel solito modo. Dato che però la soluzione non è data in forma chiusa, il calcolo del grafico è computazionalmente intensivo, e potrebbe impiegarci parecchio prima di completare il disegno. Oltre al grafico è mostrato il tempo impiegato sul mio AthlonXP64 3200, con sistema operativo WinXP Pro a 32 bit: particularsol = 8x@tD, y@tD< ê. %@@1DD ê. 8C@1D → 1, C@2D → 2<; Plot@Evaluate@particularsolD, 8t, 0, Pi<, PlotStyle → 88Blue<, 8Dashing@80.02<D, Orange<<D êê Timing 0.5 1 1.5 2 2.5 3 -500000 6 -1× 10 6 -1.5× 10 6 -2× 10 8599.25 Second, Graphics < Cavolo!!! Dopo aver aspettato tutto questo tempo, poteva almeno uscire qualcosa di più carino... Ma non ho intenzione di rifare i calcoli un'altra volta. Accontentatevi di questo. Adesso, comunque, dopo questo 'grande' finale, possiamo passare alla sezione successiva. PDE La caratteristica che differenzia le PDE dalle ODE è il fatto che in questo caso abbiamo più variabili indipendenti, invece di una soltanto. Questo complica parecchio la risoluzione, specialmente nei casi non lineari dove, a parte pochi casi, possiamo tranquillamente metterci le mani in testa e passare alla sezione successiva dove si parla della risoluzione numerica. Ma per ora tralasciamola... Vediamo subito un esempio, come il seguente: pde = D@u@x, yD, xD + y D@u@x, yD, yD y uH0,1L @x, yD + uH1,0L @x, yD Printed by Mathematica for Students y Sin@xD 505 y Sin@xD Daniele Lupo Appunti di Mathematica pdesol = DSolve@pde, u, 8x, y<D 99u → FunctionA8x, y<, 1 H−y Cos@xD + y Sin@xD + 2 C@1D@ 2 −x yDLE== Vediamo come sia stato necessario definire mediante l'operatore D la derivata perchè, come avevamo detto all'inizio, in questo caso dobbiamo specificare rispetto a quale variabile è calcolata la derivata, specialmente nelle derivate miste. Anche in questo caso, come nel caso unidimensionale, Mathematica lo interpreta tramite il comando Derivative: InputForm@pdeD y*Derivative[0, 1][u][x, y] + Derivative[1, 0][u][x, y] == y*Sin[x] Ricordate? L'abbiamo visto a proposito della risoluzione manuale delle ODE a coefficienti costanti, quando siamo andati a trovarci l'equazione caratteristica... Vi ricordo anche che, in questo caso, abbiamo le C[n] che adesso diventano delle funzioni arbitrarie, anche se l'argomento è specificato. Nel nostro caso abbiamo 2 C@1D@‰-x yD, dove l'argomento è unico, ma è dato come funzione delle due variabili indimententi x, y. Se vogliamo trovarci una soluzione particolare, dobbiamo sostituire a C@nD al posto di un valore, l'head di una funzione: pdeparticolare@x_, y_D := u@x, yD ê. pdesol@@1DD ê. C@1D → Cos In questo modo abbiamo creato la funzione. Notate come il coseno abbia lo stesso argomento di C@1D, cosa che fra l'altro era scontata. Adesso, tanto per farvi contenti, plottiamo questqa funzione: Printed by Mathematica for Students 506 Daniele Lupo Appunti di Mathematica DisplayTogetherArray@8 Plot3D@pdeparticolare@x, yD, 8x, −3, 3<, 8y, −5, 5<, PlotPoints → 100, Mesh → False, ViewPoint −> 83.252, −1.028, 2.732<D, DensityPlot@pdeparticolare@x, yD, 8x, −3, 3<, 8y, −7, 7<, PlotPoints → 100, Mesh → FalseD <D 6 4 2 -2 2 0 4 2 0 -2 -4 - 0 2 -4 -2 -4 4 2 0 -2 -6 -3 -2 -1 0 1 2 3 GraphicsArray Le equazioni PDE sono di solito classificate in tre tipi: lineati, quasilineari e non lineari. Per cominciare, vediamo un esempio di equazione omogenea a coefficienti costanti; omo = 2 ∗ D@u@x, yD, xD + 9 ∗ D@u@x, yD, yD + u@x, yD u@x, yD + 9 uH0,1L @x, yD + 2 uH1,0L @x, yD 0 0 sol = DSolve@omo, u, 8x, y<D 99u → FunctionA8x, y<, −xê2 C@1DA 1 H−9 x + 2 yLEE== 2 Non è veramente una delle equazioni più complicate che abbiamo visto finora... Vediamo delle soluzioni che si ottengono andando a sostituire alla funzione arbitraria C@1D delle funzioni specifiche: omop = u@x, yD ê. sol@@1DD ê. 88C@1D → Sin<, 8C@1D → ArcTan<, 8C@1D → H# &L<< 9 −xê2 SinA 1 H−9 x + 2 yLE, 2 −xê2 ArcTanA Andiamo a plottare adesso queste funzioni: Printed by Mathematica for Students 507 1 1 H−9 x + 2 yLE, 2 2 −xê2 H−9 x + 2 yL= Daniele Lupo Appunti di Mathematica DisplayTogetherArray@ Plot3D@#, 8x, −4, 4<, 8y, −4, 4<, PlotPoints → 40D & ê@ omop D 4 2 0 -2 -4 -4 4 4 2 -2 0 2 0 -2 4 -4 10 5 0 -4 4 4 2 -2 0 2 0 -2 4 -4 150 100 50 0 -4 4 4 2 -2 0 2 0 -2 4 -4 GraphicsArray Come potete vedere la funzione arbitrarià da una maggiore libertà alla funzione, rispetto a dei coefficienti indeterminati... Un esempio interessante di PDE è l'equazione di trasporto: ∑u ∑u ÅÅÅÅ ÅÅ + k ÅÅÅÅ ÅÅ ∑x ∑y 0 dove k rappresenta una quantità costante. Come potete vedere è un'equazione lineare, ed è una delle più semplici PDE che si possano incontrare (forse la più semplice?). trasporto = DSolve@D@u@x, yD, xD + κ D@u@x, yD, yD 88u@x, yD → C@1D@y − x κD<< 0, u@x, yD, 8x, y<D Come potete vedere è una soluzione abbastanza semplice, ed è praticamente data da qualsiasi funzione avente l'argomento specificato: trabel = u@x, yD ê. trasporto@@1DD ê. C@1D → HBesselJ@3, #D &L BesselJ@3, y − x κD Printed by Mathematica for Students 508 Daniele Lupo Appunti di Mathematica DisplayTogetherArray@8 DensityPlot@trabel ê. κ → 8, 8x, 0, 8<, 8y, 0, 8<, PlotPoints → 100, Mesh → FalseD, Plot3D@trabel ê. κ → 8, 8x, 0, 8<, 8y, 0, 8<, PlotPoints → 870, 40<D <D 8 6 0.4 0.2 0 -0.2 -0.4 0 4 8 6 4 2 2 2 4 6 80 0 0 2 4 6 8 GraphicsArray Adesso andiamo a vedere un semplice esempio di equazione differenziale PDE inomogenea, con l'aggiunta di coefficienti non costanti, data dalla seguente equazione: inom = D@u@x, yD, xD + x D@u@x, yD, yD x uH0,1L @x, yD + uH1,0L @x, yD Sin@xD + x x + Sin@xD solino = DSolve@inom, u, 8x, y<D 99u → FunctionA8x, y<, 1 1 Jx2 − 2 Cos@xD + 2 C@1DA H−x2 + 2 yLENE== 2 2 Anche in questo caso la soluzione è abbastanza semplice, seppure possa essere resa complicata a piacere scegliendo la funzione C@1D: inopar = u@x, yD ê. solino@@1DD ê. 88C@1D → Sin<, 8C@1D → ArcTan<, 8C@1D → H# &L<< 9 1 1 Jx2 − 2 Cos@xD + 2 SinA H−x2 + 2 yLEN, 2 2 1 1 1 Jx2 + 2 ArcTanA H−x2 + 2 yLE − 2 Cos@xDN, H2 y − 2 Cos@xDL= 2 2 2 Printed by Mathematica for Students 509 Daniele Lupo Appunti di Mathematica DisplayTogetherArray@8 ContourPlot@#, 8x, −6, 6<, 8y, −6, 6<, PlotPoints → 50D & ê@ inopar, Plot3D@#, 8x, −8, 8<, 8y, −8, 8<, PlotPoints → 30D & ê@ inopar <D 6 6 6 4 4 4 2 2 2 0 0 0 -2 -2 -2 -4 -4 -4 -6 -6 -4 -2 0 2 4 30 20 10 0 6 -6 -6 -4 -2 0 5 30 20 10 0 0 -5 0 2 4 -6 -6 -4 -2 0 5 5 0 -5 0 -5 0 -5 6 -5 4 6 5 0 -5 0 5 5 2 -5 5 GraphicsArray I gradi di liberta, nelle PDE, è generalmene maggiore, come possiamo vedere da questi esempi. La caterogira successiva è data dalle equazioni quasi lineari, che ha termini lineari nelle derivate prime rispetto ad x ed ad y, ma che possono contenere termini non lineari dati da queste derivate e della funzione non derivata, come in questo esempio: ql = D@u@x, yD, xD + x ^ 3 D@u@x, yD, yD x3 uH0,1L @x, yD + uH1,0L @x, yD u@x, yD ^ 2 + x x + u@x, yD2 solql = DSolve@ql, u, 8x, y<D êê FullSimplify 99u → FunctionA8x, y<, J−2 x3ê2 BesselJA− 2 2 x3ê2 2 2 x3ê2 , E + 2 x3ê2 BesselJA , E 3 3 3 3 1 2 x3ê2 1 E+ C@1DA H−x4 + 4 yLEN í J2 Jx BesselJA , 3 3 4 1 2 x3ê2 1 x BesselJA− , E C@1DA H−x4 + 4 yLENNE== 3 4 3 Printed by Mathematica for Students 510 Daniele Lupo Appunti di Mathematica Possiamo vedere come anche in questo caso cambia radicalmente la funzione se andiamo a sostituire nel nostro caso due funzioni distinte, ad esempio il seno e l'identità: DisplayTogetherArray@ DensityPlot@#, 8x, 0.1, 4<, 8y, −15, 20<, PlotPoints → 400, Mesh → FalseD & ê@ Evaluate@u@x, yD ê. solql@@1DD ê. 88C@1D → Sin<, 8C@1D → H# &L<<D D êê Timing 20 20 15 15 10 10 5 5 0 0 -5 -5 -10 -10 -15 -15 1 833.797 Second, 2 3 4 1 2 3 4 GraphicsArray < Quello che possiamo notare, come principale differenza dalle lineari alle quasilineari, è la presenza di discontinuità nella soluzione, che si possono vedere come passaggi bruschi dal bianco al nero, che sono caratteristiche di quest'ultimo tipo di equazioni, sebbene il metodo risolutivo sia per certi versi simile al caso lineare. Le cose diventano decisamente più toste quando abbiamo a che fare con equazioni PDE non lineari, nel senso che la funzione F(u,p,q) 0 Dove p, q rappresentano le derivate prime rispetto alle due variabili (sempre considerando equazioni differenziali del primo ordine ed a due variabili), ed F in questo caso non è lineare. Possiamo definire per semplicità, le nostre derivate, e la funzione, in questa maniera: z := u@x, yD; p := D@u@x, yD, xD; q := D@u@x, yD, yD; Printed by Mathematica for Students 511 Daniele Lupo Appunti di Mathematica In questa maniera non dobbiamo per forza riscrivere tutto quanto quando abbiamo bisogno di usare le derivate nelle nostre equazioni. Effettivamente potevamo farlo dall'inizio, ma ci ho pensato solamente adesso... :-( Spero che voi l'abbiate fatto e che siate stati più intelligenti di me!!! Comunque, vediamo subito un veloce esempio: nonlin = p q z^2 uH0,1L @x, yD uH1,0L @x, yD u@x, yD2 DSolve@nonlin, u, 8x, y<D — DSolve::nlpde : Solution requested to nonlinear partial differential equation. Trying to build a complete integral. 99u → FunctionA8x, y<, − y è!!!!!!!!!!!!! C@1D −x è!!!!!!!!!!!! C@1D − C@2D è!!!!!!!!!!!!! C@1D E== In questo caso Mathematica ci avverte con gentilezza e garbo che abbiamo a che fare con PDE non lineari, prima di creare l'integrale completo, che serve per calcolare simbolicamente la PDE, quando obbiamente ci riesce. In questo caso, inoltre, abbiamo ottenuto dei valori di C@1D, C@2D che non sono delle funzioni, ma delle costanti. Possiamo anche parametrizzare il tutto, in questo esempio, ponendo una costante in funzione dell'altra. In questo modo si ottiengono delle superfici dipendenti da un parametro, ovvero una famiglia di superfici, il cui inviluppo è pure soluzione della PDE, che è dato dall'integrale completo. Quest'ultimo include qualsiasi elemento della famiglia di curve ottenuta variando arbitrariamente i paramentri, l'inviluppo che abbiamo appenda detto, l'inviluppo dell'intera famiglia (integrale singolare), e altri integrali completi della PDE si possono dal processo di formazione degli inviluppi. Come potete vedere, quindi, la teoria delle PDE è alquanto corposa e cummattusa (i siciliani mi capiranno). Tuttavia ci possono essere casi in cui possono essere comunque risolte. Un caso utile si ha quando possiamo applicare la separazione delle variabili, in modo da integrare in maniera indipendente le due parti dell'equazione ottenuta: sepa = DSolve@Hy p − x qL ^ 2 + Hx p + y qL — 6, z, 8x, y<D DSolve::nlpde : Solution requested to nonlinear partial differential equation. Trying to build a complete integral. 99u@x, yD → − ArcTanA Printed by Mathematica for Students y è!!!!!!!!!!! è!!!!!!!!!!!!!!! è!!!!!!!!!!!!!!! E C@1D + C@2D + 6 LogA x2 + y2 E + C@1D LogA x2 + y2 E== x 512 Daniele Lupo Appunti di Mathematica Dato che la funzione che abbiamo ottenuto è complessa, andiamo a visualizzare la parte reale e quella immaginaria, rispettivamente, fissando dei parametri per permetterci di avere una risoluzione numerica DisplayTogetherArray@ Plot3D@#, 8x, −4, 4<, 8y, −4, 4<D & ê@ HH8Re@#D, Im@#D< &L ê@ Evaluate@u@x, yD ê. sepa@@1DD ê. 8C@1D → 1, C@2D → 2<DL D 10 5 0 -4 4 4 2 0 -2 1 0 -1 2 0 -4 4 -2 -2 0 4 2 -2 0 2 4 -4 4 -4 GraphicsArray Notate la discontinuità tipica dell' ArcTan nel grafico della fase... Il passo successivo, adesso, è quello delle PDE del secondo ordine. Come prima, partiamo dal caso più semplice, che è dato da quello lineare. In questo caso, essendo la funzione a più variabili (nel nostro caso due), dovranno comparire anche le derivate miste: ∑ u ∑ u ∑u ∑u ∑ u a ÅÅÅÅ ÅÅÅÅÅ + b ÅÅÅÅÅÅÅÅ ÅÅÅÅÅ + c ÅÅÅÅ ÅÅÅÅÅ + d ÅÅÅÅ ÅÅ + e ÅÅÅÅ ÅÅ + f u ∑x2 ∑x ∑y ∑y2 ∑x ∑y 2 2 2 g Possiamo distunguere, qua come in tutti li altri casi, le equazioni omogenee da quelle inomogenee, a seconda se g è uguale o diverso da zero. Richiamando velocemente le proprietà delle PDE lineari del secondo ordine, possiamo dividerle in tre categorie principali, che dipendono dalle derivate di ordine due, in particolar modo dai loro coefficienti. Se abbiamo b2 - 4 a c < 0, allora l'integrale è detto ellittico. Come caso particolare abbiamo l'equazione di Laplace, che si ottiene per a = 1, b = 0, c = 1. Poi abbiamo le equazioni parabolice, nel caso in cui b2 - 4 a c = 0, l'equazione è detta parabolica, come l'equazione per la diffuzione del calore. Infine, per b2 - 4 a c > 0 l'equazione viene chiamata iperbolica, come l'equazione d'onda. Come potete vedere, nelle PDE del secondo ordine, anche solo fermandoci alle lineari, sorgono molti casi importanti. La soluzione simbolica può essere trovata in non tutti i casi, per esempio si possono trovare nei casi in cui l'equazione è omogenea, e i coefficienti delle derivate prime è nullo, Printed by Mathematica for Students 513 Daniele Lupo Appunti di Mathematica come appunto l'equazione di Laplace; in questo caso la soluzione è quella generica, dato che manca ogni condizione al contorno: laplace = D@u@x, yD, 8x, 2<D + D@u@x, yD, 8y, 2<D 0; DSolve@laplace, u, 8x, y<D 88u → Function@8x, y<, C@1D@ x + yD + C@2D@− x + yDD<< Possiamo vedere un altro esempio di PDE dove compaiono solamente le derivate di ordine 2: pde2 = 8 D@u@x, yD, 8x, 2<D + 9 D@u@x, yD, x, yD + 9 D@u@x, yD, 8y, 2<D 9 uH0,2L @x, yD + 9 uH1,1L @x, yD + 8 uH2,0L @x, yD 0 0 solpde2 = DSolve@pde2, u, 8x, y<D 99u → FunctionA8x, y<, 1 1 è!!!!!! C@1DA H−9 + 3 23 L x + yE + C@2DA H−9 − 3 16 16 è!!!!!! 23 L x + yEE== pde2 ê. solpde2 êê Simplify 8True< Come possiamo vedere, le soluzioni sono rispettate. Anche in questo caso la funzione che abbiamo ottenuto ha valori sia reali che immaginari: Printed by Mathematica for Students 514 Daniele Lupo Appunti di Mathematica DisplayTogetherArray@ Plot3D@#, 8x, −6, 6<, 8y, −6, 6<, PlotPoints → 45, PlotRange → AllD & ê@ HH8Re@#D, Im@#D< &L ê@ Evaluate@u@x, yD ê. solpde2@@1DD ê. 8C@1D → Sin, C@2D → Log<DL D 100 50 0 -50 -100 -5 -2.5 5 2.5 0 -2.5 0 2.5 5 5 2.5 0 -2.5 -5 -2.5 0 2.5 -5 5 GraphicsArray Printed by Mathematica for Students 100 50 0 -50 -100 515 -5 Daniele Lupo Appunti di Mathematica DAE Molti sistemi naturali possono essere modellati tramite un insieme di equazioni che sono sia differenziali che algebriche: come non citare i circuiti elettronici, per chi fa ingegneria elettronica come me? Dalle leggi di Kirkoff, che sono algebriche, alle cariche e scariche dei condensatori ed induttori, che sono equazioni differenziali, ottenendo una bellissima ed affascinante DAE... ...................... Meglio se me ne vado a letto, per stasera....... .....zzzzzzzzzzzzzzzzzzzzzzzz...... ....................... Arieccoci, qua, per cominciare di nuovo. Di cosa parlavamo? delle letterine?!?!? Ah, ora ricordo... Bene, rientiramo nel nostro mondo fatto di matematica e di tante figurine colorate, che non sono quelle dei calciatori... Abbiamo appena definito le DAE, facendo un piccolo esempio su un sistema che le sfruttasse. Come problema analitico sono abbastanza complicate, specialmente se andiamo a considerare i casi non lineari, che possono essere veramente delle bruttissime bestie. Un caso in cui le DAE possono essere risolte è il caso dei sistemi lineari, scritti nella seguente maniera: A.x£ HtL + B.xHtL F Le matrici sono funzioni del parametro t, come pure il vettore F. Si distguono anche qua, tanto per cambiare, le equazioni omogenee da quelle inomogenee. Se A e B sono costanti, DSolve può essere in grado di trovare la soluzione. Inoltre, in questo tipo di problema bisogna porre particolare attenzione alle condizioni iniziali del problema, per assicurarci la soluzione. Possiamo vedere qua un semplicissimo esempio di DAE: dae = 8 x '@tD − 3 y@tD 0, x@tD − y@tD 0 <; Printed by Mathematica for Students 516 Daniele Lupo Appunti di Mathematica daesol = DSolve@dae, 8x, y<, tD 99x → FunctionA8t<, 1 4 3t C@1DE, y → FunctionA8t<, 1 4 3t C@1DE== Possiamo scrivere la corrispondente soluzione equazione in foma non omogenea, per compilcarci un tantinello la vita: inodae = 8 x '@tD + y@tD t, x@tD − y@tD Cos@tD <; inodaesol = DSolve@inodae, 8x, y<, tD 99x → FunctionA8t<, 1 H 4 −t FunctionA8t<, −Cos@tD + C@1D + 2 H−2 + 2 t + Cos@tD + Sin@tDLLE, y → 1 H 4 −t C@1D + 2 H−2 + 2 t + Cos@tD + Sin@tDLLE== Come possiamo vedere il risultato è dato dalla soluzione omogenea, più una soluzione particolare di quella inomogenea, come dice la teoria delle equazioni differenziali riguardo, appunto, alle eq. inomogenee. Possiamo anche imporre delle opportune condizioni iniziali per le nostre funzioni: inodaesolpar = DSolve@Flatten@8inodae, x@3D 99x → FunctionA8t<, H−4 t +2 −t t − 3 Cos@3D + t Cos@tD − 1 −t H4 3 + 2 t − 2 t t + y → FunctionA8t<, − 2 t Cos@tD + 3 Sin@3D − t Sin@tDLE== 3 −2 1 2 0<D, 8x, y<, tD êê FullSimplify t In questo modo possiamo plottare le funzioni ottenute: Printed by Mathematica for Students 517 3 3 Sin@3D + Cos@3D + t Sin@tDLE, Daniele Lupo Appunti di Mathematica Plot@Evaluate@8x@tD, y@tD, x@tD − y@tD< ê. inodaesolparD, 8t, 0, 5 Pi<, PlotRange → 8−5, 15<, PlotStyle −> 88Red<, 8Blue<, 8Thickness@0.01D, Green<<D 15 12.5 10 7.5 5 2.5 2.5 5 7.5 10 12.5 15 -2.5 -5 Graphics Possiamo vedere come la differenza delle due funzioni sia effettivamente il coseno che avevamo imposto nell'equazione algebrica della DAE. Possiamo trattare le DAE anche in generale, volendo. Per esempio, potremmo volere che la funzione della parte inomogenea sia indefinita, come nel caso che segue: indef = 8 x ''@tD + y@tD == f@tD, 2 x@tD == h@tD <; DSolve@indef, 8x, y<, tD 99x → FunctionA8t<, h @tD h@tD E== E, y → FunctionA8t<, f@tD − 2 2 In questo caso possiamo vedere che nella soluzione compaiono le funzioni arbitrarie. Possiamo notare anche un'altra cosa, però, cioè la mancanza di costanti arbitrarie, il che ci fa capire che nella soluzione del nostro problema non abbiamo gradi di libertà. Possiamo tentare di risolvere il caso in cui le funzioni diventano particolari: def = 8 x ''@tD + y@tD == Cos@tD, 2 x@tD == BesselJ@3, tD <; Printed by Mathematica for Students 518 Daniele Lupo Appunti di Mathematica sol = DSolve@def, 8x, y<, tD 99x → FunctionA8t<, 1 BesselJ@3, tDE, 2 1 1 y → FunctionA8t<, J H−BesselJ@1, tD + BesselJ@3, tDL + 4 2 1 HBesselJ@3, tD − BesselJ@5, tDLN + Cos@tDE== 2 Possiamo notare la corrispondenza fra le fue soluzioni che abbiamo ottenuto, nel caso generale ed in quello particolare. Possiamo plottare il grafico parametrico nello stesso modo di come abbiamo fatto prima tr@t_D = Evaluate@8x@tD, y@tD< ê. solD@@1DD; Module@ 8inizio = 0, fine = 100, passo = 1 ê 300, scalahue = 200<, Show@Graphics@ 8Table@8Thickness@0.003D, Hue@t ê scalahueD, Line@8tr@t − passoD, tr@tD<D<, 8t, inizio + passo, fine, passo<D<, AspectRatio → 1, Axes → True, Background → RGBColor@1, 1, .8DD D D 1 0.5 -0.15 -0.1 -0.05 0.05 -0.5 -1 Graphics Printed by Mathematica for Students 519 0.1 0.15 0.2 Daniele Lupo Appunti di Mathematica Possiamo fare adesso un passetto avanti, ed andare a considerare, per esempio, un problema DAE del terzo ordine, imponendo pure le condizioni iniziali, date sulla x: Clear @x, y, zD terzo = 8 x '''@tD − y@tD + z@tD x '@tD + z@tD − t 0, z@tD − y@tD 0, x@0D 1, x '@0D 4, x ''@0D 2 <; 0, solterzo = DSolve@terzo, 8x, y, z<, tD 88x → Function@8t<, 1 + 4 t + t2 D, y → Function@8t<, −4 − tD, z → Function@8t<, −4 − tD<< Printed by Mathematica for Students 520 Daniele Lupo Appunti di Mathematica ParametricPlot3D@Evaluate@ 8x@tD, y@tD, x@tD − y@tD, 8BlueViolet, Thickness@0.01D<< ê. solterzoD, 8t, −5, 5 Pi<, BoxRatios → 81, 1, 1<, ViewPoint −> 80.279, −3.548, 2.851<D 0 -5 -10 -15 5 -20 300 200 100 0 100 0 200 300 Graphics3D Il problema diventa invece di più difficile soluzione, diciamo impossibilile nella maggior parte dei casi, quando i coefficienti non sono costanti. In questo caso si rendono necessari metodi di risoluzione numerica. Printed by Mathematica for Students 521 Daniele Lupo Appunti di Mathematica Problema del valore al contorno Fino ad adesso abbiamo cercato, in generale, soluzioni generiche per le nostre equazioni differenziali. Queste danno idea della struttura della soluzione, ma in pratica si è interessati ad una particolare soluzione, che soddisfi la condizione al contorno specifica per il nostro problema. Possiamo imporre due tipi di condizioni: le condizioni iniziali e le condizioni al contorno. Per il primo tipo imponiamo che ad un determinato valore iniziale del parametro t (consideriamo ad esempio l'istante iniziale dell'evoluzione di un sistema), siano date delle specifiche condizioni sul valore della funzione, e sulle due derivate. Invece, per le condizioni al contorno, specifichiamo il valore che deve assumere la funzione in particolari punti del dominio della funzione stessa, per esempio il valore che deve assumere all'inizio ed alla fine dell'intervallo di spazio che prendiamo in considerazione. Analiticamente possiamo trattare allo stesso modo entrambi i tipi di restrizioni della soluzione, che chiamiamo quindi valori al contorno. La soluzione nei casi lineari è abbastanza semplice, mentre non lo è altrettanto nei casi non lineari, per i quali possono comparire più rami, per i quali soltanto alcuni sono in grado di soddisfare precise condizioni al contorno che andiamo a specificare nel nostro problema. Prendendo come sempio un caso lineare, possiamo avere un'equazione del seguente tipo: lineare = y '@tD + t ^ 2 y@tD t2 y@tD + y @tD t t La soluzione generale sarà data dal risultato dato da DSolve, senza specificare nessuna condizione al contorno: generale = DSolve@lineare, y, tD 99y → FunctionA8t<, − t3 3 C@1D − − t3 3 2 t3 3 , − 3 31ê3 H−t3 L2ê3 t2 GammaA E E== Possiamo notare come, essendo del primo ordine, compaia soltanto un coefficiente di indeterminazione. Supponiamo, adesso, di volere una precisa condizione al contorno, per esempio specifichiamo il valore che deve assumere la funzione quando si ha t = 0, per esempio 3: Printed by Mathematica for Students 522 Daniele Lupo Appunti di Mathematica particolare = DSolve@8lineare, y@1D 99y → FunctionA8t<, t3 1 J − 3 J18 1ê3 t − 3 6t 2 32ê3 H−t3 L 1ê3 3<, y, tD 2 1 2 1 , − E − 32ê3 t GammaA , − E + 3 3 3 3 31ê6 t GammaA GammaA t3 2 ,− ENNE== 3 3 y@1D ê. particolare@@1DD êê Simplify 3 Abbiamo trovato il caso particolare. Tuttavia, possiamo anche creare delle condizioni al contorno simboliche; possiamo ad esempio supporre che per t = 0 la la funzione assuma un generico valore a: particolareα = DSolve@8lineare, y@1D 99y → FunctionA8t<, t3 1 J − 3 J6 1ê3 t α − 3 6t 2 32ê3 H−t3 L 1ê3 α<, y, tD 31ê6 t GammaA GammaA 2 1 2 1 , − E − 32ê3 t GammaA , − E + 3 3 3 3 2 t3 ,− ENNE== 3 3 Possiamo notare che, non essendo il coefficiente semplicemente additivo, non basta andare a sostituire la condizione iniziale al coefficiente per ottenere il risultato, ma abbiamo bisogno di ulteriori elaborazioni: f@t_D = y@tD ê. particolareα@@1DD 1 J 6t − t3 3 J6 1ê3 tα−3 32ê3 t GammaA 31ê6 t GammaA 2 1 , − E− 3 3 1ê3 2 1 2 t3 , − E + 2 32ê3 H−t3 L GammaA , − ENN 3 3 3 3 f@1D êê Simplify α Come vedete, tutto persegue il grande cerchio della vita... Printed by Mathematica for Students 523 Daniele Lupo Appunti di Mathematica Plot@ Evaluate@Table@f@tD ê. α → αvalue, 8αvalue, −1, 5<DD , 8t, 0, 5<, PlotStyle → 8Blue, Red, Green, Cerulean, Coral, BlueViolet, Banana<, Background → LightBeigeD 6 4 2 1 2 3 4 5 -2 Graphics Possiamo vedere un altro esempio, questa volta con un'equazione del secondo oridine; in questa maniera compariranno due costanti di indeterminazione: lineare2 = y ''@tD + 3 y '@tD − y@tD −y@tD + 3 y @tD + y @tD t t generale = DSolve@lineare2, y, tD 99y → FunctionA8t<, −3 − t + I− 3 2 − è!!!! 13!!! 2 Mt C@1D + I− 3 2 + è!!!! 13!!! 2 Mt C@2DE== Come possiamo vedere sono spuntati, in questo caso, sia C@1D che C@2D; di conseguenza abbiamo bisogno di due condizioni al contorno per poter risolvere questa equazione differenziale per un caso particolare. Possiamo dare le due condizioni nei due modi citati: possiamo specificare il valore che deve assumere la funzione in due punti distinti, oppure specificare il valore che deve assumere in un punto, ed il valore della derivata in un punto, che possono pure coincidere. Un esempio del primo caso può essere: Printed by Mathematica for Students 524 Daniele Lupo Appunti di Mathematica DSolve@8lineare2, y@0D 0, y@2D 99y → −1 + FunctionA8t<, 3 2 è!!!!!! 13 +I− 3 2 − è!!!!!!! 13 2 è!!!!!! 2 13 1 Mt +8 3<, y, tD J3 − 3 è!!!!!! 13 2 è!!!!!! 3+ 13 +I− 3 2 + è!!!!!!! 13 2 I− −3 Mt 3 2 + +t− è!!!!!!! 13 2 2 Mt è!!!!!! 13 −8 è!!!!!! 3+ 13 +I− 3 2 − è!!!!!!! 13 2 Mt tNE== Mentre per il secondo caso possiamo avere, ad esempio: DSolve@8lineare2, y@0D 99y → FunctionA8t<, 15 è!!!!!! 13 I− 3 2 − 0, y '@0D 1 J−78 + 39 26 è!!!! 13!!! 2 Mt + 39 I− 3 2 2<, y, tD I− + 3 2 è!!!! 13!!! 2 − è!!!! 13!!! 2 Mt Mt + 15 − è!!!!!! 13 I− 3 2 + è!!!! 13!!! 2 Mt − 26 tNE== Il punto dove specifichiamo la derivata può essere anche diverso dal primo: DSolve@8lineare2, y@0D 0, y '@4D 99y → FunctionA8t<, è!!!!!! è!!!!!! è!!!!!! J−9 − 3 13 + 9 4 13 − 3 13 6 è!!!!!! 6+2 13 +I− è!!!!!! 6+2 13 +I− è!!!!!! I3 + 13 − 3 6 Printed by Mathematica for Students 3 2 − 3 2 + è!!!!!!! 13 2 è!!!!!!! 13 2 è!!!!!! 4 13 Mt −9 4 2<, y, tD è!!!!!! 13 +9 è!!!!!!! è!!!!!! 4 13 +I− 32 − 13 2 Mt è!!!!!! Mt − 3 t − 13 t + 3 è!!!!!! + 13 I− è!!!!!! 4 13 525 ME== 3 2 + è!!!!!!! 13 2 +3 è!!!!!! 4 13 Mt +3 è!!!!!! 13 è!!!!!! t − 13 è!!!!!! 13 è!!!!!! 4 13 +I− è!!!!!! 4 13 I− 3 2 − + è!!!!!!! 13 2 è!!!!!!! 13 2 Mt 3 2 tN í Mt + − + Daniele Lupo Appunti di Mathematica Plot@Evaluate@y@tD ê. 88%%%@@1DD<, 8%%@@1DD<, 8%@@1DD<<D, 8t, −6, 6<, PlotStyle → 8Gainsboro, Firebrick, IvoryBlack<, Background → Linen, PlotRange → 8−20, 20< D 20 15 10 5 -6 -4 -2 2 4 6 -5 -10 -15 -20 Graphics Come possiamo vedere tutte e tre le soluzioni passano per l'origine, come specificato dalle nostre condizioni al contorno. Possiamo risolvere un altro tipo di equazione differenziale: problemacontorno = y ''@xD + y@xD y@xD + y @xD E^x x prgen = DSolve@problemacontorno, y, xD 99y → FunctionA8x<, C@1D Cos@xD + C@2D Sin@xD + 1 2 x HCos@xD2 + Sin@xD2 LE== Proprio non riuscite a capire perchè l'ho chiamata così questa equazione, vero? Vedete un poco adesso quello che succede. In generale, posso porre delle condizioni al contorno per avere il risultato voluto, come in questo caso: DSolve@8problemacontorno, y@0D 99y → FunctionA8x<, 1, y@1D 1 ê 2<, y, xD 1 HCos@xD + x Cos@xD2 − Cot@1D Sin@xD + 2 Csc@1D Sin@xD − Csc@1D Sin@xD + x Sin@xD2 LE== Printed by Mathematica for Students 526 Daniele Lupo Appunti di Mathematica Guardate con più attenzione il risultato, e noterete che in effetti non possiamo trovare tutte le soluzioni che vogliamo. y@0D ê. prgen 9 1 + C@1D= 2 Se imponiamo y@0D ã 1, allora per forza di cose C@1D deve essere uguale ad 1 ê 2. In questo caso nel resto possiamo giostrarci solamente con C@2D. Se però ci mettiamo a calcolare la funzione in p: y@πD ê. prgen 9 π 2 − C@1D= Notiamo come in questo caso il valore dipenda anch'esso esclusivamente da C@1D. In questo caso, fissato y@0D, non possiamo fissare con altrettanta arbitrarietà il valore della funzione in p, perchè avrei delle condizioni inconsistenti. Abbiamo la soluzione se imponiamo i valori che abbiamo trovato: DSolve@8problemacontorno, y@0D 99y → FunctionA8x<, 1 HCos@xD + 2 1, y@πD x E ^ π ê 2 − 1 ê 2<, y, xD Cos@xD2 + 2 C@2D Sin@xD + x Sin@xD2 LE== Se invece andiamo a inserire un valore diverso nel punto p otteniamo un errore, perchè non è obiettivamente possibile andare a trovare un valore di C@1D che soddisfi entrambe le condizioni: DSolve@8problemacontorno, y@0D — 1, y@πD 1<, y, xD DSolve::bvnul : For some branches of the general solution, the given boundary conditions lead to an empty solution. More… 8< non esistono: nell'ultimo caso C@2D dovrebbe essere pari ad 1 ê 2 nel punto x = 0, mentre dovrebbe In questo caso otteniamo un errore, e non viene restituita nessuna soluzione, per il semplice fatto che essere pari a 1 - ‰p ê 2 nel punto x = p, cosa che ovviamente non è possibile in quanto C@1D rappresenta una quantità costante. Naturalmente, con le opportune condizioni al contorno, possiamo trovare soluzioni al contorno di equazioni non lineari, di grado più elevato e così via, a patto comunque, come abbiamo appena visto, che le condizioni siano consistenti: Printed by Mathematica for Students 527 Daniele Lupo Appunti di Mathematica quarto = y ''''@xD + 2 y ''@xD − 2 y@xD −2 y@xD + 2 y @xD + yH4L @xD Cos@xD Cos@xD DSolve@quarto, y, xD "#########è!!!! ######### −1+ 3 x 99y → FunctionA8x<, C@3D + "#########è!!!! ######### − −1+ 3 x "########è!!! ####### C@4D + C@1D CosA 1 + 3 xE + 1 i ######## ######## "########è!!! ####### è!!! è!!! j j− "################################ 2 H−1 + 3 L H1 + 3 L Cos@xD − C@2D SinA 1 + 3 xE + 12 j k "################################ ######## ######## "########è!!! ####### è!!! è!!! 2 H−1 + 3 L H1 + 3 L Cos@xD CosA 1 + 3 xE − 2 2 y "################################ ######## ######## "########è!!! ####### è!!! è!!! z 2 H−1 + 3 L H1 + 3 L Cos@xD SinA 1 + 3 xE z zE== { Un caso particolare si ha per: DSolve@8quarto, y@0D 0, y '@0D 99y → FunctionA8x<, i j j j j j k 1 "########è!!! ####### 24 1 + 3 "#########è!!!! ######### − −1+ 3 x 3 è!!! 6 4 "########è!!! ####### 1+ 3 "#########è!!!! ######### −1+ 3 x "########è!!! ####### 4 1+ 3 72 2 "#########è!!!! ######### −1+ 3 x +2 −3, y ''@0D i j "########è!!! ####### è!!! è!!! è!!! j j j j−9 2 + 3 6 + 2 1 + 3 + 9 2 k "########è!!! ####### 1+ 3 "#########è!!!! ######### −1+ 3 x "#########è!!!! ######### −1+ 3 x "#########è!!!! ######### −1+ 3 x "########è!!! ####### 1+ 3 Printed by Mathematica for Students − "########è!!! ####### Cos@xD + 4 1 + 3 2 "#########è!!!! ######### −1+ 3 x "########è!!! ####### Cos@xD CosA 1 + 3 xE − è!!! 2 H−1 + 3 L "########è!!! ####### %%%%%%%% SinA 1 + 3 xE − 18 $%%%%%%%%%%%%%%%% è!!! %%%%%% 1+ 3 "#########è!!!! ######### −1+ 3 x 9<, y, xD 2 "#########è!!!! ######### −1+ 3 x − "########è!!! ####### CosA 1 + 3 xE − 2 è!!! 6 H−1 + 3 L "########è!!! ####### %%%%%%%%%%%%%% SinA 1 + 3 xE + 6 $%%%%%%%%%%%%%%%% è!!! 1+ 3 4 0, y '''@0D "#########è!!!! ######### −1+ 3 x "#########è!!!! ######### −1+ 3 x "########è!!! ####### SinA 1 + 3 xE − 2 yy zz "########è!!! ####### z z zE== Cos@xD SinA 1 + 3 xE z z z zz {{ 528 Daniele Lupo Appunti di Mathematica Plot@Evaluate@y@xD ê. %D, 8x, −7, 7<, PlotStyle → DarkKhaki, Background → AzureD 10 5 -6 -4 -2 2 4 6 -5 -10 Graphics Analogamente alle equazioni possiamo trattare anche i sistemi. A volte è conveniente scrivere separatamente equazioni e condizioni iniziali, e congiungerli nel comando DSolve: A=8 80, −2, 0<, 82, 1, 0<, 81, 0, 1< <; Y@t_D := 8x@tD, y@tD, z@tD<; sislin = H#1 8x @tD #2L & ∼ MapThread ∼ 8Y '@tD, A.Y@tD< −2 y@tD, y @tD valiniziali = 8x@0D 2 x@tD + y@tD, z @tD 0, y@0D 3, z@0D x@tD + z@tD< −3<; DSolve@Join@sislin, valinizialiD, 8x, y, z<, tD 3 99x → FunctionA8t<, −4 $%%%%%% 5 y → FunctionA8t<, FunctionA8t<, − 1 5 tê2 1 10 tê2 tê2 SinA è!!!!!! 15 t EE, 2 è!!!!!! è!!!!!! i 15 t 15 t y è!!!!!! j j15 CosA zE, z → E + 15 SinA Ez 2 2 k { è!!!!!! è!!!!!! i y 15 t 15 t z è!!!!!! j j45 tê2 − 15 CosA E − 15 SinA EzE== 2 2 k { Visualizziamo separatemente, prima, e con un grafico parametrico dopo: Printed by Mathematica for Students 529 Daniele Lupo Appunti di Mathematica DisplayTogetherArray@8 Plot@#, 8t, −3, 3<, PlotStyle → 8Green, Blue, BurntUmber<D, ParametricPlot3D@Join@#, 88Red, Thickness@0.01D<<D, 8t, 0, 5<, BoxRatios −> 81, 1, .5<D< & @@ 88x@tD, y@tD, z@tD< ê. %@@1DD<, Background −> Linen D 10 5 -3 -2 -1 1 2 3 -5 0 -200 -400 -600 -10 0 -20 -20 -10 0 -15 10 GraphicsArray parecchi casi reali. Un esempio è l'equazione logistica y£ HtL r I1 - ÅÅÅÅKÅÅÅÅÅ M yHtL , detta anche di Molte volte abbiamo a che fare con equazioni non lineari, invece che lineari, che sono utili in yHtL Verhlust, che altro non è se un'equazione di Riccati: eqlog = y '@tD y @tD r H1 − Hy@tD ê KLL ∗ y@tD r y@tD J1 − y@tD N K DSolve@eqlog, y, tD 99y → FunctionA8t<, r t+K C@1D −1 + K r t+K C@1D E== Quest'equazione serve per descrivere l'andamento di una popolazione. Il tasso di crescita è dato da r, mentre K rappresenta la saturazione della popolazione: avremo una crestica più o meno veloce, a seconda della popolazione iniziale, che viene definita nelle condizioni iniziali: andamento = DSolve@8eqlog, y@0D — α<, y, tD Solve::ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. More… 99y → FunctionA8t<, rt K−α+ Kα rt α E== Il warning ottenuto è dato dall'inversione dell'esponenziale, ma siccoma sappiamo essere nel nostro caso una funzione monotona, possiamo ignorarlo. In questo modo possiamo plottare il grafico per Printed by Mathematica for Students 530 Daniele Lupo Appunti di Mathematica diversi valori iniziali della popolazione. Se normalizziamo il tutto, quindi ponendo K = 1, allora possiamo disegnare diversi grafici, che rappresentano l'andamento della popolazione che varia a seconda del tasso di crescita e della popolazione iniziale. Ho dovuto lavorare un poco per ottenere un array bidimensionale, ma è sempre meglio di dover scrivere i vari plot singolarmente: DisplayTogetherArray@ 8Take@#, 3D, Take@#, −3D< & @ HDisplayTogether ê@ Table@ Plot@ Evaluate@ y@tD ê. andamento@@1DD ê. 8α → a, r → b, K → 1< D, 8t, 0, 5<, PlotRange → 80, 1<, PlotStyle → Hue@aD, PlotLabel −> "Tasso crescita r = " <> ToString@bD D, 8b, −2, 3<, 8a, 0, 1, .03< DL D Tasso crescita r = −2 Tasso crescita r = −1 Tasso crescita r = 0 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 Tasso crescita r = 1 1 0.8 0.6 0.4 0.2 Tasso crescita r = 2 1 0.8 0.6 0.4 0.2 1 2 3 4 5 1 2 3 4 5 Tasso crescita r = 3 1 0.8 0.6 0.4 0.2 1 2 3 4 5 GraphicsArray Possiamo vedere che all'aumentare oppure al diminuire del tasso di crescita, il raggiungimento della saturazione (o l'estinzione, nel caso di tasso di crescita negativo), sia più o meno rapido. Nelle nostre condizioni al contorno, possiamo anche specificare l'infinito, come nel caso seguente di un'equazione non lineare del secondo ordine: nonlinord2 = y ''@xD ê 3 y @xD 3 y@xD ^ 3 − y@xD −y@xD + y@xD3 Printed by Mathematica for Students 531 Daniele Lupo Appunti di Mathematica DSolve@nonlinord2, y, xD — Solve::ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. More… 99y → FunctionA8x<, 1 è!!! 3 $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%!!!!!! %%%%%%%% JacobiSNA è!!! è!!!!!!!!!!!!!!!! −3 + 3 3 − 2 C@1D 1 , è!!! è!!!!!!!!!!!!!!!!!!!!!! è!!! è!!!!!!!!!!!!!!!!!!!!!! − è!!! I H3 x2 + 3 x2 3 − 2 C@1D + 6 x C@2D + 2 3 x 3 − 2 C@1D C@2D + 2 è!!!!!!!!!!!!!!!!!!!!!! 3 − 9 − 6 C@1D è!!! è!!!!!!!!!!!!!!!!!!!!!! 3 C@2D2 + 3 3 − 2 C@1D C@2D2 LM, è!!!!!!!!!!!!!!!!!!!!!! E − 3 + 9 − 6 C@1D 1 è!!!!!!!!!!!!!!!!!!!!!! $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%!!!!!! %%%%%%%% 3 − 2 C@1D JacobiSNA è!!! è!!!!!!!!!!!!!!!! −3 + 3 3 − 2 C@1D 1 , è!!! è!!!!!!!!!!!!!!!!!!!!!! è!!! è!!!!!!!!!!!!!!!!!!!!!! − è!!! I H3 x2 + 3 x2 3 − 2 C@1D + 6 x C@2D + 2 3 x 3 − 2 C@1D C@2D + 2 è!!!!!!!!!!!!!!!!!!!!!! 3 − 9 − 6 C@1D è!!! è!!!!!!!!!!!!!!!!!!!!!! 3 C@2D2 + 3 3 − 2 C@1D C@2D2 LM, è!!!!!!!!!!!!!!!!!!!!!! EE=, 3 + 9 − 6 C@1D 9y → FunctionA8x<, 1 è!!! 3 $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%!!!!!! %%%%%%%% JacobiSNA è!!! è!!!!!!!!!!!!!!!! −3 + 3 3 − 2 C@1D 1 , è!!! 2 è!!!!!!!!!!!!!!!!!!!!!! è!!! è!!!!!!!!!!!!!!!!!!!!!! 2 3 − 2 C@1D + 6 x C@2D + 2 3 x 3 − 2 C@1D C@2D + è!!! I H3 x + 3 x 2 è!!!!!!!!!!!!!!!!!!!!!! 3 − 9 − 6 C@1D è!!! è!!!!!!!!!!!!!!!!!!!!!! 3 C@2D2 + 3 3 − 2 C@1D C@2D2 LM, è!!!!!!!!!!!!!!!!!!!!!! E − 3 + 9 − 6 C@1D 1 è!!!!!!!!!!!!!!!!!!!!!! $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%!!!!!! %%%%%%%% 3 − 2 C@1D JacobiSNA è!!! è!!!!!!!!!!!!!!!! −3 + 3 3 − 2 C@1D 1 , è!!! 2 è!!!!!!!!!!!!!!!!!!!!!! è!!! è!!!!!!!!!!!!!!!!!!!!!! 2 3 − 2 C@1D + 6 x C@2D + 2 3 x 3 − 2 C@1D C@2D + è!!! I H3 x + 3 x 2 è!!!!!!!!!!!!!!!!!!!!!! 3 − 9 − 6 C@1D è!!! è!!!!!!!!!!!!!!!!!!!!!! 3 C@2D2 + 3 3 − 2 C@1D C@2D2 LM, è!!!!!!!!!!!!!!!!!!!!!! EE== 3 + 9 − 6 C@1D Possiamo andare a trovarci una soluzione particolare, andando ad invocare la derivata calcolata all'infinito, il che permette di dire che la funzione deve avere un asintoto orizzontale: Printed by Mathematica for Students 532 Daniele Lupo Appunti di Mathematica DSolve@8nonlinord2, y@0D 0, y '@∞D 0<, y, xD — Solve::ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. More… — DSolve::bvlim : For some branches of the general solution, unable to compute the limit at the given points. Some of the solutions may be lost. More… — DSolve::bvlim : For some branches of the general solution, unable to compute the limit at the given points. Some of the solutions may be lost. More… 3 è!!!!!! 99y → FunctionA8x<, −TanhA$%%%%%% x2 EE=, 2 3 è!!!!!! 9y → FunctionA8x<, TanhA$%%%%%% x2 EE== 2 Nonostante i warning, che avvertono che non possiamo trovare soluzioni per tutti i rami, abbiamo ottenuto due soluzioni che soddisfano la nostra condizione al contorno Plot@8y@xD ê. %@@1DD, y@xD ê. %@@2DD, 1, −1<, 8x, −2, 2<, PlotRange → All, PlotStyle → 88Red<, 8Green<, 8Dashing@80.03<D, Blue<, 8Dashing@80.03<D, Violet<<D 1 0.5 -2 -1 1 2 -0.5 -1 Graphics Una caratteristica molto importante, spuntata dalla versione 5.0 di Mathematica, è la caratteristica di poter utilizzare come coefficienti non costanti funzioni definite a tratti, tramite Piecewise. Possiamo quindi definire delle funzioni definite a tratti, come la seguente: Printed by Mathematica for Students 533 Daniele Lupo Appunti di Mathematica tratto = PiecewiseExpand@UnitStep@xD − UnitStep@x − 2D + Max@x, Sin@xDDD Ø ≤ ≤ ≤ ≤ ≤ ∞ ≤ ≤ ≤ ≤ ≤ ± Hx ≥ 2 && x − Sin@xD ≥ 0L »» Hx < 0 && x − Sin@xD ≥ 0L 0 ≤ x < 2 && x − Sin@xD ≥ 0 Hx ≥ 2 && x − Sin@xD < 0L »» Hx < 0 && x − Sin@xD < 0L True x 1+x Sin@xD 1 + Sin@xD Come possiamovedere otteniamo una funzione definita a tratti. Plot@tratto, 8x, −4, 5<, PlotStyle → 8Red, Thickness@0.01D<, Background → GainsboroD 5 4 3 2 1 -4 -2 2 4 -1 Graphics Possiamo utilizzare quindi anche questi coefficienti. Visualizziamo il risultato generico (stavolta in forma tradizionale. Troppa Mathematica d'estate mi fa male...): eqtratto = y '@xD Max@x, x ^ 2D == y@xD + x; DSolve@eqtratto, y, xD êê TraditionalForm ::y Ø FunctionB8x<, Ø - ÅÅ1xÅÅ ≤ ≤ ≤ ≤ ≤ logHxL ∞ ≤ ≤ ≤ ≤ ≤ 1 ‰±1- ÅÅÅÅx x§0 0<x§1 True Ø - ÅÅ1xÅÅ ≤ ≤ ≤ ≤ ≤ logHxL ∞ ≤ ≤ ≤ ≤ ≤1- ÅÅ1ÅÅ c1 + ‰± x x§0 0<x§1 True 1 ijØ ≤ -EiH ÅÅÅÅx L jj≤ ≤ jj≤ ≤ logHxL jjj∞ jj≤ ≤ jj≤ ≤ EiH1L EiI ÅÅ1xÅÅ M ≤ ÅÅÅÅ - ÅÅÅÅÅÅÅÅ ÅÅÅÅÅ k± ÅÅÅÅÅÅÅÅ ‰ ‰ yz zz z 0 < x § 1 zzzzF>> zz zz True { x§0 Anche in questo caso possiamo imporre delle condizioni iniziali opportune per il nostro problema: Printed by Mathematica for Students 534 Daniele Lupo Appunti di Mathematica DSolve@8eqtratto, y@1D 0<, y, xD 99y → FunctionA8x<, Ø − 1x ≤ ≤ ≤ ≤ ∞Log@xD ≤ ≤ ≤ 1 ≤ ±1− x x≤0 0<x≤1 True Ø −ExpIntegralEi@ 1x D i ≤ j ≤ j ≤ j j j≤ ∞ Log@xD j j j ≤ j≤ ≤ ≤ ExpIntegralEi@1D − ExpIntegralEi@ 1x D k± x≤0 y z z z 0 < x ≤ 1z z E== z z z z True { Possiamo plottare adesso la funzione, assieme alla sua derivata, per far vedere meglio dove sono presenti le discontinuità nella derivata, che diventano cuspidi nella funzione Plot@ Evaluate@8y@xD, D@y@xD, xD< ê. %D, 8x, −3, 3<, PlotStyle → 88Blue, Thickness@0.01D<, 8Red, Thickness@0.005D, Dashing@80.02<D<<, Background → Linen D — Power::infy : Infinite expression 1 encountered. 0 More… — Power::infy : Infinite expression 1 encountered. 0 More… 1 -3 -2 -1 1 2 3 -1 -2 Graphics Come abbiamo visto, abbiamo il plottaggio della soluzione, a parte i warning dovuti alla singolarità nel punto 0. Con questo, credo di aver concluso la panoramica sulla risoluzione simbolica delle equazioni differenziali. Credo che adesso possiamo passare alla risoluzione numerica, che dal punto di vista teorico è più difficile, dato che fino ad adesso ha fatto tutto DSolve... Printed by Mathematica for Students 535 Daniele Lupo Appunti di Mathematica Risoluzione numerica Quando non possiamo trovare la soluzione esatta di un'equazione differenziale, magari perchè non si riesce a trovarla, l'equazione è troppo difficile da risolvere etc etc etc, allora si rende necessaria la risoluzione numerica. ü NDSolve Il comando di Mathematica dedito a questo fine è NDSolve. Il suo utilizzo per equazioni abbastanza elementari è abbastanza semplice. Basta infatti specificare l'equazione differenziale (oppure il sistema), come abbiamo fatto per DSolve, con un paio di modifiche: prima di tutto, siccome la soluzione è numerica e non simbolica, NDSolve non utilizza i parametri indeterminati C@nD, per cui abbiamo sempre bisogno di specificare le condizioni al contorno del nostro problema, per trovare una soluzione esclusivamente numerica. Inoltre, è necessario anche specificare il range entro cui vogliamo calcolare la funzione. Infatti Mathematica crea una funzione interpolata all'interno dell'intervallo definito: al di fuori dell'intervallo dobbiamo utilizzare delle estrapolazioni, che rendono impreciso il risultato all'aumentare della distanza dal punto. In teoria basterebbe estendere il range ad un intervallo più grande di quello che ci interessa, per evitarci problemi. Tuttavia questo può essere eccessivamente dispendioso in termini di tempo di calcolo della nostra funzione interpolata, oltre che eccessivo in termini di memoria e di pesantezza di gestione della funzione. Occorre sempre studiarsi prima il problema per cercare l'intervallo che ci interessa, senza sprecare tempo prezioso. Se questo appare abbastanza superfluo per problemi semplici, può diventare determinante quando i calcoli si fanno pesanti, come per esempio si ha nelle equazioi di Navier-Stokes. Inoltre, a seconda della natura delle equazioni, sono adatti certi algoritmi di risoluzione rispetto ad altri. Trovare il giusto algoritmo è sempre stato un problema; NDSolve effettua ogni volta un pre-processing della nostra equazione, cercando di scoprire quale algoritmo meglio si presta alla risoluzione del nostro problema. Comunque in casi particolarmente ostici, oppure quando semplicemente vogliamo maggior controllo su quello che facciamo, possiamo selezionare manualmente le varie opzioni, per rendere più specifico il risultato e forzare NDSolve a lavorare come vogliamo noi. Data l'estrema generalizzazione della risoluzione numerica, possiamo trovare soluzioni per tutti i tipi di equazioni differenziali, tempo di calcolo permettendo; quindi lavora tranquillamente anche con tutti i tipi di equazioni differenziali finora elencate nella risoluzione simbolica: ODE, PDE, DAE. Vediamo un esempio per capire il funzionamento base di NDSolve: Printed by Mathematica for Students 536 Daniele Lupo Appunti di Mathematica esnum = NDSolve@8x '@tD Sin@x@tD + Cos@t ^ 2DD, x@0D 88x → InterpolatingFunction@880., 9.42478<<, <>D<< 3<, x, 8t, 0, 3 Pi<D Possiamo notare come abbiamo utilizzato per forza anche le condizioni iniziali, e come assieme alla variabile indipendente sia stato specificato l'intervallo di integrazione dell'equazione differenziale. Mathematica naturalmante restituisce un errore se mettiamo un numero insufficiente di condizioni iniziali, oppure se non specifichiamo l'intervallo: NDSolve@x '@tD — NDSolve::ndnco : The number of constraints H0L Hinitial conditionsL is not equal to the total differential order of the system H1L. More… NDSolve@x @tD NDSolve@8x '@tD — Sin@x@tD + Cos@t ^ 2DD, x, 8t, 0, 3 Pi<D Sin@Cos@t2 D + x@tDD, x, 8t, 0, 3 π<D Sin@x@tD + Cos@t ^ 2DD, x@0D 3<, x, tD NDSolve::ndlim : Range specification t is not of the form 8x, xend< or 8x, xmin, xmax<. NDSolve@8x @tD Sin@Cos@t2 D + x@tDD, x@0D More… 3<, x, tD Una volta ottenuta la soluzione, possiamo utilizzarla come abbiamo fatto per il risultato di DSolve. Inatti anche in questo caso il risultato è formito come regola di sostituzione da utilizzare dove serve: x@4D ê. esnum 83.13899< Possiamo utilizzarla ovunque possiamo mettere una regola, ergo anche nel plottaggio, naturalmente: Printed by Mathematica for Students 537 Daniele Lupo Appunti di Mathematica Plot@x@tD ê. esnum, 8t, 0, 3 Pi<, PlotRange → All, PlotStyle → Green, Background → LinenD 3.2 3.1 2 4 6 8 2.9 2.8 2.7 2.6 Graphics Come potete notare la funzione viene smorzata, fino a stabilizzarsi attorno ad un valore limite. Tuttavia, se proviamo a calcolare la funzione in punti che cadono al di fuori dell'intervallo di definizione, viene utilizzata l'estrapolazione, che porta ad un grande errore, che cresce, in generale, man mano che ci si allontana dagli estremi dell'intervallo. Mathematica segna questo problema con un warning, da prendere abbastanza seriamente, anche se esegue comunque l'estrapolazione: Printed by Mathematica for Students 538 Daniele Lupo Appunti di Mathematica Plot@x@tD ê. esnum, 8t, 0, 4 Pi<, PlotRange → All, PlotStyle → Green, Background → LinenD — — — — InterpolatingFunction::dmval : Input value 89.47541< lies outside the range of data in the interpolating function. Extrapolation will be used. More… InterpolatingFunction::dmval : Input value 89.99568< lies outside the range of data in the interpolating function. Extrapolation will be used. More… InterpolatingFunction::dmval : Input value 89.72325< lies outside the range of data in the interpolating function. Extrapolation will be used. More… General::stop : Further output of InterpolatingFunction::dmval will be suppressed during this calculation. More… 1500 1250 1000 750 500 250 2 4 6 8 10 12 Graphics Come potete vedere, al di fuori dell'intervallo praticamente la funzione interpolata non serve a niente; questo è vero in particolar modo per le funzioni oscillanti come questa. Analogamente a questo, possiamo anche risolvere sistemi di equazioni differenziali, come in questo esempio: sol = NDSolve@8 x ''@tD y '@tD ê x@tD, y '@tD x@tD ê Hy@tD + 30 x@tD ^ Sin@tDL, x@0D 1, x '@0D .1, y@0D 4 <, 8x, y<, 8t, 0, 40<D 88x → InterpolatingFunction@880., 40.<<, <>D, y → InterpolatingFunction@880., 40.<<, <>D<< Printed by Mathematica for Students 539 Daniele Lupo Appunti di Mathematica Plot@Evaluate@8x@tD, y@tD< ê. %D, 8t, 0, 40<, PlotRange → All, PlotPoints → 200, PlotStyle → 8Red, Blue<D; 30 25 20 15 10 5 10 20 30 40 Il sistema può contenere un numero arbitrario di incognite, ovviamente. L'esempio serve solo per farvi capire come funziona il tutto. Possiamo risolvere numericamente anche delle PDE, come nel caso seguente, dove consideriamo l'evoluzione di un'equazione differenziale che dipende da due coordinate spaziali e da una temporale; questo esempio può prendere un po' di tempo in più rispetto alle equazioni banali. Potremmo dire che è un esempio quasi serio... sol = NDSolve@8 D@u@t, x, yD, t, tD D@u@t, x, yD, x, xD + D@u@t, x, yD, y, yD − Sin@u@t, x, yDD, u@0, x, yD Exp@−Hx2 + y2 LD, Derivative@1, 0, 0D@uD@0, x, yD 0, u@t, −5, yD u@t, 5, yD Sin@t y ê 2D ê 2, u@t, x, −5D u@t, x, 5D 0 <, u, 8t, 0, 25<, 8x, −5, 5<, 8y, −5, 5<D — — NDSolve::ibcinc : Warning: Boundary and initial conditions are inconsistent. More… NDSolve::eerr : Warning: Scaled local spatial error estimate of 508.4402959177667` at t = 25.` in the direction of independent variable x is much greater than prescribed error tolerance. Grid spacing with 91 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or you may want to specify a smaller grid spacing using the MaxStepSize or MinPoints options. More… 88u → InterpolatingFunction@880., 25.<, 8−5., 5.<, 8−5., 5.<<, <>D<< Adesso possiamo vedere il risultato tramite un'animazione che mostra l'evoluzione temporale della soluzione: Printed by Mathematica for Students 540 Daniele Lupo Appunti di Mathematica animazione = MoviePlot3D@First@u@t, x, yD ê. solD, 8x, −5, 5<, 8y, −5, 5<, 8t, 0, 25<, PlotRange → 8−1, 1<, PlotPoints → 40, Frames → 110, Mesh → False, LightSources → 8882, 2, 2<, White<, 88−2, 2, 2<, Blue<<D 1 0.5 4 0 -0.5 -1 2 0 -4 -2 -2 0 2 -4 4 Se procedete con l'animazione, vedrete nel vostro schermo una simpatica cosuccia... Quello che fa NDSolve è considerare ina sequenza di passi nella variabile indipendente t, e cercare il metodo migliore per la risoluzione dell'equazione differenziale. Inoltre, adatta automaticamente lo step size della variabile, in maniera da rispettare la precisione voluta da PrecisionGoal, che rappreseta l'errore relativo, ed AccuracyGoal, che rappresenta quello assoluto. Nel nostro caso abbiamo una discontinuità ai vertici della regione d'integrazione, perchè in due lati la condizione al contorno è 0, mentre negli altri 2 è Sin@tD: allora ai vertici ho contemporaneamente 0 e Sin@tD, e Mathematica ci avvisa di questo, non riuscendo fra l'altro ad imporre un giusto step size per l'algoritmo, in quanto le discontinuità non sono ben digerite. In ogni caso otteniamo un risultato che giudico personalmente più che soddisfacente, almeno per poter fare un esempio pratico. In altri casi si arriva alla produzione di artefatti, che possono falsare il risultato. Possiamo anche aumentare la precisione della soluzione numerica, quando ci serve: Printed by Mathematica for Students 541 Daniele Lupo Appunti di Mathematica NDSolve@8 x ''@tD x@tD, x@0D 1, x '@0D 0 <, x, 8t, 3<, AccuracyGoal → 20, PrecisionGoal → 20, WorkingPrecision → 25 D 88x → InterpolatingFunction@880, 3.000000000000000000000000<<, <>D<< Una volta ottenuta la soluzione con la precisione voluta, possiamo calcolare i risultati: x@2D ê. % 83.76219569108363660919707< Come possiamo vedere, abbiamo ottenuto una precisione superiore a quella di macchina, che avremmo ottenuto sfruttando le impostazioni di default. Moltissimo c'è da dire ancora sulla risoluzione numerica delle equazioni differenziali, ma purtroppo il tempo mi manca.... Se avete dubbi leggete il manuale, sperimentate e non esistate a cercare nuove vie e nuovi comandi!!!!! Printed by Mathematica for Students 542