Solving the 1D heat equation for a symmetric rod

We’ve derived the 1D heat diffusion equation, so now let’s take a look at one of the solutions. This solution is for a symmetric rod. Simply, imagine a rod of some length $L$. This rod is insulated across its length so that no heat is lost to the atmosphere across that length. The only place where the rod is not insulated is at the ends which will change temperature. We want to calculate how the temperature changes across the length of the rod in time when the temperature at the ends of the rod change.

Let’s recall the heat diffusion equation that we derived:

We want to find a solution for the temperature profile of the rod in time \textit{and} space. This would take the form of an equation:

We want to find an equation for $T(x,t)$ where we can plug in the location within the rod and the time after the temperature change to find the temperature.

We will assume that the entire rod starts off at a uniform temperature $T_i$ at $t=0$, that is to say: $T_i=T(x,0)$. Then at an infinitesimal time after $t=0$, the edges ends of the rod immediately change temperature to $T_f$. Therefore $T_f = T(0,t>0) = T(L,t>0)$. One can imagine that after a long time the rod will reach one uniform temperature, because the side of the rod is insulated and no heat can be lost, the entire rod must reach thermal equilibrium with the rod ends. Therefore we can also say: $T_f = T(x,\infty)$

Notice we have just defined our boundary conditions. Let’s take another look:

Our temporal (initial) condition is the temperature of the ends of the rod at $t=0$ And the two spacial boundary conditions refer to the two ends of the rod: .

We can rewrite these another way which may make more sense:

  • At $t=0$, $T=T_i$ (temporal)
  • At $x=0$, $T=T_f$ (spacial)
  • At $x=L$, $T=T_f$ (spacial)

In order to solve the differential equation for this problem it helps to non-dimensionalise our variables, which are $x$, time and temperature. We define our non-dimensional variables as follows:

I will not discuss how to non-dimensionalise the variables here, instead you can see \textbf{this} article.

The PDE now looks like this:

Why did we just do this? It makes solving the PDE mathematically easier because now the initial temperature is $\hat{T_i} = 1$ and the final temperature is $\hat{T_f} = 0$, the initial length is $\hat{x} = 0$, the final length is $\hat{x}=1$ and the initial time is $\hat{t} = 0$.

Let’s rewrite our boundary conditions in non-dimensional form:

  • At $\hat{t}=0$, $\hat{T}=1$ (temporal)
  • At $\hat{x}=0$, $T=0$ (spacial)
  • At $\hat{x}=1$, $T=0$ (spacial)

We will now use a technique called \textbf{separation of variables} to solve this PDE in non-dimensional form. This means that we are going to assume that there exists a solution to the differential which is the product of two functions which each are a function only of either space or time.

Let’s unpack that. Recall that we want a solution which describes temperature in space and time like so: $\hat{T} = \hat{T}(\hat{x},\hat{t})$. We are going to assume that we can write $\hat{T}(\hat{x},\hat{t})$ (which is a two-variable function) as the product of two one-variable functions like so:

Where $\beta (\hat{x})$ is some function of $\hat{x}$ alone, and $\gamma(\hat{t})$ is a function only of $\hat{t}$. Multiplying these two functions together gives our two-variable function which is our solution.

We are now going to look at the left hand side of our PDE: $\frac{\partial^2\hat{T}}{\partial \hat{x}^2}$, notice that this is what we get if we differentiate $\hat{T}(\hat{x},\hat{t})$ twice with respect to $\hat{x}$ whilst keeping $\hat{t}$ constant:

Because $\gamma(\hat{t})$ is being treated as a constant, we can easily apply the product rule when differentiating, and so we get

Where $\beta’’(\hat{x}) = \frac{d^2 \beta(\hat{x})}{d\hat{x}^2}$. Notice it’s a normal derivative because $\beta$ can only be differentiated with respect to $\hat{x}$ because it is a function of $\hat{x}$ only.

Moving on to the right hand side of the PDE we see that

Again applying the product rule we get

Substituting both of these results into the original PDE we get an interesting equation composed only of single-variable functions:

Which can be re-written as:

Why did I equate this equation to a constant $-\lambda^2$? Well we know that each side of the equation is actually a constant because the left hand side of the equation is a function of $\hat{t}$ only, meaning that changing that side cannot affect the right hand side of the equation which is a function only of $\hat{x}$ and vice versa. Okay then, well why is the constant negative? Why is it squared? You will see further on why choosing a negative constant is the only one which can give you an answer which makes sense in the real world, and squaring it simply makes the maths more convenient.

Anyway - we now have two equations which are easy to solve independently:

  • $\frac{\gamma'(\hat{t})}{\gamma(\hat{t})} = -\lambda^2$
  • $\frac{\beta''(\hat{x})}{\beta (\hat{x})} = -\lambda^2$

Let’s solve the first by re-arranging and integrating directly:

Combining the constant and exponentiating:

Integrating the second function is a little bit more involved. Let’s rearrange it:

This is a \textbf{homogeneous 2nd order ordinary differential equation with constant coefficients}, this means we can guess a solution of the form $\beta = e^{k\hat{x}}$. Let’s plug this solution into the ODE:

Notice how the exponential terms drop out and we are left with:


We can now take linear combinations of the solution as we have two of them (the $\pm i \lambda$)

We can write

Recall Euler’s identity which states that $e^{i\theta}=\cos(\theta)+i\sin(\theta)$

We can substitute this into our general solution for $\beta$ and we get the following:

Noting sin and cos properties:

Collecting terms:

And collecting constants into new ones:

We now have a general solution for $\hat{T}$!

Now we must substitute in our boundary conditions to find the final solution.

Let’s apply the boundary conditions: At $\hat{x}=0$, $\hat{T}=0$ (spacial)

For this to be satisfied, $C_2 = 0$. Therefore we now have

$\hat{x}=1$, $\hat{T}=0$ (spacial)

Therefore $C_3\sin(\lambda) = 0$ because the exponential cannot equal 0. This is true when $\lambda=n\pi$. Taking linear combinations of solutions we get

And finally: At $\hat{t}=0$, $\hat{T}=1$ (temporal)

Therefore $C_1 = 1$. We get a final solution of:

Now we must find the the constants $C_n$, we will use the property that:

Now we apply the temporal (initial) boundary condition:

Multiplying each side of this equation by $\sin(n\pi \hat{x})$ and integrating between limits yields:

Notice that the only terms in this infinite sum that will not be zero integrated are when $n=\hat{x}$, therefore we get a formula for the unknown coefficients:

For the denominator recall that $\sin^2(x)=\frac{1}{2}[1-\cos(2x)]$ and therefore

If you evaluate the right hand integral you find that:


So it can be seen that the coefficients are:

Note that $C_1$ is not the same constant mentioned earlier.

Substituting these coefficients into the final solution we get our (more) final solution: