Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda amatrix » 19/07/2019, 17:09

Immagine
Sto cercando di rappresentare in Matlab in un grafico loglog che per i vari metodi LMF (https://en.wikipedia.org/wiki/Linear_multistep_method) risolutivi di equazioni differenziali ordinarie l'errore $ e_n=y(t_n) -y_n $ (ovvero la differenza tra la soluzione reale e quella approssimata) tende a zero con ordini di convergenza diversa in base all'ordine del metodo.
Per fare ciò ho considerato l'equazione test $y'(t)=-y(t), \quad y(0)=1$ con soluzione esatta ugaule a $y(t)=e^{-t}$ e applicando il seguente algoritmo mi aspettavo che le rette che mi rappresentano l'errore di ciascun metodo avessero pendenze diverse in funzione dell'ordine del metodo (ho considerato metodi di ordine 1,2,3 e 4).

Vi allego sia l'algoritmo che l'immagine che ne esce fuori. Qualcuno sa spiegarmi perchè sembrerebbe che tutti i metodi abbiano lo stesso ordine di convergenza?

Codice:
function []=OrdineErrore(t0,T,nmin,nmax)
l=-1;
y0=1;
n=linspace(nmin,nmax,50);
a=length(n);
for i=1:a
    h(a+1-i)=(T-t0)/n(i);
    q(i)=h(i)*l;
end
for i=1:a
%Costruzione soluzione esatta
    t=linspace(t0,T,n(i));
    solt=y0*exp(l*t);
    for j=2:n(i)
%Costruzione Eulero Esplicito (ordine 1, una condizione iniziale)
        y(1)=y0;
        y(j)=y(j-1)/(1-q(i));
%Costruzione trapezi (ordine 2, una condizione iniziale)
        p(1)=y0;
        p(j)=((1+q(i)/2)/(1-q(i)/2))*p(j-1);
%Costruzione Eulero Implicito (ordine 1, una condizione iniziale)
        b(1)=y0;
        b(j)=(1/(1-q(i)))*b(j-1);     
    end
    for r=3:n(i)
%Costruzione BashFort (ordine 2, 2 condizioni iniziali)
        v(1)=y0;
        v(2)=solt(2);
        v(r)=v(r-1)+q(i)*(3*v(r-1)-v(r-2))/2;
%Costruzione BDF2 (ordine 2, 2 condizioni iniziali)
        z(1)=y0;
        z(2)=solt(2);
        z(r)=((4/3)*z(r-1)-(1/3)*z(r-2))/(1-(2/3)*q(i));
%Costruzione AM-3 (ordine 3, 2 condizioni iniziali)
        w(1)=y0;
        w(2)=solt(2);
        w(r)=w(r-1)*((1+8*q(i)/12)/(1-5*q(i)/12))-w(r-2)*((q(i)/12)/(1-5*q(i)/12));
%Costruzione Simpson (ordine 4, 2 condizioni iniziali)
        d(1)=y0;
        d(2)=solt(2);
        d(r)=(((4/3)*q(i))/(1-q(i)/3))*d(r-1)+((1+q(i)/3)/(1-q(i)/3))*d(r-2);
%Costruzione MidPoint (ordine 2, 2 condizioni iniziali)
        mi(1)=y0;
        mi(2)=solt(2);
        mi(r)= 2*q(i)*mi(r-2)+y(r-1);
    end
    for k=4:n(i)
%Costruzione AM-4 (ordine 4, 3 condizioni iniziali)
        x(1)=y0;
        x(2)=solt(2);
        x(3)=solt(3);
        x(k)=x(k-1)*(1+19*q(i)/24)/(1-q(i)*9/24)-x(k-2)*(q(i)*5/24)/(1-q(i)*9/24)+x(k-3)*(q(i)/24)/(1-q(i)*9/24);
    end
end
%errori
for i=1:a
    eE(i)=abs(y(i)-solt(i)); %Eulero Esplicito (1)
    eI(i)=abs(b(i)-solt(i)); %Eulero Implicito (1)
    eV(i)=abs(v(i)-solt(i)); %Adams-Bashfort 2
    eT(i)=abs(z(i)-solt(i)); %BDF 2
    eP(i)=abs(p(i)-solt(i)); %trapezi (2)
    eA(i)=abs(w(i)-solt(i)); %Adams-Moulton 3
    eX(i)=abs(x(i)-solt(i)); %Adams-Moulton 4
    eD(i)=abs(d(i)-solt(i)); %Simpson (4)
    eM(i)=abs(mi(i)-solt(i)); %Midpoint (2)
   
end
%grafici
hold on

loglog(eE,'b') %EE
loglog(eI,'m') %EI
loglog(eP,'c') %trapezi
loglog(eT,'b') %BDF 2
loglog(eV,'g') %AB 2
loglog(eX,'k') %AM 4
loglog(eA,'r') %AM 3
loglog(eD,'--') %Simpson (4)
loglog(eM,'c') %MidPoint (2)
end
amatrix
Starting Member
Starting Member
 
Messaggio: 2 di 6
Iscritto il: 11/07/2019, 16:02

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda gugo82 » 19/07/2019, 18:04

Parlo da profano, però è possibile che non ti sia messo in un caso “brutto” per nessuno dei metodi utilizzati.

Prova a cercare qualcosa sul tuo testo (o in rete) circa i casi peggiori.
Sono sempre stato, e mi ritengo ancora un dilettante. Cioè una persona che si diletta, che cerca sempre di provare piacere e di regalare il piacere agli altri, che scopre ogni volta quello che fa come se fosse la prima volta. (Freak Antoni)
Avatar utente
gugo82
Moderatore globale
Moderatore globale
 
Messaggio: 22021 di 22141
Iscritto il: 13/10/2007, 00:58
Località: Napoli

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda arnett » 19/07/2019, 19:03

Il codice è un disastro. Ti è chiaro cos'è l'ordine di convergenza di un metodo?
Tu stai raffigurando l'errore punto per punto in funzione dell'indice del nodo spaziale (!), una roba che non c'entra nulla, e solo dell'ultima iterazione che fai (cioè della soluzione approssimata con griglia spaziale più fitta), perché anche se risolvi il problema a volte, salvi solo l'ultima. Tu devi raffigurare l'errore in funzione di $h$.
"ci scruta poi gira se ne va"
arnett
Senior Member
Senior Member
 
Messaggio: 1001 di 1036
Iscritto il: 18/07/2018, 09:08

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda feddy » 19/07/2019, 20:53

gugo ha ragione, per alcuni dei metodi espliciti che hai proposto hai delle restrizioni sul passo (tipo per Eulero esplicito). Assicurati prima di fare i dovuti test che siano rispettati.

Un esempio pratico con Eulero implicito dovrebbe chiarirti definitivamente le idee.
Codice:
clear all
close all
t0=0;
y0=1;
tf=1;
tsrange=[10:10:100];
count=0;
for ts=tsrange
   count=count+1;
   t=linspace(t0,tf,ts+1);
   k=(tf-t0)/ts;
   uex=exp(-t);
   y=y0;
   for n=1:ts
      y=y/(1+k);
   end
   
   err(count)=abs(y-uex(end));
end
figure
loglog(tsrange,err,'*',tsrange,err(end)*(tsrange/tsrange(end)).^(-1),'-')
xlabel('time steps')
ylabel('error wrt analytical solution')
legend('err','O(k)')
Avatar utente
feddy
Advanced Member
Advanced Member
 
Messaggio: 2540 di 2554
Iscritto il: 26/06/2016, 01:25
Località: Italia

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda arnett » 19/07/2019, 21:45

@feddy: calcoli l'errore direttamente nell'ultimo nodo perché sai a priori che l'errore è crescente o perché l'ordine può essere verificato in un nodo qualsiasi fissato? Io sono stato abituato in questo contesto ad usare $\max |u_h -u_{ex}|$.

Edit: no, ho fatto qualche prova e un nodo fissato qualsiasi non basta, mi dava stime fasulle di tipo $p=1.5$.
Ultima modifica di arnett il 20/07/2019, 09:22, modificato 1 volta in totale.
"ci scruta poi gira se ne va"
arnett
Senior Member
Senior Member
 
Messaggio: 1002 di 1036
Iscritto il: 18/07/2018, 09:08

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda amatrix » 20/07/2019, 00:24

Grazie arnett, ho capito l'errore!

E grazie feddy, sto prendendo spunto dal tuo algoritmo.
Ho notato che se utilizzo "hold on" per graficare più metodi insieme mi scompare la scala logaritmica sugli assi. Sapete come ovviare a questo problema?
amatrix
Starting Member
Starting Member
 
Messaggio: 3 di 6
Iscritto il: 11/07/2019, 16:02

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda feddy » 20/07/2019, 12:14

@arnett sì, un nodo qualsiasi non va bene!

La definizione di errore che ho in testa io è che se risolvo un sistema differenziale fino al tempo $\bar{t}$ con due discretizzazioni temporali a passo costante $k_1=\frac{\bar{t}}{m_1}, k_2=\frac{\bar{t}}{m_2}$, avrò come errori finali (supponendo le costanti "uguali") \[ \boldsymbol{e}_{m_1,k_1}=C k_1^p\] e \[ \boldsymbol{e}_{m_2,k_2}=C k_2^p\]

Dunque \[ \frac{\boldsymbol{e}_{m_1,k_1}}{\boldsymbol{e}_{m_2,k_2}}=(\frac{k_2}{k_1})^p \]

da cui \[ \log \boldsymbol{e}_{m_2,k_2} - \log \boldsymbol{e}_{m_1,k_1}= -p(\log m_2 - \log m_1) \]
Avatar utente
feddy
Advanced Member
Advanced Member
 
Messaggio: 2541 di 2554
Iscritto il: 26/06/2016, 01:25
Località: Italia

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda amatrix » 22/07/2019, 13:52

Ho utilizzato la bozza di algoritmo di feddy per sviluppare il seguente algoritmo.
(Ad ogni passo ho usato l'ultimo valore del vettore errore perché, in questo caso, il vettore errore è monotono)

Codice:
function []=ErrOrdini
t0=0;
y0=1;
tf=1;
tsrange=[10:1000:100010];
count=0;
for ts=tsrange
   count=count+1;
   t=linspace(t0,tf,ts);
   h=(tf-t0)/ts;
   uex=exp(-t);
   y=y0;
   z=y0;
   for n=1:ts
       y=y/(1+h); %Eulero Implicito
       z=z*(1-h/2)/(1+h/2) ; %trapezi
   end
   w(1)=y0;
   w(2)=uex(2);
   w(3)=uex(3);
   for i=4:ts %Adams Bashfort ordine 3
       w(i)=w(i-1)*(1-23*h/12)+w(i-2)*(16*h/12)-w(i-3)*(5*h/12);
   end
   er1(count)=abs(y-uex(end));
   er2(count)=abs(z-uex(end))
   er3(count)=abs(w(end)-uex(end));
end
loglog(tsrange,er1,'b');
loglog(tsrange,er2,'r');
loglog(tsrange,er3,'g');
end



La soluzione che ottengo mi soddisfa solo in parte:
infatti per i metodi di Eulero Implicito e Trapezi, detto h il passo, ho rispettivamente convergenza di ordine O(h) e O(h^2) come si può notare in figura (rispettivamente blu e rosso).
Tuttavia per quanto riguarda Adams Bashfort di ordine 3 (e ho provato anche con altri metodi) non riesco a raggiungere una convergenza di ordine O(h^3).
Dove sbaglio? Cosa posso usare per avere una convergenza di ordine maggiore?

Immagine
amatrix
Starting Member
Starting Member
 
Messaggio: 4 di 6
Iscritto il: 11/07/2019, 16:02

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda feddy » 22/07/2019, 15:39

Il for di AB non è scritto giusto... fatti due iterate a mano su un foglio e scova l'errore
Avatar utente
feddy
Advanced Member
Advanced Member
 
Messaggio: 2543 di 2554
Iscritto il: 26/06/2016, 01:25
Località: Italia

Re: Ordine di convergenza di metodi risolutivi di ODE a *k* passi

Messaggioda amatrix » 22/07/2019, 16:50

Considerando che per il metodo di AB3 abbiamo (https://en.wikipedia.org/wiki/Linear_multistep_method)

$y_{n+3}-y_{n+2} = h ( \frac{23}{12} f_{n+2} - \frac{16}{12} f_{n+1} + \frac{5}{12} f_n) $

e che nel nostro caso $ f_n = - y_n $ ho rifatto i conti su carta e mi sembrano corretti :?
amatrix
Starting Member
Starting Member
 
Messaggio: 5 di 6
Iscritto il: 11/07/2019, 16:02

Prossimo

Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite

cron