Pagina 1 di 2

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

MessaggioInviato: 19/07/2019, 17:09
da amatrix
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

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

MessaggioInviato: 19/07/2019, 18:04
da gugo82
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.

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

MessaggioInviato: 19/07/2019, 19:03
da arnett
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$.

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

MessaggioInviato: 19/07/2019, 20:53
da feddy
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)')

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

MessaggioInviato: 19/07/2019, 21:45
da arnett
@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$.

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

MessaggioInviato: 20/07/2019, 00:24
da amatrix
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?

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

MessaggioInviato: 20/07/2019, 12:14
da feddy
@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) \]

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

MessaggioInviato: 22/07/2019, 13:52
da amatrix
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

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

MessaggioInviato: 22/07/2019, 15:39
da feddy
Il for di AB non è scritto giusto... fatti due iterate a mano su un foglio e scova l'errore

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

MessaggioInviato: 22/07/2019, 16:50
da amatrix
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 :?