Sign in to follow this  
DarwinNE

Maxima, il nipotino del grande MACSYMA

Recommended Posts

Macsyma è stato il primo programma di manipolazione algebrica sviluppato al MIT fra la fine degli anni sessanta e gli inizi degli anni ottanta. Programmi tipo Mathematica devono molto a questi pioneristici sforzi. Ne parlo qui, perché so che ci sono molti scienziati qui nel forum e spero di fare cosa gradita:

http://en.wikipedia.org/wiki/Macsyma

Esiste una versione OpenSource rilasciata secondo la licenza GPL e chiamata Maxima ed adattata a Common Lisp.

http://en.wikipedia.org/wiki/Maxima

Ieri sera, sono riuscito ad installare senza nessuna difficoltà Maxima sotto MacOSX direttamente dai sorgenti. Ad ogni modo, esiste anche un port sotto Darwinports, quindi la vita è ancora semplificata:

http://darwinports.opendarwin.org/ports/?by=name&substr=maxima

Come tanti programmi OpenSource puri e duri (come Gnuplot e Octave), si usa da terminale, a me la cosa non dispiace affatto perché sono abituato ad usare strumenti del genere e quasi tutti i miei file di lavoro (testo puro) sono cosi' compatibili con qualunque ambiente e sistema operativo.

Non l'ho usato ancora molto, ma ecco un esempio di quello che Maxima puo' fare:

Ecco il fattoriale di 300:

(%i37) 300!;(%o37) 30605751221644063603537046129726862938858880417357699941677674125947653\317671686746551529142247757334993914788870172636886426390775900315422684292790\697455984122547693027195460400801221577625217685425596535690350678872526432189\626429936520457644883038890975394348962543605322598077652127082243763944912012\867867536830571229368194364995646049816645022771650018517654646934011222603472\972406633325858350687015016979416885035375213755491028912640715715483028228493\795263658014523523315693648223343679925459409527682060806223281238738388081704\9600000000000000000000000000000000000000000000000000000000000000000000000000(%i38) 

Mi ha calcolato senza fare una piega il fattoriale di 1000 e di 10000, solo che prendono troppo spazio ed ho pensato che fosse inutile ingombrare inutilmente il database del forum ;)

Derivata di sin(x)^5*exp(x) rispetto ad x:

(%i38) diff(sin(x)^5*%e^x,x);                        x    5          x           4(%o38)                %e  sin (x) + 5 %e  cos(x) sin (x)

Integrale indefinito di sin(x)^5 rispetto ad x:

(%i39) integrate(sin(x)^5,x);                             5           3                          cos (x)   2 cos (x)(%o39)                  - ------- + --------- - cos(x)                             5          3(%i40) 

Non so a voi, ma a me già piace!!!

Share this post


Link to post
Share on other sites

Aggiungo qualche commento e qualche esempio. Ho installato Maxima sull'iMac G5 procedendo dai sorgenti con la procedura seguente:

-ho caricato ed installato clisp tramite DarwinPorts (sudo port install clisp, ovviamente DarwinPorts deve essere installato)

-ho scaricato i sorgenti di Maxima

-ho fatto l'installazione classica:

.configure

make all

sudo make install

-credo che ci sia bisogno di Gnuplot, sulla mia macchina, l'avevo già installato, se no, c'è il pacchetto con DarwinPorts. Ho anche installato il terminale 'Aquaterm', che posso usare da Gnuplot.

-per far girare Maxima, basta a questo punto battere "maxima" da terminale

Sul MacBook, l'ho invece installato direttamente con DarwinPorts, ma la gestione della linea di comando non è perfetta: sull'iMac G5 con i tasti cursore posso infatti navigare nello storico dei comandi digitati in precedenza, mentre sul MacBook mi compaiono solo dei codici strani.

Una volta lanciato, Maxima si presenta come segue:

[email protected]]$ maximaMaxima 5.11.0 http://maxima.sourceforge.netUsing Lisp CLISP 2.38 (2006-01-24)Distributed under the GNU Public License. See the file COPYING.Dedicated to the memory of William Schelter.This is a development version of Maxima. The function bug_report()provides bug reporting information.(%i1) 
A questo punto, si possono introdurre i comandi. Iniziamo con il disegno di una funzione semplice, tipo sin(x)/x:

(%i1) f: sin(x)/x;                                    sin(x)(%o1)                               ------                                      x(%i2) plot2d(f,[x,-3*%pi,3*%pi], [gnuplot_term, aqua]);Output file "/Users/davidebucci/maxplot.aqua".(%o2) (%i3)

Ed ecco il risultato:

Posted Image

Naturalmente, le ultime versioni di Gnuplot permettono di fare grafici piuttosto complessi, ma bisogna conoscere questo (splendido) programma e si tratta di un impresa un po' ostica all'inizio.

Ecco un esempio di come specificare dei comandi Gnuplot all'interno di Maxima per ottenere un grafico tridimensionale:

(%i19) f: sin(sqrt(x^2+y^2));                                        2    2(%o19)                        sin(sqrt(y  + x ))(%i20) plot3d(f,[x,-3*%pi,3*%pi],[y,-3*%pi,3*%pi],[gnuplot_term,aqua], [gnuplot_pm3d, true], \[gnuplot_preamble,"set pm3d at s;unset surface;unset key"]);Output file "/Users/davidebucci/maxplot.aqua".(%o20) (%i21) 

Ed ecco il risultato. Almeno nella mia installazione di Gnuplot, fintantoché non mi decido di configurarlo come si deve, il programma parte senza avere un terminale grafico definito. Bisogna quindi specificarlo ogni volta con l'opzione [gnuplot_term,aqua].

Posted Image

Altrimenti, si puo' benissimo usare una rappresentazione testuale che, seppure poco attrattiva, conserva un certo suo fascino vintage:

(%i21) f: sin(x)/x;                                    sin(x)(%o21)                              ------                                      x(%i22) plot2d(f,[x,-3*%pi,3*%pi], [gnuplot_term, dumb]);    1 ++-----+------+------+------+----$$$$$----+------+------+------+-----++      +      +      +      +      +   $$ + $$   +      +    sin(x)/x $$$$$$ +      |                              $$     $$                              |  0.8 ++                            $$       $$                            ++      |                             $         $                             |  0.6 ++                           $           $                           ++      |                           $$           $$                           |      |                           $             $                           |  0.4 ++                         $$             $$                         ++      |                          $               $                          |  0.2 ++                        $                 $                        ++      |     $$$$$              $$                 $$              $$$$$     |      |  $$$$   $$$           $$                   $$           $$$   $$$$  |    0 ++$$         $$        $$                     $$        $$         $$++      |             $$$     $$                       $$     $$$             | -0.2 ++              $$$  $$                         $$  $$$              ++      |                 $$$                             $$$                 |      +      +      +      +      +      +      +      +      +      +      + -0.4 ++-----+------+------+------+------+------+------+------+------+-----++     -10    -8     -6     -4     -2      0      2      4      6      8      10Output file "/Users/davidebucci/maxplot.dumb".(%o22) (%i23) 

E dopotutto, il comportamento della funzione è interpretabile...

Quando avro' un po' di tempo continuero', se la cosa interessa.

Share this post


Link to post
Share on other sites

continua continua, mi interessa molto! Anche se prima devo risolvere i problemi che ho con darwinport e fink (entrambi, forse perché li ho trascurati da tempo, hanno smesso di funzionare correttamente..)

se riesci metti anche tutti i file di configurazione per definire il terminale grafico.

Presto io tenterò l'esperimento di installare Koctave su Mac (interfaccia grafica di Octave per KDE), installando prima KDE e poi il porting tramite Portage.. solo che devo prima liberare un po' di spazio su HD.

Share this post


Link to post
Share on other sites

se riesci metti anche tutti i file di configurazione per definire il terminale grafico.

In realtà, non ho configurato niente: specifico ogni volta a Gnuplot di utilizzare il terminale Aqua.

Ecco qualche esempio di calcolo con espressioni trigonometriche. Io le ho sempre trovate molto rognose perché non mi sono mai ricordato a memoria le varie formule (sarà per quello che ho poi trovato molto comoda la notazione esponenziale).

Iniziamo con il definire l'espressione di partenza: sin(x+y)^3, una possibile esigenza puo' esser quella di scrivere tutto utilizzando dei termini in sin() e cos(), con argomento semplice composto solo da una variabile. Il comando da utilizzare è trigexpand:

(%i25) ex: sin(x+y)^3;                                     3(%o25)                            sin (y + x)(%i26) trigexpand(ex);                                                      3(%o26)                 (cos(x) sin(y) + sin(x) cos(y))

Un'altra esigenza puo' essere quella di convertire le potenze di funzioni trigonometriche utilizzando le note proprietà, di modo da ottenere una somma di termini che contengono solo un sin() o cos(), pero' con argomenti più complessi. La funzione da utilizzare è trigreduce:

(%i27) trigreduce(ex);                         3 sin(y + x) - sin(3 (y + x))(%o27)                   -----------------------------                                       4(%i28) 

Una cosa che capita alle volte è avere delle funzioni tipo tangente, cotangente ed altro che si desidera riportare in una espressione compatta utilizzando solo seni e coseni. La funzione è trigsimp. Ecco un esempio in cui si vedono anche gli effetti dei comandi visti in precedenza.

(%i29) ex1: sin(x)*tan(2*x)^2;                                         2(%o29)                         sin(x) tan (2 x)(%i30) trigexpand(ex1);                                           2                               4 sin(x) tan (x)(%o30)                         ----------------                                        2    2                                (1 - tan (x))(%i31) trigreduce(ex1);                                         2(%o31)                         sin(x) tan (2 x)(%i32) trigsimp(ex1);                                         2                               sin(x) sin (2 x)(%o32)                         ----------------                                     2                                  cos (2 x)(%i33) 

Quanto è noioso calcolare gli integrali definiti o indefiniti di funzioni trigonometriche... Integrare un seno alla quinta non presenta nessuna difficoltà concettuale, a parte lo svolgere i calcoli. Che sono noiosi. Adatti quindi ad essere eseguiti in maniera automatica:

(%i33) f: sin(x)^5;                                       5(%o33)                              sin (x)(%i34) integrate(f);Wrong number of arguments to integrate -- an error.  To debug this try debugmode(true);(%i35) integrate(f,x);                             5           3                          cos (x)   2 cos (x)(%o35)                  - ------- + --------- - cos(x)                             5          3(%i36) trigreduce(%);       - cos(5 x) - 5 cos(3 x) - 10 cos(x)   cos(3 x) + 3 cos(x)(%o36) ----------------------------------- + ------------------- - cos(x)                       80                             6(%i37) 
Si noti come nell'input 33 ho dimenticato il punto e virgola. Maxima non continua fintantoché non l'inserisco, per esempio nella riga successiva.

Nell'input 34 ho invece dimenticato la variabile secondo cui derivare.

Il comando alla linea 36 come visto sopra indica invece di eliminare le potenze ed ottenere delle somme in cui le funzioni trigonometriche compaiono solo con potenza 1. Il simbolo di percento % indica di utilizzare nel processo l'ultimo risultato ottenuto.

Ovviamente, si puo' derivare, con il comando diff(). Riprendo qui l'espressione iniziale, il seno alla quinta:

(%i37) diff(f);                                        4(%o37)                      5 cos(x) sin (x) del(x)
In questo caso, non avendo specificato la variabile secondo cui differenziare, Maxima resta sul generale e la forma del(x) indica appunto il differenziale di x.

Volendo dire di differenziare rispetto ad x, non c'è più quest'ambiguità:

(%i48) diff(sin(x)^5,x);                                           4(%o48)                         5 cos(x) sin (x)(%i49) 

Share this post


Link to post
Share on other sites

Dato che questa sera a quanto pare non riesco a combinare nulla di buono per il mio lavoro, continuo questa esplorazione di Maxima, tanto per invogliare qualcuno ad usare questo bello strumento.

Vediamo qualche manipolazione algebrica su polinomi.

Definiamo un polinomio per esempio di grado due, di cui possiamo estrarre facilmente le due radici (reali):

(%i54) p: 3*x^2-5*x+2;                                   2(%o54)                          3 x  - 5 x + 2(%i55) sol: solve(p=0, x);                                            2(%o55)                          [x = 1, x = -]                                            3(%i56) 

Se le radici non sono reali, questo non costituisce un gran problema per Maxima, dato che naturalmente si puo' lavorare in campo complesso:

(%i65) p: 3*x^2+2*x+2;                                   2(%o65)                          3 x  + 2 x + 2(%i66) sol: solve(p=0, x);                         sqrt(5) %i + 1      sqrt(5) %i - 1(%o66)            [x = - --------------, x = --------------]                               3                   3(%i67) solve(p=0, x), numer;`rat' replaced 4.47213595499958 by 24476/5473 = 4.472135940069432`rat' replaced 4.47213595499958 by 24476/5473 = 4.472135940069432(%o67) [x = - 6.090504902856446E-5 (12238 %i + 5473),                                     x = 6.090504902856446E-5 (12238 %i - 5473)](%i68) expand(%);(%o68) [x = - .7453559900115719 %i - .3333333333333333,                                   x = .7453559900115719 %i - .3333333333333333](%i69) 

Come ci aspettavamo, dato che il discriminante è negativo, si ottengono le due radici complesse e coniugate. Maxima fornisce un'espressione in forma compatta resultante dalle manipolazioni analitiche. Delle volte, fa comodo avere i risultati numerici e quindi ho utilizzato nella linea %i67 il modificatore numer, che indica al programma di presentare i risultati in forma numerica, che poi ho fatto scrivere in forma estesa utilizzando expand;

Fino a qui, niente di trascendente. Diventa più interessante calcolare tutte le radici di un polinomio di grado alto. Non è un problema che nel caso generale è risolvibile analiticamente, quindi anche Maxima propone la funzione allroots(), che permette di fare i calcoli numericamente (ottenendo se si ha fortuna circa 16 cifre significative, se il problema non è mal condizionato).

Ecco un esempio di grado 12 (che ho preso da un mio libro di elettronica, per verificare il risultato):

%i71) p: 1+2*x^3+2*x^6+2*x^9+x^12;                          12      9      6      3(%o71)                   x   + 2 x  + 2 x  + 2 x  + 1(%i72) allroots(p=0);(%o72) [x = .8660254151728854 %i + .5000000152932945, x = .5000000152932945 - .8660254151728854 %i, x = .4999999999999991 %i - .8660254037844395, x = - .4999999999999991 %i - .8660254037844395, x = .9999999999999994 %i + 5.866392267360735E-16, x = 5.866392267360735E-16 - .9999999999999994 %i, x = .8660253923959913 %i + .4999999847067038, x = .4999999847067038 - .8660253923959913 %i, x = .5000000000000003 %i + .8660254037844387, x = .8660254037844387 - .5000000000000003 %i, x = - .9999999683898616, x = - 1.000000031610135](%i73) 
Confrontando i risultati, si vede che la precisione che si ottiene è di circa 7-8 cifre (è banale osservare che x=-1 è soluzione). Problemi come questi non sono banali ed è sempre buona norma darsi gli strumenti per controllare la loro efficacia, per esempio provando altri metodi numerici e confrontando i risultati.

I polinomi sono spesso utilizzati in molti campi e le loro proprietà dipendono in molti contesti dalla posizione delle radici sul piano complesso. Purtroppo, come abbiamo visto, fare calcoli con polinomi di grado alto è difficile e procedendo solo per via numerica puo' capitare che le approssimazioni si facciano sentire in modo preoccupante. Spesso pero', quello che serve è solamente osservare la parte reale delle radici, per vedere se è positiva, nulla o negativa, per delle ragioni legate alla stabilità di un sistema. Per un polinomio con soluzioni ideali, basta la nota regoletta dei segni di Cartesio. Quando pero' si prendono in considerazione soluzioni in campo complesso questa regola non funziona più e bisogna ricorrere ad una prova matematica più complessa denominata prova di Hurwitz.

La prova di Hurwitz è un algoritmo puramente algebrico che consente di osservare se un polinomio ha tutte le radici con parte reale non positiva (polinomio hurwitziano in senso lato), oppure negativa (polinomio strettamente hurwitziano). Vediamo qui come effettuare la prova di Hurwitz utilizzando Maxima. Ecco il polinomio di partenza, preso sempre dal mio libro di teoria delle reti elettriche (uso p per indicare la pulsazione complessa, occhio che non deve essere stata definita prima):

(%i95) kill(p);(%o95)                               done(%i96) M(p):=56*p^5+7*p^4+161*p^3+18*p^2+98*p+8;                           5      4        3       2(%o96)         M(p) := 56 p  + 7 p  + 161 p  + 18 p  + 98 p + 8(%i97) 
Si noti che ho definito M(p) come una funzione dipendente dalla variabile p, grazie all'operatore :=

La parte pari di un polinomio è il polinomio formato da tutti i termini con potenza di p pari mentre la parte dispari è ovviamente quella che ha tutti i termini con potenza di p dispari.

Bisogna procedere per divisioni successive dividendo la parte pari per la parte dispari del polinomio o viceversa, a seconda di quella che ha il grado più elevato. Ovviamente, si ottiene un termine di solito di grado 1, più un resto.

Con Maxima è semplice estrarre parte pari e dispari di un polinomio, basta ricordare alcune proprietà ed in particolare le seguenti:

(%i105) M(p):=56*p^5+7*p^4+161*p^3+18*p^2+98*p+8;                           5      4        3       2(%o105)        M(p) := 56 p  + 7 p  + 161 p  + 18 p  + 98 p + 8(%i106) Pa(p):=1/2*(M(p)+M(-p));                                   1(%o106)                   Pa(p) := - (M(p) + M(- p))                                   2(%i107) Di(p):=1/2*(M(p)-M(-p));                                   1(%o107)                   Di(p) := - (M(p) - M(- p))                                   2(%i108) expand(Pa(p));                                  4       2(%o108)                        7 p  + 18 p  + 8(%i109) expand(Di(p));                                 5        3(%o109)                      56 p  + 161 p  + 98 p(%i110) 

Bene. Incominciamo con le derivate successive. Dividiamo Di(p) per Pa(p), perché la prima ha grado più alto. La funzione da utilizzare è divide, che vuole il numeratore, il denominatore e la variabile principale del polinomio (dato che possono essercene diverse):

(%i110) divide(Di(p),Pa(p),p);                                        3(%o110)                       [8 p, 17 p  + 34 p](%i111) 

Il risultato è scritto sotto forma di lista, ovvero una serie di elementi racchiusi fra parentesi quadre. Il primo elemento è il risultato, il secondo è il resto. La prova di Hurwitz consiste a continuare le divisioni successive, questa volta fra la parte pari ed il resto (che è dispari). Si accede al resto prendendo il secondo elemento della lista e facendo la divisione:

(%i128) Di1(p):=second(%o110);(%o128)                     Di1(p) := second(%o110)(%i129) Di1(p);                                     3(%o129)                          17 p  + 34 p(%i130) divide(Pa(p), Di1(p),p);                                 7 p     2(%o130)                         [---, 4 p  + 8]                                 17(%i131)

Si continua cosi' fino a che non rimane più nulla da dividere:

(%i131) Pa1(p):=second(%o130);(%o131)                     Pa1(p) := second(%o130)(%i132) Pa1(p);                                      2(%o132)                            4 p  + 8(%i133) divide(Di1(p),Pa1(p),p);                                    17 p(%o133)                            [----, 0]                                     4(%i134) 

Si sono ottenuti quindi come quozienti di primo grado 8p, 7/17p e 17/4p. Sono tutti positivi, ma sono solo tre invece di cinque (tali quanti il grado). Il polinomio dato è hurwitziano in senso lato.

Il testo da cui ho preso gli esempi è "Sintesi dei circuiti passivi", di Claudio Beccari edito da CLUT.

Share this post


Link to post
Share on other sites

Grazie Jeby! Nei prossimi giorni saro' un po' impegnato, ma spero di continuare qui l'esplorazione di questo strumento appassionante. Cerco nel mio piccolo di dare un contributo al livello tecnico delle discussioni del forum.

Tanto per restare in tema, l'operatore bfloat() permette di definire tutti gli elementi dell'espressione che è stata fornita come dei big float, ovvero dei numeri in virgola mobile anche molto lunghi. Il numero di cifre è stabilito dal valore della variabile fpprec.

Ecco come stampare le prime cinquantamila cifre di pi greco (ci vuole un po' di tempo):

(%i143) fpprec:50000;(%o143)                              50000(%i144) bfloat(%pi);(%o144) 3.1415926535897...Qui taglio, senno' riempio io da solo il database del forum...712307568610450045489703600795698276263923441071465848957802414081584052295\369374997106655948944592462866199635563506526234053394391421112718106910522900\24657423488b0

Per la cronaca, ho controllato con questo sito, che riporta i primi 4 milioni di decimali ed il risultato è corretto nell'ambito della precisione richiesta:

http://zenwerx.com/pi.php

Alla prossima puntata. Nel frattempo, rimango in ascolto delle vostre domande, esperienze, suggerimenti, correzioni. Cerchiamo di far vedere come il Mac sia uno strumento sempre migliore per la ricerca scientifica!

Share this post


Link to post
Share on other sites

Visto che stanotte il sonno è ormai una cosa piuttosto remota, continuo le esplorazioni.

In questa puntata, vedremo brevemente come fare qualche calcolo con le matrici e come utilizzare TeX per mostrare i risultati.

Ecco uno dei modi per definire una matrice: il comando entermatrix, che vuole come argomento le dimensioni:

(%i154) a: entermatrix(4,4);Is the matrix  1. Diagonal  2. Symmetric  3. Antisymmetric  4. GeneralAnswer 1, 2, 3 or 4 : 4;Row 1 Column 1: k^2-1;Row 1 Column 2: 3;Row 1 Column 3: 1;Row 1 Column 4: 1;Row 2 Column 1: 5;Row 2 Column 2: 3;Row 2 Column 3: 2;Row 2 Column 4: 1;Row 3 Column 1: k;Row 3 Column 2: 2;Row 3 Column 3: 1;Row 3 Column 4: -1;Row 4 Column 1: k;Row 4 Column 2: -1;Row 4 Column 3: 2;Row 4 Column 4: 0;Matrix entered.                            [  2                  ]                            [ k  - 1   3   1   1  ]                            [                     ](%o154)                     [   5      3   2   1  ]                            [                     ]                            [   k      2   1  - 1 ]                            [                     ]                            [   k     - 1  2   0  ](%i155) 
Dapprima viene chiesto il tipo di matrice e poi si procede elemento per elemento. Da notare che non si è limitati ad input numerici, ma si possono anche utilizzare delle variabili, come per esempio k.

Rappresentiamo con il TeX il risultato:

(%i156) tex(a);$$\pmatrix{k^2-1&3&1&1\cr 5&3&2&1\cr k&2&1&-1\cr k&-1&2&0\cr }$$
Il codice racchiuso fra il doppio dollaro $$ è l'espressione scritta nel formato TeX. Per chi non conoscesse TeX e LaTeX, si veda qui:

http://it.wikipedia.org/wiki/TeX

http://it.wikipedia.org/wiki/LaTeX

Io ho una versione di teTeX installata a cui ho aggiunto alcuni pacchetti. Per comodità, alle volte utilizzo un softwarino molto carino, LaTeXiT:

http://ktd.club.fr/programmation/latexit_en.php

Con cui si puo' fare molto comodamente copia/incolla in presentazioni, testi o salvare facilmente l'immagine.

Il solo problema che ho incontrato è che il codice generato da Maxima segue le vecchie convenzioni TeX puro ed ho per esempio dovuto disattivare il pacchetto amsmath modificando il preambolo standard usato da LaTeXiT, escludendo questo pacchetto.

Ad ogni modo, ecco il risultato.

Posted Image

Già: TeX e LaTeX sono considerati attualmente i migliori sistemi al mondo per la composizione automatica delle equazioni. In ogni istante si puo' ottenere il codice TeX ed ottenere l'immagine dell'espressione senza fare nessuno sforzo.

Proseguiamo. Ecco il calcolo della matrice inversa con il comando invert; puo' risultare comodo tenere a fattore comune il determinante con l'opzione detout:

(%i160) invert(a);                 [                 13                  ]                 [ ----------------------------------- ]                 [      2                              ]                 [ 13 (k  - 1) + 3 k + 3 (k - 10) - 30 ]                 [                                     ]                 [               k - 10                ]                 [ ----------------------------------- ]                 [      2                              ]                 [ 13 (k  - 1) + 3 k + 3 (k - 10) - 30 ](%o160)  Col 1 = [                                     ]                 [              - 6 k - 5              ]                 [ ----------------------------------- ]                 [      2                              ]                 [ 13 (k  - 1) + 3 k + 3 (k - 10) - 30 ]                 [                                     ]                 [              9 k - 25               ]                 [ ----------------------------------- ]                 [      2                              ]                 [ 13 (k  - 1) + 3 k + 3 (k - 10) - 30 ]         [                   12                  ]         [ - ----------------------------------- ]         [        2                              ]         [   13 (k  - 1) + 3 k + 3 (k - 10) - 30 ]         [                                       ]         [                  2                    ]         [              2 (k  - 1)               ]         [  -----------------------------------  ]         [       2                               ]         [  13 (k  - 1) + 3 k + 3 (k - 10) - 30  ] Col 2 = [                                       ]         [              2                        ]         [             k  + 6 k - 1              ]         [  -----------------------------------  ]         [       2                               ]         [  13 (k  - 1) + 3 k + 3 (k - 10) - 30  ]         [                                       ]         [               2                       ]         [           5 (k  - 1) - 6 k            ]         [  -----------------------------------  ]         [       2                               ]         [  13 (k  - 1) + 3 k + 3 (k - 10) - 30  ]         [                   1                   ]         [  -----------------------------------  ]         [       2                               ]         [  13 (k  - 1) + 3 k + 3 (k - 10) - 30  ]         [                                       ]         [              2                        ]         [          2 (k  - 1) + k - 10          ]         [  -----------------------------------  ]         [       2                               ]         [  13 (k  - 1) + 3 k + 3 (k - 10) - 30  ] Col 3 = [                                       ]         [                 2                     ]         [                k  - 6                 ]         [  -----------------------------------  ]         [       2                               ]         [  13 (k  - 1) + 3 k + 3 (k - 10) - 30  ]         [                                       ]         [       2                               ]         [ - 8 (k  - 1) + 3 k + 3 (10 - 2 k) + 5 ]         [ ------------------------------------- ]         [       2                               ]         [  13 (k  - 1) + 3 k + 3 (k - 10) - 30  ]         [                  5                  ]         [ ----------------------------------- ]         [      2                              ]         [ 13 (k  - 1) + 3 k + 3 (k - 10) - 30 ]         [                                     ]         [              2                      ]         [        - 3 (k  - 1) - k + 10        ]         [ ----------------------------------- ]         [      2                              ]         [ 13 (k  - 1) + 3 k + 3 (k - 10) - 30 ] Col 4 = [                                     ]         [     2                               ]         [ 5 (k  - 1) + 3 k + 3 (- k - 5) - 10 ]         [ ----------------------------------- ]         [      2                              ]         [ 13 (k  - 1) + 3 k + 3 (k - 10) - 30 ]         [                                     ]         [       2                             ]         [    - k  - 3 k - 3 (5 - 2 k) + 11    ]         [ ----------------------------------- ]         [      2                              ]         [ 13 (k  - 1) + 3 k + 3 (k - 10) - 30 ](%i161) 

Si noti come il risultato venga fornito colonna per colonna.

Facciamo la prova che la matrice per la sua inversa dia la matrice identità:

(%i168) a . b;        [  2                  ]        [ k  - 1   3   1   1  ]        [                     ](%o168) [   5      3   2   1  ] . (matrix([13, - 12, 1, 5],         [                     ]        [   k      2   1  - 1 ]        [                     ]        [   k     - 1  2   0  ]             2           2                      2[k - 10, 2 (k  - 1), 2 (k  - 1) + k - 10, - 3 (k  - 1) - k + 10],              2             2          2[- 6 k - 5, k  + 6 k - 1, k  - 6, 5 (k  - 1) + 3 k + 3 (- k - 5) - 10],                2                   2[9 k - 25, 5 (k  - 1) - 6 k, - 8 (k  - 1) + 3 k + 3 (10 - 2 k) + 5,    2                                  2- k  - 3 k - 3 (5 - 2 k) + 11])/(13 (k  - 1) + 3 k + 3 (k - 10) - 30))(%i169) c:expand(a . b);                 [          2                                             ]                 [      13 k                6 k                 73        ]                 [ ---------------- + ---------------- - ---------------- ]                 [     2                  2                  2            ]                 [ 13 k  + 6 k - 73   13 k  + 6 k - 73   13 k  + 6 k - 73 ](%o169)  Col 1 = [                                                        ]                 [                           0                            ]                 [                                                        ]                 [                           0                            ]                 [                                                        ]                 [                           0                            ]         [                           0                            ]         [                                                        ]         [          2                                             ]         [      13 k                6 k                 73        ]         [ ---------------- + ---------------- - ---------------- ] Col 2 = [     2                  2                  2            ]         [ 13 k  + 6 k - 73   13 k  + 6 k - 73   13 k  + 6 k - 73 ]         [                                                        ]         [                           0                            ]         [                                                        ]         [                           0                            ]         [                           0                            ]         [                                                        ]         [                           0                            ]         [                                                        ]         [          2                                             ] Col 3 = [      13 k                6 k                 73        ]         [ ---------------- + ---------------- - ---------------- ]         [     2                  2                  2            ]         [ 13 k  + 6 k - 73   13 k  + 6 k - 73   13 k  + 6 k - 73 ]         [                                                        ]         [                           0                            ]         [                           0                            ]         [                                                        ]         [                           0                            ]         [                                                        ]         [                           0                            ] Col 4 = [                                                        ]         [          2                                             ]         [      13 k                6 k                 73        ]         [ ---------------- + ---------------- - ---------------- ]         [     2                  2                  2            ]         [ 13 k  + 6 k - 73   13 k  + 6 k - 73   13 k  + 6 k - 73 ](%i170) 
L'operatore prodotto fra matrici è il punto. Non si ottiene la matrice identità (o almeno cosi' non sembra) perché i termini non sono fattorizzati, problema a cui si rimedia con il comando factor:

(%i174) factor(c);                                [ 1  0  0  0 ]                                [            ]                                [ 0  1  0  0 ](%o174)                         [            ]                                [ 0  0  1  0 ]                                [            ]                                [ 0  0  0  1 ](%i175) 

Naturalmente, si puo' benissimo calcolare autovalori (eigenvalues in inglese) o autovettori (eigenvectors) con i comandi eigenvectors(a) ed eigenvalues(a), ma non riporto l'output perché è molto lungo.

Per finire, ecco un esempio di un'altra espressione ottenuta grazie a TeX, un integrale questa volta:

(%i187) i:'integrate((x^2+2*x+1)/(x-2),x);                               /  2                               [ x  + 2 x + 1(%o187)                        I ------------ dx                               ]    x - 2                               /(%i188) o:integrate((x^2+2*x+1)/(x-2),x);                                            2                                           x  + 8 x(%o188)                     9 log(x - 2) + --------                                              2(%i189) tex(i=o);$$\int {{{x^2+2\,x+1}\over{x-2}}}{\;dx}=9\,\log \left(x-2\right)+{{x^2 +8\,x}\over{2}}$$(%o189)                              false(%i190) 

Qui ho usato il comando integrate due volte: la prima per calcolare la primitiva e la seconda per farmi scrivere l'espressione. Si noti che ho usato l'operatore ' (apice) che impedisce a Maxima di svolgere l'integrazione, ma la lascia non elaborata, di modo da ottenere il codice TeX corrispondente. Ecco il risultato:

Posted Image

(ah... l'uguaglianza è vera, ma nel caso generale ci vuole una costante arbitraria nella primitiva).

Alla prossima...

P.S. Nell'ultima immagine ottenuta con TeX, vedo un piccolo bug. Secondo le norme internazionali, la 'd' del differenziale non deve essere in corsivo, ma in tondo. Comunque, si tratta di un errore che ho trovato anche in diverse pubblicazioni. Basta sostituire a mano la 'd' nell'espressione TeX con {\rm d} per risolvere il problema.

Share this post


Link to post
Share on other sites

Continuiamo con qualche esempio di analisi.

Abbiamo già visto che Maxima non si fa problemi a calcolare primitive e derivate:

(%i194) f:sin(x)/(atan(x)+x^2);                                    sin(x)(%o194)                          ------------                                            2                                 atan(x) + x(%i195) g:diff(f,x);                                       1                                    (------ + 2 x) sin(x)                                      2                        cos(x)       x  + 1(%o195)              ------------ - ---------------------                                2                  2 2                     atan(x) + x       (atan(x) + x )(%i196) 'integrate(g,x)=integrate(g,x);                             1                          (------ + 2 x) sin(x)        /                   2        [     cos(x)       x  + 1(%o196) I (------------ - ---------------------) dx =         ]             2                  2 2        /  atan(x) + x       (atan(x) + x )                              cos(x) sin(2 x) - sin(x) cos(2 x) + sin(x)                        -------------------------------------------------------                                        2     2                      2     2                        (2 atan(x) + 2 x ) sin (x) + (2 atan(x) + 2 x ) cos (x)(%i197) 'integrate(g,x)=trigreduce(integrate(g,x));                               1                            (------ + 2 x) sin(x)          /                   2          [     cos(x)       x  + 1                        sin(x)(%o197)   I (------------ - ---------------------) dx = ------------          ]             2                  2 2                     2          /  atan(x) + x       (atan(x) + x )           atan(x) + x(%i198) tex(%);$$\int {{{\cos x}\over{\arctan x+x^2}}-{{\left({{1}\over{x^2+1}}+2\,x \right)\,\sin x}\over{\left(\arctan x+x^2\right)^2}}}{\;dx}={{\sin x }\over{\arctan x+x^2}}$$(%o198)                              false(%i199) 

Ho estratto il codice TeX alla linea %o198, che produce come risultato:

Posted ImageNaturalmente, si puo' facilmente calcolare anche dei limiti con il comando limit; vediamone uno facile: il limite per x che tende a zero di k*sin(x)/x.

(%i202) f:k*sin(x)/x;                                   k sin(x)(%o202)                            --------                                      x(%i203) limit(f,x,0);(%o203)                                k(%i204) 

Vediamo ancora qualche esempio, per esempio con la funzione tangente:

(%i3) limit(tan(x),x,%pi/2);(%o3)                                 und(%i4) limit(tan(x),x,%pi/2,plus);(%o4)                                minf(%i5) limit(tan(x),x,%pi/2,minus);(%o5)                                 inf

Nel primo caso, si è calcolato il limite della tangente per x che tende verso pi/2. Maxima risponde correttamente che tale limite non esiste (und corrisponde a undefined). Si sono quindi calcolati i limiti destro e sinistro allo stesso punto singolare, ottenendo meno infinito (minf), oppure più infinito (inf).

Il limite è svolto utilizzando alcune sostituzioni che alla fine portano ad una libreria di limiti notevoli, oppure confrontando gli ordini di infinito/infinitesimo (il manuale di Maxima fornisce come riferimento una tesi di dottorato molto interessante: http://www.lcs.mit.edu/publications/specpub.php?id=660 La cosa divertente è che il relatore è stato Seymour Papert, il papà del LOGO).

Esiste alternativamente le funzione tlimit() che calcola il limite utilizzando uno sviluppo di Taylor.

Attenzione pero' che in alcuni casi se la funzione ha un comportamento patologico, si puo' mettere in crisi Maxima. Ecco un esempio:

(%i58) p:(log(1+%e^(1/x))*sin(x^7))/((1/(1-x^2))^(sin(x)^2)-1-x^4);                                 7        1/x                            sin(x ) log(%e    + 1)(%o58)                     ------------------------                                  1           4                           --------------- - x  - 1                                      2                                 2 sin (x)                           (1 - x )(%i59) limit(factor(p(x)),x,0,minus);^CMaxima encountered a Lisp error: EXT:GC: User breakAutomatically continuing.To reenable the Lisp debugger set *debugger-hook* to nil.(%i60) limit(factor(p(x)),x,0,plus);                           2             2                        sin (x) log(1 - x )      7        1/x                      %e                    sin(x ) log(%e    + 1)(%o60)      - limit   --------------------------------------------              x -> 0+                 2                 2                         4       2 sin (x)         2 sin (x)                        x  (1 - x )        + (1 - x )        - 1(%i61) tlimit(factor(p(x)),x,0,plus);(%o61)                                 0(%i62) tlimit(factor(p(x)),x,0,minus);(%o62)                                 0(%i63) 

Quando ho lanciato il comando limit(factor(p(x)),x,0,minus);, Maxima si è messo a lavorare intensivamente per un po', fino a che mi sono stufato e l'ho interrotto con un CTRL-C.

Il comando tlimit sembra incontrare meno problemi e fornisce zero quasi subito, ma il risultato E' SBAGLIATO!!!

Il limite per x che tende a zero più della funzione che ho fornito è infatti 6. Lo si dimostra sviluppando in serie di Taylor un po' più in profondità rispetto a quanto non si faccia di solito. Ad ogni modo, la funzione ha un comportamento un parecchio ostico attorno all'origine, come si vede dal suo grafico:

(%i80) plot2d(p(x), [x, -0.1,.1], [gnuplot_term, aqua]); 

Posted Image

Una cosa comoda in questi casi, ma sempre da usare con attenzione, è il comando taylor, che consente appunto di sviluppare secondo Taylor in un intorno di un punto singolare fino ad un certo limite.

Il manuale ci dice:

taylor (expr, x, a, n) expands the expression expr in a truncated Taylor or Laurent series in the variable x around the point a, containing terms through (x - a)^n.

Ecco un esempio applicato al nostro caso:

(%i81) p:(log(1+%e^(1/x))*sin(x^7))/((1/(1-x^2))^(sin(x)^2)-1-x^4);                                 7        1/x                            sin(x ) log(%e    + 1)(%o81)                     ------------------------                                  1           4                           --------------- - x  - 1                                      2                                 2 sin (x)                           (1 - x )(%i82) taylor(p(x),x,0,1);(%o82)/T/                          0 + . . .(%i83) taylor(p(x),x,0,2);(%o83)/T/                          0 + . . .(%i84) taylor(p(x),x,0,3);(%o84)/T/                          0 + . . .(%i85) taylor(p(x),x,0,4);(%o85)/T/                          0 + . . .(%i86) taylor(p(x),x,0,5);(%o86)/T/                          0 + . . .(%i87) taylor(p(x),x,0,6);                                   6(%o87)/T/                         x  + . . .(%i88) taylor(p(x),x,0,7);                   2          4            6              128 x    51209 x    2021759 x(%o88)/T/ 6 - ------ + -------- - ---------- + . . .                5        525         5250               3          5            7          128 x    51209 x    2021759 x             - 1/x + (6 x - ------ + -------- - ---------- + . . .) %e            5        525         5250                3          5            7            64 x    51209 x    2021759 x              - 1/x 2 + (- 3 x + ----- - -------- + ---------- + . . .) (%e     )              5       1050       10500               3          5            7          128 x    51209 x    2021759 x              - 1/x 3 + (2 x - ------ + -------- - ---------- + . . .) (%e     )            15       1575       15750                3          5            7      3 x   32 x    51209 x    2021759 x              - 1/x 4 + (- --- + ----- - -------- + ---------- + . . .) (%e     )       2      5       2100       21000               3          5            7    6 x   128 x    51209 x    2021759 x              - 1/x 5 + (--- - ------ + -------- - ---------- + . . .) (%e     )     5      25       2625       26250              3          5            7          64 x    51209 x    2021759 x              - 1/x 6 + (- x + ----- - -------- + ---------- + . . .) (%e     )           15       3150       31500               3          5            7    6 x   128 x    51209 x    2021759 x              - 1/x 7 + (--- - ------ + -------- - ---------- + . . .) (%e     )  + . . .     7      35       3675       36750(%i89) 

Si ottiene un risultato utilizzabile soltanto quando abbiamo forzato il calcolo del polinomio di Taylor almeno fino al settimo grado. Il motivo per cui vengono riportati anche dei termini esponenziali che non sono sviluppati credo sia dovuto al fatto che e^(1/x) è decisamente un oggetto da trattare con molto rispetto nell'intorno di zero e quindi il programma non tenta di svilupparlo.

Si noti che quest'espressione è uguale a 6 per x reale che tende verso zero positivo.

Ho fatto questi commenti per fare vedere come anche usando strumenti automatici, bisogna sempre fare molta attenzione ed avere ben presente la sostanza del problema prima ancora di introdurlo. Parafrasando qualcuno che conosco: per usare un programma di calcolo simbolico, bisogna saper fare i calcoli meglio di lui.

Ovviamente, rimango sempre aperto a commenti, richieste e correzioni eventuali.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
Sign in to follow this