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