APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2...
Transcript of APPUNTI DEL CORSO DI LABORATORIO DI …labcalc.fisica.uniroma1.it/docenti/limongi/public_html/...2...
APPUNTI DEL CORSO DI LABORATORIO DI CALCOLO
AVANZATO Metodi Numerici per le Equazioni Differenziali
Ordinarie
MARCO LIMONGIIstituto Nazionale di Astrofisica – Osservatorio Astronomico di Roma
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
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:
2
sostituendo questa espressione nell’ODE , per si ha:
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=+
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
1
da cui
2
Risolvendo per la derivata prima si ha:
2
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).
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.
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
7
sostituendo in questa equazione l’e pansione in serie di Taylor, si ha:
2
s
da cui:
2
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:
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
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’:
2
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:
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
2
,
,
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
o
,
2
da cui:
/ , / ,2
,
La formula di Runge, percio’, sara’:
,
2,
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
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
, , ,
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
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,
2
,
6 3 3 6
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
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 ,
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
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
, , .
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:
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
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
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
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
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
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
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
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
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.
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
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:
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
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
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
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, … ,
e
35
Percio’ anche in questo caso, analogamente a quanto fatto sopra, dovremo imporre che:
∆ 1, … ,
tano:
,
, , , … , , ∆ 1, … ,
e
∆ 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
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