# Heat Equation and Schrödinger Equation

There is an intriguing relationship between Schrödinger equation for a free particle and homogeneous heat equation.

1-dimentional Schrödinger equation for a free particle is
\begin{equation}
\label{eq:se}
i\hbar\frac{\partial\psi(x,t)}{\partial t}=-\frac{\hbar^2}{2m}\frac{\partial^2\psi(x,t)}{\partial x^2}.
\end{equation}
Take the Wick rotation $t\mapsto \tau=it$. Then the Schrödinger equation \eqref{eq:se} turns into
\begin{equation}
\label{eq:hheq}
\frac{\partial\phi(x,\tau)}{\partial\tau}=\frac{\hbar}{2m}\frac{\partial^2\phi(x,\tau)}{\partial x^2},
\end{equation}
where $\phi(x,\tau)=\psi\left(x,\frac{\tau}{i}\right)$. \eqref{eq:hheq} is a homogeneous heat equation with diffusion coefficient $\alpha^2=\frac{\hbar}{2m}$. Conversely, apply the Wick rotation $t\mapsto\tau=-it$ to the 1-dimensional homogeneous heat equation
\begin{equation}
\label{eq:hhe2}
\frac{\partial u(x,t)}{\partial t}=\alpha^2\frac{\partial^2 u(x,t)}{\partial x^2}.
\end{equation}
Then the resulting equation is
\begin{equation}
\label{eq:se2}
i\hbar\frac{\partial w(x,\tau)}{\partial t}=-\alpha^2\hbar\frac{\partial^2 w(x,\tau)}{\partial x^2},
\end{equation}
where $w(x,\tau)=u\left(x,-\frac{\tau}{i}\right)$. \eqref{eq:se2} is a Schrödinger equation for a free particle with $m=\frac{1}{2\alpha^2}$. The solution of \eqref{eq:hhe2} with homogeneous boundary conditions takes the form
$$u(x,t)=\sum_{n=0}^\infty A_ne^{-\lambda_n^2\alpha^2 t}X_n(x).$$
Its Wick rotated solution is
\begin{align*}
w(x,\tau)&=\sum_{n=0}^\infty w_n(x,\tau)\\
&=\sum_{n=0}^\infty A_ne^{-i\lambda_n^2\alpha^2\tau}X_n(x).
\end{align*}
$i\hbar\frac{\partial w(x,\tau)}{\partial\tau}=\lambda_n^2\alpha^2\hbar w(x,\tau)$. Thus for each $n=0,1,2,\cdots$, $E_n=\lambda_n^2\alpha^2\hbar$ is the energy and $\omega_n=\lambda_n^2\alpha^2$ is the frequency of the wave $w_n(x,\tau)$.

# The Semi-Homogeneous Heat Problem

In this lecture, we study how to solve the semi-homogeneous heat problem:
\begin{aligned} &u_t=\alpha^2 u_{xx}+F(x,t),\ 0<x<L,\ t>0\\ &{\rm BCs:}\ \left\{\begin{aligned} -k_1u_x(0,t)+h_1u(0,t)&=0,\ t>0\\ k_2u_x(L,t)+h_2u(L,t)&=0, \end{aligned} \right.\\ &{\rm IC:}\ u(x,0)=\phi(x),\ 0<x<L\ \ \ \ \ (1) \end{aligned}

We consider a solution of the form
$$u(x,t)=\sum_{n=0}^\infty T_n(t)X_n(x)$$
$$F(x,t)=\sum_{n=0}^\infty F_n(t)X_n(x)$$
Let $A_n:=T_n(0)$. Then
$$\phi(x)=\sum_{n=0}^\infty A_nX_n(x)$$
For each $n=0,1,2,\cdots$, $X_n(x)$ is the solution of the Sturm-Liouville problem:
\begin{aligned} &X_n^{\prime\prime}=-\lambda_n^2X_n\\ &{\rm BCs:}\ \left\{\begin{aligned} -k_1X_n'(0)+h_1X_n(0)&=0\\ k_2X_n'(L)+h_2X_n(L)&=0 \end{aligned} \right.\ \ \ \ \ (2) \end{aligned}
and $T_n(t)$ satisfies the first order ODE with initial condition:
\left\{\begin{aligned} \dot{T}_n+\lambda_n^2\alpha^2T_n=F_n\\ T_n(0)=A_n. \end{aligned} \right.\ \ \ \ \ (3)
Then $T_n(t)$ is given by
$$T_n(t)=\left[A_n+\int_0^t F_n(s)e^{\lambda_n^2\alpha^2 s}ds\right]e^{-\lambda_n^2\alpha^2 t}.$$
Finally, $F_n(t)$ and $A_n$ are to be determined by
\begin{align*}
F_n(t)&=\frac{1}{L_n}\int_0^L F(x,t)X_n(x)dx\ \ \ \ \ (3)\\
A_n&=\frac{1}{L_n}\int_0^L\phi(x)X_n(x)dx\ \ \ \ \ (4)
\end{align*}
As a summary we have the recipe for solving the semi-homogeneous heat IBVP (1):

1. Find the eigenvalue $\lambda_n$ and the corresponding eigenfunction $X_n$ by solving the Sturm-Liouville problem (2)
2. Find the coefficients $F_n(t)$ and $A_n$ using the formulae (4) and (5) respectively.
3. For each $n=0,1,2,\cdots$, find $T_n(t)$ using the formula (3)
4. Put all things together as $$u(x,t)=\sum_{n=0}^\infty T_n(t)X_n(x)$$

Example. Solve the semi-homogeneous heat IBVP:
\begin{align*}
&u_t=\alpha^2 u_{xx}+b\sin\pi x,\ 0<x<1,\ 0<t<\infty\\
&{\rm BCs:}\ \left\{\begin{aligned}
u(0,t)&=0,\ 0<t<\infty\\
u(1,t)&=0,
\end{aligned}
\right.\\
&{\rm IC:}\ u(x,0)=1,\ 0<x<1
\end{align*}

Solution:

Step 1. The eigenfunctions are
$$X_n(x)=\sin n\pi x,\ n=1,2,3,\cdot.$$

Step 2. \begin{align*}
F_n(t)&=2b\int_0^1\sin\pi x\sin n\pi xdx\\
&=\left\{\begin{array}{ccc}
b & {\rm if} & n=1\\
0 & {\rm otherwise}.\end{array}\right.\\
A_n&=2\int_0^1\sin n\pi xdx\\
&=\frac{2[1-(-1)^n]}{n\pi}
\end{align*}
For $n=1,2,3,\cdots$, $A_{2n}=0$ and $A_{2n-1}=\frac{4}{2(n-1)\pi}$

Step 3. \begin{align*}
T_1(t)&=\left[\frac{4}{\pi}+\int_0^t be^{\alpha^2\pi^2s}ds\right]e^{-\alpha^2\pi^2t}\\
&=\frac{b}{\alpha^2\pi^2}+\left(\frac{4}{\pi}-\frac{b}{a^2\pi^2}\right)e^{-\alpha^2\pi^2t}
\end{align*}
For $n=2,3,\cdots,$
$$T_{2n-1}(t)=\frac{4}{(2n-1)\pi}e^{-\alpha^2(2n-1)^2\pi^2t}$$

Step 4. The solution is therefore
\begin{align*}
u(x,t)=\left[\frac{b}{\alpha^2\pi^2}+\left(\frac{4}{\pi}-\frac{b}{a^2\pi^2}\right)e^{-\alpha^2\pi^2t}\right]\sin\pi x\\
+\sum_{n=2}^\infty \frac{4}{(2n-1)\pi}e^{-\alpha^2(2n-1)^2\pi^2t}\sin(2n-1)\pi x
\end{align*}

The eventual temperature is
$$\lim_{t\to\infty}u(x,t)=\frac{b}{\alpha^2\pi^2}\sin\pi x$$

References:

David Betounes, Partial Differential Equations for Computational Science with Maple and Vector Analysis, TELOS, Springer-Verlag

# Solving Heat Equation with Non-Homogeneous BCs 2: Time-Dependent BCs

Consider the heat IBVP:
\begin{align*}
&u_t=\alpha^2 u_{xx},\ 0<x<L,\ t>0\\
&{\rm BCs:}\ \left\{\begin{aligned}
&u(0,t)=g_1(t),\ t>0\\
&u_x(L,t)+hu(L,t)=g_2(t),
\end{aligned}
\right.\\
&{\rm IC:}\ u(x,0)=\phi(x),\ 0<x<L.
\end{align*}
Again we seek a solution of the form:
$$u(x,t)=S(x,t)+U(x,t)$$
where $S(x,t)$ satisfies the non-homogeneous BCs:
\begin{align*}
S(0,t)&=g_1(t)\ \ \ \ \ (1)\\
S_x(L,t)+hS(L,t)&=g_2(t)
\end{align*}
Suppose that
$$S(x,t)=A(t)x+B(t)$$
From the BCs (1) we obtain
$$S(x,t)=\frac{g_2(t)-hg_1(t)}{1+hL}x+g_1(t)$$
$U(x,t)$ satisfies the semi-homogeneous heat IBVP:
\begin{align*}
&U_t=\alpha^2 U_{xx}-S_t,\ 0<x<L,\ t>0\\
&{\rm BCs:}\ \left\{\begin{aligned}
&U(0,t)=0,\ t>0\\
&U_x(L,t)+hU(L,t)=0,
\end{aligned}
\right.\\
&{\rm IC:}\ U(x,0)=\phi(x)-S(x,0),\ 0<x<L.
\end{align*}
Since $u(x,t)\approx S(x,t)$ for large time, the eventual temperature is $\lim_{t\to\infty}S(x,t)$. The homogeneous BCs are pretty much similar to the ones we studied here but the heat equation is non-homogeneous. Such a problem (non-homogeneous heat equation + homogeneous BCs) is called a semi-homogeneous IBVP. Unfortunately, we cannot solve the semi-homogenous heat IBVP yet. We will discuss how to solve semi-homogeneous IBVPs later.

Example. Solve the heat IBVP:
\begin{align*}
&u_t=\alpha^2 u_{xx}-ae^{-t}x,\ 0<x<1,\ t>0\\
&{\rm BCs:}\ \left\{\begin{aligned}
&u(0,t)=1,\ t>0\\
&u(1,t)=ae^{-t},
\end{aligned}
\right.\\
&{\rm IC:}\ u(x,0)=0,\ 0<x<1.
\end{align*}

Solution: Let $u(x,t)=S(x,t)+U(x,t)$
where $S(x,t)=A(t)x+B(t)$ satisfies the BCs. Then
$S(x,t)=(ae^{-t}-1)x+1$. $U(x,t)$ satisfies the homogeneous heat IBVP:
\begin{align*}
&U_t=\alpha^2 U_{xx},\ 0<x<1,\ t>0\\
&{\rm BCs:}\ \left\{\begin{aligned}
&U(0,t)=0,\ t>0\\
&U(1,t)=0,
\end{aligned}
\right.\\
&{\rm IC:}\ U(x,0)=-(a-1)x-1,\ 0<x<1.
\end{align*}
We know how to solve this problem and the solution $U(x,t)$ is given by
$$U(x,t)=\sum_{n=1}^\infty A_n e^{-\alpha^2n^2\pi^2t}\sin n\pi x$$
where for each $n=1,2,\cdots$
\begin{align*}
A_n&=2\int_0^1[-(a-1)x-1]\sin n\pi xdx\\
&=2\frac{a(-1)^n-1}{n\pi}
\end{align*}
That is,
$$U(x,t)=2\sum_{n=1}^\infty\frac{a(-1)^n-1}{n\pi}e^{-\alpha^2n^2\pi^2t}\sin n\pi x$$
Therefore, the solution to the original heat problem is
$$u(x,t)=(ae^{-t}-1)x+1+2\sum_{n=1}^\infty\frac{[a(-1)^n-1]}{n\pi}e^{-\alpha^2 n^2\pi^2 t}\sin n\pi x.$$
The eventual temperature is
$$\lim_{t\to\infty}S(x,t)=-x+1$$

References:

David Betounes, Partial Differential Equations for Computational Science with Maple and Vector Analysis, TELOS, Springer-Verlag

# Solving Heat Equation with Non-Homogeneous BCs 1: Time-Indepdendent BCs

Consider the following 1-dimensional heat IBVP:
\begin{align*}
&u_t=\alpha^2 u_{xx},\ 0<x<L,\ t>0\\
&{\rm BCs:}\ \left\{\begin{aligned}
u(0,t)&=c_1,\ t>0\\
u(L,t)&=c_2,
\end{aligned}
\right.\\
&{\rm IC:}\ u(x,0)=\phi(x),\ 0<x<L.
\end{align*}
Since the BCs are not homogeneous, the method of separation of variables cannot be used. Let us assume that $u(x,t)$ can be decomposed to
$$u(x,t)=S(x,t)+U(x,t)$$
where $S(x,t)$ satisfies the non-homogeneous BCs i.e.
\begin{align*}
S(0,t)&=c_1\ \ \ \ \ (1)\\
S(L,t)&=c_2
\end{align*}
This implies that $U(x,t)$ satisfies the homogeneous BCs
$$U(0,t)=U(L,t)=0$$
We don’t know the function $S(x,t)$ but we can make a guess. We can first try a possibly simplest one, a linear function in $x$
$$S(x,t)=A(t)x+B(t)$$
If this works, it is our luck day. If not, try something else that is next simplest one, say a quadratic function in $x$. It turns out that our guess works! This means that we can determine the coefficients $A(t)$ and $B(t)$ from the BCs (1) and we find that
$$A(t)=\frac{c_2-c_1}{L},\ B(t)=c_1$$
Thus
$$S(x,t)=\frac{c_2-c_1}{L}x+c_1$$
$U(x,t)$ satisfies the homogeneous IBVP:
\begin{align*}
&U_t=\alpha^2 U_{xx},\ 0<x<L,\ t>0\\
&{\rm BCs:}\ \left\{\begin{aligned}
U(0,t)&=0,\ t>0\\
U(L,t)&=0,
\end{aligned}
\right.\\
&{\rm IC:}\ U(x,0)=\phi(x)-S(x,0)=\phi(x)-\frac{c_2-c_1}{L}x-c_1,\ 0<x<L.
\end{align*}
We already know how to solve the homogeneous IBVP (see here) and
$$U(x,t)=\sum_{n=1}^\infty A_ne^{-\alpha^2\frac{n^2\pi^2}{L^2}t}\sin\frac{n\pi}{L}x$$
where
$$A_n=\frac{2}{L}\int_0^L[\phi(x)-\frac{c_2-c_1}{L}x-c_1]\sin\frac{n\pi}{L}xdx,\ n=1,2,\cdots$$

Note that $\lim_{t\to\infty}U(x,t)=0$ and so $\lim_{t\to\infty}u(x,t)=S(x)$. $U(x,t)$ is called the transient and $S(x)$ is called the steady state. The steady state is eventual state for large time.

Example. Solve the IBVP:
\begin{align*}
&u_t=\alpha^2 u_{xx},\ 0<x<1,\ t>0\\
&{\rm BCs:}\ \left\{\begin{aligned}
u(0,t)&=c_1,\ t>0\\
u(1,t)&=c_2,
\end{aligned}
\right.\\
&{\rm IC:}\ u(x,0)=0,\ 0<x<1.
\end{align*}

Solution: First we solve
\begin{align*}
&U_t=\alpha^2 U_{xx},\ 0<x<1,\ t>0\\
&{\rm BCs:}\ \left\{\begin{aligned}
U(0,t)&=0,\ t>0\\
U(1,t)&=0,
\end{aligned}
\right.\\
&{\rm IC:}\ U(x,0)=-c_1-(c_2-c_1)x,\ 0<x<1.
\end{align*}
The transient $U(x,t)$ is given by
$$U(x,t)=\sum_{n=1}^\infty A_ne^{-\alpha^2n^2\pi^2 t}\sin n\pi x$$
where
\begin{align*}
A_n&=-2\int_0^1[c_1+(c_2-c_1)]\sin n\pi xdx\\
&=\frac{2}{n\pi}[c_2(-1)^n-c_1].
\end{align*}
Hence the solution is
$$u(x,t)=(c_2-c_1)x+c_1+2\sum_{n=1}^\infty\frac{[c_2(-1)^n-c_1]}{n\pi}e^{-\alpha^2n^2\pi^2 t}\sin n\pi x$$

References:

David Betounes, Partial Differential Equations for Computational Science with Maple and Vector Analysis, TELOS, Springer-Verlag

# 1-Dimensional Heat Initial Boundary Value Problems 3: An Example of Heat IBVP with Mixed Boundary Conditions (Insulated and Specified Flux)

Let us consider the heat IBVP:
\begin{align*}
u_t=\alpha^2u_{xx},\ 0<x<1,\ t>0\\
u_x(0,t)=0,\ t>0\\
u_x(1,t)+u(1,t)=0,\ t>0\\
u(x,0)=1,\ 0<x<1.
\end{align*}
First we solve the Sturm-Liouville problem:
\begin{align*}
X^{\prime\prime}=kX\ \mbox{with BCs}\\
\left\{
\begin{aligned}
&X'(0)=0\\
&X'(1)+X(1)=0.
\end{aligned}
\right.
\end{align*}
The two cases $k=0$ and $k=\lambda^2>0$ are left to be checked for students. Note that the case $k=0$ leads to a trivial solution and for $k=\lambda^2>0$ case there is no solution. Let $k=-\lambda^2$. Then
$$X(x)=A\cos\lambda x+B\sin \lambda x$$
and
$$X'(x)=-A\lambda\sin\lambda x+B\lambda\cos\lambda x.$$
We obtain $B=0$ from the condition $X'(0)=0$. So $X'(x)=-A\lambda\sin\lambda x$. The BC $X'(1)+X(1)=0$ then implies that the eigenvalues $\lambda$ satisfies
$$\lambda=\cot\lambda\ \ \ \ \ (1)$$ The roots of (1) can be found by the Newton’s method. In Maxima, the Newton’s method can be performed as follows.

First we define $\cot(x)-x$ as a function $f(x)$:

(%i1) f(x):=cot(x)-x;

Newton’s method applied to f(x) is performed by

(%i2) newton(t):=subst(t,x,x-f(x)/diff(f(x),x));

To use Newton’s method we need the initial approximation. We can make a guess from its graph. As you can see the first zero is near 0.5, so we choose $x_0=0.5$.

If you just run

(%i3) newton(0.5);

Maxima will return the output

(%o3)                              0.74865744241706

In order to calculate next approximation $x_1$, simply run

(%i4) newton(%);

(%o4)                              0.85239848545338

For the next approximation $x_2$, again run

(%i5) newton(%);

(%o5)                              0.86029880244504

By running the same command, the next approximation $x_3$ is obtained as $x_3=0.86033358835819$. This value is accurate to eight decimal places. If you think this is good enough, you can stop.

Similarly we obtain $\lambda_n$ for $n=2,3,4,\cdots$:
$$\begin{array}{|c||c|c|} \hline n & \lambda_n & \frac{(2n-1)\pi}{2}\\ \hline 1 & .8603358 & 1.570796327\\ \hline 2 & 3.425618 & 4.712388981\\ \hline 3 & 6.437298 & 7.853981635\\ \hline 4 & 9.529334 & 10.99557429\\ \hline 5 & 12.645287 & 14.13716694\\ \hline \cdots & \cdots & \cdots\\ \hline \end{array}$$
The table shows that as $n$ increases $\frac{(2n-1)\pi}{2}$ gets closer to $\lambda_n$. So for a practical purpose one may use the approximation
$$\lambda_n\approx \frac{(2n-1)\pi}{2},\ \lambda_n=1,2,3,\cdots$$
The temperature distribution $u(x,t)$ is then given by
$$u(x,t)=\sum_{n=1}^\infty A_ne^{-\alpha^2\lambda_n^2t}\cos\lambda_n x$$
The $L_n$’s and $A_n$’s are computed to be:
\begin{align*}
L_n&=\int_0^1\cos^2\lambda_nxdx\\
&=\frac{1}{2}\int_0^1[\cos(2\lambda_nx)+1]dx\\
&=\frac{1}{2}\left[\frac{\sin 2\lambda_n}{2\lambda_n}+1\right]\\
&=\frac{1}{2}\left[\frac{\sin\lambda_n\cos\lambda_n}{\lambda_n}+1\right]\ (\sin2\lambda_n=2\sin\lambda_n\cos\lambda_n)\\
&=\frac{1}{2}\left[\frac{\sin\lambda_n\cos\lambda_n+\lambda_n}{\lambda_n}\right]
\end{align*}
and
\begin{align*}
A_n&=\frac{1}{L_n}\int_0^1\phi(x)\cos\lambda_nxdx\\
&=\frac{2\lambda_n}{\sin\lambda_n\cos\lambda_n+\lambda_n}\int_0^1\cos\lambda_nxdx\\
&=\frac{2\sin\lambda_n}{\sin\lambda_n\cos\lambda_n+\lambda_n}.
\end{align*}

This calculation is done by assuming that $\cos\lambda_n x$, $n=1,2,\cdots$ are orthogonal functions. In fact they are orthogonal functions if $\lambda_n=\frac{(2n-1)\pi}{2}$, $n=1,2,\cdots$.

Finally we have
$$u(x,t)=\sum_{n=1}^\infty\frac{2\sin\lambda_n}{\sin\lambda_n\cos\lambda_n+\lambda_n}e^{-\alpha^2\lambda_n^2 t}\cos\lambda_n x.$$
For an approximation one may take
$$u(x,t)\approx\sum_{n=1}^\infty\frac{4(-1)^{n-1}}{(2n-1)\pi}e^{-\alpha^2\frac{(2n-1)^2\pi^2}{4}t}\cos\frac{(2n-1)\pi x}{2}\ \ \ \ \ (2)$$

The 3D plot of $u(x,t)$ in (2) is

References:

David Betounes, Partial Differential Equations for Computational Science with Maple and Vector Analysis, TELOS, Springer-Verlag