% % 31-10-94 Tode Prima impaginazione fine % % Daniele Finocchiaro & Gianluca Marcoccia % % INTEGRA.library % % articolo di presentazione - IPISA '94 % 14.10.1994 Dany % % ********************************************************************** % L'uso e l'eventuale modifica di questo testo è riservato agli autori % del software e agli organizzatori di IPISA94. Tutti gli altri sono % caldamente invitati a rinunciarvi. % ********************************************************************** % % \nuovoarticolo{Daniele Finocchiaro e Gianluca Marcoccia}{Integra.library}{% Daniele Finocchiaro \\ Via P.~Scuderi 1 \\ I-95015 Linguaglossa CT \\ Tel.:~(095)643014 \\ \medskip \\ Gianluca Marcoccia \\ Via Garibaldi 78 \\ I-56100 Pisa PI } {Internet:~marcocci@cli.di.unipi.it~(Marcoccia)}{italian}{% \selectlanguage{english}The primary purpose of the ¶{DATA} statement is to give names to constants; instead of referring to $\pi$ as $3.141592653589793$ at every appearance, the variable ¶{PI} can be given that value with a ¶{DATA} statement and used instead of the longer form of the constant. This also simplifies modifying the program, should the value of $\pi$ change.}{FORTRAN manual for Xerox Computers} %Idea e implementazione: Daniele Finocchiaro %Supporto libraries Amiga: Gianluca Marcoccia La libreria di integrazione che presentiamo cerca di venire incontro alle esigenze degli utilizzatori più diversi, che abbiano bisogno di calcolare integrali, offrendo loro routine che costituiscono lo stato dell'arte delle conoscenze attuali. Per la risoluzione di questo problema esistono su grosse piattaforme programmi celebri quali CADRE o QUADPACK, mentre su \amiga{} non si è ancora visto nulla di simile. Il calcolo di una qualche forma di integrale è un problema che si incontra di frequente nelle applicazioni matematiche e ingegneristiche. Ad esempio, per effettuare misurazioni sismografiche, si utilizzano degli strumenti (accelerometri) che forniscono l'accelerazione $a(t)$ in ogni istante di tempo $t$; per risalire da questa allo spostamento effettivo del terreno bisogna effetture due integrazioni successive dei dati. Da un punto di vista matematico, l'integrale di una funzione $f(x)$ è definito come l'area compresa tra l'asse delle ascisse e la funzione stessa (per una definizione accurata si rimanda ovviamente a un qualunque libro di calcolo). Si parla di integrale definito tra i punti $a$ e $b$ quando quest'area venga limitata da due rette verticali passanti per $a$ e per $b$. L'integrale definito di $f(x)$ nell'intervallo $[a,b]$ si indica allora con $$ \int_a^b f(x)dx . $$ Il teorema fondamentale del calcolo integrale afferma che, se $g(x)$ è una funzione tale che $g'(x)=f(x)$, allora per l'integrale precedente vale $$ \int_a^b f(x)dx = g(b) - g(a) . $$ Poiché esistono vari metodi per risalire dall'espressione esplicita di una funzione $f(x)$ a quella della sua funzione {\em primitiva\/} $g(x)$, è spesso possibile calcolare un integrale facendo ricorso alle primitive ed ottenendo così un risultato matematicamente esatto. Tuttavia questo teorema non risolve per niente il problema del calcolo di integrali definiti. Tanto per cominciare, della funzione $f(x)$ non si ha, nelle applicazioni reali, una espressione nota, ma se ne conosce solamente il valore $f(x_i)$ assunto in alcuni punti $x_i$. In secondo luogo, anche quando si conosce esplicitamente $f(x)$, spesso non è semplice trovare un'espressione analitica per una sua primitiva. Infine, per alcune funzioni, come la celebre $e^{-x^2}$, è stata addirittura dimostrata la {\em non\/} esistenza di una primitiva esprimibile analiticamente. Si pone dunque il problema di calcolare integrali definiti senza lavorare sulle espressioni delle funzioni ma soltanto su numeri (i valori $x_i$ e $f(x_i)$) e cercando di ottenere solamente dei numeri (il valore dell'integrale). Questo è il campo dell'integrazione numerica. L'integrazione numerica introduce innanzitutto due concetti fondamentali: {\it efficienza} e {\it precisione}. L'obiettivo che ci si pone è quello di trovare il valore numerico di un integrale definito. Questa operazione deve essere eseguita velocemente, in quanto è presumibilmente solo un passo nella risoluzione di un problema più grande. Inoltre, escludendo in generale di poter dare un risultato esatto, bisogna raggiungere la massima accuratezza possibile. È evidente tuttavia che maggiore è la precisione richiesta, maggiore è il tempo di calcolo necessario per ottenerla. Molti dei metodi oggi utilizzati per calcolare integrali attraverso l'uso del computer sono noti da diversi secoli: essi venivano applicati a mano per calcolare gli integrali delle funzioni più ``difficili'', ad esempio per tabulare i valori della ``funzione di errore'' $\int e^{-x^2}$. Newton scriveva a Leibniz addirittura nel 1676 delle formule oggi ancora usate e note sotto il nome di Newton-Cotes. Esse applicavano già quella che è un'idea comune a quasi tutti i metodi. Supponiamo di conoscere una funzione solo attraverso i valori $f(x_i)$ che essa assume in determinati punti $x_i$. Come possiamo calcolarne l'integrale? Per metterci l'anima in pace, la teoria ci garantisce che esistono {\em sempre\/} funzioni tanto cattive da farci commettere un errore arbitrariamente grande. Esistono cioè funzioni $g(x)$ tali che $g(x_i)=f(x_i)$ (e quindi per noi indistinguibili da $f(x)$) ma tali che $\int g(x)dx$ sia un numero qualunque. Bisogna quindi stabilire un certo insieme di funzioni ``plausibili'' e richiedere che almeno su queste funzioni il nostro metodo sia abbastanza accurato. Per funzioni base sono stati universalmente scelti i polinomi di grado via via crescente: $1, x, x^2,\ldots$. Gran parte dei metodi si basano sul supporre che per i punti assegnati $(x_i,f(x_i))$ passi un polinomio, ed è questo che viene integrato effettivamente. Integrare un polinomio è cosa relativamente semplice, così come trovare un polinomio $p(x)$ tale che $p(x_i)=f(x_i)$ (ovvero un polinomio che {\em interpola\/} $f(x)$ nei punti $x_i$). Quindi il gioco è fatto. Si può inoltre mostrare che gran parte del lavoro necessario per implementare questo metodo può essere fatto a priori, calcolando dei coefficienti $w_i$ (detti pesi) che identificano il metodo. A quel punto, effettuare l'integrazione consiste semplicemente in un prodotto scalare: $$ I = \sum_i w_i f(x_i) . $$ Se si hanno a disposizione $k$ punti $x_1,x_2,\ldots,x_k$ dove la funzione è nota, si può interpolarla con un polinomio di grado $k-1$. Queste formule sono dette interpolatorie, e se i punti $x_i$ sono equidistanziati sono dette formule di Newton-Cotes. Il primo problema grosso che si pone è il fatto che i polinomi di grado alto sono dei cattivi interpolanti, nel senso che oscillano violentemente e, pur passando per i punti $(x_i,f(x_i))$, sono vistosamente poco ``plausibili''. È per questo motivo che in genere i punti su cui la funzione è nota vengono divisi in gruppetti, per ciascuno dei quali si fa passare un polinomio di grado basso, che si comporta ``bene'' (si usano in genere delle spline). Questo approccio dà vita alle formule composte, come la celebre formula di Cavalieri-Simpson, che, intuitivamente, fa passare una parabola per ogni gruppo di tre punti consecutivi. Quando i punti $x_i$ sono assegnati dall'esterno, e in un numero fissato, non si può fare molto di meglio. Ci sono però dei casi in cui la $f(x)$ può essere calcolata in punti arbitrari, e quindi chi definisce il metodo di integrazione è libero di piaz\-za\-re gli $x_i$ dove è più conveniente. Si possono allora trovare formule molto più precise. Esiste ad e\-sem\-pio una vasta let\-te\-ra\-tu\-ra sulle formule gaussiane, in cui i nodi $x_i$ sono le radici di opportuni polinomi ortogonali. Esistono dei problemi in cui la funzione $f(x)$ può essere valutata in un numero qualunque di punti. Questo accade ad e\-sem\-pio quando sia fornita un'espressione esplicita di $f(x)$, o quando essa sia la risposta data da uno strumento fisico a un certo input. Si sono quindi sviluppati metodi che decidono autonomamente in quanti e quali punti valutare la funzione, anziché ricevere i valori $(x_i,f(x_i))$ in ingresso. L'utilizzo massiccio di computer negli ultimi anni ha favorito lo studio di questo tipo di algoritmi nel contesto di quella che è nota come {\em quadratura automatica}. In questo ambito, quello che costa di più è valutare la funzione, e quindi gli algoritmi vengono confrontati in base a quante valutazioni richiedono per poter dare un risultato con una precisione prestabilita. Si cerca allora di studiare algoritmi che riescano a capire quanto una funzione sia ``cattiva'' e quindi quanto tempo deve essere dedicato per raggiungere una precisione ragionevole. La strategia comune a tutti è quella del {\em divide et impera\/}: se si ritiene (in base a certi calcoli) che l'approssimazione non sia buona, si divide l'intervallo da integrare in più sottointervalli, e su ciascuno di questi si calcola un integrale accurato, sommando i risultati. Ovviamente la strategia viene poi applicata ricorsivamente a tutti i sottointervalli. Negli intervalli più piccoli il calcolo dell'integrale viene infine fatto attraverso formule note, come quelle gaussiane. In questo modo si può vedere che l'integrazione di una funzione semplice, come un polinomio, richiede poche valutazioni. Una funzione ``difficile'' (come le forze impulsive che si trovano in fisica) richiede invece molte più valutazioni. All'interno della quadratura automatica si distinguono due grandi filoni, a seconda che la tecnica utilizzata sia {\em adattiva\/} oppure no. Le tecniche adattive sono più sensibili alla particolare funzione che si sta integrando, in quanto riescono a capire {\em dove\/} la funzione è più cattiva, all'interno dell'intervallo in esame, e concentrare lì il maggior numero di valutazioni. Le tecniche non adattive invece valutano la funzione uniformemente in tutto l'intervallo, anche in zone dove la funzione è più ``buona''. Un gruppo di ricerca nato a Pisa intorno al professor Francesco Romani studia da tempo algoritmi e formule per la quadratura automatica. Alcuni dei risultati maggiori sono quelli esposti negli articoli citati in bibliografia. Sono state trovate delle formule (dette RMS dalle i\-ni\-zia\-li di Ricorsive, Monotòne e Stabili) che si prestano meglio di altre al loro utilizzo all'interno di una tecnica adattiva. L'algoritmo ottenuto inserendo queste formule nelle funzioni della libreria QUADPACK, battezzato ``Algorithm 691'', è attualmente il migliore integratore nel suo genere. Dalle prove empiriche infatti risulta che esso riesce ad ottenere, a parità di valutazioni della funzione, una precisione molto migliore degli altri algoritmi noti. Il punto di forza consiste nel fatto che le formule RMS consentono di applicare una strategia adattiva riutilizzando tutti i valori di $f(x)$ già noti, laddove l'utilizzo di altre formule costringe l'algoritmo a valutare $f(x)$ in nuovi punti. Le funzioni che formano la libreria che presentiamo cercano di fornire buoni algoritmi per la risoluzione dei diversi tipi di problemi di cui abbiamo parlato. Infatti la funzione può essere ({\em a\/}) nota in $k$ punti equidistanti; ({\em b\/}) nota in $k$ punti qualsiasi; ({\em c\/}) calcolabile ovunque attraverso una chiamata di funzione; ({\em d\/}) nota esplicitamente. Nei casi ({\em a\/}) e ({\em b\/}) si utilizzano le formule di Newton-Cotes o di Cavalieri-Simpson (i due casi sono distinti per ragioni di efficienza, in quanto il caso di punti equidistanti è di trattazione molto più agile). Questi casi si applicano quando si siano effettuate delle registazioni attraverso apparecchiature fisiche, e se ne deve effettuare uno studio a posteriori. Nei casi ({\em c\/}) e ({\em d\/}) vengono invece utilizzate le formule RMS all'interno di una tecnica adattiva. Il caso ({\em c\/}) è utile quando il computer è direttamente interfacciato a un dispositivo che può calcolare la $f(x)$, o quando questo valore può essere ottenuto attraverso elaborati calcoli, o infine quando esso possa essere chiesto da tastiera. Infine, il caso ({\em d\/}) si applica quando della funzione abbiamo una espressione esplicita, sotto forma di stringa. In questa prima versione della libreria si è cercato di favorire la qualità, mettendola davanti all'efficienza. Nella speranza che la libreria trovi degli utilizzatori in\-te\-res\-sa\-ti, è intenzione degli autori affinare ulteriormente le tecniche utilizzate, rendendole al contempo più efficienti. È importante, a questo proposito, il feedback degli utilizzatori. Infatti ogni algoritmo va in un certo senso ``sintonizzato'' sulle particolari classi di funzioni che deve integrare, così come una macchina di formula 1 va regolata in base alla pista su cui deve correre. \begin{thebibliography}{0} \bibitem[1]{} P.~J.~Davis, P.~Rabinowitz. \newblock {\em Methods of Numerical Integration}. \newblock Academic Press, New York, 1975. \bibitem[2]{} R.~Bevilacqua, D.~Bini, M.~Capovani, O.~Menchi. \newblock {\em Metodi Numerici}. \newblock Zanichelli, Bologna, 1992. \bibitem[3]{} P.~Favati, G.~Lotti, F.~Romani. \newblock {\em Interpolatory Integration Formulas for Optimal Composition}. \newblock {\em ACM Transactions on Mathematical Software}. \newblock Vol.~17, N.~2, Giugno 1991, pp.~207--217. \bibitem[4]{} P.~Favati, G.~Lotti, F.~Romani. {\em Algorithm 691. Improving QUADPACK Integration Routines}. \newblock {\em ACM Transactions on Mathematical Software}. \newblock Vol.~17, N.~2, Giugno 1991, pp.~218--232. \end{thebibliography} \finearticolo