dif_finitas_tridi.m

% #######################  PARTICULA CATALITICA ESFERICA ###################
%
%
% TRABALHO DE MODELAGEM E SIMULAO - 2006
% 
% 
% Problema:
% - partcula cataltica esfrica (s=2 )
% - reao-difuso em estado estacionrio
% 
% Soluo via Diferenas Finitas usando o solve tridi.dll
%
% -> Exata 1 ordem
%
% Giovani Tonel ([e-mail]) - Fevereiro de 2007


clear all
clc





disp(' ');
disp('========================================================================================================');
disp(' 	                         PARTICULA CATALITICA ESFERICA (estacionario e m=1)');
disp('========================================================================================================');
disp(' ');





% ######################################
op=menu('Escolha o numero de pontos (n):', ...
   'n= 10 pontos',...
   'n= 20 pontos',...
   'n= 25 pontos',...
   'n= 50 pontos');

switch op
    case 1
       n=10;
    case 2
       n=20;
    case 3
       n=25;   
    case 4
       n=50;
    end

% ######################################

run_sim	 	= 1;
phi			=	4;
N				=	n+2;
x0				=	0;
xf				=	1;
h				=	(xf-x0)/(n+1);
x				=	[x0:h:xf]';

ordem		=	1;




x(1)		=	1e-5;
xi= [x(1) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0];
xexat	= 	x;
yexat	= 	sinh(phi*xexat)./(xexat*sinh(phi));
y0			= 	yexat;
[tam]	=	size(xexat);



% ########################## Construindo a Matriz Tridiagonal ############################

% x		= 	0.1:0.1:0.9;
x=x(2:n+1);


for i=1:n,
   if ordem==1,
      me(i)	=		(x(i)/h -1);
      mc(i)	= -	(2*x(i)/h + phi^2.*x(i)*h);
      md(i)	=		(x(i)/h + 1);
   end
end


a	=	[ me 0];
b	=	[1 mc 1];
c	=	[-1 md ];

d		=	[zeros(1,n+1) 1];



% Tridiagonal matrix.
% T = tridiag(a, b, c, n) returns an n by n matrix that has 
% a, b, c as the inferior diagonal, main diagonal, and superior diagonal 
% entries in the matrix.
% M = diag(a,-1) + diag(b,0) + diag(c,1)
% M*y=d

%inf 	= 	diag(ones(n,1),-1); 		% inferior diagonal matrix
%sup	=	diag(ones(n,1), 1);		% superior diagonal matrix

% M 		= a* inf + b*eye(N)  + c* sup;

ysol	= tridi(a, b, c, d);

ysol	=	ysol';

yint = (xexat).^2.*(ysol).^ordem; 		% x^s*[y(x)]^m


op 	= 	3; % ordem do polinmio

inty	= quad8('func',x(1), 1,[],[], xexat,yint,op);


etai	 	= 	(3)*inty
etad 	=	(3/phi^2)*(ysol(end)-ysol(end-1))/(xexat(end,1)-xexat(end-1,1))



ysol2		=	interp1(xexat ,ysol, xi');
 
sol_inter		=	[xi' ysol2]	;
 
 
figure(1)
plot(xexat, yexat,'r' ,xexat, ysol,'-sb');
xlabel('Valor de x')
ylabel('Valor de y')
title('Anlise Comparativa - Diferenas Finitas (m=1)');
legend('Soluo Analtica', 'Soluo via tridi.dll')
axis([0 1 0 1])
   

% ##########################################################################################