• La parte izquierda de la varilla tiene una longitud d, está hecha de un material de calor específico c₁, densidad ρ₁ y conductividad térmica K₁

• La parte derecha de la varilla tiene una longitud L-d, está hecha de un material de calor específico c₂, densidad ρ₂ y conductividad térmica K₂

Estudiaremos una situación en la que los extremos de la varilla se mantienen a temperatura constante T(0,t)=T₀ y T(L, t)=T_L.

La temperatura inicial de la varilla es T(x ,0)

En un instante t, la distribución de temperaturas en la parte izquierda es T₁(x, t) y en la parte derecha, T₂(x, t). Soluciones de las ecuaciones diferenciales en derivadas parciales

${\displaystyle \begin{array}{l}\frac{\partial T_1(x,t)}{\partial t}=\frac{K_1}{{\rho }_1{c}_1}\frac{{\partial }^2{T}_1(x,t)}{\partial x^2},0\le x<d\\ \frac{\partial T_2(x,t)}{\partial t}=\frac{K_2}{{\rho }_2{c}_2}\frac{{\partial }^2{T}_2(x,t)}{\partial x^2},d\le x\le L\end{array}}$

Con las condiciones iniciales y de controno especificadas. En la posición d de contacto de los dos tramos se ha de cumplir

${\displaystyle \begin{array}{l}{T}_1(d,t)=T_2(d,t)\\ K_1{\frac{\partial T_1(x,t)}{\partial x}|}_{x=d}=K_2{\frac{\partial T_2(x,t)}{\partial x}|}_{x=d}\end{array}}$

La función T(x,t) es continua y el flujo de calor a través de la sección de la varilla en la posición de contacto d, a uno y otro lado, es el mismo

Una situación similar la encontramos en la propagación de ondas en una cuerda no homogénea

Estado estacionario

El estado estacionario está caracterizado por una temperatura no uniforme en la varilla

Como en una varilla con temperaturas fijas en los extremos, la temperatura varía linealmente entre T₀ y T_d en la parte izquierda y entre T_d y T_L en la parte derecha. La ecuación de las rectas es

T₁(x,∞)=m₁x+b₁, parte izquierda
T₂(x,∞)=m₂x+b₂, parte derecha

${\displaystyle \begin{array}{l}{T}_1(x,\infty )=\frac{T_d-T_0}{d}x+T_0\\ T_2(x,\infty )=\frac{T_L-T_d}{L-d}x+\frac{T_{d}L-T_{L}d}{L-d}=\frac{T_d(L-x)+T_L(x-d)}{L-d}\end{array}}$

Calculamos la temperatura T_d en la posición de contacto, d, mediante la relación

${\displaystyle \begin{array}{l}{K}_1{\frac{d{T}_1(x,\infty )}{dx}|}_{x=d}=K_2{\frac{d{T}_2(x,\infty )}{dx}|}_{x=d}\\ K_1{m}_1=K_2{m}_2\\ T_d=\frac{K_{2}d{T}_L+K_1(L-d)T_0}{K_{2}d+K_1(L-d)}\end{array}}$

Ejemplo

Sea la varilla

• Tramo izquierdo, el material tiene las siguientes caraterísticas: densidad ρ₁=1, calor específico, c₁=1, conductividad térmica, K₁=1
• Tramo derecho, el material tiene las siguientes caraterísticas: densidad ρ₂=2, calor específico, c₂=2, conductividad térmica, K₂=2
• La longitud de la varilla es L=1, la del tramo izquierdo es d=L/3
• La temperatura en los extremos es T₀=10 y T_L=0

k1=1; %conductividad térmica
c1=1; %calor específico
r1=1; %densidad
k2=2;
c2=2;
r2=2;

T0=10; %temperatura, extremo izquierdo
TL=0; %extremo derecho
L=1; %longitud varilla
d=L/3; %tramo izquierdo

%temperatura en la posición de contacto
T_d=(TL/(L-d)+k1*T0/(d*k2))/(k1/(k2*d)+1/(L-d)); 
m1=(T_d-T0)/d; %pendiente
b1=T0; %ordenada en el origen
m2=(TL-T_d)/(L-d);
b2=(T_d*L-TL*d)/(L-d);
line([0,d],[T0,T_d])
line([d,L],[T_d,TL])
line([d,d],[0,T_d],'lineStyle','--','color','k')
xlabel('x')
ylabel('T')
grid on
title('Estado estacionario')

Solución de la ecuación de la conducción del calor

Definimos las funciones: u₁(x,t)=T₁(x,t)-T₁(x,∞) y u₂(x,t)=T₂(x,t)-T₂(x,∞). En términos de estas nuevas funciones, las ecuaciones de la conducción del calor se esciben,

${\displaystyle \begin{array}{l}\frac{\partial u_1(x,t)}{\partial t}=\frac{K_1}{{\rho }_1{c}_1}\frac{{\partial }^2{u}_1(x,t)}{\partial x^2},0\le x<d\\ \frac{\partial u_2(x,t)}{\partial t}=\frac{K_2}{{\rho }_2{c}_2}\frac{{\partial }^2{u}_2(x,t)}{\partial x^2},d\le x\le L\end{array}}$

• Condiciones iniciales

${\displaystyle \begin{array}{l}{u}_1(x\mathrm{,0})=T_1(x\mathrm{,0})-\frac{T_d-T_0}{d}x+T_0\\ u_2(x\mathrm{,0})=T_2(x\mathrm{,0})-\frac{T_d(L-x)+T_L(x-d)}{L-d}\end{array}}$

• Condiciones de contorno

${\displaystyle \begin{array}{l}{u}_1(\mathrm{0,}t)=T_1(\mathrm{0,}t)-T_1(\mathrm{0,}\infty )=T_0-T_0=0\\ u_2(\mathrm{0,}t)=T_2(\mathrm{0,}t)-T_2(\mathrm{0,}\infty )=T_L-T_L=0\end{array}}$

• Posición d de contacto

${\displaystyle \begin{array}{l}{u}_1(d,t)=u_2(d,t)\\ K_1{\frac{\partial u_1(x,t)}{\partial x}|}_{x=d}=K_2{\frac{\partial u_2(x,t)}{\partial x}|}_{x=d}\end{array}}$

Variables separadas

Buscamos una solución de la forma u(x,t)=F(x)·G(t), variables separadas

${\displaystyle \begin{array}{l}\frac{\partial u(x,t)}{\partial t}=\alpha \frac{{\partial }^{2}u(x,t)}{\partial x^2}\\ F(x)\frac{dG(t)}{dt}=\alpha G(t)\frac{d^{2}F(x)}{d{x}^2}\\ \frac{1}{G(t)}\frac{dG(t)}{dt}=\frac{\alpha }{F(x)}\frac{d^{2}F(x)}{d{x}^2}=-{\omega }^2\end{array}}$

Integramos la primera ecuación diferencial

${\displaystyle \begin{array}{l}\frac{dG(t)}{dt}+{\omega }^{2}G(t)=0\hfill \\ G(t)=G(0)·\exp (-{\omega }^{2}t)\hfill \end{array}}$

Integramos la segunda ecuación diferencial

${\displaystyle \begin{array}{l}\frac{d^{2}F(x)}{d{x}^2}+\frac{{\omega }^2}{\alpha }F(x)=0\\ \frac{d^{2}F(x)}{d{x}^2}+{\mu }^2{\omega }^{2}F(x)=\mathrm{0,}{\mu }^2=\frac{1}{\alpha }\\ F(x)=A\sin \left(\mu \omega t\right)+B\cos \left(\mu \omega t\right)\end{array}}$

Las soluciones de la segunda, para cada una de las tramos de la varilla, son

${\displaystyle \begin{array}{l}{F}_1(x)=A_1\sin ({\mu }_1\omega x)+B_1\cos ({\mu }_1\omega x),0\le x<d\\ F_2(x)=A_2\sin ({\mu }_2\omega x)+B_2\cos ({\mu }_2\omega x),d\le x\le L\end{array}}$

• La condición de contorno en el extremo izquierdo, x=0, u(0,t)=0, F₁(0)=0, lo que implica que, B₁=0

• La condición de contorno en el extremo derecho, x=L, u(L, t)=0, F₂(L)=0, lo que implica que, $A_2\sin ({\mu }_2\omega L)+B_2\cos ({\mu }_2\omega L)=0$

• En la posición d de contacto, se cumplirá que

${\displaystyle \begin{array}{l}{F}_1(d)=F_2(d)\\ K_1{\frac{d{F}_1(x)}{dx}|}_{x=d}=K_2{\frac{d{F}_2(x)}{dx}|}_{x=d}\\ \{\begin{array}{l}{A}_1\sin ({\mu }_1\omega d)=A_2\sin ({\mu }_2\omega d)+B_2\cos ({\mu }_2\omega d)\\ {\mu }_1\omega A_1\cos ({\mu }_1\omega d)={\mu }_2\omega \left(A_2\cos ({\mu }_2\omega d)-B_2\sin ({\mu }_2\omega d)\right)\end{array}\end{array}}$

Tenemos un sistema homogéneo de tres ecuaciones con tres incógnitas A₁ ,A₂ y B₂. Despejamos A₂ y B₂ en función de A₁

${\displaystyle \begin{array}{l}{A}_2=-A_1\frac{\sin ({\mu }_1\omega d)·\cos ({\mu }_2\omega L)}{\sin \left({\mu }_2\omega (L-d)\right)}\\ B_2=A_1\frac{\sin ({\mu }_1\omega d)·\sin ({\mu }_2\omega L)}{\sin \left({\mu }_2\omega (L-d)\right)}\end{array}}$

Las funciones F₁(x) y F₂(x) se expresan

${\displaystyle \begin{array}{l}\{\begin{array}{l}{F}_1(x)=A_1\sin ({\mu }_1\omega x)\\ F_2(x)=-A_1\frac{\sin ({\mu }_1\omega d)·\cos ({\mu }_2\omega L)}{\sin \left({\mu }_2\omega (L-d)\right)}\sin ({\mu }_2\omega x)+A_1\frac{\sin ({\mu }_1\omega d)·\sin ({\mu }_2\omega L)}{\sin \left({\mu }_2\omega (L-d)\right)}\cos ({\mu }_2\omega x)\end{array}\\ \{\begin{array}{l}{F}_1(x)=A_1\sin ({\mu }_1\omega x),0\le x<d\\ F_2(x)=A_1\frac{\sin ({\mu }_1\omega d)}{\sin \left({\mu }_2\omega (L-d)\right)}\sin \left({\mu }_2\omega (L-x)\right),d\le x\le L\end{array}\end{array}}$

Además, obtenemos una ecuación transcendente (el determinante de los coeficientes del sistema homogéneo ha de ser nulo)

${\displaystyle K_1{\mu }_1\cot ({\mu }_1\omega d)+K_2{\mu }_2\cot \left({\mu }_2\omega (L-d)\right)=0}$

Representamos la función f(ω)

k1=1; %conductividad térmica
c1=1; %calor específico
r1=1; %densidad
k2=2;
c2=2;
r2=2;

T0=10; %temperatura, extremo izquierdo
TL=0; %extremo derecho
L=1; %longitud varilla
d=L/3; %tramo izquierdo
a1=sqrt(c1*r1/k1);
a2=sqrt(c2*r2/k2);

%f=@(w) cos(a1*w*d).*sin(a2*w*(L-d))+(k2*a2/(k1*a1))*cos(a2*w*(L-d)).
*sin(a1*w*d);
f=@(x) cot(a1*x*d)+(k2*a2/(k1*a1))*cot(a2*x*(L-d));
fplot(f,[0,20])
ylim([-10,10])
xlabel('\omega')
ylabel('f(\omega)')
grid on
title('Raíces de la ecuación transcendente')

Para calcular las raíces, es más conveniente expresar esta ecuación de la forma

${\displaystyle K_1{\mu }_1\cos ({\mu }_1\omega d)\sin \left({\mu }_2\omega (L-d)\right)+K_2{\mu }_2\cos \left({\mu }_2\omega (L-d)\right)\sin ({\mu }_1\omega d)=0}$

La primera raíz es ω=0

Utilizamos la función raices que llama al procedimiento fzero de MATLAB para calcular las raíces múltiples de la ecuación transcendente.

...
w=linspace(0,20,20);
f=@(w) cos(a1*w*d).*sin(a2*w*(L-d))+(k2*a2/(k1*a1))*cos(a2*w*(L-d)).
*sin(a1*w*d);
raiz=raices(f,w); 
...

>> raiz
raiz =    2.0901    4.9664    7.7754    9.7136   11.9808   14.8978   17.6080   
19.4447

La representación gráfica nos ayuda a establecer el paso de exploración de la función y el intervalo o el número de raíces que deseamos obtener. Es conveniente, comprobar que las raíces obtenidas numéricamente se corresponden con los ceros de la función f(ω)

Solución completa

La solución completa u(x,t) es la superposición de los productos de las funciones F_n(x)·G_n(x) para cada ω_n

${\displaystyle \begin{array}{l}{u}_1(x,t)={\displaystyle \sum _{n=1}^{\infty }{A}_n\sin ({\mu }_1{\omega }_{n}x)\exp (-{\omega }_n^{2}t)},0\le x<d\\ u_2(x,t)={\displaystyle \sum _{n=1}^{\infty }{A}_n\frac{\sin ({\mu }_1{\omega }_{n}d)}{\sin \left({\mu }_2{\omega }_n(L-d)\right)}\sin \left({\mu }_2{\omega }_n(L-x)\right)\exp (-{\omega }_n^{2}t)},d\le x\le L\end{array}}$

La primera raíz ω=0 no contribuye a la suma

De forma abreviada expresamos la superposición

${\displaystyle \begin{array}{l}{u}_1(x,t)={\displaystyle \sum _{n=1}^{\infty }{A}_n{X}_{1n}(x)\exp (-{\omega }_n^{2}t)},0\le x<d\\ u_2(x,t)={\displaystyle \sum _{n=1}^{\infty }{A}_n{X}_{2n}(x)\exp (-{\omega }_n^{2}t)},d\le x\le L\end{array}}$

Relación de ortogonalidad

Se puede demostar que

${\displaystyle c_1{\rho }_1{\displaystyle \underset{0}{\overset{d}{\int }}{X}_{1n}(x)}{X}_{1m}(x)dx+c_2{\rho }_2{\displaystyle \underset{d}{\overset{L}{\int }}{X}_{2n}(x)}{X}_{2m}(x)dx=0}$

Evaluamos esta relación, utilizando procedimientos numéricos

k1=1; %conductividad térmica
c1=1; %calor específico
r1=1; %densidad
k2=2;
c2=2;
r2=2;

T0=10; %temperatura, extremo izquierdo
TL=0; %extremo derecho
L=1; %longitud varilla
d=L/3; %tramo izquierdo
a1=sqrt(c1*r1/k1);
a2=sqrt(c2*r2/k2);

w=linspace(0,20,20);
f=@(w) cos(a1*w*d).*sin(a2*w*(L-d))+(k2*a2/(k1*a1))*cos(a2*w*(L-d)).
*sin(a1*w*d);
raiz=raices(f,w); 
w1=raiz(4); 
w2=raiz(5);
X1_1=@(x) sin(a1*w1*x);
X1_2=@(x) sin(a1*w2*x);
X2_1=@(x) sin(a1*w1*d)*sin(a2*w1*(L-x))/sin(a2*w1*(L-d));
X2_2=@(x) sin(a1*w2*d)*sin(a2*w2*(L-x))/sin(a2*w2*(L-d));
Z1=@(x) X1_1(x).*X1_2(x);
Z2=@(x) X2_1(x).*X2_2(x);
res=c1*r1*integral(Z1,0,d)+c2*r2*integral(Z2,d,L);
disp(res)

   1.4710e-15

Coeficientes A_n

Los coeficientes A_n se calculan a partir de las distribución inicial de temperatura u(x,0)

${\displaystyle u(x\mathrm{,0})\{\begin{array}{l}{u}_1(x\mathrm{,0})={\displaystyle \sum _{n=1}^{\infty }{A}_n{X}_{1n}},0\le x<d\\ u_2(x\mathrm{,0})={\displaystyle \sum _{n=1}^{\infty }{A}_n{X}_{2n}},d\le x\le L\end{array}}$

Multiplicamos la primera ecuación por X_1m y la segunda por X_2m e integramos entre 0 y d, la primera y entre d y L, la segunda

${\displaystyle \begin{array}{l}\{\begin{array}{l}{\displaystyle \underset{0}{\overset{d}{\int }}{u}_1(x\mathrm{,0})X_{1m}dx}={\displaystyle \underset{0}{\overset{d}{\int }}\left({\displaystyle \sum _{n=1}^{\infty }{A}_n{X}_{1n}{X}_{1m}}\right)}dx\\ {\displaystyle \underset{d}{\overset{L}{\int }}{u}_2(x\mathrm{,0})X_{2m}dx}={\displaystyle \underset{d}{\overset{L}{\int }}\left({\displaystyle \sum _{n=1}^{\infty }{A}_n{X}_{2n}{X}_{2m}}\right)}dx\end{array}\\ \{\begin{array}{l}{\displaystyle \underset{0}{\overset{d}{\int }}{u}_1(x\mathrm{,0})X_{1m}dx}={\displaystyle \sum _{n=1}^{\infty }{A}_n}·{\displaystyle \underset{0}{\overset{d}{\int }}\left(X_{1n}{X}_{1m}dx\right)}dx\\ {\displaystyle \underset{d}{\overset{L}{\int }}{u}_2(x\mathrm{,0})X_{2m}dx}={\displaystyle \sum _{n=1}^{\infty }{A}_n}·{\displaystyle \underset{d}{\overset{L}{\int }}\left(X_{2n}{X}_{2m}dx\right)}dx\end{array}\end{array}}$

Sumamos

${\displaystyle \begin{array}{l}{\rho }_1{c}_1{\displaystyle \underset{0}{\overset{d}{\int }}{u}_1(x\mathrm{,0})X_{1m}dx}+{\rho }_2{c}_2{\displaystyle \underset{d}{\overset{L}{\int }}{u}_2(x\mathrm{,0})X_{2m}dx}={\displaystyle \sum _{n=1}^{\infty }{A}_n}·\left\{{\rho }_1{c}_1{\displaystyle \underset{0}{\overset{d}{\int }}\left(X_{1n}{X}_{1m}dx\right)}dx+{\rho }_2{c}_2{\displaystyle \underset{d}{\overset{L}{\int }}\left(X_{2n}{X}_{2m}dx\right)}dx\right\}\\ {\rho }_1{c}_1{\displaystyle \underset{0}{\overset{d}{\int }}{u}_1(x\mathrm{,0})X_{1m}dx}+{\rho }_2{c}_2{\displaystyle \underset{d}{\overset{L}{\int }}{u}_2(x\mathrm{,0})X_{2m}dx}=A_m\left\{{\rho }_1{c}_1{\displaystyle \underset{0}{\overset{d}{\int }}{X}_{1m}^2}dx+{\rho }_2{c}_2{\displaystyle \underset{d}{\overset{L}{\int }}{X}_{2m}^2}dx\right\}\end{array}}$

Los coeficientes A_n se pueden calcular de forma analítica, basta resolver integrales por partes, pero es muy tedioso. Por lo que utilizaremos procedimientos numéricos.

Una vez calculados los coeficientes A_n, se ha completado la función u(x, t). Ahora queda la temperatura de la varilla, T(x, t)=u(x, t)+T(x, ∞)

${\displaystyle \begin{array}{l}{T}_1(x,t)={\displaystyle \sum _{n=1}^{\infty }{A}_n\sin ({\mu }_1{\omega }_{n}x)\exp (-{\omega }_n^{2}t)}+\frac{T_d-T_0}{d}x+T_0,0\le x<d\\ T_2(x,t)={\displaystyle \sum _{n=1}^{\infty }{A}_n\frac{\sin ({\mu }_1{\omega }_{n}d)}{\sin \left({\mu }_2{\omega }_n(L-d)\right)}\sin \left({\mu }_2{\omega }_n(L-x)\right)\exp (-{\omega }_n^{2}t)}+\frac{T_d(L-x)+T_L(x-d)}{L-d},d\le x\le L\end{array}}$

Dado que la exponencial disminuye rápidamente con el tiempo y por otra parte, los valores de las raíces ω_n se hacen grandes, basta sumar unos pocos términos a la serie para obtener buenos resultados.

k1=1; %conductividad térmica
c1=1; %calor específico
r1=1; %densidad
k2=2;
c2=2;
r2=2;

T0=10; %temperatura, extremo izquierdo
TL=0; %extremo derecho
Tini=0; %temperatura inicial
L=1; %longitud varilla
d=L/3; %tramo izquierdo
a1=sqrt(c1*r1/k1);
a2=sqrt(c2*r2/k2);

%estado estacionario
T_d=(TL/(L-d)+k1*T0/(d*k2))/(k1/(k2*d)+1/(L-d));
m1=(T_d-T0)/d;
b1=T0;
m2=(TL-T_d)/(L-d);
b2=(T_d*L-TL*d)/(L-d);

%raíces wn
w=linspace(0,40,40);
f=@(w) cos(a1*w*d).*sin(a2*w*(L-d))+(k2*a2/(k1*a1))*cos(a2*w*(L-d)).*
sin(a1*w*d);
raiz=raices(f,w); 

%amplitudes An
A=zeros(1,length(raiz));
for n=1:length(raiz)
    w=raiz(n);
    X1=@(x) sin(a1*w*x);
    X2=@(x) sin(a1*w*d)*sin(a2*w*(L-x))/sin(a2*w*(L-d));
    Z1=@(x) X1(x).^2;
    Z2=@(x) X2(x).^2;
    norm=integral(Z1,0,d)+(r2*c2/(r1*c1))*integral(Z2,d,L);
    Y1=@(x) X1(x).*(Tini-m1*x-b1);
    Y2=@(x) X2(x).*(Tini-m2*x-b2);
    A(n)=(integral(Y1,0,d)+(r2*c2/(r1*c1))*integral(Y2,d,L))/norm; 
end

hold on
x1=linspace(0,d,50);
x2=linspace(d,L,100);
 for t=[0,0.05,0.5,1]
    u1=zeros(1,length(x1));
    u2=zeros(1,length(x2));
    for n=1:length(raiz)
        w=raiz(n);
        u1=u1+A(n)*sin(a1*w*x1)*exp(-w^2*t);
        u2=u2+A(n)*sin(a1*w*d)*sin(a2*w*(L-x2))*exp(-w^2*t)/sin(a2*w*(L-d));
    end
    T=[u1+m1*x1+b1,u2+m2*x2+b2];
    x=[x1,x2];
    plot(x,T, 'displayName',num2str(t))
end
hold off
xlabel('x')
ylabel('T')
grid on
legend('-DynamicLegend','location','best')
title('Conducción del calor, dos barras')

La distribución inicial de temperatura T(x,0)=0 en el instante t=0, requiere de muchos sumandos para que aparezca casi horizontal y luego, vertical, T(0, t)=T₀=10. Véase, desarrollo en serie de Fourier

Observamos que en el instante t=1, se ha alcanzado aproximadamente el estado estacionario

Consideremos ahora, una situación más realista que utiliza para la varilla dos sustancias, acero y aluminio que difieren en su conductividad térmica

• Tramo izquierdo, hecho de aluminio, tiene las siguientes caraterísticas: densidad ρ₁=2700, calor específico, c₁=880, conductividad térmica, K₁=209.3
• Tramo derecho, hecho de acero, tiene las siguientes caraterísticas: densidad ρ₂=7800, calor específico, c₂=460, conductividad térmica, K₂=45. Todas las magnitudes en unidades S.I.
• La longitud de la varilla es L=1, la del tramo izquierdo es d=L/3
• La temperatura en los extremos es T₀=100 y T_L=0

En primer lugar, representamos la función f(ω), para hacernos una idea de los valores de las raíces ω_n

k1=209.3; %aluminio
c1=880;
r1=2700;
k2=45; %acero
c2=460;
r2=7800;
a1=sqrt(c1*r1/k1);
a2=sqrt(c2*r2/k2);

L=1; %longitud de la varilla
d=L/3; %posición de contacto

f=@(w) cos(a1*w*d).*sin(a2*w*(L-d))+(k2*a2/(k1*a1))*cos(a2*w*(L-d)).*
sin(a1*w*d);
w=linspace(0,0.2,50);
raiz=raices(f,w); %calcula las raíces
fplot(f,[0,0.2])
xlabel('\omega')
ylabel('f(\omega)')
grid on
title('Raíces de la ecuación transcendente')

raiz =    0.0150    0.0293    0.0423    0.0551    0.0691    0.0839    
0.0990    0.1136    0.1270    0.1396    0.1532    0.1679    0.1829    0.1977

Comprobamos que las raíces calculadas (14) se corresponden con los ceros de la función f(ω) en el intervalo seleccionado (0, 0.2).

Completamos el programa para representar la distribución de temperaturas en los instantes t: 50, 500, 1000, 10000 s. Las líneas a trazos corresponden al estado estacionario. La parte derecha de la varilla, está hecha de acero y la parte izquierda, de aluminio de conductividad térmica K mucho mayor

k1=209.3; %aluminio
c1=880;
r1=2700;
k2=45; %acero
c2=460;
r2=7800;
a1=sqrt(c1*r1/k1);
a2=sqrt(c2*r2/k2);

T0=100; %temperaturas en los extremos
TL=0;
Tini=0; %temperatura inicial
L=1; %longitud varilla
d=L/3; %posición contacto

%estado estacionario
T_d=(TL/(L-d)+k1*T0/(d*k2))/(k1/(k2*d)+1/(L-d));
m1=(T_d-T0)/d;
b1=T0;
m2=(TL-T_d)/(L-d);
b2=(T_d*L-TL*d)/(L-d);

%raíces wn
w=linspace(0,0.2,50);
f=@(w) cos(a1*w*d).*sin(a2*w*(L-d))+(k2*a2/(k1*a1))*cos(a2*w*(L-d)).
*sin(a1*w*d);
raiz=raices(f,w); 

%amplitudes An
A=zeros(1,length(raiz));
for n=1:length(raiz)
    w=raiz(n);
     X1=@(x) sin(a1*w*x);
     X2=@(x) sin(a1*w*d)*sin(a2*w*(L-x))/sin(a2*w*(L-d));
    Z1=@(x) X1(x).^2;
    Z2=@(x) X2(x).^2;
    norm=integral(Z1,0,d)+(r2*c2/(r1*c1))*integral(Z2,d,L);
    Y1=@(x) X1(x).*(Tini-m1*x-b1);
    Y2=@(x) X2(x).*(Tini-m2*x-b2);
    A(n)=(integral(Y1,0,d)+(r2*c2/(r1*c1))*integral(Y2,d,L))/norm; 
end

hold on
%estado estacionario
line([0,d],[T0,T_d], 'lineStyle','--','color','k')
line([d,L],[T_d,TL],'lineStyle','--','color','k')
line([d,d],[-20,T_d],'lineStyle','--','color','k')
x1=linspace(0,d,50);
x2=linspace(d,L,100);
%estado transitorio
for t=[50,500,1000, 10000]
    u1=zeros(1,length(x1));
    u2=zeros(1,length(x2));
    for n=1:length(raiz)
        w=raiz(n);
        u1=u1+A(n)*sin(a1*w*x1)*exp(-w^2*t);
        u2=u2+A(n)*sin(a1*w*d)*sin(a2*w*(L-x2))*exp(-w^2*t)/sin(a2*w*(L-d));
    end
    T=[u1+m1*x1+b1,u2+m2*x2+b2];
    x=[x1,x2];
    plot(x,T)
end
hold off
xlabel('x')
ylabel('T')
grid on
title('Conducción del calor, dos barras')

Referencias

M. M. Smirnov. problemas de ecuaciones de la física matemática. Editorial Mir Moscú. 1976. Problema 161, enunciado, pág. 49, solución, págs. 131-132