19/07/2019, 16:09
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
19/07/2019, 17:04
19/07/2019, 19:53
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)')
19/07/2019, 23:24
20/07/2019, 11:14
22/07/2019, 12:52
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
22/07/2019, 14:39
22/07/2019, 15:50
22/07/2019, 19:19
%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
%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
Skuola.net News è una testata giornalistica iscritta al Registro degli Operatori della Comunicazione.
Registrazione: n° 20792 del 23/12/2010.
©2000—
Skuola Network s.r.l. Tutti i diritti riservati. — P.I. 10404470014.
Powered by phpBB © phpBB Group - Privacy policy - Cookie privacy
phpBB Mobile / SEO by Artodia.