Diferencia entre revisiones de «Ecuación de Calor: "Calentamiento Varilla"»

De MateWiki
Saltar a: navegación, buscar
(→‎Variación de la temperatura para tiempos grandes)
 
(No se muestran 13 ediciones intermedias del mismo usuario)
Línea 169: Línea 169:
 
[[Archivo:EULEREXPL_SURF.jpg|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma ]]
 
[[Archivo:EULEREXPL_SURF.jpg|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma ]]
  
−
EULER IMPLÍCITO
+
A continuación resolvemos el problema por el método de Euler implícito:
  
 
{{matlab|codigo=
 
{{matlab|codigo=
Línea 235: Línea 235:
  
 
En este caso, al estar estudiando el problema a través de un método implícito no tenemos problemas de inestabilidad, por lo que no ha sido necesario modificar el tamaño de paso.
 
En este caso, al estar estudiando el problema a través de un método implícito no tenemos problemas de inestabilidad, por lo que no ha sido necesario modificar el tamaño de paso.
−
RUNGE KUTTA
+
 
 +
Por último, empleamos el método de Runge-Kutta para resolver el problema.
  
 
{{matlab|codigo=
 
{{matlab|codigo=
Línea 245: Línea 246:
 
tM=10;
 
tM=10;
 
q=2;
 
q=2;
 +
%-------------------------------------------------
 
%Datos de entrada
 
%Datos de entrada
 
%ca: Funcion f(t) de la condicion de contorno en x=a (t vector)
 
%ca: Funcion f(t) de la condicion de contorno en x=a (t vector)
Línea 254: Línea 256:
 
%g: funcion término independiente (x será vector, t escalar)
 
%g: funcion término independiente (x será vector, t escalar)
 
g=@(x,t)0*x;
 
g=@(x,t)0*x;
 +
%--------------------------------------------------------------------------
 
Discretización del intervalo [a,b]
 
Discretización del intervalo [a,b]
 
h=0.1;
 
h=0.1;
Línea 260: Línea 263:
 
x=x';
 
x=x';
 
xx=x(2:end-1);
 
xx=x(2:end-1);
−
 
+
%-----------------------------------------------------------------------
 
%Matriz K (de N-1 x N-1)
 
%Matriz K (de N-1 x N-1)
 
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
 
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
−
   
+
  %-------------------------------------------------
 
%valor inicial  
 
%valor inicial  
 
U0=u0(xx);
 
U0=u0(xx);
−
   
+
  %-------------------------------------------------
 
%discretización del tiempo
 
%discretización del tiempo
 
k=0.001;%paso en t
 
k=0.001;%paso en t
 
M=round((tM-t0)/k);
 
M=round((tM-t0)/k);
 
t=linspace(t0,tM,M+1);
 
t=linspace(t0,tM,M+1);
−
+
%-------------------------------------------------
 
sol=zeros(N-1,M+1);
 
sol=zeros(N-1,M+1);
 
sol(:,1)=U0;
 
sol(:,1)=U0;
Línea 297: Línea 300:
 
    
 
    
 
end  
 
end  
−
   
+
  %---------------------------------------------------------
 
%añadimos las condiciones Dirichlet
 
%añadimos las condiciones Dirichlet
 
Ua=ca(t);
 
Ua=ca(t);
 
Ub=cb(t);
 
Ub=cb(t);
 
sol=[Ua;sol;Ub]
 
sol=[Ua;sol;Ub]
−
 
+
%-------------------------------------------------
 
%gráfico 3D
 
%gráfico 3D
 
[Mt,Mx]=meshgrid(t,x);
 
[Mt,Mx]=meshgrid(t,x);
Línea 308: Línea 311:
 
xlabel('t');
 
xlabel('t');
 
ylabel('x');
 
ylabel('x');
−
 
+
%-------------------------------------------------
 
%resultados
 
%resultados
 
xx=2;
 
xx=2;
Línea 327: Línea 330:
  
 
Al emplear el tamaño de paso k=0.1 hemos detectado inestabilidad en la solución, por lo que hemos tratado de solventarla utilizando un tamaño de paso menor, siendo este k=0.001. El código adjunto es este último.
 
Al emplear el tamaño de paso k=0.1 hemos detectado inestabilidad en la solución, por lo que hemos tratado de solventarla utilizando un tamaño de paso menor, siendo este k=0.001. El código adjunto es este último.
 +
  
  
 
== Variación de la temperatura para tiempos grandes ==
 
== Variación de la temperatura para tiempos grandes ==
 +
 +
A continuación comprobamos gráficamente como evoluciona la temperatura en el tiempo. Se observa que a medida que avanza el tiempo la solución tiende a parecerse a la obtenida en el caso estacionario. Para que sea fácil de apreciar, se representa la distancia entre la solución estacionaria y la solución u(x,t) para distintos tiempos.
  
 
{{matlab|codigo=
 
{{matlab|codigo=
Línea 465: Línea 471:
 
}}  
 
}}  
  
−
[[Archivo:Distancia5.jpeg|marco|centro| ]]
+
[[Archivo:Dist apart 5.JPG|marco|centro| ]]
−
 
+
  
 +
La ecuación que debería verificar el estado estacionario es una recta, la cual podemos obtener una vez que hemos despreciado el término Ut. Esta recta es \[\frac{5x}{4} +5=0\]. Como se observa en las gráficas, la función a la que debería parecerse la solución u(x,t) para tiempos grandes es precisamente la recta correspondiente al estado estacionario.
  
 +
== Extremo con pieza aislante ==
  
 +
Se coloca en el extremo izquierdo una pieza aislante que hace que no haya pérdida de calor en ese extremo, por lo que el flujo de temperatura es nulo. Esto nos cambia las condiciones de contorno del problema, pasando de tener dos condiciones Dirichlet a tener una condición Neumann y otra Dirichlet.
  
 +
{{matlab|codigo=
  
 +
%DATOS DEL PROBLEMA
 +
a=0;
 +
b=4;
 +
t0=0;
 +
tM=10;
 +
q=2;
 +
%----------------------------------------------------------------
 +
%Discretización
 +
h=0.1;
 +
N=round((b-a)/h);
 +
x=linspace(a,b,N+1);
 +
x=x';
 +
xx=x(1:end-1); %tomamos los elementos del primero al penúltimo
 +
%----------------------------------------------------------------
 +
%Funciones.Datos de entrada.
 +
%ca:Función f(t) de la condición de contorno en x=a(t vector)
 +
ca=0; %Condición Neumann
 +
%cb:Función f(t) de la condición de contorno en x=b(t vector)
 +
cb=0;%Condición Dirichlet
 +
%u0Forma analitica f(x) funcion valor inicial(x vector)
 +
u0=exp(-3*((xx-3).^2))+5+xx;
 +
%Matriz K
 +
K=q/h^2*(2*diag(ones(1,N))-diag(ones(1,N-1),-1) -diag(ones(1,N-1),1));
 +
K(1,2)=-4/h^2;
 +
%----------------------------------------------------------------
 +
%Discretización del tiempo
 +
k=0.1; %Anchura de paso en el tiempo(k=0.5*h^2)
 +
M=round((tM-t0)/k);%numero de subintervalos en t
 +
t=linspace(t0,tM,M+1);
 +
%sol matriz que tendrá la solución para cada punto e instante
 +
sol=zeros(N,M+1);
 +
sol(:,1)=u0;
 +
%----------------------------------------------------------------
 +
%Términos independiente
 +
G=0*xx;
 +
G(1)=G(1)-2*q*ca/h;
 +
G(end)=G(end)+q*cb/h^2;
  
 +
for i=1:M
 +
   
 +
    sol(:,1+i)=(eye(size(K))+k*K)\(sol(:,i)+k*G);
 +
   
 +
end
 +
%----------------------------------------------------------------
 +
%Añado frontera(Solo Dirichlet)
 +
Ub=q*cb*ones(1,length(t));
 +
sol=[sol;Ub];
 +
%grafico
 +
[Mt,Mx]=meshgrid(t,x);
 +
surf(Mt,Mx,sol);
 +
xlabel('t');
 +
ylabel('x');
 +
zlabel('Temperatura')
 +
}}
 +
[[Archivo:Aislante.jpg|marco|centro| ]]
  
  
−
FOURIER (apartado 7)
+
== Método de Fourier ==
  
 +
Partiendo de una solución analítica previa obtenemos la solución en forma de sumatorio de infinitos términos, reduciendo este número de términos numéricamente calcularemos una solución aproximada.
 +
Mediante el siguiente código de Matlab resolvemos el problema numéricamente para 1,3,5,10 y 20 términos de la serie de Fourier.
  
 
{{matlab|codigo=
 
{{matlab|codigo=
Línea 560: Línea 625:
 
[[Archivo:PLOTFOURER.jpg|marco|centro| ]]
 
[[Archivo:PLOTFOURER.jpg|marco|centro| ]]
  
−
Con este gráfico queremos representar una vez pasadas 0.5 unidades de tiempo la temperatura a lo largo de la barra. Hemos hecho la representación usando Q= 1,3,5,10 y 20 como términos de la Serie de Fourier. La aproximación más exacta la obtenemos con Q=20 y todas son proporcionales decreciendo cumpliendo así las condiciones de contorno impuestas en el enunciado u(0,t)=5 y u(4,t)=0.   
+
Con este gráfico queremos representar una vez pasadas 0.5 unidades de tiempo la temperatura a lo largo de la barra. La aproximación más exacta la obtenemos con Q=20 y todas son proporcionales decreciendo cumpliendo así las condiciones de contorno impuestas en el enunciado u(0,t)=5 y u(4,t)=0.   
  
 +
== Pérdida de calor en la varilla ==
  
 +
Suponemos ahora una pérdida de calor de la varilla a través del aire que tiene una temperatura constante de 16 grados. Con la siguiente ecuación \[\frac{\partial U}{\partial t}-\frac{\partial ^2U}{\partial x^2}+ (u-9)=0\] y manteniendo las condiciones de contorno establecidas en el primer apartado, calculamos el estado estacionario y el tiempo que se tarda en alcanzar el estado estacionario con un error de 10-³.
  
 +
{{matlab|codigo=
  
 +
%CONDICIONES DIRICHLET & DIRICHLET
 +
%-------------------------------------
 +
%concidiones iniciales
 +
a=0;%extremo izquierdo
 +
b=4;%extremo derecho
 +
t0=0;%instante inicial
 +
tM=10;%instante final
 +
q=1;%difusividad termica
 +
%-------------------------------------------------
 +
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
 +
cb=@(t)0*t;%función f(t) de la condición de contorno en x=b
 +
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
 +
%-----------------------------------------------
 +
%discretización intervalo[a,b]
 +
h=0.1;
 +
N=round((b-a)/h);%subintervalos en x
 +
x=linspace(a,b,N+1);
 +
x=x';%el vector x tiene N+1 elementos
 +
xx=x(2:end-1);%tomamos los nodos intermedios
 +
%----------------------------------------------------------
 +
%discretización del tiempo
 +
k=0.1;%tamaño de paso del tiempo
 +
M=round((tM-t0)/k);%como cuando hACEMOS N=round((b-a)/h)
 +
t=linspace(t0,tM,M+1);%vector con los tiempos
 +
%----------------------------------------------------
 +
%MATRIZ K ( de N-1 x N-1 )
 +
Ka=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
 +
I=eye(size(Ka));%matriz para sumar los terminos u1,u2,u3 a la matriz K
 +
K=Ka+I;
 +
%-----------------------------------------------------------
 +
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
 +
sol=zeros(N-1,M+1);
 +
U0=u0(xx);%valores iniciales
 +
sol(:,1)=U0;
 +
%----------------------------------------------------------------
 +
%APLICAMOS METODO DEL TRAPECIO
 +
for i=1:M
 +
F1=9*ones(N-1,1);%matriz términos independientes instante i
 +
F1(1)=F1(1)+(ca(t(i))/(h^2)); %no es ca(i) sino ca( t(i)  ), lo mismo en los siguientes
 +
F1(end)=F1(end)+(cb(t(i))/(h^2));
 +
%------------------------------------------------------------------
 +
for i=1:M
 +
F2=9*ones(N-1,1);%matriz términos independientes instante i+1
 +
F2(1)=F2(1)+(ca(t(i+1))/(h^2));
 +
F2(end)=F2(end)+(cb(t(i+1))/(h^2));
 +
%-----------------------------------------------------------------
 +
B= eye(size(K))+((k/2)*K);%matriz resultante de despejar el metodo implicito
 +
B=inv(B);
 +
%k minúscula es tamaño de paso del tiempo
 +
sol(:,i+1)=B*(sol(:,i)+((k/2)*(F1-K*(sol(:,i))))+(k/2)*F2));
 +
end
 +
end
 +
%------------------------------------------------------------
 +
Ua=ca(t);
 +
Ub=cb(t);
 +
sol=[Ua;sol;Ub];
 +
%---------------------------------
 +
[Mt,Mx]=meshgrid(t,x);
 +
subplot(1,2,1);
 +
surf(Mt,Mx,sol);
 +
title('Gráfico de la solucion global')
 +
subplot(1,2,2);
 +
plot(x',sol(:,end)) %grafico de la solucion estacionaria
 +
xlabel('Longitud de la barra');
 +
ylabel('Temperatura');
 +
title('Gráfico de la solucion estacionaria')
 +
%------------------------------------------------------------
 +
%EL ESTADO ESTACIONARIO PODEMOS VISUALIZARLO DIRECTAMENTE DESDE MESHGRID
 +
%-----------------------------------------------------------
 +
%hallamos instante en que la T pasa a estacionaria
 +
for i=1:length(t)
 +
urt=sol(i,:);%sol en instante i
 +
ut=sol(N,:);%penultima fila de la matriz sol, que ya esta en estado estacionario
 +
if max(urt-ut)<=10^-3
 +
inst=t(i);
 +
break
 +
end
 +
end
 +
inst
 +
}}
  
 +
[[Archivo:EJ8PLOT.jpg|marco|centro| ]]
  
 +
El estado estacionario se alcanza en t=3,9 con un error máximo de 10-³.
  
 +
== Cambio de condiciones de contorno ==
  
 +
Ahora tratamos de resolver la ecuación anterior, pero con unas nuevas condiciones de contorno
 +
\[ \left \{\begin{matrix} Ux(0,t)=3\\ U(1,t)=-1,5 \end{matrix} \right . \]
 +
Estas condiciones implican que en el punto de la varilla x=0 hay una pieza aislante que mantiene la temperatura del mismo a 3 grados y evita que haya flujo de temperatura, mientras que el punto x=1 se encuentra situado sobre un objeto que mantiene una temperatura constante de -1,5 grados.
  
  
 +
{{matlab|codigo=
 +
%-------------------------------------
 +
%concidiones iniciales
 +
a=0;%extremo izquierdo
 +
b=1;%extremo derecho
 +
t0=0;%instante inicial
 +
tM=10;%instante final
 +
q=1;%difusividad termica
 +
%-------------------------------------------------
 +
ca=@(t)(0*t)+3;%Condición Neumann
 +
cb=@(t)(0*t)-1.5;%Condición dirichlet
 +
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
 +
%-----------------------------------------------
 +
%discretización intervalo[a,b]
 +
h=0.1;
 +
N=round((b-a)/h);%subintervalos en x
 +
x=linspace(a,b,N+1);
 +
x=x';%el vector x tiene N+1 elementos
 +
xx=x(1:end-1);%tomamos los nodos intermedios
 +
%----------------------------------------------------------
 +
%discretización del tiempo
 +
k=0.1;%tamaño de paso del tiempo
 +
M=round((tM-t0)/k);%como cuando hACEMOS N=round((b-a)/h)
 +
t=linspace(t0,tM,M+1);%vector con los tiempos
 +
%----------------------------------------------------
 +
%MATRIZ K ( de N-1 x N-1 )
 +
Ka=(q/h^2)*(2*diag(ones(1,N))-diag(ones(1,N-1),-1)-diag(ones(1,N-1),1));
 +
Ka(1,2)=-2*q/h^2;
 +
I=eye(size(Ka));%matriz para sumar los términos u1,u2,u3 a la matriz K
 +
K=Ka+I;
 +
%K(1,2)=-2*q/h^2;
 +
%-----------------------------------------------------------
 +
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
 +
sol=zeros(N,M+1);
 +
U0=u0(xx);%valores iniciales
 +
sol(:,1)=U0;
 +
%----------------------------------------------------------------
 +
%APLICAMOS METODO DEL TRAPECIO
 +
for i=1:M
 +
F1=9*ones(N,1);%matriz términos independientes instante i
 +
F1(1)=F1(1)-2*q*ca(t(i))/(h); %no es ca(i) sino ca( t(i)  ), lo mismo en los siguientes
 +
F1(end)=F1(end)+(cb(t(i))/(h^2));
 +
%------------------------------------------------------------------
 +
for i=1:M
 +
F2=9*ones(N,1);%matriz términos independientes instante i+1
 +
F2(1)=F2(1)-2*q*ca(t(i+1))/(h);
 +
F2(end)=F2(end)+(cb(t(i+1))/(h^2));
 +
%-----------------------------------------------------------------
 +
B= eye(size(K))+((k/2)*K);%matriz resultante de despejar el metodo implicito
 +
B=inv(B);
 +
%k minúscula es tamaño de paso del tiempo
 +
sol(:,i+1)=B*(sol(:,i)+((k/2)*(F1-K*(sol(:,i))))+(k/2)*F2);
 +
%en la ecuación del método pone h
 +
%nuestra h es el tamaño de paso en la x
 +
%el tamaño de paso en el tiempo es k, por eso ponemos k
 +
end
 +
end
 +
%------------------------------------------------------------
 +
Ub=cb(t);
 +
sol=[sol;Ub];
 +
%---------------------------------
 +
[Mt,Mx]=meshgrid(t,x);
 +
surf(Mt,Mx,sol);
 +
xlabel('t');
 +
ylabel('x');
 +
zlabel('Temperatura')
 +
}}
  
−
== Conclusiones ==
+
[[Archivo:Graf9.jpg|marco|centro| ]]
−
 
+
−
== Anejos ==
+
  
  

Revisión actual del 17:23 28 abr 2017

Trabajo realizado por estudiantes
Título Ecuación del Calor: "Calentamiento de una varilla" (Grupo 5)
Asignatura Ecuaciones Diferenciales
Curso Curso 2016-17
Autores César Blanco Posadas

Sara González Bravo

Lucía Granados Casado

Patricia del Pozo García

Marta Nogal Prata

Christian Balic Stefanovic

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Descripción del problema a estudiar

La ecuación del calor fue estudiada y propuesta por ‎Jean-Baptiste Joseph Fourier en 1807, en su memoria sobre la propagación del calor en os cuerpos sólidos. En ella propone ademas el gérmen de lo que pasaría a ser la Teoría de las Series de Fourier. La ecuación del calor es un modelo matemático (quizás el más sencillo) que trata de describir la evolución de la temperatura en un cuerpo sólido. Su expresión matemática aplicada a nuestro objeto (varilla), en este caso un objeto unidimensional, tratando como x la variable de la longitud, sería: \[\frac{\partial U}{\partial t}-\frac{\partial ^2U}{\partial x^2}=0\]

Siendo U(x,t) la temperatura de cada punto de la varilla en cada instante t, consideramos una varilla de longitud L=4 y de un cierto material de espesor constante. La varilla es conductora de calor y los extremos x=0 y x=4 están colocados sobre objetos que mantienen una temperatura constante de 5 y 0 grados respectivamente. En el instante inicial la temperatura de la varilla viene determinada por la siguiente ecuación: \[u(x,0)=e^{-3(x-3)^2}+5+x\]

Las condiciones físicas de la varilla son densidad=1 y conductividad térmica (k)=2. Suponiendo que la temperatura u(x,t) de la varilla satisface la ecuación del calor Ut-Uxx=0, el sistema completo de ecuaciones que satisface u(x,t) será: \[ \left \{\begin{matrix} Ut – 2Uxx = 0 ; t>0, x € (0,4) \\ U(0,t) = 5 ; U(4,t) = 0 ; t>0 \\ U(x,0) = e^{-3(x-3)^2}+5+x ; t>0, x € (0,4) \end{matrix} \right . \]


2 Resolución por distintos métodos

En un primer lugar estudiaremos el problema por el método de diferencias finitas aplicando el método del trapecio en este caso. Con un tamaño de paso h=0.1 para la longitud de la varilla y k=0.1 para la discretización del tiempo construimos el siguiente código: 

%concidiones iniciales
a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=2;%difusividad térmica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%funcion f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
g=@(x,t)0*x;%términos independientes
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);
x=linspace(a,b,N+1);
x=x'
xx=x(2:end-1);
%----------------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
K=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
%valor inicial
U0=u0(xx);
%--------------------------------------
%discretización del tiempo
k=0.1;%tamaño de paso del tiempo
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);
%---------------------------------------------------- 
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS MÉTODO DEL TRAPECIO
for i=1:M
G1=g(xx,t(i));%Matriz términos independientes instante i
G2=g(xx,t(i+1));% matriz términos independientes instante i+1
G1(1)=G1(1)+q*ca(t(i))/h^2;
G2(1)=G2(1)+q*ca(t(i+1))/h^2;
G1(end)=G1(end)+q*cb(t(i))/h^2;
G2(end)=G2(end)+q*cb(t(i+1))/h^2;
B= eye(size(K))+((k/2)*K);
B=inv(B);
sol(:,i+1)=B*(sol(:,i)+((k/2)*(G1-K*(sol(:,i))))+(k/2)*G2);
end
%----------------------------------------------------------
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];

%grafico
figure(1)
[Mt,Mx]=meshgrid(t,x)
surf(Mt,Mx,sol);
xlabel('t');
ylabel('x');
figure(2)
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))


Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.

Primeramente lo que obtenemos es una representación en 3D de la evolución de la temperatura en la varilla a lo largo de tiempo. En esta podemos ver como debido a las condiciones iniciales, efectivamente los dos extremos de la varilla mantienen constante su temperatura a lo largo del tiempo (siendo 5 para x=0 y 0 para x=4). Con respecto al resto de la varilla vemos las distintas variaciones que experimentan distintas secciones de ella.


Evolución de la temperatura en el punto medio de la varilla a lo largo del tiempo

En esta gráfica podemos observar como al principio en el punto medio de la varilla se alcanza la temperatura máxima y a medida que discurre el tiempo esta va disminuyendo hasta hacerse nula.


A continuación vamos a estudiar el mismo problema por el método de Euler explícito, Euler implícito y Runge-Kutta de orden 4. Con todo esto compararemos las distintas soluciones y comentaremos las divergencias que podamos encontrar en los resultados.

Primeramente adjuntamos el código correspondiente al método de Euler explícito: 

%concidiones iniciales
a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=2;%difusividad térmica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%función f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
g=@(x,t)0*x;%términos independientes
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);
x=linspace(a,b,N+1);
x=x'
xx=x(2:end-1);
%----------------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
K=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
%valor inicial
U0=u0(xx);
%--------------------------------------
%discretización del tiempo
k=0.001;%tamaño de paso del tiempo
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);
%---------------------------------------------------- 
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS MÉTODO DE EULER EXPLÍCITO
for i=1:M
Gi=g(xx,t(i));%Matriz términos independientes instante i
Gi(1)=Gi(1)+q*ca(t(i))/h^2;
Gi(end)=Gi(end)+q*cb(t(i))/h^2;

%k minúscula es tamaño de paso del tiempo
  sol(:,i+1)=sol(:,i)+k*(-K*sol(:,i)+Gi);
end
%------------------------------------------------------------
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];
%resultados
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))


Evolución de la temperatura en el punto medio de la varilla frente al tiempo (Inestabilidad)
Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma (Inestabilidad)

En estas dos primeras gráficas debido al tamaño de paso utilizado, excesivamente grande (k=0.1), obtenemos unos resultados los cuales denotan inestabilidad. Esto lo solucionamos, como se puede ver en las dos siguientes representaciones, reduciendo dicho tamaño de paso a k=0.001 .

Evolución de la temperatura en el punto medio de la varilla frente al tiempo
Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma

A continuación resolvemos el problema por el método de Euler implícito: 

%concidiones iniciales
a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=2;%difusividad térmica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%función f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
g=@(x,t)0*x;%términos independientes
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);
x=linspace(a,b,N+1);
x=x'
xx=x(2:end-1);
%----------------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
K=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
%valor inicial
U0=u0(xx);
%--------------------------------------
%discretización del tiempo
k=0.1;%tamaño de paso del tiempo
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);
%---------------------------------------------------- 
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS MÉTODO DE EULER IMPLÍCITO
for i=1:M

G2=g(xx,t(i+1));% matriz terminos independientes instante i+1
G2(1)=G2(1)+q*ca(t(i+1))/h^2;
G2(end)=G2(end)+q*cb(t(i+1))/h^2;
B= eye(size(K))+((k)*K);
B=inv(B);

sol(:,i+1)=B*(sol(:,i)+(k*(G2)));

end

Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];
%resultados
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))


Evolución de la temperatura en el punto medio de la varilla frente al tiempo
Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.

En este caso, al estar estudiando el problema a través de un método implícito no tenemos problemas de inestabilidad, por lo que no ha sido necesario modificar el tamaño de paso.

Por último, empleamos el método de Runge-Kutta para resolver el problema. 

a=0;
b=4;
t0=0;
tM=10;
q=2;
%-------------------------------------------------
%Datos de entrada
%ca: Funcion f(t) de la condicion de contorno en x=a (t vector)
ca=@(t)0*t+5;
%cb: Funcion f(t) de la condicion de contorno en x=b (t vector)
cb=@(t)0*t;
%u0: forma analítica f(x) funcion valor inicial (x vector)
u0=@(x)(exp(-3*(x-3).^2))+5+x;
%g: funcion término independiente (x será vector, t escalar)
g=@(x,t)0*x;
%--------------------------------------------------------------------------
Discretización del intervalo [a,b]
h=0.1;
N=round((b-a)/h);
x=linspace(a,b,N+1);
x=x';
xx=x(2:end-1);
%-----------------------------------------------------------------------
%Matriz K (de N-1 x N-1)
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
 %-------------------------------------------------
%valor inicial 
U0=u0(xx);
 %-------------------------------------------------
%discretización del tiempo
k=0.001;%paso en t
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);
%------------------------------------------------- 
sol=zeros(N-1,M+1);
sol(:,1)=U0;
 
for i=1:M
   
   G1=g(xx,t(i)); 
   G1(1)=G1(1)+q*ca(t(i))/h^2;
   G1(end)=G1(end)+q*cb(t(i))/h^2;
   
   G2=g(xx,t(i)+k/2);
   G2(1)=G2(1)+q*ca(t(i)+k/2)/h^2;
   G2(end)=G2(end)+q*cb(t(i)+k/2)/h^2;
   
   G4=g(xx,t(i+1));
   G4(1)=G4(1)+q*ca(t(i+1))/h^2;
   G4(end)=G4(end)+q*cb(t(i+1))/h^2;
   
   K1=-K*sol(:,i)+G1;
   K2=-K*(sol(:,i)+1/2*K1*k)+G2;
   K3=-K*(sol(:,i)+1/2*K2*k)+G2;
   K4=-K*(sol(:,i)+K3*k)+G4;
   sol(:,i+1)=sol(:,i)+k/6*(K1+2*K2+2*K3+K4);
   
   
end 
 %---------------------------------------------------------
%añadimos las condiciones Dirichlet
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub]
%-------------------------------------------------
%gráfico 3D
[Mt,Mx]=meshgrid(t,x);
surf(Mt,Mx,sol);
xlabel('t');
ylabel('x');
%-------------------------------------------------
%resultados
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))


Evolución de la temperatura en el punto medio de la varilla frente al tiempo (Inestabilidad)
Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma (Inestabilidad)
Evolución de la temperatura en el punto medio de la varilla frente al tiempo
Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma

Al emplear el tamaño de paso k=0.1 hemos detectado inestabilidad en la solución, por lo que hemos tratado de solventarla utilizando un tamaño de paso menor, siendo este k=0.001. El código adjunto es este último.


3 Variación de la temperatura para tiempos grandes

A continuación comprobamos gráficamente como evoluciona la temperatura en el tiempo. Se observa que a medida que avanza el tiempo la solución tiende a parecerse a la obtenida en el caso estacionario. Para que sea fácil de apreciar, se representa la distancia entre la solución estacionaria y la solución u(x,t) para distintos tiempos. 

a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=2;%difusividad termica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%función f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
g=@(x,t)0*x;%términos indepencientes
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);
x=linspace(a,b,N+1);
x=x'
xx=x(2:end-1);
%----------------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
K=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
%valor inicial
U0=u0(xx);
%--------------------------------------
%discretización del tiempo
k=0.1;%tamaño de paso del tiempo
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);
%---------------------------------------------------- 
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS METODO DEL TRAPECIO
for i=1:M

G1=g(xx,t(i));%Matriz terminos independientes instante i
G2=g(xx,t(i+1));% matriz terminos independientes instante i+1
G1(1)=G1(1)+q*ca(t(i))/h^2;
G2(1)=G2(1)+q*ca(t(i+1))/h^2;
G1(end)=G1(end)+q*cb(t(i))/h^2;
G2(end)=G2(end)+q*cb(t(i+1))/h^2;
B= eye(size(K))+((k/2)*K);%matriz resultante de despejar el metodo implicito
B=inv(B);

sol(:,i+1)=B*(sol(:,i)+((k/2)*(G1-K*(sol(:,i))))+(k/2)*G2);

end
%------------------------------------------------------------
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];
%------------------------------------------------------------------------
U=@(x,t)((-5/4)*x)+5;% f aproximada para valores estacionarios
%------------------------------------------------------------
%hallamos instante en que la T pasa a estacionaria
for i=1:length(t)
urt=sol(i,:);
ut=U(x,t(i));
if max(urt-ut)<=0.05
inst=t(i);
break
end
end
inst %instante en el que la temperatura permanece estacionaria
%CALCULAMOS VARIACIONES ENTRE LA SOLUCIÓN POR TRAPECIO Y LA APROXIMADA
%instantes
ta=0;
tb=1;
tc=2;
td=10;
%números de las columnas correspondientes a nuestros instantes de tiempo
cola=round((ta-t0)/k+1);
colb=round((tb-t0)/k+1);
colc=round((tc-t0)/k+1);
cold=round((td-t0)/k+1);
ut=U(x);%vector soluciones estacionarias
solu=[sol(:,cola),sol(:,colb),sol(:,colc),sol(:,cold),ut];
%  sol2 = [Tªs en t=0, Tªs en t=1,........., Tªs ESTACIONARIAS]
%DISTANCIA  para t=0
dista=(solu(:,1)-solu(:,5))';
%DISTANCIA  para t=1
distb=(solu(:,2)-solu(:,5))';
%DISTANCIA  para t=2
distc=(solu(:,3)-solu(:,5))';
%DISTANCIA  para t=10
distd=(solu(:,4)-solu(:,5))';
%diagramas----------------------------------------------------------------
subplot(2,2,1);
hold on
plot(x,ut,'r')
plot(x,solu(:,1)','b')
xlabel('longitud de la barra')
ylabel('temperatura')
title('distancia U(x,t) y solución estacionaria en t=0')
legend('solución estacionaria','solución numérica)');
hold off
%%--------------------------------------------------------------
subplot(2,2,2);
hold on
plot(x,ut,'r')
plot(x,solu(:,2)','b')
xlabel('longitud de la barra')
ylabel('temperatura')
title('distancia U(x,t) y solución estacionaria en t=1')
legend('solución estacionaria','solución numérica)');
hold off
%%--------------------------------------------------------------
subplot(2,2,3);
hold on
plot(x,ut,'r')
plot(x,solu(:,3)','b')
xlabel('longitud de la barra')
ylabel('temperatura')
title('distancia U(x,t) y solución estacionaria en t=2')
legend('solución estacionaria','solución numérica)');
hold off
%%--------------------------------------------------------------
subplot(2,2,4);
hold on
plot(x,ut,'r')
plot(x,solu(:,4)','b')
xlabel('longitud de la barra')
ylabel('temperatura')
title('distancia U(x,t) y solución estacionaria en t=10')
legend('solución estacionaria','solución numérica)');
hold off
%%-----------------------------------------------------------


Dist apart 5.JPG

La ecuación que debería verificar el estado estacionario es una recta, la cual podemos obtener una vez que hemos despreciado el término Ut. Esta recta es \[\frac{5x}{4} +5=0\]. Como se observa en las gráficas, la función a la que debería parecerse la solución u(x,t) para tiempos grandes es precisamente la recta correspondiente al estado estacionario.

4 Extremo con pieza aislante

Se coloca en el extremo izquierdo una pieza aislante que hace que no haya pérdida de calor en ese extremo, por lo que el flujo de temperatura es nulo. Esto nos cambia las condiciones de contorno del problema, pasando de tener dos condiciones Dirichlet a tener una condición Neumann y otra Dirichlet. 

%DATOS DEL PROBLEMA
a=0;
b=4;
t0=0;
tM=10;
q=2;
%----------------------------------------------------------------
%Discretización
h=0.1;
N=round((b-a)/h);
x=linspace(a,b,N+1);
x=x';
xx=x(1:end-1); %tomamos los elementos del primero al penúltimo
%----------------------------------------------------------------
%Funciones.Datos de entrada.
%ca:Función f(t) de la condición de contorno en x=a(t vector)
ca=0; %Condición Neumann
%cb:Función f(t) de la condición de contorno en x=b(t vector)
cb=0;%Condición Dirichlet
%u0Forma analitica f(x) funcion valor inicial(x vector)
u0=exp(-3*((xx-3).^2))+5+xx;
%Matriz K
K=q/h^2*(2*diag(ones(1,N))-diag(ones(1,N-1),-1) -diag(ones(1,N-1),1));
K(1,2)=-4/h^2;
%----------------------------------------------------------------
%Discretización del tiempo
k=0.1; %Anchura de paso en el tiempo(k=0.5*h^2)
M=round((tM-t0)/k);%numero de subintervalos en t
t=linspace(t0,tM,M+1);
%sol matriz que tendrá la solución para cada punto e instante
sol=zeros(N,M+1);
sol(:,1)=u0;
%----------------------------------------------------------------
%Términos independiente
G=0*xx;
G(1)=G(1)-2*q*ca/h;
G(end)=G(end)+q*cb/h^2;

for i=1:M
    
    sol(:,1+i)=(eye(size(K))+k*K)\(sol(:,i)+k*G);
    
end
%----------------------------------------------------------------
%Añado frontera(Solo Dirichlet)
Ub=q*cb*ones(1,length(t));
sol=[sol;Ub];
%grafico
[Mt,Mx]=meshgrid(t,x);
surf(Mt,Mx,sol);
xlabel('t');
ylabel('x');
zlabel('Temperatura')
Aislante.jpg


5 Método de Fourier

Partiendo de una solución analítica previa obtenemos la solución en forma de sumatorio de infinitos términos, reduciendo este número de términos numéricamente calcularemos una solución aproximada. Mediante el siguiente código de Matlab resolvemos el problema numéricamente para 1,3,5,10 y 20 términos de la serie de Fourier. 

a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
h=0.1;k=0.1;
N=round((b-a)/h);
M=round((tM-t0)/k);

x=linspace(a,b,N+1);
t=linspace(t0,tM,M+1);

u0=(exp(-3*((x-3).^2)))+5+x;
[xx,tt]=meshgrid(x,t);

%-------------------------------------------------------
U=0;
Q=1
for j=1:Q
    
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));

hold on

Q=3;
for j=1:Q
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));

Q=5;
for j=1:Q
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));


Q=10;
for j=1:Q
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));

Q=20;
for j=1:Q
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));

legend('Q=1','Q=3','Q=5','Q=10','Q=20')
title('Serie de Fourier con 1,3,5,10,20 términos de la serie')
hold off


PLOTFOURER.jpg

Con este gráfico queremos representar una vez pasadas 0.5 unidades de tiempo la temperatura a lo largo de la barra. La aproximación más exacta la obtenemos con Q=20 y todas son proporcionales decreciendo cumpliendo así las condiciones de contorno impuestas en el enunciado u(0,t)=5 y u(4,t)=0.

6 Pérdida de calor en la varilla

Suponemos ahora una pérdida de calor de la varilla a través del aire que tiene una temperatura constante de 16 grados. Con la siguiente ecuación \[\frac{\partial U}{\partial t}-\frac{\partial ^2U}{\partial x^2}+ (u-9)=0\] y manteniendo las condiciones de contorno establecidas en el primer apartado, calculamos el estado estacionario y el tiempo que se tarda en alcanzar el estado estacionario con un error de 10-³. 

%CONDICIONES DIRICHLET & DIRICHLET
%-------------------------------------
%concidiones iniciales
a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=1;%difusividad termica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%función f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);%subintervalos en x
x=linspace(a,b,N+1);
x=x';%el vector x tiene N+1 elementos
xx=x(2:end-1);%tomamos los nodos intermedios
%----------------------------------------------------------
%discretización del tiempo
k=0.1;%tamaño de paso del tiempo
M=round((tM-t0)/k);%como cuando hACEMOS N=round((b-a)/h)
t=linspace(t0,tM,M+1);%vector con los tiempos
%----------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
Ka=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
I=eye(size(Ka));%matriz para sumar los terminos u1,u2,u3 a la matriz K
K=Ka+I;
%-----------------------------------------------------------
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
U0=u0(xx);%valores iniciales
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS METODO DEL TRAPECIO
for i=1:M
F1=9*ones(N-1,1);%matriz términos independientes instante i
F1(1)=F1(1)+(ca(t(i))/(h^2)); %no es ca(i) sino ca( t(i)  ), lo mismo en los siguientes
F1(end)=F1(end)+(cb(t(i))/(h^2));
%------------------------------------------------------------------
for i=1:M
F2=9*ones(N-1,1);%matriz términos independientes instante i+1
F2(1)=F2(1)+(ca(t(i+1))/(h^2));
F2(end)=F2(end)+(cb(t(i+1))/(h^2));
%-----------------------------------------------------------------
B= eye(size(K))+((k/2)*K);%matriz resultante de despejar el metodo implicito
B=inv(B);
%k minúscula es tamaño de paso del tiempo
sol(:,i+1)=B*(sol(:,i)+((k/2)*(F1-K*(sol(:,i))))+(k/2)*F2));
end
end
%------------------------------------------------------------
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];
%---------------------------------
[Mt,Mx]=meshgrid(t,x);
subplot(1,2,1);
surf(Mt,Mx,sol);
title('Gráfico de la solucion global')
subplot(1,2,2);
plot(x',sol(:,end)) %grafico de la solucion estacionaria
xlabel('Longitud de la barra');
ylabel('Temperatura');
title('Gráfico de la solucion estacionaria')
%------------------------------------------------------------
%EL ESTADO ESTACIONARIO PODEMOS VISUALIZARLO DIRECTAMENTE DESDE MESHGRID
%-----------------------------------------------------------
%hallamos instante en que la T pasa a estacionaria
for i=1:length(t)
urt=sol(i,:);%sol en instante i
ut=sol(N,:);%penultima fila de la matriz sol, que ya esta en estado estacionario
if max(urt-ut)<=10^-3
inst=t(i);
break
end
end
inst


EJ8PLOT.jpg

El estado estacionario se alcanza en t=3,9 con un error máximo de 10-³.

7 Cambio de condiciones de contorno

Ahora tratamos de resolver la ecuación anterior, pero con unas nuevas condiciones de contorno \[ \left \{\begin{matrix} Ux(0,t)=3\\ U(1,t)=-1,5 \end{matrix} \right . \] Estas condiciones implican que en el punto de la varilla x=0 hay una pieza aislante que mantiene la temperatura del mismo a 3 grados y evita que haya flujo de temperatura, mientras que el punto x=1 se encuentra situado sobre un objeto que mantiene una temperatura constante de -1,5 grados. 


%-------------------------------------
%concidiones iniciales
a=0;%extremo izquierdo
b=1;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=1;%difusividad termica
%-------------------------------------------------
ca=@(t)(0*t)+3;%Condición Neumann
cb=@(t)(0*t)-1.5;%Condición dirichlet
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);%subintervalos en x
x=linspace(a,b,N+1);
x=x';%el vector x tiene N+1 elementos
xx=x(1:end-1);%tomamos los nodos intermedios
%----------------------------------------------------------
%discretización del tiempo
k=0.1;%tamaño de paso del tiempo
M=round((tM-t0)/k);%como cuando hACEMOS N=round((b-a)/h)
t=linspace(t0,tM,M+1);%vector con los tiempos
%----------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
Ka=(q/h^2)*(2*diag(ones(1,N))-diag(ones(1,N-1),-1)-diag(ones(1,N-1),1));
Ka(1,2)=-2*q/h^2;
I=eye(size(Ka));%matriz para sumar los términos u1,u2,u3 a la matriz K
K=Ka+I;
%K(1,2)=-2*q/h^2;
%-----------------------------------------------------------
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N,M+1);
U0=u0(xx);%valores iniciales
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS METODO DEL TRAPECIO
for i=1:M
F1=9*ones(N,1);%matriz términos independientes instante i
F1(1)=F1(1)-2*q*ca(t(i))/(h); %no es ca(i) sino ca( t(i)  ), lo mismo en los siguientes
F1(end)=F1(end)+(cb(t(i))/(h^2));
%------------------------------------------------------------------
for i=1:M
F2=9*ones(N,1);%matriz términos independientes instante i+1
F2(1)=F2(1)-2*q*ca(t(i+1))/(h);
F2(end)=F2(end)+(cb(t(i+1))/(h^2));
%-----------------------------------------------------------------
B= eye(size(K))+((k/2)*K);%matriz resultante de despejar el metodo implicito
B=inv(B);
%k minúscula es tamaño de paso del tiempo
sol(:,i+1)=B*(sol(:,i)+((k/2)*(F1-K*(sol(:,i))))+(k/2)*F2);
%en la ecuación del método pone h
%nuestra h es el tamaño de paso en la x
%el tamaño de paso en el tiempo es k, por eso ponemos k
end
end
%------------------------------------------------------------
Ub=cb(t);
sol=[sol;Ub];
%---------------------------------
[Mt,Mx]=meshgrid(t,x);
surf(Mt,Mx,sol);
xlabel('t');
ylabel('x');
zlabel('Temperatura')


Graf9.jpg
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_de_Calor:_%22Calentamiento_Varilla%22&oldid=39391»