Práctica 7. La ecuación del calor


La ecuación del calor [Graphics:images/eqcalor_gr_1.gif] modela la distribución de la temperatura [Graphics:images/eqcalor_gr_2.gif] en una barra delgada de longitud fija L. En esta práctica, estudiaremos dicha ecuación con condiciones de contorno homogéneas y no homogénas. En ambos casos, la distribución inicial de la temperatura vendrá dada por una función [Graphics:images/eqcalor_gr_3.gif].

La ecuación del calor con condiciones homogéneas

El método de resolución es separación de variables, en el que la solución [Graphics:images/eqcalor_gr_4.gif] se expresa como producto de dos funciones, una de la variable x y otra de la variable t. La ecuación se transforma entonces en un par de ecuaciones diferenciales ordinarias que, en la mayoría de los casos, se resuelven sin dificultad.

Ejemplo 1

En este ejemplo, analizaremos el caso de temperatura cero en los extremos de una varilla con [Graphics:images/eqcalor_gr_5.gif] y [Graphics:images/eqcalor_gr_6.gif]. Este problema de contorno viene dado por:

[Graphics:images/eqcalor_gr_7.gif]

Resolvemos el problema para la distribución inicial de temperatura dada por la  función

f[x_]:=12 x - 26 x^2 + 20 x^3 - 5 x^4

La gráfica de esta función es la siguiente:

Plot[f[x],{x,0,2},PlotStyle->RGBColor[1,0,0]];

[Graphics:images/eqcalor_gr_8.gif]

Sabemos que la solución es una suma infinita (en el índice n) de términos de la forma:

[Graphics:images/eqcalor_gr_9.gif]

donde los coeficientes [Graphics:images/eqcalor_gr_10.gif] están dados por la siguiente fórmula integral. Observemos que la integral se calcula numéricamente, con el comando NIntegrate. El comando Chop se usa para redondear a cero cantidades muy pequeñas.

b[n_]:= NIntegrate[f[x] Sin[n Pi x/2],{x,0,2}]//Chop;

Por ejemplo, para obtener el primer coeficiente basta pedir:

b[1]
     1.73685

A continuación, generamos una tabla con los seis primeros coeficientes [Graphics:images/eqcalor_gr_11.gif].

coefi=Table[b[n],{n,1,6}]
     {1.73685, 0, 0.890548, 0, 0.206635, 0}

Ahora presentamos estos coeficientes, junto con su índice, en forma de tabla.

Table[{n,coefi[[n]]},{n,1,6}]//TableForm
     1   1.73685
     
     2   0
     
     3   0.890548
     
     4   0
     
     5   0.206635
     
     6   0

El término n-ésimo del desarrollo en serie de la solución viene dado por:

u[x_,t_,n_]:=coefi[[n]] Sin[n Pi x/2] Exp[-t (n Pi/2)^2]

Por tanto, la aproximación dada por los seis primeros términos del desarrollo es

aprou[x_,t_]=Sum[u[x,t,j],{j,1,6}]
                 Pi x                 3 Pi x
     1.73685 Sin[----]   0.890548 Sin[------]
                  2                     2
     ----------------- + -------------------- +
            2                      2
         (Pi  t)/4            (9 Pi  t)/4
        E                    E
      
                    5 Pi x
       0.206635 Sin[------]
                      2
       --------------------
                 2
           (25 Pi  t)/4
          E

Para dibujar esta solución aproximada, construiremos una tabla de gráficos. Estos gráficos representan la solución para los valores del tiempo comprendidos entre [Graphics:images/eqcalor_gr_12.gif] y  [Graphics:images/eqcalor_gr_13.gif], usando un incremento de 0.7/11. Esto produce exactamente doce gráficos de la solución. Observemos que cuando t crece, la temperatura [Graphics:images/eqcalor_gr_14.gif] se aproxima a cero.

grafico=Table[Plot[aprou[x,t],{x,0,2},PlotRange->{0,2},
              DisplayFunction->Identity], {t,0,.7,.7/11}];
gratabla=Partition[grafico,3];
Show[GraphicsArray[gratabla]];

[Graphics:images/eqcalor_gr_15.gif]

Ejemplo 2

Analizamos el caso de extremos aislados, es decir, el problema de contorno dado por:

[Graphics:images/eqcalor_gr_16.gif]

Sean [Graphics:images/eqcalor_gr_17.gif] y [Graphics:images/eqcalor_gr_18.gif]. Consideremos la misma distribución inicial de temperatura que en el ejemplo anterior. Entonces

[Graphics:images/eqcalor_gr_19.gif]

Los coeficientes [Graphics:images/eqcalor_gr_20.gif] están dados por la siguiente fórmula integral.

a[0]=(1/2)NIntegrate[f[x],{x,0,2}]//Chop
     1.33333

a[n_]:=NIntegrate[f[x] Cos[n Pi x/2],{x,0,2}]//Chop
coeficos=Table[a[n],{n,1,6}]
     {0, 0.0321273, 0, -0.453937, 0, -0.239772}

El término n-ésimo del desarrollo en serie de la solución viene dado por

u[x_,t_,n_]:=coeficos[[n]] Cos[n Pi x/2] Exp[-t(n Pi/2)^2]

La aproximación de la solución mediante los seis primeros términos del desarrollo es

aprou[x_,t_]=a[0]+Sum[u[x,t,j],{j,1,6}]
               0.0321273 Cos[Pi x]   0.453937 Cos[2 Pi x]
     1.33333 + ------------------- - -------------------- -
                        2                       2
                      Pi  t                 4 Pi  t
                     E                     E
      
       0.239772 Cos[3 Pi x]
       --------------------
                  2
              9 Pi  t
             E
grafico2=Table[Plot[aprou[x,t],{x,0,2},PlotRange->{0,2},
              DisplayFunction->Identity], {t,0,0.1,.1/8}];             
gratabla2=Partition[grafico2,3];
Show[GraphicsArray[gratabla2]];

[Graphics:images/eqcalor_gr_21.gif]

La ecuación del calor con condiciones no homogéneas

Este problema de contorno viene dado por:

[Graphics:images/eqcalor_gr_22.gif]

siendo [Graphics:images/eqcalor_gr_23.gif] y [Graphics:images/eqcalor_gr_24.gif] las temperaturas de los extremos de la varilla.

El método de separación de variables falla si alguna de las condiciones de contorno no es homogénea. En el caso de la ecuación del calor podemos hacer alguna modificación para seguir aplicando este método. Observemos que cuando t crece la temperatura se aproxima a un estado estacionario que sólo depende de x, digamos [Graphics:images/eqcalor_gr_25.gif].
Entonces tenemos que [Graphics:images/eqcalor_gr_26.gif] donde
[Graphics:images/eqcalor_gr_27.gif]
y [Graphics:images/eqcalor_gr_28.gif] verifica el siguiente problema de contorno:

[Graphics:images/eqcalor_gr_29.gif]

Ejercicio

Analizar  el problema de contorno propuesto, para los valores [Graphics:images/eqcalor_gr_30.gif]
y la distribución inicial de temperatura dada por [Graphics:images/eqcalor_gr_31.gif] Entonces,

s[x_]:=x/2; f[x_]:=x^2;
Plot[{f[x],s[x]},{x,0,1}, PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}];

[Graphics:images/eqcalor_gr_32.gif]

A continuación, definimos los coeficientes   [Graphics:images/eqcalor_gr_33.gif] usando NIntegrate y Chop, con

a[n_]:=2 NIntegrate[(f[x]-s[x]) Sin[n Pi x], {x,0,1}]//Chop;

Los seis primeros coeficientes son:

coe=Table[a[n],{n,1,6}]
     {0.0602976, -0.159155, 0.0965473, -0.0795775, 0.0615979, 
      
       -0.0530516}

Ahora definimos el término n-ésimo de la parte transitoria de la temperatura en [Graphics:images/eqcalor_gr_34.gif],

v[x_,t_,n_]:=coe[[n]] Sin[n Pi x] Exp[- t(n Pi)^2]

La suma de los seis primeros términos se obtiene con

vapro[x_,t_]=Sum[v[x,t,j],{j,1,6}]
     0.0602976 Sin[Pi x]   0.159155 Sin[2 Pi x]
     ------------------- - -------------------- +
              2                       2
            Pi  t                 4 Pi  t
           E                     E
      
       0.0965473 Sin[3 Pi x]   0.0795775 Sin[4 Pi x]
       --------------------- - --------------------- +
                  2                        2
              9 Pi  t                 16 Pi  t
             E                       E
      
       0.0615979 Sin[5 Pi x]   0.0530516 Sin[6 Pi x]
       --------------------- - ---------------------
                   2                       2
              25 Pi  t                36 Pi  t
             E                       E

La solución aproximada de la ecuación del calor con condiciones no homogéneas es la suma de la solución estacionaria y la solución transitoria. Se define con uapro.

uapro[x_,t_]=s[x]+vapro[x,t]
     x   0.0602976 Sin[Pi x]   0.159155 Sin[2 Pi x]
     - + ------------------- - -------------------- +
     2            2                       2
                Pi  t                 4 Pi  t
               E                     E
      
       0.0965473 Sin[3 Pi x]   0.0795775 Sin[4 Pi x]
       --------------------- - --------------------- +
                  2                        2
              9 Pi  t                 16 Pi  t
             E                       E
      
       0.0615979 Sin[5 Pi x]   0.0530516 Sin[6 Pi x]
       --------------------- - ---------------------
                   2                       2
              25 Pi  t                36 Pi  t
             E                       E

Por último, dibujaremos la solución aproximada para valores del tiempo comprendidos entre [Graphics:images/eqcalor_gr_35.gif] y  [Graphics:images/eqcalor_gr_36.gif], usando un incremento de [Graphics:images/eqcalor_gr_37.gif]. Observemos que cuando t crece la solución se aproxima a la función que describe el estado estacionario de temperatura.

grafico3=Table[Plot[uapro[x,t],{x,0,1},
              PlotRange->{0,1}, PlotStyle->RGBColor[0,0,1],
              DisplayFunction->Identity], {t,0,.1,.1/8}];
gratabla3=Partition[grafico3,3];
Show[GraphicsArray[gratabla3]];

[Graphics:images/eqcalor_gr_38.gif]

La gráfica de la solución aproximada   [Graphics:images/eqcalor_gr_39.gif] se obtiene con

Plot3D[uapro[x,t],{x,0,1},{t,0,.1},
Boxed->True,ViewPoint->{-1,3,1}];

[Graphics:images/eqcalor_gr_40.gif]


Creada por MBilbao Lab     27 de Septiembre de 2001