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

[amatrix1]

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?

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 

[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.

[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.

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)')

[amatrix1]

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?

[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) \]

[amatrix1]

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)

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?

[feddy]

Il for di AB non è scritto giusto... fatti due iterate a mano su un foglio e scova l'errore

[amatrix1]

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

[amatrix1]

Ho provato a inserire anche i seguenti metodi di ordine 4 (AM e AB) e ricontrollato svariate volte i coefficienti ma dal grafico convergono come O(h)...

%AB4
for i=5:ts
u(1)=y0;
u(2)=uex(1);
u(3)=uex(2);
u(4)=uex(3);
u(i)=u(i-1)*(1-h*55/24)+u(i-2)*(h*59/24)-u(i-3)*h*37/24+u(i-4)*h*9/24 ;
end

e

%AM3
for i=4:ts
b(1)=y0;
b(2)=uex(1);
b(3)=uex(2);
b(i)=b(i-1)*(1-19*h/24)/(1+h*9/24)-b(i-2)*(-h*5/24)/(1+h*9/24)+b(i-3)*(-h/24)/(1+h*9/24);
end

E' possibile che definendoli attraverso il ciclo for incombo automaticamente in degli errori o che magari questi dipendano dalle condizioni iniziali benché stia inserendo come dati iniziali direttamente i valori della funzione?