Heat Equation

Heat Equation with homogeneous boundary conditions.

Consider the heat equation

a D[u[x,t],x,x] == D[u[x,t],t]; 0 < x < 20; t > 0;

subject to the boundary conditions

u[ 0,t] == 0;
u[20,t] == 0;  t > 0

and the initial condition

u[x,0] == 100; 0 < x < 20

Use Fourier series to find an expression for u[x,t]. For each of the following three materials, a) silver (a==1.71), b) aluminum (a==0.86), and c) iron (a==0.12), use the first five non-zero terms of the series to do the following:
1. Plot u[x,t] as a function of x and t
2. Plot u[x,c] as a function of x for various choices of t==c. Compare these plots to similar plots using only the first term of the series. When do the plots start to coincide?
3. Plot the temperature function of the center of the bar, u[10,t], as a function of t. Determine how long it takes the center of the bar to cool to 25. Compare this answer to that obtained using only one term of the series.

Solution

Let's make the constant a and the number of terms used in the series approximation parameters to the function u:

u[x_,t_,a_,k_] =
Sum[b[n] Sin[n Pi x/20] E^(-(n Pi/20)^2 a t),{n,1,k}]

The coefficients are:

b[n_] = 2/20 Integrate[100 Sin[n Pi x/20],{x,0,20}]//Together

Since the even coefficients are 0, the first five non-zero terms are given by k=9.

u[x,t,a,9]

Let's set the constants used for a:

silv = 1.71;
alum = 0.89;
iron = 0.12;

a) Silver

Plot3D[u[x,t,silv,9],{x,0,20},{t,0,30}]

An animation of the temperature distribution as a function of time:

Do[
    Plot[u[x,t,silv,9], {x,0,20}, PlotRange->{0,110}],
{t,0,30,3}]

After 10 seconds, the first term of the series (red curve) gives a good approximation to first five non-zero terms of the series (black curve):

Do[
    Plot[{u[x,t,silv,1],u[x,t,silv,9]}, {x,0,20},
    PlotStyle->{RGBColor[1,0,0], RGBColor[0,0,0]},
    PlotRange->{0,110}],
{t,0,10,2}]

u[10,10,silv,1]//N
u[10,10,silv,9]//N

(u[10,10,silv,1]-u[10,10,silv,9])/u[10,10,silv,9]//N

A plot of the temperature at the center of the bar as a function of time:

Plot[u[10,t,silv,9], {t,0,30}, PlotRange->{0,105}]

To reach 25:

FindRoot[u[10,t,silv,9] == 25, {t,20}]

We get almost the same answer using only the first term of the series:

FindRoot[u[10,t,silv,1] == 25, {t,20}]

b) Aluminum

Plot3D[u[x,t,alum,9],{x,0,20},{t,0,30}]

It takes longer, about 20 seconds, for the first term of the series (red curve) to give a good approximation to the first five non-zero terms of the series (black curve):

Do[
    Plot[{u[x,t,alum,1],u[x,t,alum,9]}, {x,0,20},
    PlotStyle->{RGBColor[1,0,0], RGBColor[0,0,0]},
    PlotRange->{0,110}],
{t,0,20,4}]

u[10,20,alum,1]//N
u[10,20,alum,9]//N

(u[10,20,alum,1]-u[10,20,alum,9])/u[10,20,alum,9]//N

The temperature at the center of the bar drops off more slowly:

Plot[u[10,t,alum,9], {t,0,30}, PlotRange->{0,105}]

To reach 25:

FindRoot[u[10,t,alum,9] == 25, {t,20}]

We get the same answer using only the first term of the series:

FindRoot[u[10,t,alum,1] == 25, {t,20}]

c) Iron

Plot3D[u[x,t,iron,9],{x,0,20},{t,0,30}]

It takes even longer, around 140 seconds, for the first term of the series (red curve) to give a good approximation to the first five non-zero terms of the series (black curve):

Do[
    Plot[{u[x,t,iron,1],u[x,t,iron,9]}, {x,0,20},
    PlotStyle->{RGBColor[1,0,0], RGBColor[0,0,0]},
    PlotRange->{0,110}],
{t,0,100,20}]

u[10,140,iron,1]//N
u[10,140,iron,9]//N

(u[10,140,iron,1]-u[10,140,iron,9])/u[10,140,iron,9]//N

Plot[u[10,t,iron,9], {t,0,140}, PlotRange->{0,105}]

To reach 25:

FindRoot[u[10,t,iron,9] == 25, {t,20}]

We get almost the same answer using only the first term of the series:

FindRoot[u[10,t,iron,1] == 25, {t,20}]

Heat Equation with non-homogeneous boundary conditions.

Consider the heat equation

a D[u[x,t],x,x] == D[u[x,t],t]; 0 < x < 20; t > 0;

subject to the boundary conditions

u[ 0,t] == 10;
u[20,t] == 110;  t > 0

and the initial condition

u[x,0] == 50; 0 < x < 20

Use Fourier series plus a steady state term to find an expression for u[x,t]. Assuming the material is aluminum (a=0.86), use the first ten non-zero terms of the series to do the following:
1. Plot u[x,t] as a function of x and t
2. Plot u[x,c] as a function of x for various choices of t==c. Compare these plots to similar plots using only the first term of the series. When do the plots start to coincide?
3. Plot u[5,t] as a function of t. Why does the temperature drop below the steady state temperature 35 at x==5 before climbing back up to this temperature? Determine how long it takes the point x==5 to reach the its steady state temperature of 35. Compare this answer to that obtained using only one term of the series.

Solution

The steady state temperature is a straight line:

v[x_] = (110-10)/20 x + 10

With non-homogeneous boundary conditions, we look for a solution of the form
    u[x_,t_] = v[x] + w[x,t]
where w[x,t] satisfies the same heat equation with homogeneous boundary conditions and initial temperature distribution
    w[x,0] = 50 - v[x]

As before, we define the series approximation

w[x_,t_,a_,k_] =
Sum[b[n] Sin[n Pi x/20] E^(-(n Pi/20)^2 a t),{n,1,k}]

where the coefficients are:

b[n_] = 2/20 Integrate[(50-v[x])Sin[n Pi x/20],{x,0,20}]/.
         Sin[n Pi]->0//Together

The series approxiamtion to the original problem is then

u[x_,t_,a_,k_] = v[x] + w[x,t,a,k]

Since none of the coefficients are 0, the first ten non-zero terms of the series are given by k=10.

u[x,t,a,10]

alum = 0.86;

Plot3D[u[x,t,alum,10],{x,0,20},{t,0,30},BoxRatios->{1,1,1}]

After 50 seconds, the first non-zero term of the series (red curve) gives a good approximation to first five non-zero terms of the series (black curve):

Do[
    Plot[{u[x,t,alum,1],u[x,t,alum,10]}, {x,0,20},
    PlotStyle->{RGBColor[1,0,0], RGBColor[0,0,0]},
    PlotRange->{0,115}],
{t,0,50,5}]

u[5,50,alum,1]//N
u[5,50,alum,10]//N

(u[5,50,alum,1]-u[5,50,alum,10])/u[5,50,alum,10]//N

Plot[u[5,t,alum,10], {t,0,100}, PlotRange->{20,60}]

The temperature dips below the steady state temperature, 35, before asymptotically approaching it from below. The reason appears to be that the starting temperature and position (x==5) is nearer the cold end and so cools off before the heat from the warm end arrives and brings the temperature back up to its limiting value.

We can use FindRoot[] to determine the time it takes for the point x==5 to arrive within 1% of its steady state temperature 35:

FindRoot[u[5,t,alum,10] == 35-(0.35), {t,100}]

This is essentially the same answer we would get using only one term of the series:

FindRoot[u[5,t,alum,1] == 35-(0.35), {t,100}]