The One-Dimensional Heat Equation

Part 1. Equal Boundary Conditions

Enter the next line to define and see a graph of the initial temperature function f.  We plot two copies, one in red for immediate display, and one in blue for future comparision with our Fourier approximation.  Since we will be plotting graphs with similar names, we suppress the spelling warnings.

Off[General::spell1];

f[x_] :=  
        Which[x <= 0.5, 10*x, x > 0.5  &&  x < 1.5, 5 - 7*(x - 0.5), 
          x >= 1.5, -2 + 4*(x - 1.5)];
graph1a  =  
        Plot[f[x], {x, 0, 2}, 
          PlotStyle -> RGBColor[1, 0, 0]];
graph1b  =  
        Plot[f[x], {x, 0, 2}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];

The next lines define and calculate the first 15 coefficients of the Fourier sine expansion for f.  Note the we have expressed f in terms of UnitStep.  Convince yourself that this expression is the same as f defined earlier.  This is necessary so that Mathematica will carry out the integration symbolically.  We will take this approach for all subsequent Fourier coefficient calculations.

f[x_] :=  
        10*x*UnitStep[1/2 - x] + (5 - 7*(x - 1/2))*
            UnitStep[-(x - 1/2) (x - 3/2)] + (-2 + 
                4*(x - 3/2))*UnitStep[x - 3/2];

Integrate[f[x]*Sin[k*(Pi/2)*x], {x, 0, 2}];
 b[k_] = %

Table[b[k], {k, 15}]

Compare the Fourier approximations with the original function.

s[x_] :=  Sum[b[k]*Sin[k*(Pi/2)*x], {k, n}];

n  =  1;   
graph2a  =  
        Plot[s[x], {x, 0, 2}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
Show[graph2a, graph1b, 
        DisplayFunction -> $DisplayFunction];

We define the approximation to the solution of the initial/boundary value problem as a function u[x] with t as a parameter.

Clear[n];
u[x_] :=  
        Sum[b[k]*Sin[k*(Pi/2)*x]*Exp[-k^2*(Pi/2)^2*t], {k, 
            n}];

Increase t from t = 0.001 to t = 2 to see snapshots in time of the temperature evolution.

n  =  15;
t  =  0.001;
graph3  =  
        Plot[u[x], {x, 0, 2}, 
          PlotStyle -> RGBColor[1, 0, 0]];

Part 2. Changing the Diffusivity Constant

Copy command lines from Part 1 into the following input cell, and modify them as needed to carry out steps 1 and 2.

Answer to step 3:

Part 3. An Initial Condition Not Compatible with the Boundary Conditions

We look at an initial condition u[x,0] = g[x] such that g[0] and g[2] do not match up with the boundary conditions.

g[x_] :=  
        Which[x <= 1, 1 - 3*x, x > 1, -2 + 5*(x - 1)];

graph4a  =  
        Plot[g[x], {x, 0, 2}, 
          PlotStyle -> RGBColor[1, 0, 0]];
graph4b  =  
        Plot[g[x], {x, 0, 2}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];

Again we calculate the first 15 terms in the half-range sine expansion.

g[
          x_] :=  (1 - 3*x)*
            UnitStep[1 - x] + (-2 + 5*(x - 1))*
            UnitStep[x - 1];

Integrate[g[x]*Sin[k*(Pi/2)*x], {x, 0, 2}];
b[k_] = %

Table[b[k], {k, 15}]

Compare the Fourier approximations with the original function.

s[x_] :=  Sum[b[k]*Sin[k*(Pi/2)*x], {k, n}];

\(n\  = \ 1; \  (*\ An\ Example\ *) 

    graph5a\  = \ 
      Plot[s[x], {x, 0, 2}, PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity];\)
Show[graph5a, graph4b, 
        DisplayFunction -> $DisplayFunction];

Clear[n];
u[x_] :=  
        Sum[b[k]*Sin[k*(Pi/2)*x]*Exp[-k^2*(Pi/2)^2*t], {k, 
            n}];

Now we look at time snapshots of the evolution of the temperature distribution for this model.

n  =  15;
t  =  0.001;
graph6  =  
        Plot[u[x], {x, 0, 2}, 
          PlotStyle -> RGBColor[1, 0, 0]];

Part 4. Unequal Boundary Conditions

The next line defines and graphs the initial condition function h.

h[x_] :=  
        Which[x <= 0.5, 1 + 8*x, x > 0.5  &&  x < 1.5, 5 - 8*(x - 0.5), 
          x >= 1.5, -3 + 2*(x - 1.5)];

graph4a  =  
        Plot[h[x], {x, 0, 2}, 
          PlotStyle -> RGBColor[1, 0, 0]];
graph4b  =  
        Plot[h[x], {x, 0, 2}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];

Define w and graph h[x]-w[x].

w[x_] :=  1 - (3/2)*x;

graph5a  =  
        Plot[h[x] - w[x], {x, 0, 2}, 
          PlotStyle -> RGBColor[1, 0, 0]];
graph5b := 
        Plot[h[x] - w[x], {x, 0, 2}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];

The next lines define and calculate the first 15 coefficients for the half-range sine expansion of h-w.

h[
          x_] :=  (1 + 8*x)*
            UnitStep[1/2 - x] + (5 - 8*(x - 1/2))*
            UnitStep[-(x - 1/2) (x - 3/2)] + (-3 + 
                2*(x - 3/2))*UnitStep[x - 3/2];

Integrate[(h[x] - w[x])*Sin[k*(Pi/2)*x], {x, 0, 
          2}];
B[k_] = %

Table[B[k], {k, 15}]

Increase n from 1 to 15 and compare the Fourier sine approximation with the graph of h[x]-w[x].

s[x_] :=  Sum[B[k]*Sin[k*(Pi/2)*x], {k, n}];

n  =  1;   
graph6  =  
        Plot[s[x], {x, 0, 2}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
Show[graph5b, graph6, 
        DisplayFunction -> $DisplayFunction];

The next lines define the approximation U, again as a function of x with t as a parameter.  We will use a fifteen-term approximation.

U[x_] :=  
        Sum[B[k]*Sin[k*(Pi/2)*x]*Exp[-k^2*t], {k, n}] + 
          w[x];

n  =  15;

Increase t from 0.001 to 3, and examine the "time snapshots."

\(t\  = \ 0.001; \  (*\ An\ Example\ *) 

    graph7\  = \ 
      Plot[U[x], {x, 0, 2}, PlotStyle -> RGBColor[1, 0, 0]];\)

Now we consider unequal boundary conditions with an initial condition that is not compatible.  Our new initial condition is u[x,0] = h2[x], where h2 has the graph given below.

h2[x_] :=  
        Which[x <= 0.5, -4*x, 
          x > 0.5  &&  x < 1.5, -2 + 6*(x - 0.5), x >= 1.5, 
          4 - 6*(x - 1.5)];

graph8  =  
        Plot[h2[x], {x, 0, 2}, 
          PlotStyle -> RGBColor[1, 0, 0]];
graph8b  =  
        Plot[h2[x], {x, 0, 2}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];

We compare the initial condition with the steady-state.

graph9  =  
        Plot[w[x], {x, 0, 2}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[graph8, graph9, 
        DisplayFunction -> $DisplayFunction];

Now we calculate the approximate solution as we did above.  First, the coefficients for the Fourier sine approximations:

h2[
          x_] :=  (-4*x)*
            UnitStep[1/2 - x] + (-2 + 6*(x - 1/2))*
            UnitStep[-(x - 1/2) (x - 3/2)] + (4 - 
                6*(x - 3/2))*UnitStep[x - 3/2];

Integrate[(h2[x] - w[x])*Sin[k*(Pi/2)*x], {x, 0, 
          2}];
B[k_] = %

Table[B[k], {k, 15}]

Now we define the approximate solution function U[x,t].

U[x_] :=  
        Sum[B[k]*Sin[k*(Pi/2)*x]*Exp[-k^2*t], {k, n}] + 
          w[x];

Compare the time snapshots of the approximate solution for small time with the initial condition.

n  =  15;

t  =  0.001;
graph10  =  
        Plot[U[x], {x, 0, 2}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
Show[graph8b, graph10, 
        DisplayFunction -> $DisplayFunction];

Now examine the time snapshots of the approximate solution for larger times t.

\(t\  = \ 0.1; \  (*\ An\ Example\ *) 

    graph11\  = \ 
      Plot[U[x], {x, 0, 2}, PlotStyle -> RGBColor[1, 0, 0]];\)

Part 5. Summary

Answer to the questions in Part 5: