APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2...

37
APPUNTI DEL CORSO DI LABORATORIO DI CALCOLO AVANZATO Metodi Numerici per le Equazioni Differenziali Ordinarie MARCO LIMONGI Istituto Nazionale di Astrofisica Osservatorio Astronomico di Roma

Transcript of APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2...

Page 1: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

  

APPUNTI DEL CORSO DI LABORATORIO DI CALCOLO 

AVANZATO Metodi Numerici per le Equazioni Differenziali 

Ordinarie    

MARCO LIMONGIIstituto Nazionale di Astrofisica – Osservatorio Astronomico di Roma

 

 

 

 

 

 

 

 

 

 

Page 2: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

1  

1. EQUAZIONI DIFFERENZIALI ORDINARIE 

I problemi fisici che andremo ad affrontare in questo corso implicano la soluzione di di  equazioni o  sistemi di  equazioni differenziali ordinarie  (ODE). Nel  caso  in  cui  il grado  dell’ODE  sia maggiore  di  uno,  questa  si  puo’  in  generale  ricondurre  ad  un sistema di ODE del primo ordine. Per esempio l’ODE del secondo ordine: 

 

 

 

si puo’ riscrivere come l’insieme delle due equazioni del primo ordine 

 

 

 

 

dove  z(x)  e’  una  nuova  variabile.    Percio’  il  generico  problema  che  implica  la soluzione di una ODE e’ in generale ricondotto allo studio di un sistema di N ODE del primo ordine della forma: 

 

, , … .                   1, … ,  

dove le funzioni  , , … .  sono note. 

Un  problema  che  implica  la  soluzione  di  una  ODE  non  e’  completamente 

 

determinato  fino  a  quando  non  vengano  definite  le  condizioni  al  contorno.  Le condizioni al contorno sono condizioni algebriche sui valori delle funzioni   in punti specifici  del  dominio  di  integrazione.  Tipicamente  e’  la  natura  delle  condizioni  al 

Page 3: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

2  

contorno che determina  il metodo numerico di  integrazione piu’ consono al tipo di problema. Le condizioni al contorno dividono i problemi in due grandi categorie. 

Problemi  ai Valori  Iniziali.  In  questo  genere  di  problemi  le  condizioni  al  contorno sono date ad un estremo del dominio di  integrazione. Tipicamente  le   sono date nel punto iniziale   del dominio di integrazione e la soluzione e’ rappresentata dalla determinazione dei valori delle funzioni   nel punto finale del dominio   . 

Problemi ai Valori al Contorno.  In questo tipo di problemi  le condizioni al contorno sono date agli estremi del dominio di integrazione. Tipicamente alcune sono fornite nel punto iniziale  , mente le restanti altre nel punto finale  . 

 

1.1. PROBLEMI AI VALORI INIZIALI.  

1.1.1. Metodo di Eulero Esplicito (o “in avanti”): 

Il metodo  di  Eulero  rappresenta  l’idea  di  base  di  tutti  i metodi  numerici  per  la risoluzione dei problemi ai valori iniziali per ODE: 

,  

Il metodo di Eulero esplicito si basa sull’approssimazione della derivata prima della 

12

funzione   mediante la formula di Taylor (al secondo ordine) applicata al punto  : 

 

la quantita’    e’  il  resto della  formula di  Taylor  con   un punto 

ell’intorno di  . D

2

opportuno n efinendo   si ha: 

 

da cui si ricava il valore della derivata: 

sostituendo questa espressione nell’ODE   ,   per   si ha: 

Page 4: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

3  

2,  

Trascurando ora  il  termine del  resto di Taylor       non abbiamo piu’  i valori 

della soluzione esatta, ma otterremo i valori della soluzione approssimata.  

  e   si ha: 

Chiamando    la  soluzione  approssimata  e ponendo        e   

,  

Da cui la formula di Eulero Esplicito: 

,  

Questa  formula  avanza  la  soluzione  da    a  .  Il  metodo  si  dice  “esplicito” poiche’ il calcolo della funzione in   si basa su quantita’ calcolate in   che sono 

valonza e’ data da  ,  

Figura 1: Nel metodo di Eulero 

tutte note (esplicite).  

Da un punto di vista geometrico il  re    e’ approssimato utilizzando il valore della retta la cui pende

 

si approssima, ad ogni passo,  il valore della soluzione con il valore 

che assume la retta tangente per lo stesso nodo infatti la formula è

)( ixy

),( ii yxfh ix 1 iiy +y=+  

 

 

Page 5: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

4  

.1.2. Metodo di Eulero Implicito (o “all’indietro”): 

me fatto in precedenza, ma nel  viceversa:

2

 

1

Applichiamo la formula di Taylor non al punto     copunto di arrivo  , cioe’ scambiando    con    e  

da cui 

Risolvendo per la derivata prima si ha: 

Sostituendo nella ODE di partenza, quindi, si ha: 

2,  

Al solito, trascurando il resto avremo una soluzione approssimata. Facendo uso della notazione definita sopra, avremo: 

,  

 

La differenza rispetto al metodo esplicito e’ che la   non e’ piu’ valutata in    dove ra nota, ma al contrario  in   dove e’ ancora  incognita, poiche’ e’  incognita  la 

,di dipe

 

e. Poiche’  in generale  la   dipendera’ da     a secondo membro comparira’  

l’incognita  . Percio’  il calcolo      nde  implicitamente da    Per questo il metodo si dice “implicito”. Se la dipendenza di   da   e’ lineare, la valutaz e di

. Se di e  nu

ri di 

stesso. 

ione di   si ricondurra’ alla risoluzion  una equazione di primo grado e sara’ triviale  invece la dipendenza   da   sara’ non lineare la valutazione di 

 sara’ piu’  laboriosa e richiedera’  l’applicazione di un m todo merico per  la ricerca degli ze una funzione (es. metodo di Newton‐Raphson). 

Page 6: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

5  

1.1.3

Partiamo dall’ODE 

,  

e integriamo ambo i membri nell’intervallo  , : 

 

. Metodo di Crank‐Nicolson: 

,  

,  

do l’integrale con il metodo dei trapezi ed utilizzando la solita notazione per la soluzione approssimata si ha: 

, ,2

da cui: 

Risolven

 

2

cioe’: 

, ,  

 

La stessa formula si poteva ottenere prendendo la media aritmetica delle formule di Eulero implicito ed esplicito: 

,  

 e dividendo per due, otteniamo: 

2

,  

sommando le due equazioni e

, ,  

Anche  in  questo  caso  il metodo  e’  implicito  poiche’  richiede  la  conoscenza  della funzione nel punto di arrivo che e’ ancora incognita. 

Page 7: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

6  

 

1.1.4. Accuratezza. 

a alcune definizioni preliminari, estendibili, con opprtune modifiche, anche a ad altri metodi. 

Errore  globale  di  discretizzazione:  Si  dice  errore  globale  di  discretizzazione  nel enza tra  la soluzione vera e quella calcolata a partire dal punto 

  dice  errore  locale  di  troncamento,  l’errore ttuare il singolo passo   

 

rdine del metodo: Si dice ordine del metodo, il piu’ grande intero positivo p per il 

Convergenza: Un metodo si dice convergente se: 

o equivalentemente se 

lim 0 

 

Sulla base di tali definizioni andiamo ad analizzare i vari metodi fin qui descritti. 

Metodo di Eulero esplicito

Con riferimento a metodi di questo tipo, daremo or

punto    la differiniziale del dominio di integrazione: 

 

 

Errore  locale  di  troncamento:  Siintrodotto dal metodo numerico nell’effe

Oquale risulti: 

 

lim  

: In questo caso l’errore locale di troncamento e’ dato da: 

,  

icordando  che  i  valori  veri  e  quelli  approssimati  coincidono  per  definizione  nel rpunto iniziale, si ha che  ,  , da cui 

 

Page 8: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

7  

sostituendo in questa equazione l’e pansione in serie di Taylor, si ha: 

2

s

 

da cui: 

Questa espressione, inoltre, implica che 

 

percio’ il metodo e’ del primo ordine e l’errore va a zero come  . 

Metodo di Eulero  plicitoim : In questo caso l’errore locale di troncamento e’ dato da: 

,  

Utilizzando l’equazione (vedi sopra) 

2,  

2

e sostituendo, si ha: 

 

Anche in questo caso, quindi, come nel caso esplicito, si ha: 

 

el primo ordine e l’errore va a zero come  . 

Metodo di Crank‐Nicolson

percio’ il metodo e’ d

: In questo caso l’errore locale di troncamento e’ dato da: 

2, ,  

Ricordando che 

,  

e che applicando il metodo dei trapezi (con l’errore) abbiamo: 

Page 9: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

8  

,, ,

2 12 

do, semplificando e ricordando che   e che quindi    

si ha: 

Quindi, sostituen

12 

In questo caso quindi, 

 

 zero 

1.1.5. Analisi di Stabilita’. 

Un  algoritmo  e’  sta ile  se  non  amplifica  gli  errori  che  compaiono  nel  corso  del udieremo  ora  la  stabilita’  dei metodi  di  Eulero  esplicito  e 

implicito e del metodo di Crank‐Nicolson. Per fare questo studiamo il caso specifico dell’equazione: 

percio’ il metodo e’ del secondo ordine e l’errore va a come  . 

 

bprocesso  risolutivo.  St

 

0 .  La  soluzione  esatta  di  questa  equazione  e’ . 

con  condizione  iniziale 

Metodo  di  Eulero  esplicito:  Dobbiamo  provare  che  l’errore  si mantiene  limitato. Applicando il metodo all’equazione  specifica: 

1 1  

Poiche’ la soluzione esatta tende a zero per valori crescenti di x, vogliamo che anche la soluzio e numerica tenda a zero per valori crescenti di n, in modo da mantenere 

one numen, solo se: 

nlimitato l’errore. La soluzi rica, fissato h, tende a zero per valori crestenti di 

|1 | 1 

cioe’ per  

Page 10: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

9  

1 1 1     0 

 Eulero esplicito e’ stabile sotto condizione. 

Metodo  di    implicito

2

essndo h positivo per definizione la condizione di stabilita’ del metodo e’: 

Il metodo di

Eulero :  In  questo  caso  la  soluzione  numerica  dell’equazione specifica e’: 

11

11

 

iasi valore di h. Il metodo di Eulero implicito e’ incondizionatamente stabile. In questo caso la soluzione tende a zero per valori crescenti di n per quals

Metodo  di  Crank‐Nicolson:  In  questo  caso  la  soluzionspecifica e’: 

e  numerica  dell’equazione 

222

22

 

In questo caso la soluzione tende a zero per valori crescenti di n per qualsiasi valore di h. Il metodo di Crank‐Nicolson e’ incondizionatamente stabile. 

1.1.6. Metodi Runge‐Kutta: 

ge‐Kutta sono metodi  in cui  l’approssimazione della soluzione  in    dipende  solo  dai  valori  della  funzione  in    .  Sono  quindi metodi  espliciti  ad  un passo.  Il  vantaggio  rispetto  al  metodo  esplicito  di  Eulero  e’  che  questi  metodi simulano gli effetti delle derivate di ordine superiore valutando  la  funzione    

etodo di Heun: 

semplice dei metodi Runge‐Kutta ed e’ un analogo del metodo  di  Crank‐Nicolson,  ma  rimane  esplicito  poiche’  usa  Eulero  esplicito  per 

ula

 

I metodi Run

in piu’ punti dell’intervallo  ,   in modo da  avere un errore di  troncamento 

locale che sia di ordine superiore a uno. 

M

Il metodo di Heun e’  il piu’ 

predire il valore di    da usare in  , . Consideriamo, quindi, la form  di Crank‐Nicolson: 

Page 11: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

10  

2, ,  

Approssimiamo ora  il valore  che si trova  in  ,    effettuando un passo del metodo di Eulero esplicito. Abb

 di   iamo quindi: 

2, , ,  

,  

,   , , ,  

Si ottiene la formula di Heun: 

Ponendo  

,  

,

 di Heun e’ un analogo del metodo di Crank‐Nicolson, per questo varrano tutte le considerazioni fatte in precedenza per il metodo di Crank‐Nicolson. Percio’, tra le altre cose, il metodo di Heun e’ un metodo del  secondo ordine. 

Metodo di Runge: 

passo  con  il metodo di Eulero esplicito  fino a meta’ dell’intervallo  / /2    per  calcolare  il  valore  della  funzione  /  . Effettuiamo poi un pass completo  in cui  la derivata e’ calcolata  in  / ,  /  

(cioe’ a meta’ dell’intervallo). Si ottiene cosi’ il metodo di Runge: 

/ 2

 

Poiche’  il metodo

Supponiamo di effettuare un 

,  

2

da cui: 

/ , / ,2

,  

La formula di Runge, percio’, sara’: 

Page 12: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

 

,  

2,

 

 

Figura 2: Nel metodo di Eulero  la derivata nel punto  iniziale di ogni  intervallo e' estrapolata per trovare  il  successivo della funzione. Questo metodo ha un'accuratezza del primo ordine 

 

 

Figura 3: Nel metodo di Runge  l'accuratezza del secondo ordine viene ottenuta utilizzando  il valore della 

e  il  valore  nel  punto nale.  Nella  figura  i  cerchi  pieni  rappresentano  i  valori  finali  della  funzione  mentre  i  cerchi  aperti 

  i valor della  funzione che servono per calcolare  la derivata da utilizzare nel singolo passo. Questi valori vengono scartati alla fine del passo. 

valore

derivata nel punto iniziale per determinare il valore della funzione nel punto intermedio dell'intervallo epoi utilizzando  questo  valore  per  determinare  la  derivata  da  utilizzare  per  determinarfirappresentano

11  

Page 13: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

12  

1.1.7. Sui metodi di Runge‐Kutta: 

I  due metodi  che  abbiamo  fino  ad  ora  (Heun  e  Runge)  si  possono  vedere  piu’  in generale nella forma: 

 

con 

,  

,  

icienti ripartiamo dallo sviluppo  in serie di Taylor della funzione  : 

12

I coefficienti  ,   , ,   sono scelti in modo che l’errore locale di troncamento sia  

 

cioe’ i metodi siano del secondo ordine. 

Per ricavare  i valori di questi coeff

16

 

dove al solito: 

16

 

e’ il resto della formula di Taylor con   un punto opportuno nell’intorno di  . 

Ricordando che: 

,  

,  

e che, 

, , , ,  

si ottiene 

,2

, , ,  

Page 14: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

13  

Trascuarando il resto  si ha la formula approssimata: 

,2

, , ,        1  

 confrontare questa espressione con quella generale data all’inizio  .  Prima  di  fare  questo,  pero’,  conviene espandere   in se e di Taylor. Ricordando che in generale la formula di Taylor per una funzione a due variabili  ,   nell’intorno di  ,  e’ data da: 

∆ , ∆ , , ∆ , ∆ ∆ , ∆  

si puo’ espandere    nell’espressione: 

, , ,  

 che  ,   si ha: 

, , , ,  

Sostituendo questa espressione in quella generale si ha: 

  ,, , , ,  

  ,

1       12

A questo punto vorremo

ri

ricordando

Raccogliendo a fattor comune si ha: 

, , ,  

Confrontando questa espressione con  lo sviluppo  in serie di Taylor approssimato si hanno le seuqnti identita’: 

       12 

Questo  sistema di  tre  equazioni  in  quattro  incognite non  ammette  una  soluzione unica. Cio’ significa che non esiste un solo metodo del secondo ordine, ma un’intera 

Nel caso del metodo di Heun si ha: 

1

famiglia. 

2      

21        1        1 

Page 15: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

14  

Nel caso del metodo di Runge si ha: 

0       1       12       

12 

Le  considerazioni  fin  qui  fatte dimostrano  che  i metodi  di Heun  e  di Runge  sono entrambe metodi del secondo ordine, poiche’    in entrambe  i casi nell’epansione in 

iti poiche’ per procedere dal punto  inizale     al punto finale        si  richiede  solo  la  conoscenza  della  soluzione  nel  punto  iniziale     percio’  il  passo  dovra’  soddisfare  la  condizione  di  stabilita’.    Evidentemente  si possono costruire metodi di Runge‐Kutta di qualsiasi ordine, troncando  lo sviluppo 

nta sempre piu’ difficile da ricavare.  

u’  usato  e’  quello  del  4°  ordine,  definito  dalle seguenti relazoni: 

 

,  

2

 serie  di  Taylor  della  funzione  si  e’  trascurato  il  resto  .  Entrambe  i metodi  sono metodi esplic

in serie di Taylor a termini via via di ordine piu’ elevato. All’aumentare dell’ordine il procedimento diventa via via piu’ complitato e allo stesso tempo la stima degli errori locali dive

Tra  i metodi  di  Runge‐Kutta,  il  pi

,2

 

2,

,  

6 3 3 6 

 

Page 16: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

15  

 

Figura 4: Nel metodo di Runge‐Kutta del quarto ordine  la derivata viene calcolata quattro volte per ogni passo: una volta nel punto  iniziale (1), due volte per un valore di prova nel punto  intermedio (2,3), e una volta per un valore di prova nel punto finale (4). Utilizzando queste derivate si calcola il valore finale della nzione. 

 

Dato che fino all’ordine 4 sono richieste tante valutazioni del sistema ODE quanto e’ l’ordine  stesso,  mentre  per  ordini  maggiori  occorre  effettuare  un  numero  di valutazioni  via  via  sempre  piu’  grande  dell’ordine  del metodo,  il  Runge‐Kutta  del quarto  ordine  rappresente  un  ottimo  compromesso  tra  precisione,  sforzo computazionale e stabilita’. 

 

I pro e i contro dei metodi Runge‐Kutta sono i seguenti: 

fu

PRO: 

1. Sono efficiente anche quando  la  soluzione non viene ben approssimata  con polinomi; 

2. Sono  di  solito  poco  sensibili  a  eventuali  discontinuita’  delle  funzioni  del 

rse  di  calcolo,  il  tempo  complessivo  di  integrazione risulta basso. 

 

sistema; 3. Richiedono poca memoria; 4. E’ facile cambiare il passo di integrazioni in un momento qualsiasi; 5. Il tempo di calcolo dell’algoritmo non e’ elevato. Quindi per sistemi ODE che 

richiedano  poche  riso

Page 17: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

16  

RO:CONT  

altre alternative; 

cale non è facile da calcolare; problemi in cui siano richiesti valori delle variabili 

  inferiori  a 

odi Multistep: 

recedenti ed eventualmente nel punto  finale.  In i se     dipende     ma non dai valori 

erico metodo multistep a  1  passi 

 

dove  , . Se  0    il metodo e’ ESPLICITO. Se, al contrario, 

0    il metodo e’ IMPLICITO. Alcuni esempi di metodi multistep sono: 

 Nel metodo del punto medio la derivata  ,  si assume 

 

Da cui per i valori approssimati avremo: 

2  

1. Il  numero  di  calcoli  del  sistema  ODE  e’  generalmente  maggiore  rispetto aquello di 

2. Non sono adatti a sistemi stiff; 3. L’errore lo4. Non si prestano a risolvere 

dipendenti molto  ravvicinati,  ossia  con  valori  di  h  decisamentequelli richiesti dalla precisione dell’algoritmo; 

5. Non  e’  possibile  risolvere  sistemi  in  cui  le  derivate  siano  presenti  in  forma implicita. 

 

1.1.8. Met

I  metodi  Multistep,  utilizzano  i  valori  della  variabile  dipendente    e  della  sua derivata   accumulati nei punti pgenerale un metodo si dice a     pass

precedenti     1 . Pertanto un gene’ della forma: 

Metodo del punto medio:costante nell’intervallo  ,   ed uguale al valore in  ,  , . Pertanto si avra’: 

, 2 ,  

Page 18: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

17  

al  solito,  .  Confrongenerale  si o  del secondo ordine a due passi, esplicito, in cui 

0,   

Metodo di Simpson: Nel metodo di Simpson,  l’integrale viene approssimato con  la formula di Simpson. Percio’ in questo caso si avra’: 

,

2

dove,  , tando  questa  espressione  con  quella   deduce  che  1  ,    e  0.  Pertanto  questo  e’  un metod

1, 2,    0 

 

 

64  

imati  o: 

3

Da cui per i valori appross avrem

14  

Anche  in questo caso, confrontando questa  formula con quella generale si deduce he  1  e  1/3.  Pertanto,  questo  e’  un metodo  del  quarto  ordine  a  due 

  1/3 

 

Metodi di Adams: Nei metodi di Adams la funzione  ,   viene sostituita con un polinomio interpolatore di Lagrange di grado    ,   ,che interpola la funzione in 

1 punti. Cioe’: 

  ,    

di di  ms‐Bashfe  ,     nell’integrale  con  un  polinomio  interpolatore  di 

Lagrange di grado   che interpola la funzione  ,   nei  1  nodi  ,   , . 

cpassi, implicito, in cui: 

0,    1, 4/3,  

 

Meto Ada ort: Un metodo di Adams‐Bashfort a  1  passi e’ ottenuto sostituendo  la  funzion

Page 19: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

18  

 

Figura 5: Nel metodo di Adams‐Bashfort a q+1 passi la funzione y’(t)=f(t,y(t)) viene sostituita con un polinomio di Lagrange di grado q che interpola f(t,y(t)) nei nodi tn‐q,...,tn. 

 

Un metodo di Adams‐Bashfort a   passi e’ di ordine  . 

Metodi di Adams‐Moulton: Un metodo di Adams‐Moulton a  1  passi e’ ottenuto sostituendo  la  funzione  ,     nell’integrale  con  un  polinomio  interpolatore  di Lagrange  di  grado  1  che  interpola  la  funzione  ,     nei  2    nodi 

,   , . 

 

 

 

 

 

 

 

 

 

 

Page 20: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

 

19  

igura 6: Nel metodo di Adams‐Moulton a  q+1 passi la funzione y’(t)=f(t,y(t)) viene sostituita con un olinomio di Lagrange di grado q+1 che interpola f(t,y(t)) nei q+2 nodi tn‐q,...,tn+1. 

Un metodo di Adams‐Moulton a   passi e’ di ordine  1. 

 

In  generale un  Polinomio  interpolatore  di  Lagrange  che passa per  i    1   punti 

, , , , ,  , … , , , ,  cioe’ 

, , ,  , … , , ,  e’  un  polinomio  di  grado    ed  e’  dato 

dall’espressione: 

 

dove  

Fp

 

Scritto esplicitamente diventa: 

Page 21: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

20  

 

 

  

 

 l’epressione generale per i metodi di Adams‐Bashfort a  1  passi sara’: 

 

L’espressione generale per i metodi di Adams‐Moulton a  1  passi sara’: 

 

 

Facciamo ora alcuni esempi dei metodi di Adams‐Bashfort: 

• Metodo  a  un  passo,  0:  in  questo  caso  abbiamo  solo  un  nodo  , , percio’ il metodo sara’ del primo ordine: 

 da cui: 

   

        Eulero Esplicito 

biam , , , , percio’ il metodo sara’ del secondo ordine: 

Quindi

• Metodo  a  due  passi,  1:  in  questo  caso  ab o  due  nodi 

 

da cui: 

23  

 

Page 22: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

21  

• Metodo  a tre  passi,  2:  in  questo  caso  abbiamo  tre  nodi  , , ,  , , percio’ il metodo sara’ del terzo ordine: 

12

 ,

23 16 5  

 

• Metodo a quattro passi, 3: in questo caso abbiamo quattro nodi  , , , ,  , ,  , ,  percio’  il metodo  sara’  del  quarto 

ordine: 

 

2455 59 37 9  

 

Faccia alcuni esempi dei metodi di Adams‐Moulton: 

• Metodo  a  un  passo,  0:  in  questo  caso  abbiamo  due  nodi  , , , , percio’ il metodo sara’ del secondo ordine: 

mo ora 

 

2

da cui: 

        Crank‐Nicolson 

 

, ,  il metodo

12

• Metodo  a  due  passi,  1:  in  questo  caso  abbiamo  tre  nodi  , , ,   percio’  sara’ del terzo ordine: 

5 8     

 

• Metodo a tre passi,  2: in questo caso abbiamo quattro nodi  , , sara’ del quarto ordine: , ,  , ,  , , percio’ il metodo 

249 19 5     

 

 

Page 23: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

22  

I pro e i contro dei metodi multistep sono i seguenti: 

PRO: 

1. Forniscono senza difficolta’ una stima dell’errore locale; egra

CONTRO:

2. Richiedono un minor numero di calcoli di funzione per passo di int zione 

 

1. Tutti i metodi Multistep espliciti sono instabili, quindi occorre utilizzare quelli impliciti; 

2. Gli algoritmi multistep non sono autosufficienti. Occorre cioe’ disporre di altri algoritmi per avviare l’integrazione e calcolare i punti coinvolti nella formula; 

i mopasso integ  nuovamente dei punti necessari alla formula stessa. 

 

.1.9. Metodi Predictor Corrector. 

metodi  implici  (qualunque  essi  siano)  hanno  una  difficolta’  aggiuntiva  legata  al generale  in  forma  implicita.  In 

generale quindi, ogni metodo implicito, porta ad una equazione del tipo: 

Ψ  

 un’equazione di questo tipo, quindi, deve essere effettuata nella maggior parte dei 

ricamente. Tra  le varie possibilita’  (metodo di Newton‐Raphson, metodo 

prova di   con un metodo esplicito. Correggere  il valore con un 

 

3. Le  formule  sono  adatte  ad  un  passo  costante.  Ogniqualvolta  s difica  il  di  razione occorre reinizializzare  l’integrazione per poter disporre

1

I fatto che  l’equazione per  il nuovo valore   e’  in 

dove  la  funzione  Ψ  puo’  essere  anche  molto  complicata.  La  risoluzione  di

casi numedel  punto  fisso)  c’e’  quella  del  Predictor‐Corrector.  L’idea  e’  quella  di  fornire  un primo valore di metodo  implicito  ed  eventualmente  iterare  il  metodo  implicito  fino  alla convergenza. 

Piu’ specificatamente si ha la seguente procedura.  

Dato un generico metodo esplicito 

Page 24: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

23  

, p

 

nel nuovo punto  ,   ), 

e  ,  senza ricalcolarlo. 

Quindi lo schema generale base e’ del tipo: 

Predizione:   ∑ ∑   (metodo esplicito) 

  (metodo implicito)

Valutazione (Evaluation):   ,  

 

Tuttavia si possono utilizzare per uno schema predictor‐corrector piu’ varianti: 

• PEC 

• P(EC)k 

 

mediante  esso  si  determina  un  valore  approssimato  .  Per  migliorare 

l’approssimazione si usa questo valore,   er calcolare  ,  in 

modo da correggere   con un metodo implicito 

Percio’  il  metodo  esplicito  predice  una  approssimazione  e  quello  implicito  la corregge. La formula implicita puo’ essere ripetuta un certo numero di volte 

 

Dopo   correzioni si puo’ ricalcolare  il valore di   

oppure utilizzare il valore precedent

Valutazione (Evaluation):   ,  

Correzione:   ∑ ∑  

• PECE (quello base appena descritto) 

• P(EC)kE 

Page 25: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

24  

a predictor e’ di ordine   e  lo  schema corrector  e’  di  ordine    ,  allora  il metodo  predictor‐corrector  PECE  si  comporta 

todo di ordine min  , 1 . 

corrector sono i seguenti: 

ctor‐Corrector EE+CN: 

predictor      ,            (Eulero Esplicito, ordine  1) 

corrector     

Ordine del predictor‐corrector: Se  lo  schem

come un me

Esempi di schemi predictor‐

Metodo Predi

, ,      (Crank‐Nicolson,  2) 

In questo esempio il metodo predictor‐corrector e’ di ordine 2. 

Metodo Predictor‐Corrector AB2+AM3: 

predictor        3 , ,      (Adams‐Bashfort a due 

e  2) 

corrector       

passi , ordin

5 , 8 , ,      

‐Moulton a tre passi,  3) 

quando si tratta con ODE di ordine superiore al primo. Un sistema di ODE del primo ordine si dice stiff quando esistono due o piu’ 

i della variabile indipend ile dipendente. Per esempio, consideriamo il seguente set di ODE: 

998 1998  

999 1999  

con condizioni iniziali  0 1, 0 0. Effettuando la trasformazione: 

2  

 

si ha la soluzione analitica: 

(Adams

In questo esempio il metodo predictor‐corrector e’ di ordine 3. 

 

1.1.10. Sistemi Stiff: 

I sistemi stiff si possono presentare 

scale different ente sulle quali cambia la varib

Page 26: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

25  

2  

 

ione: 

 

0  

e sia del tipo: 

 

0  

con  0  la condizione di stabilita’ per  il metodo di eulero esplicito e’  2/ . Percio’ e’ chiaro che se si integrasse il sistema di partenza con un metodo esplicito, questo  rich derebbe un passo  1/1000 per  garantire  la  stabilita’. D’altronde, 

  e    il termine   e’ trascurabile non appena ci si allontana dall’origine. Quindi  il problema generico  legato all’applicazione di metodi espliciti  per  equazioni  stiff  e’  quello  di  dover  seguire  la  soluzione  sulle  scale  piu’ piccole, per garantire  la stabilita’ del metodo di  integrazione, anche se  i requisiti di accuratezza della soluzione permettano un passo molto maggiore. L’unico modo di 

e gia’ visto, il metodo di Eulero implicito e’ un metodo del primo ordine. Se si vuole ottenere una maggior accuratezza per sistemi stiff, si puo’ ricorrere a metodi  impliciti di ordine superiore  l primo.  

Nell’analisi di stablita’ abbiamo gia’ visto che nel caso dell’equaz

 

con  soluzione  esatta    i  metodi  impliciti  erano  incondizionatamente stabili, mentre  il metodo di eulero esplicito era stabile solo per  2. Ora, si puo’ facilmente dimostrare che nel caso in cui l’equazion

 

ieper determinare i valori di 

ovviare a questo problema e’ quello di utilizzare metodi impliciti. Com

e a

Page 27: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

26  

stante questo, pero’, essa determina la stabilita’ della soluzione (linea tratteggiata e puntinata) anche per scale piu’ grandi 

 

 

1.2. PROBLEMI AI VALORI AL CONTORNO.   

1.2.1. Introduzione 

Quando una ODE deve soddisfare delle condizioni al contorno che corrispondono a piu’ di un valore per la variabile dipendente si ha un problema ai valori al contorno. Tipicamente,  le  condizioni  che  devono  essere  soddisfatte  vengono  speficicate  nel punto  iniziale  e  nel  punto  finale  dell’intervallo  di  integrazione.  La  distinzione principale  tra problemi ai valori  iniziali e problemi ai valori al  contorno e’  che nel 

al o  caso, 

to  iniziale non determinano una soluzione nica con cui partire, cioe’ una soluzione che soddisfi  le condizioni al contorno nel 

punto  iniziale non  verifichera’  in  generale  le  rimanenti  condizioni  al  contorno nel unto finale. Nel caso dei problemi ai valori al contorno, quindi, sara’ necessaria una 

procedura iterativa attraverso la quale tra tutte le soluzioni possibili che soddisfano 

 

Figura 7: Esempio di instabilita’ incontrata integrando un’equazione stiff. Qui si suppone che l’equazione abbia due soluzioni (linea solida e linea tratteggiata). La seconda, va rapidamente a zero non appena ci si discosta dall’origine. Nono

primo caso la soluzione, nota nel punto iniziale, si aggiorna passo dopo passo fino punto  iniziale  attraverso  integrazioni  numeriche  successive;  nel  secondinvece,  le condizioni al contorno nel punu

p

Page 28: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

27  

lezioni” quella che soddisfi anche le condizioni  al  contorno  nel  punto  finale.  Questo,  in  generale,  richiedera’ l’integrazione delle ODE su tutto l’intervallo di integrazione, ovvero l’applicazione di 

lassamento” (vedi sotto), per certo un numero di iterazioni. Per 

le condizioni al contorno nel punto iniziale si “se

una procedura di “riquesto motivo la soluzione dei problemi ai valori al contorno richiedono uno sforzo maggiore  sia  in  termini di programmazione dell’algoritmo,  sia  in  termini di  tempo calcolo.  

Un problema ai valori al contorno standard ha  la seguente  forma: si vuole che un dato  set  di    equazione  differenziali  ordinarie  del  primo  ordine,  soddisfi   condizioni  al  contorno  nel  punto  iniziale    e  un  rimenente  set  di   condizioni al contorno nel punto finale  . 

Percio’ potremo scrivere il set di ODE come: 

, , , … ,           1,2, … ,  

Nel punto iniziale  , la soluzione dovra’ soddisfare le   condizioni al contorno: 

, , , … , 0                  1,     

mentre nel punto finale   dovra’ soddisfare le rimanenti   condizioni  

, , , … , 0                1,   

Esistono  due  classi  distinte  di metodi  numerici  per  risolvere  problemi  ai  valori  al 

 un  le’e le  equazion

  i metodi numerici per problemi ai valori valori  di  tutte  le  variabili  dipendenti 

nell’altro  estremo  dell’intervallo  di  integrazione.  In  generale  questi  valori  non saranno  consistenti  con  le  condizioni  al  contorno  per  questo  estremo,  percio’  si dovra’ operare  in modo tale da aggiustare opportuname te  i valori forniti al punto 

i  al contorno volute. 

 

 

contorno. Nel metodo  di  shooting  vengono  scelti  dei  valori  per  tutte  le  variabili dipendenti in  estremo; tali valori ovviamente dovranno essere consistenti con  condizioni  al  contorno  relative  a  quell stremo.  Fatto  questo,  i differenziali vengono  integrate con uno deiniziali.  In  questo  modo  si  determinano  i 

niniziale  per  ottenere  nel  punto  finale  dei  valori  consistenti  con  le  condizion

Page 29: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

28  

hem   prova  

na variante di questo metodo e’  il metodo del  fitting.  In questo caso si scelgono dei  valori  di  prova  per  tutte  le  variabili  dipendenti  in  entrambe  gli  estremi  di integrazione  (e che siano consistenti con  le condizioni al contorno ai due estremi). Fatto questo si integrano le equazioni partendo dall’interno e partendo dall’esterno, con  uno  dei metodi  per  i  problemi  ai  valori  iniziali,  fino  ad  un  punto  intermedio detto punto di  fitting. Le discrepanze  tra  i valori delle varibiali dipendenti ottenuti con  l’integrazione  dall’interno  e  con  quella  dall’esterno  vengono  utilizzate  per correggere i valori di prova  in modo da azzerare le discrepanze stesse. 

Il metodo  del  rilassamento,  invece,  utilizza  un  approccio  differente.  Le  equazioni differenziali  sono  rimpiazzate  da  equazioni  alle  differenze  finite  su  una  griglia  di punti (mesh points) che copre tutto l’intervallo di integrazione. Si fornisce quindi una 

rilassamento, consiste nell’aggiustare tutti  i valori di tutte  le variabili dipendenti  in odo  tale ottenere  la  verifica delle equazioni  in  tutti  i mesh, entro un  tolleranza 

 

Figura  8:  Sc atizzazione  del metodo  di shooting.  I  valori  di    della  soluzione  che soddisfano  le condizioni al contorno nel punto iniziale vengono “lanciati” al punto finale. Le discrepanze tra i valori cosi’ ottenuti e  le condizioni al contorno nel punto finale vengono utilizzate per aggiustare  i valori di prova nel punto iniziale. Il processo iterativo avra’ termine solo quando le discrepanze saranno inferiori ad una certa tolleranza fissata a priori. 

 

U

prima  soluzione  di  prova  che  consiste  nell’insieme  dei  valori  di  tutte  le  variabili dipendenti  corrispondenti ad ogni mesh point. Questi  valori non  soddisferanno  in generale  le equazioni alle differenze finite per ogni coppia di mesh, ne’ tantomeno le  condizioni  al  contorno  ai  due  estremi.  L’iterazione,  detta  comunemente 

mdata, nonche’ delle condizioni al contorno ai due estremi. 

Page 30: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

 

29  

c

 metodo  del  rilassamento  e’  piu’  efficiente  del  metodo  di  shooting  quando  le equazioni  sono  tali  che  i  valori  della  soluzione  agli  estremi  dell’intervallo  di integrazione sono molto sensibili. Percio’, piccole variazioni dei valori di prova della soluzione  ad  un  estremo  determinano  grandi  variazioni  dei  valori  ottenuti all’estremo opposto. In questo caso, quindi, e’ difficile modificare i valori iniziali per ottenere i valori finali compatibili con le condizioni al contorno volute. 

Essendo,  in  generale,  il  metodo  del  rilassamento  piu’  robusto  del  metodo  di shooting, questo e’ preferito quando le funzioni non variano lentamente ma hanno variazioni ripide in funzione della variabile indipendente. 

 

va  iniziali  di  tutte  le  variabili  dipendenti    nel  punto  iniziale intervallo  di  integrazione    compatibili  con  le    condizioni  al  contorno. 

Pertanto  ci  saranno      valori  iniziali  di  prova  che  potranno  essere 

 

Figura 9: Schematizzazione del metodo del rilassamento. Inizialmente si fornisce una soluzione di prova che soddisfa in modo approssimato sia le equazioni  he le condizioni al contorno. In un processo iterativo, poi, la soluzione di prova viene via via modificata in modo da ottenere la miglior verifica delle equazioni e delle condizioni al contorno. 

 

Il

1.2.2. Il Metodo di Shooting 

Il primo passo che occorre effettuare nel metodo di shooting e’ quello di specificare gli    valori  di  prodell’

assegnati liberamente senza condizioni. Supponiamo che questi valori “liberi” siano contenuti  in  un  vettore    di  dimensioni    .  Quindi,  ad  ogni  particolare    

Page 31: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

30  

.  In generale le no il

,                 1, … ,  

 trovare un vettore di valori   che porti a zero tutti  i 

er  g

corrispondera’  un  particolare    ,  cioe’  il  vettore  dei  valori  di  prova  iniziali. L’integrazione del sistema di   ODE con uno dei metodi numerici per  i problemi ai valori  iniziali  portera’  ad  un  particolare  ,  cioe’  il  vettore  dei  valori  della soluzione  corrispondenti al punto  finale dell’intervallo di  integrazione

    non  verifichera’      condizioni  al  contor   date  per    punto  finale. Supponiamo che le discrepanze tra le condizioni al contorno date e quelle ottenute con  l’integrazione  appena  effettuata  siano  contenute  in  un  vettore    ,  anch’esso ovviamente di dimensioni  , percio’ risultera’: 

Il problema quindi e’ quello divalori  del  vettore  .  Dobbiamo  quindi  studiare  come  varia    al  variare  di  .  Il diff enziale totale della generica componente di   sara’ dato in enerale da: 

 

Percio’ per azzerare   e’ su iciente variare     modo che: 

          1, … ,  

Cioe’: 

ff in 

       1, … ,  

che scritto in forma compatta risulta: 

·

essendo   la matrice Jacobiana   la cui generica componente e’ data da:  

 

Una volta  risolto questo sistema avremo  le correzioni   da apportare ai valori di 

prova di  , che chiameremo ora  . Percio’ avremo che  i nuovi valori ”liberi “ di prova saranno: 

Page 32: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

31  

Ora  se  le  equazioni  differenziali  sono  lineari,  allora questa  procedura  porta direttamente  alla  soluzione  finale.  Se  invece  le  equazioni  sono  non  lineari,  allora 

ecessario  iterare  l’intera procedura  fino a  che  le discrepanze  saranno  tutte inferiori ad una certa tolleranza scelta a priori. 

E’  importante notare  che  in generale non  sara’  facile  calcolare  le derivate parziali che compaiono nella matrice Jacobiana in modo analitico. Si dovra’ quindi ricorrere 

merica  effettuando    integrazioni  separate delle di   ODE, in cui ogni valore del vettore   viene incrementato di una quantita’ scelta a priori, seguite dal calcolo della derivata parziale dato da: 

sara’ n

alla  valutazione  di  queste  per  via  nu

, … , ∆ , … ,… , , …∆

 

 

1.2.3. Metodo del Fitting 

oting assume implicitamente che i “lanci” siano sempre in grado di attraversare  il  dominio  di  integrazione,  dal  punto  iniziale  al  punto  finale,  e convergere  ad  una  soluzione  corretta,  anche  nelle  prima  fase  del  processo  di 

  m

fanno  semplicemente  “morire”  il programma.  Inoltre  in  alcuni  casi  l’alta  non  linearita’  delle  equazioni  implica  una 

 dei 

  fitting  e  tale  che 

,  e poi una  integrazione dell’esterno, da    a  .  Se,  come  abbiamo 

detto in precedenza, il numero di condizioni al contorno da imporre in   sono   ed il numero da  imporre  in   sono  , allora ci saranno  valori “liberi”  in   ed   alori “liberi”  in   delle variabili dipendenti da assegnare. Supponiamo che  i valori liberi in   siano contenuti in un vettore   di dimensioni   e che i valori liberi in   

    di  dimensioni    .  Allora  ad  ogni  scelta  di   corrispondera’ un  vettore di  valori  iniziali  in  ,      ed  analogamente  ad ogni 

Il metodo di sho

convergenza.  In  realta’,  in  alcuni  casi,  se i  valori  di  prova  sono  olto  lontani  dai valori  corretti,  la  soluzione  non  e’  in  grado  di  essere  “lanciata”  da    a    senza incontrare  espressioni  incalcolabili  o  catastrofiche,  come,  ad  esempio,  radici  di numeri  negativi  o  divisioni  per  zero,  che 

diffolta’ di convergenza dovuta all’elevata sensibilita’ delle condizioni in   da quelle in  .  In  questi  casi  si  puo’  superare  questa  difficolta’  attraverso  il metodo  dello shooting to a fitting point o piu’   semplicemente  il metodo l fitting. L’idea e’ che invece  d integrare  dal  punto    al  punto  ,  si  effettua  prima  una  integrazione dall’interno,  da    fino  ad  un  certo  punto  ,  detto  punto  di

 v

siano  contenuti  in  un  vettore

Page 33: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

32  

i

an , m

e i v

differenziale 

cioe’: 

scelta di   corrispondera’ un vettore di valori  iniziali  in  ,    . L’integrazione del sistema d    ODE da   a  , con uno dei metodi numerici per i problemi ai valori 

iniziali, permettera’ di definire tutti i valori di tutte le variabili indipendente al punto di  fitting,  .  Analogamente  l’integrazione  dall’esterno  da    a    portera’  a 

definire il vettore  . In generale questi due vettori non coincider no a anzi 

si potra’ un vettore di discrepanze al punto di fitting, di dimensioni   

 

dove la generica componente e’: 

          1, … ,  

Queste discrepanze dipenderanno in general dai valori di   e  , percio’ si dovra’ operare su quest alori affinche’  le discrepanze al punto di fitting vadano a zero. Il 

totale di   sara’ dato in generale da: 

          1, … ,  

co  e    sa

          1, … ,  

Percio’ la  ndizione sulle variazioni dei valori “liberi” in   in ra’: 

          1, … ,  

La  soluzione  di  ques ema    fornira’,  quindi,  le  correzioni  e    da apportare ai valori di prova di   e di  , percio’ avremo che i nuovi valori 

to  sist”liberi “ di 

prova saranno: 

 

 

 

Come nel caso precedente, se  le equazioni differenziali sono  lineari, allora, questa procedura porta direttamente alla soluzione finale. Se invece le equazioni sono non 

Page 34: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

33  

ch

 occorrera’,  alla valutazione numercia delle derivate parziali come descritto in precedenza.  

 1.2.4. Metodo del Rilassamento 

Il metodo del rilassamento ha un approccio al problema completamente differente trovare  la 

combinazione di valori “liberi” per  le variabili dipendenti agli estremi dell’intervallo di integrazione, che minimizza gli scarti (o al punto di fitting o a uno dei due estremi di  integrazione),  in  questo  caso  l’idea  e’  quello  di  fornire  una  intera  soluzione  di prova su tutto l’intervallo di integrazione, suddiviso in una griglia di mesh points,  e 

soluzi r o nerrt o

glio il metodo. 

Consideriamo il generico sistema di    ODEs: 

lineari,  allora  sara’ necessario  iterare  l’intera procedura  fino  a  e  le discrepanze saranno tutte inferiori ad una certa tolleranza scelta a priori. 

Anche  in questo caso   in generale, procedere

 

da  quello  dei  metodi  di  shooting  e  di  fitting.  Infatti,  invece  di 

variare  in maniera opportuna  le  oni di prova nei singoli mesh pe tte e  la verifica  delle  equazioni  di  pa enza,  pportunamente  linarizzate.  Andiamo  ora  a spiegare piu’ in detta

, , … ,                  1, … ,  

  primo  passo  suddividiamo  l’intervallo  di  integrazione  in    punti  che Comechiameremo mesh points (o piu’ semplicemente mesh) e trasformiamo le equazioni differenziali  in  equazioni  algebriche  alle  differenze  finite.  Percio’  per  la  generica coppia di mesh  , 1   e per la generica equazione  ‐esima risultera’: 

, , … ,      , ,

, ,, , … ,  

/ 2

/ , / , /

dove 

 

era soluzione di prova che consiste negli   

valori  , , … , , , 1, … , . In generale le equazioni linearizzate  

Supponiamo di fornire ora un’int

Page 35: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

34  

,, ,

, ,/ , , / , … , , /       1, . . ; 1, … , 1 

non saranno verificate, cioe’ risultera’: 

, , 0      1, . . ; 1, … , 1 

dove  ,  e’ la “non” verifica dell’equazione  ‐esima per la coppia di mesh  , 1 . 

Poiche’  ,  dipendera’ dalla soluzione di prova, studiando da  la sua dipendenza da 

questa, si potra’ operare sui valori di prova per minimizzare gli scarti. Il differenziale totale di   e’ dato da: ,

,,

,,

,

,,  

e  

,

Ovviamente  ,  dipendera’ dalle  variabili  in      in  1. Percio’ per azzerare gli errori  si  dovra’  imporre  che  le  correzioni  ai  valori  iniziali  di  prova  soddisfino  il seguente sistema di equazioni: 

,,

, ,

,,

,,       1, . .   1, … , 1 

Si  ha  quindi  un  sistema  di  1   equazioni  nelle    incognite 

, … , , , 1, … , ,  che  rappresentano  le  correzioni  da  apportare  ai 

valori di prova per azzerare le “non” verifiche. 

Le  restanti    equa oni  derivano  dalle    condizioni  al  contorno.  Supponiamo  al : 

, , ,

e che nel secondo estremo siano date le rimanenti   condizioni 

, , , … , , 0       1, … ,  

 queste equazioni linearizzate non saranno in generale verificate, percio’ anche in questo  analogamente a sopra avremo: 

,

zisolito che nel primo estremo siano date   condizioni al contorno nell forma

… , , 0       1, … ,  

Dati i valori di prova della soluzione anche caso,

, , , … , , ∆      1, … ,  

Page 36: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

35  

Percio’ anche in questo caso, analogamente a quanto fatto sopra, dovremo imporre che: 

∆      1, … ,  

tano: 

,

, , , … , , ∆      1, … ,  

∆      1, … ,  

che scritte esplicitamente diven

,,

∆                 1, … ,  

,,

,

∆                 1, … ,  

 

che  sono  appunto  le    che mancano per portare  il  sistema  a   equazioni 

nelle   incognite  , , … , , , 1, … ,  e che rendono percio’ il sistema 

in linea di principio risolubile. 

i mesh  , de  dalle variabili  in  , 1 ,  la matrice  dei  coefficienti  di  questo  grande  sistema  sara’ evidentemente  una matrice  con  una  ba da  sulla  diagonale  principale  (vedi  figura 

i sistemi di equazioni, soprattutto attraverso tecniche specifiche per  la risoluzione di matrici sparse. 

 

 

 

Poiche’  ogni  equazione  per  la  generica  coppia  d 1   dipen

n10).  La  soluzione  di  questo  sistema,  quindi,  potra’  essere  effettuata  oltre  che attraverso  le tecniche standard (gauss‐jordan, decomposizione LU) di risoluzione d

 

 

Page 37: APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2 Ü æ Ü Ù æ Ù @ U : T ; @ T contorno che determina il metodo numerico di integrazione

 

36  

te  , cioe’ le correzioni da apportare ai valori di prova di tutte le variabili in tutti i mesh, mentre  le “B”  indicano  i  termini noti, cioe’  le “non” verifiche delle varie equazioni. Gli spazi vuoti ppresentano gli zeri. 

ome di consueto, se le equazioni sono non lineari, questa procedura dovra’ essere   fino  a  raggiungere  l’azzeramento  delle  “non”  verifiche  e  percio’  la 

 verso la soluzione. 

 

Figura 10: Struttura della matrice dei coefficienti del sistema di   equazioni differenziali trasormate in equazioni alle differenze finite, con condizioni al contorno ai due estremi dell’intervallo di integrazione. Le “X” rappresentano  le derivate della generica equazione   (righe) rispetto alla varibile   (colonne). Le “V” rappresentano le incogni

ra

 

Citerataconvergenza