Método de Crank-Nicholson
% Programa que permite resolver la transferencia de calor en estado
transitorio mediante el Método de Crank-Nicholson
clear;
% Parámetros para definir la ecuación de calor y el rango en espacio y
tiempo
L = 1.; % Longitud del alambre
T =1.; % Tiempo final
% Parámetros necesarios para resolver la ecuación dentro del método de
Crank-Nicholson
maxk = 2500; % Pasos de tiempo
dt = T/maxk;
n = 50.; % Número de pasos de espacio
dx = L/n;
cond = 1/2; % Conductividad
b = cond*dt/(dx*dx); % Parametro del metodo
% Temperatura inicial del alambre: un seno.
for i = 1:n+1
x(i) =(i-1)*dx;
u(i,1) =sin(pi*x(i));
end
% Temperatura en el límite (T = 0)
for k=1:maxk+1
u(1,k) = 0.;
u(n+1,k) = 0.;
time(k) = (k-1)*dt;
end
% Definiendo las matrices M_right y M_left en el método.
aal(1:n-2)=-b;
bbl(1:n-1)=2.+2.*b;
ccl(1:n-2)=-b;
MMl=diag(bbl,0)+diag(aal,-1)+diag(ccl,1);
aar(1:n-2)=b;
bbr(1:n-1)=2.-2.*b;
ccr(1:n-2)=b;
MMr=diag(bbr,0)+diag(aar,-1)+diag(ccr,1);
% Implementacion del metodo de Crank-Nicholson
for k=2:maxk % bucle de tiempo
uu=u(2:n,k-1);
u(2:n,k)=inv(MMl)*MMr*uu;
end
% Representación gráfica de la temperatura en diferentes tiempos
seleccionados.
figure(1)
plot(x,u(:,1),'-',x,u(:,100),'-',x,u(:,300),'-',x,u(:,600),'-')
title('Temperature within the Crank-Nicholson method')
xlabel('X')
ylabel('T')
figure(2)
mesh(x,time,u')
title('Temperature within the Crank-Nicholson method')
xlabel('X')
ylabel('Temperature')
function x=gseidel(A,b,tol,imax)
n=max(size(A)); % Obtenemos de la matriz A
% el número de ecuaiones
% y por ende de incognitas.
c=diag(A); % Almacenamos en 'c' los elementos
% de la diagonal principal.
d=1;
i=1;
while (d~=0)&&(i<n) % revisamos si en la diagonal principal 'c'
d=c(i,1)*c(i+1,1); % existe algún elemento igual a 0.
i=i+1; % En caso de que lo haya
end % salir del bucle y d=0.
while d==0 % Si se ha detectado algún cero en la
B(1,:)=A(n,:); % diagonal principal se entra en este bucle
U(1,1)=b(n,1); % para así cambiar las filas. Se alternarán
for i=2:n % las filas hasta que todos los elementos
B(i,:)=A(i-1,:);% de la diagonal principal
U(i,1)=b(i-1,1);% sean distintos de cero.
end % Todo esto para evitar una división
A=B; % por cero. Aunque se recomienda
b=U; % siempre ordenar las ecuaciones apriori.
c=diag(A);
d=1;i=1;
while (d~=0)&&(i<n)
d=c(i,1)*c(i+1,1);
i=i+1;
end
end
% Se selecciona un vector X(n,1) sólo para evaluar
% el error en la primera iteración.
for i=1:n
X(i,1)=b(i,1)/A(i,i);
end
% se especifica el vector x(n,1) con valores iguales a cero
% para iniciar las iteraciones.
x=zeros(n,1);
err=tol+1;
iter=1;
% Se inician las iteraciones
while (tol<err)&&(iter<imax)
for i=1:n
if i==1
sum=0;
for j=i+1:n
sum=A(i,j)*x(j,1)+sum;
end
x(i,1)=(b(i,1)-sum)/A(i,i);
end
if i==n
sum=0;
for j=i-1:-1:1
sum=A(i,j)*x(j,1)+sum;
end
x(i,1)=(b(i,1)-sum)/A(i,i);
end
if 2<=i<=n-1
sum1=0;
for j=i+1:n
sum1=A(i,j)*x(j,1)+sum1;
end
sum2=0;
for j=i-1:-1:1
sum2=A(i,j)*x(j,1)+sum2;
end
x(i,1)=(b(i,1)-sum1-sum2)/A(i,i);
end
end
err=abs((x-X)./x);
err=max(err);
X=x;
iter=iter+1;
end
Solución :
Las temperaturas correspondientes son
EJEMPLO DE ECUACION DEL CALOR BIDIMIENSIONAL
Encontrar la solución de transmisión de calor de una placa cuadrada, cuyas dimensiones
son 1m por cada lado, está sometida en sus bordes a una temperatura de 0°, excepto en el
tramo CB donde está sometida a una temperatura
% Programa para graficar la solución del ejercicio de ecuación bidimensional del calor de
una placa Cuadrada
x=0:0.01:1;
y=x;
[X,Y]=meshgrid(x,y);
T=(1/sinh(pi))*sin(pi*X).*sinh(pi*Y);
surf(X,Y,T)
surf(X,Y,T);view(2)
surf(X,Y,T);view(2);colorbar
title('Solucion de la transferencia de Calor en una Placa vista como un Plano')
xlabel('coordenada x')
ylabel('coordenada y')
