PDE solution by separation of variables

Assume a product solution

First, we assume a product solution of the form u_p(t,x) = T(t)X(x) where T and X are unknown functions on t > 0 and x ∈ [0,L].

Separate PDE

Second, we substitute the product solution into [@eq:1D_heat_ex_pde] and separate variables: $$\begin{aligned} T' X &= k T X'' \Rightarrow \\ \frac{T'} {k T} &= \frac{X''} {X}. \end{aligned}$$ So it is separable! Note that we chose to group k with T, which was arbitrary but conventional.

Set equal to a constant

Since these two sides depend on different independent variables (t and x), they must equal the same constant we call  − λ, so we have two ODEs: $$\begin{aligned} \frac{T'} {k T} &= -\lambda \quad \Rightarrow T' + \lambda k T = 0 \\ \frac{X''} {X} &= -\lambda \quad \Rightarrow X'' + \lambda X = 0. \end{aligned}$$

Solve the boundary value problem

The latter of these equations with the boundary conditions [@eq:1D_heat_ex_bcs] is precisely the same sturm-liouville boundary value problem from [@ex:sturm_liouville1], which had eigenfunctions

with corresponding (positive) eigenvalues $$\begin{aligned} \lambda_n &= \left(\frac{n \pi} {L}\right)^2. \end{aligned}$$

Solve the time variable ODE

The time variable ODE is homogeneous and has the familiar general solution $$\begin{aligned} T(t) = c e^{-k \lambda t} \end{aligned}$$ with real constant c. However, the boundary value problem restricted values of λ to λ_n, so $$\begin{aligned} T_n(t) = c e^{-k \left(n \pi/L\right)^2 t}. \end{aligned}$$

Construct the product solution

The product solution is $$\begin{aligned} u_p(t,x) &= T_n(t) X_n(x) \nonumber \\ &= c e^{-k \left(n \pi/L\right)^2 t} \sin\left(\frac{n \pi} {L} x\right). \end{aligned}$$ This is a family of solutions that each satisfy only exotically specific initial conditions.

Apply the initial condition

The initial condition is u(0,x) = f(x). The eigenfunctions of the boundary value problem form a fourier series that satisfies the initial condition on the interval [0,L] if we extend f to be periodic and odd over x [@kreyszig2011 p. 550]; we call the extension f^*. The odd series synthesis can be written $$\begin{aligned} f^*(x) = \sum_{n=1}^\infty b_n \sin\left(\frac{n \pi} {L} x\right) \end{aligned}$$ where the fourier analysis gives $$\begin{aligned} \label{eq:1D_heat_ex_fourier} b_n &= \frac{2} {L} \int_0^L f^*(\chi) \sin\left(\frac{n \pi} {L} \chi\right). \end{aligned}$$ So the complete solution is $$\begin{aligned} u(t,x) &= \sum_{n=1}^\infty b_n e^{-k \left(n \pi/L\right)^2 t} \sin\left(\frac{n \pi} {L} x\right). \end{aligned}$$ Notice this satisfies the PDE, the boundary conditions, and the initial condition!

Plotting solutions

If we want to plot solutions, we need to specify an initial condition u(0,x) = f^*(x) over [0,L]. We can choose anything piecewise continuous, but for simplicity let’s let The odd periodic extension is an odd square wave. The integral [@eq:1D_heat_ex_fourier] gives $$\begin{aligned} b_n &= \frac{4} {n \pi} (1 - \cos(n\pi)) \nonumber \\ &= \begin{cases} 0 & n\text{ even} \\ \frac{4} {n\pi} & n\text{ odd}. \end{cases} \end{aligned}$$ Now we can write the solution as $$\begin{aligned} u(t,x) &= \sum_{n=1,\, n \text{ odd}}^\infty \frac{4} {n\pi} e^{-k \left(n \pi/L\right)^2 t} \sin\left(\frac{n \pi} {L} x\right). \end{aligned}$$

Plotting in Python

First, load some Python packages.

import numpy as np
import matplotlib.pyplot as plt

Set k = L = 1 and sum values for the first N terms of the solution.

L = 1
k = 1
N = 100
x = np.linspace(0,L,300)
t = np.linspace(0,2*(L/np.pi)**2,100)
u_n = np.zeros([len(t),len(x)])
for n in range(N):
  n = n+1 # because index starts at 0
  if n % 2 == 0: # even
    pass # already initialized to zeros
  else: # odd
    u_n += 4/(n*np.pi)*np.outer(
      np.exp(-k*(n*np.pi/L)**2*t),
      np.sin(n*np.pi/L*x)
    )

Let’s first plot the initial condition.

fig, ax = plt.subplots()
ax.plot(x,u_n[0,:])
plt.xlabel('space $x$')
plt.ylabel('$u(0,x)$')
plt.draw()

Now we plot the entire response.

fig, ax = plt.subplots()
plt.contourf(t,x,u_n.T)
c = plt.colorbar()
c.set_label('$u(t,x)$')
plt.xlabel('time $t$')
plt.ylabel('space $x$')
plt.show()

We see the diffusive action proceeds as we expected.