03/06/2018, 21:29
04/06/2018, 10:59
D(1,1:2)=[1 - (2/h^2),2/h^2]
b(1)=2/h
A(m,m-1:m)=[0,0]; B(m,m-1:m)=[0,0]
b(m)=1;
clc
clear all
close all
mrif=2^10+1;
m=mrif;
x=linspace(0,pi/2,m)';
h=(pi/2)/ (m-1);
%discretizzazione derivate
A = toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
B = toeplitz(sparse(1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m));
I=speye(m);
%% CONDIZIONI AL BORDO
%DIRICHLET al nodo a destra.
A(m,m-1:m)=[0,0];
B(m,m-1)=0; % AZZERO PER FAR SOPRAVVIVERE SOLO LA MATRICE IDENTITÀ NELL'ULTIMA RIGA
D=A+B+I;
% NEUMANN, risolte con nodo fantasma
D(1,1:2) = [1-2/h^2,2/h^2];
b=cos(x); %termine noto
b(1)=2/h;
b(m)=1;
F=@(u) D*u - b; %funzione da azzerare
u=D\b;
figure
plot(x,u,'b',x,sin(x),'r-o')
axis([0, pi/2,0,1])
legend('sol numerica','sol analtica')
%% Costruzione soluzione di riferimento
Urif=u;
mrange=2.^(4:8)+1;
count=0;
for m=mrange
count=count+1;
x=linspace(0,pi/2,m)';
h=(pi/2)/ (m-1);
A = toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
B = toeplitz(sparse(1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m));
I=speye(m);
%CONDIZIONI AL BORDO
A(m,m-1:m)=[0,0];
B(m,m-1)=0; % AZZERO PER FAR SOPRAVVIVERE SOLO LA MATRICE IDENTITÀ NELL'ULTIMA RIGA
D=A+B+I;
D(1,1:2) = [1-2/h^2,2/h^2];
b=cos(x);
b(1)=2/h;
b(m)=1;
u=D\b;
err(count) = norm(u-Urif(1:(mrif-1)/(m-1):mrif),inf)
end
figure
loglog(mrange,err,'o',mrange,mrange.^-2*(err(end)/mrange(end)^-2))
05/06/2018, 17:18
05/06/2018, 17:33
05/06/2018, 18:04
s-o-p-h-i ha scritto:non ho capito il passaggio per cui b=1.
s-o-p-h-i ha scritto:Poi avrei un altro problema da proporti
b=ones(m,1);
b(1)=0;b(m)=0;
A(1,1:2)=[1,0]; B(1,2)=0;
M=-diag(1+x)*D2-D1;
u=M\b
count ++
count = count +1
clear all
close all
mrange=[10,20,100,200];
count=0;
for m=mrange+1
count++;
h=1/(m-1);
x=linspace(0,1,m)';
A = toeplitz(sparse([1,2],[1,1],[-2,1]/(h^2),m,1));
B = toeplitz(sparse(2,1,-1/(2*h),m,1), sparse(2,1,1/(2*h),m,1));
%boundary cond.
A(1,1:2)=[1,0];
A(m,m-1:m)=[0,1];
B(1,1:2)=[0,0];
B(m,m-1:m)=[0,0];
b=ones(m,1);
b(1)=0;
b(m)=b(1);
M=-diag(1+x)*A-B; %matrice del problema
u=M\b;
s=@(x) -x + (log(1 + x))/(log(2)); %sol analitica
err(count)=norm(u-s(x),inf);
end
figure(1)
plot(x,s(x),'k',x,u,'o')
axis([0,1,0,0.09])
legend('sol analitica','sol numerica')
figure(2)
loglog(mrange,err,'*',mrange,(err(end))*(mrange/mrange(end)).^-2)
legend('errore','ordine 2')
06/06/2018, 13:24
06/06/2018, 18:23
s-o-p-h-i ha scritto:Non ho capito come hai fatto a trovare la soluzione analitica,che formula ci sta dietro?
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.