%{
Metodo de diferencias Finitas una Dimension para resolver el problema
con condiones Dirichlet
p(x)Uxx + q(x)Ux + r(x)U = f(x), con solucion analitica s(x)
asi, para el problema particular
-Uxx + Vux = 0 en 0 < x < 1 con U(0) = 0 y U(1) = 1
Escribir en consola
v=@(x) 100.0; % velocidad, posibles valore 1,25,100, etc
p=@(x) -1.0;
q=@(x) v(x);
r=@(x) 0;
f=@(x) 0;
s=@(x) (exp(v(x)*x)-1.0)/(exp(v(x))-1.0);
[error,A,b,u,x,V] = fdm1d_DD(p,q,r,f,0,1,0,1,11,1,s,1);
asi, para el problema particular
Uxx = -pi*pi*cos(pi*x) en 0 < x < 1 con U(0)=1 y U(1)=-1
Escribir en consola
p=@(x) 1;
q=@(x) 0;
r=@(x) 0;
f=@(x) -pi*pi*cos(pi*x);
s=@(x) cos(pi*x);
[error,A,b,u,x,V] = fdm1d_DD(p,q,r,f,0,1,1,-1,40,1,s,1);
asi, para el problema particular
Uxx = -pi*pi*cos(pi*x) en -1 < x < 2 con U(-1)=-1 y U(2)=1
Escribir en consola
p=@(x) 1;
q=@(x) 0;
r=@(x) 0;
f=@(x) -pi*pi*cos(pi*x);
s=@(x) cos(pi*x);
[error,A,b,u,x,V] = fdm1d_DD(p,q,r,f,-1,2,-1,1,11,1,s,1);
%}
%% [error,A,b,u,x,V] = fdm1d_DD(p,q,r,f,s,xi,xf,vi,vf,n,grf,s,ssw)
%{
Definicion del problema
a(x)Uxx + b(x)Ux + c(x)U = f(x)
p=@(x) ; Coeficiente para Uxx
q=@(x) ; Coeficiente para Uxx
r=@(x) ; Coeficiente para Uxx
f=@(x) ; Lado derecho de la ecuacion
s=@(x) ; Solucion analitica
a Inicio del dominio
b Fin del dominio
vi Valor de la condicion inicial
vf Valor de la condicion final
n Numero de nodos en la particion
grf Muestra las graficas con (1), otro valor no las muesta
sws Se proporciona solucion analitica o no
%}
function [error,A,b,u,x,V] = fdm1d_DD(p,q,r,f,xi,xf,vi,vf,n,grf,s,sws)
if n < 3
return
end
% Numero de incognitas del problema
tm = n -2;
% Declaracion de la matriz y vectores de trabajo
%A = sparse(tm,tm);
A = zeros(tm,tm); % Matriz de carga
b = zeros(tm,1); % Vector de carga
u = zeros(tm,1); % Vector de solucion
x = zeros(n,1); % Vector de coordenadas de la particion
V = zeros(n,1) ; % Vector solucion
h = (xf-xi)/(n-1);
h1 = h*h;
% Llenado de los puntos de la malla
for i = 1: n,
x(i) = xi + (i-1)*h;
end
% Llenado de la matriz y vector
b(1) = f(xi) - p(xi)*(vi/h1);
A(1,2) = p(x(1))/h1 - q(x(1))/(2.0*h);
A(1,1) = ((-2.0 * p(x(1))) / h1) + r(x(1));
for i=2:tm-1,
A(i,i-1) = p(x(i))/h1 - q(x(i))/(2.0*h);
A(i,i) = ((-2.0 * p(x(i))) / h1) + r(x(i));
A(i,i+1) = p(x(i))/h1 + q(x(i))/(2.0*h);
b(i) = f(x(i));
end
A(tm,tm-1) = p(x(tm))/h1 - q(x(tm))/(2.0*h);
A(tm,tm) = ((-2.0 * p(x(tm))) / h1) + r(x(tm));
b(tm) = f(x(tm+1))-p(x(tm+1))*(vf/h1);
% Soluciona el sistema
u = inv(A)*b;
% Copia la solucion obtenida del sistema lineal al vector solucion
V(1) = vi;
for i=1:tm,
V(i+1) = u(i);
end
V(n) = vf;
% Encuentra el error en norma infinita usando una particion de PAR = 10
error = 0;
if sws == 1
par = 10;
for i = 1: n-1,
inc = (x(i+1)-x(i))/par;
for j = 1:par+1,
px = x(i)+inc*(j-1);
e = abs(s(px)-l(px,x(i),V(i),x(i+1),V(i+1)));
if e > error
error = e;
end
end
end
end
% Revisa si se graficara
if grf == 1
if sws == 1
% Calcula la solucion analitica en la particion de calculo
ua = zeros(n,1);
for i = 1: n,
ua(i) = s(x(i));
end
end
% Graficar la solucion numerica
plot(x,V,'o');
hold
% Grafica la solucion analitica en una particion de tamano xPart
if sws == 1
xPart = 1000;
h = (xf-xi)/(xPart-1);
xx = zeros(xPart,1);
xa = zeros(xPart,1);
for i = 1: xPart,
xx(i) = xi + (i-1)*h;
xa(i) = s(xx(i));
end
plot(xx,xa);
% Grafica el error
figure(2);
plot(x,V-ua);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evalua el punto x en la recta dada por los puntos (x1,y1) y (x2,y2), se usa para el calculo
% de la norma infinito
function y = l(x,x1,y1,x2,y2)
y = y1+((y2-y1)/(x2-x1))*(x-x1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%