Three-dimensional heat conduction boundary value problem

Let us consider the three-dimensional nonsymmetric heat conduction problem in a ball of radius $R$. The heat conduction equation (1.8$\displaystyle \frac{\partial T(M, t)}{\partial t}$ $\displaystyle = a^2 \cdot \Delta T(M, t)$ $\displaystyle + f(M, t),$) in spherical coordinates takes the form

where $r$ $\in (0,R)$, $\theta$ $\in (0,\pi)$, $\phi$ $\in (0,2\pi)$, $R$ is the radius, $T(r, \theta, \phi, t)$ is the temperature, $a^2$ is the thermal diffusivity, $f(r, \theta, \phi, t)$ is the heat source density function.

The Dirichlet boundary conditions (1.13$\displaystyle T(M, t)$ $\displaystyle = \Phi(M, t), \quad M$ $\displaystyle \in S.$) in the three-dimensional case have the form

where $T_0(r, \theta, \phi)$ is the initial condition, $\Phi(\theta, \phi, t)$ is the boundary condition, which is, of course, periodic in $\phi$; at the poles ($\theta$ $= 0$ and $\theta$ $= \pi$) the boundary conditions are replaced by the requirement that the solution be bounded.

To get rid of the inhomogeneity, we introduce $T(r, \theta, \phi, t)$ $= \widehat{T}(r, \theta, \phi, t)$ $+ U(\theta, \phi, t)$, let us choose the function $U(\theta, \phi, t)$ $= \Phi(\theta, \phi, t)$, in accordance with (2.2$\displaystyle \widehat{T}(M, t)$ $\displaystyle + U(M, t)$ $\displaystyle = \Phi(M, t), \quad \widehat{T}(M, t)$ $\displaystyle = 0 \;$ $\displaystyle \Rightarrow\; U(M, t)$ $\displaystyle = \Phi(M, t), \quad M$ $\displaystyle \in S.$). Then equation (2.71$\displaystyle \frac{\partial T(r, \theta, \phi, t)}{\partial t}$ $\displaystyle = a^2 \cdot \frac{\partial^2 T(r, \theta, \phi, t)}{\partial r^2}$ $\displaystyle + \frac{2 \cdot a^2}{r} \cdot \frac{\partial T(r, \theta, \phi, t)}{\partial r}$ $\displaystyle + \frac{a^2}{r^2 \cdot \sin(\theta)} \cdot \frac{\partial}{\partial \theta} \left( \sin(\theta) \cdot \frac{\partial T(r, \theta, \phi, t)}{\partial \theta} \right)$ $\displaystyle + \frac{a^2}{r^2 \cdot \sin^2 \theta} \cdot \frac{\partial^2 T(r, \theta, \phi, t)}{\partial \phi^2}$ $\displaystyle + f(r, \theta, \phi, t),$) and the boundary conditions (2.72$\displaystyle T(r, \theta, \phi, 0)$ $\displaystyle = T_0(r, \theta, \phi),$$\displaystyle T(R, \theta, \phi, t)$ $\displaystyle = \Phi(\theta, \phi, t),$$\displaystyle |T(r, 0, \phi, t)|$ $\displaystyle < \infty, \quad |T(r, \pi, \phi, t)|$ $\displaystyle < \infty,$$\displaystyle T(r, \theta, \phi, t)$ $\displaystyle = T(r, \theta, \phi + 2\pi, t),$) take the form

where

The Neumann boundary conditions (1.14$\displaystyle \lambda \cdot \frac{\partial T(M, t)}{\partial \vec{n}}$ $\displaystyle = \Phi(M, t), \quad M$ $\displaystyle \in S,$) in the three-dimensional case have the form

where $\frac{\displaystyle \partial T(r, \theta, \phi, t)}{\displaystyle \partial \mathbf{n}} \bigg|_{r=R}$ is the normal derivative of the temperature on the sphere of radius $R$; the direction of the normal coincides with the direction of the coordinate $r$.

To get rid of the inhomogeneity, we introduce $T(r, \theta, \phi, t)$ $= \widehat{T}(r, \theta, \phi, t)$ $+ U(r, \theta, \phi, t)$, let us choose the function $U(r, \theta, \phi, t)$ $= \Phi(\theta, \phi, t) \cdot r$ — it satisfies condition (2.3$\displaystyle \frac{\partial \widehat{T}(M, t)}{\partial \vec{n}}$ $\displaystyle + \frac{\partial U(M, t)}{\partial \vec{n}}$ $\displaystyle = \Phi(M, t), \quad \frac{\partial \widehat{T}(M, t)}{\partial \vec{n}}$ $\displaystyle = 0 \;$ $\displaystyle \Rightarrow\; \frac{\partial U(M, t)}{\partial \vec{n}}$ $\displaystyle = \Phi(M, t), \quad M$ $\displaystyle \in S.$): $\frac{\displaystyle \partial U(r, \theta, \phi, t)}{\displaystyle \partial \mathbf{n}} \bigg|_{r=R}$ $= \Phi(\theta, \phi, t)$. Then equation (2.71$\displaystyle \frac{\partial T(r, \theta, \phi, t)}{\partial t}$ $\displaystyle = a^2 \cdot \frac{\partial^2 T(r, \theta, \phi, t)}{\partial r^2}$ $\displaystyle + \frac{2 \cdot a^2}{r} \cdot \frac{\partial T(r, \theta, \phi, t)}{\partial r}$ $\displaystyle + \frac{a^2}{r^2 \cdot \sin(\theta)} \cdot \frac{\partial}{\partial \theta} \left( \sin(\theta) \cdot \frac{\partial T(r, \theta, \phi, t)}{\partial \theta} \right)$ $\displaystyle + \frac{a^2}{r^2 \cdot \sin^2 \theta} \cdot \frac{\partial^2 T(r, \theta, \phi, t)}{\partial \phi^2}$ $\displaystyle + f(r, \theta, \phi, t),$) and the boundary conditions (2.74$\displaystyle T(r, \theta, \phi, 0)$ $\displaystyle = T_0(r, \theta, \phi),$$\displaystyle \frac{\displaystyle \partial T(r, \theta, \phi, t)}{\displaystyle \partial \mathbf{n}} \bigg|_{r=R}$ $\displaystyle = \Phi(\theta, \phi, t),$$\displaystyle |T(r, 0, \phi, t)|$ $\displaystyle < \infty, \quad |T(r, \pi, \phi, t)|$ $\displaystyle < \infty,$$\displaystyle T(r, \theta, \phi, t)$ $\displaystyle = T(r, \theta, \phi + 2\pi, t),$) take the form

where

Equation (2.8$\displaystyle \Delta \Psi(M)$ $\displaystyle + \gamma^2 \cdot \Psi(M)$ $\displaystyle = 0.$) takes the form

Let us represent the function $\Psi(r, \theta, \phi)$ as a product of three functions $\Psi(r, \theta, \phi)$ $= R(r) \cdot \Theta(\theta) \cdot \Upsilon(\phi)$, substitute it into equation (2.76$\displaystyle \frac{\displaystyle \partial^2 \Psi(r, \theta, \phi)}{\displaystyle \partial r^2}$ $\displaystyle + \frac{\displaystyle 2}{\displaystyle r} \cdot \frac{\displaystyle \partial \Psi(r, \theta, \phi)}{\displaystyle \partial r}$ $\displaystyle + \frac{\displaystyle 1}{\displaystyle r^2 \cdot \sin(\theta)} \cdot \frac{\displaystyle \partial}{\displaystyle \partial \theta} \left( \sin(\theta) \cdot \frac{\displaystyle \partial \Psi(r, \theta, \phi)}{\displaystyle \partial \theta} \right)$ $\displaystyle + \frac{\displaystyle 1}{\displaystyle r^2 \cdot \sin^2 \theta} \cdot \frac{\displaystyle \partial^2 \Psi(r, \theta, \phi)}{\displaystyle \partial \phi^2}$ $\displaystyle + \gamma^2 \cdot \Psi(r, \theta, \phi)$ $\displaystyle = 0.$) and carry out the transformations

The expression in the first bracket depends only on $r$, the one in the second bracket only on the angles; therefore, dividing by $R(r) \cdot \Theta(\theta) \cdot \Upsilon(\phi)$, both sides can be set equal to a constant $\gamma_1^2$ — by analogy with the separation of variables for time and geometry. The equation for the radius takes the form

and continue the transformations

Here the expression in brackets depends only on $\theta$, while the last term depends only on $\phi$; therefore, multiplying by $\frac{\displaystyle \sin^2(\theta)}{\displaystyle \Theta(\theta) \cdot \Upsilon(\phi)}$, we separate the variables again with the constant $m^2$ and obtain the equations for the polar and azimuthal angles

The solution of equation (2.79$\displaystyle \frac{\displaystyle d^2 \Upsilon(\phi)}{\displaystyle d\phi^2}$ $\displaystyle + m^2 \cdot \Upsilon(\phi)$ $\displaystyle = 0.$) we already found in the two-dimensional case, so we refer to (2.54$\displaystyle \Upsilon_{1m}(\theta)$ $\displaystyle = C_{1m} \cdot \sin(m \cdot \theta), m$ $\displaystyle \in (1..\infty),$$\displaystyle \Upsilon_{2m}(\theta)$ $\displaystyle = C_{2m} \cdot \cos(m \cdot \theta), m$ $\displaystyle \in (0..\infty).$). And to solve equation (2.78$\displaystyle \frac{\displaystyle 1}{\displaystyle \sin(\theta)} \cdot \frac{\displaystyle d}{\displaystyle d\theta} \left( \sin(\theta) \cdot \frac{\displaystyle d \Theta(\theta)}{\displaystyle d\theta} \right)$ $\displaystyle + \left[ \gamma_1^2 - \frac{\displaystyle m^2}{\displaystyle \sin^2(\theta)} \right] \cdot \Theta(\theta)$ $\displaystyle = 0,$) let us introduce the substitution $z$ $= \cos(\theta)$: for $\theta$ $\in (0, \pi)$ we have $z$ $\in (-1, 1)$ and $dz$ $= - \sin(\theta) \cdot d\theta$. Then the equation takes the form

This equation is analyzed in the appendix “The Legendre equation”; its solution has the form (E.6$\displaystyle \Theta_{km}(z)$ $\displaystyle = P_k^{(m)}(z)$ $\displaystyle = (1 - z^2)^{\frac{\displaystyle m}{\displaystyle 2}} \cdot \frac{\displaystyle d^m}{\displaystyle dz^m} P_k(z),$)

where $\gamma_{1k}^2$ $= k \cdot (k + 1)$ are the eigenvalues, $P_k^{(m)}(z)$ are the associated Legendre polynomials, and nontrivial solutions exist only for $m$ $\leq k$. The norm of the associated Legendre polynomial has the form (F.6$\displaystyle \int_{-1}^1 \left[ P_k^{(m)}(z) \right]^2 \,dz$ $\displaystyle = \frac{\displaystyle 2}{\displaystyle 2 \cdot k + 1} \cdot \frac{\displaystyle (k+m)!}{\displaystyle (k-m)!}, \quad m$ $\displaystyle \leq k.$). Now let us return to the angle $\theta$, making the inverse substitution

Generally speaking, this solution is defined up to a constant factor — let us set it equal to one: it is absorbed into the expansion coefficients anyway. Now let us return to equation (2.77$\displaystyle r^2 \cdot \frac{\displaystyle d^2 R(r)}{\displaystyle dr^2}$ $\displaystyle + 2 \cdot r \cdot \frac{\displaystyle d R(r)}{\displaystyle dr}$ $\displaystyle + \left( r^2 \cdot \gamma^2 - \gamma_1^2 \right) \cdot R(r)$ $\displaystyle = 0,$) and rewrite it taking into account $\gamma_1^2$ $= k \cdot (k + 1)$

This equation is analyzed in the appendix “The Bessel equation in spherical coordinates”: the substitution $R(r)$ $= \frac{\displaystyle \widehat{R}(r)}{\displaystyle \sqrt{r}}$ reduces it to the Bessel equation with half-integer index, and the solution has the form (G.2$\displaystyle R(r)$ $\displaystyle = \frac{\displaystyle 1}{\displaystyle \sqrt r} \cdot J_{k + 1/2}(\gamma \cdot r),$). It should be noted that the solution is bounded at the point $r$ $= 0$: the Bessel function $J_{k + 1/2}(\gamma \cdot r)$ tends to zero no slower than $\sqrt{r}$.

Thus, the solutions of equation (2.76$\displaystyle \frac{\displaystyle \partial^2 \Psi(r, \theta, \phi)}{\displaystyle \partial r^2}$ $\displaystyle + \frac{\displaystyle 2}{\displaystyle r} \cdot \frac{\displaystyle \partial \Psi(r, \theta, \phi)}{\displaystyle \partial r}$ $\displaystyle + \frac{\displaystyle 1}{\displaystyle r^2 \cdot \sin(\theta)} \cdot \frac{\displaystyle \partial}{\displaystyle \partial \theta} \left( \sin(\theta) \cdot \frac{\displaystyle \partial \Psi(r, \theta, \phi)}{\displaystyle \partial \theta} \right)$ $\displaystyle + \frac{\displaystyle 1}{\displaystyle r^2 \cdot \sin^2 \theta} \cdot \frac{\displaystyle \partial^2 \Psi(r, \theta, \phi)}{\displaystyle \partial \phi^2}$ $\displaystyle + \gamma^2 \cdot \Psi(r, \theta, \phi)$ $\displaystyle = 0.$) have the form

Let us substitute the solutions (2.82$\displaystyle \begin{aligned} &\Psi_{1km}(r, \theta, \phi) = A_{km} \cdot \frac{\displaystyle 1}{\displaystyle \sqrt{r}} \cdot J_{k + 1/2}(\gamma \cdot r) \cdot P_k^{(m)}(\cos(\theta)) \cdot \sin(m \cdot \phi),\\ &1 \leq m \leq k, \quad k \geq 1, \end{aligned}$$\displaystyle \begin{aligned} &\Psi_{2km}(r, \theta, \phi) = A_{km} \cdot \frac{\displaystyle 1}{\displaystyle \sqrt{r}} \cdot J_{k + 1/2}(\gamma \cdot r) \cdot P_k^{(m)}(\cos(\theta)) \cdot \cos(m \cdot \phi),\\ &0 \leq m \leq k, \quad k \geq 0, \end{aligned}$) into the boundary conditions (2.73$\displaystyle \frac{\displaystyle \partial \widehat{T}(r, \theta, \phi, t)}{\displaystyle \partial t}$ $\displaystyle = a^2 \cdot \frac{\displaystyle \partial^2 \widehat{T}(r, \theta, \phi, t)}{\displaystyle \partial r^2}$ $\displaystyle + \frac{2 \cdot a^2}{r} \cdot \frac{\displaystyle \partial \widehat{T}(r, \theta, \phi, t)}{\displaystyle \partial r}$ $\displaystyle + \frac{a^2}{r^2 \cdot \sin(\theta)} \cdot \frac{\displaystyle \partial}{\displaystyle \partial \theta} \left( \sin(\theta) \cdot \frac{\displaystyle \partial \widehat{T}(r, \theta, \phi, t)}{\displaystyle \partial \theta} \right)$ $\displaystyle + \frac{a^2}{r^2 \cdot \sin^2 \theta} \cdot \frac{\displaystyle \partial^2 \widehat{T}(r, \theta, \phi, t)}{\displaystyle \partial \phi^2}$ $\displaystyle + \widehat{f}(r, \theta, \phi, t),$$\displaystyle \widehat{T}(r, \theta, \phi, 0)$ $\displaystyle = T_0(r, \theta, \phi)$ $\displaystyle - U(\theta, \phi, 0),$$\displaystyle \widehat{T}(R, \theta, \phi, t)$ $\displaystyle = 0,$$\displaystyle |\widehat{T}(r, 0, \phi, t)|$ $\displaystyle < \infty, \quad |\widehat{T}(r, \pi, \phi, t)|$ $\displaystyle < \infty,$$\displaystyle \widehat{T}(r, \theta, \phi, t)$ $\displaystyle = \widehat{T}(r, \theta, \phi + 2\pi, t),$) and obtain

Clearly, the solution makes sense only when $A_{km}$ $\neq 0$, which is possible only when $J_{k + 1/2}(\gamma \cdot R)$ $= 0$. Let us denote the roots of the equation $J_{k + 1/2}(\mu)$ $= 0$ by $\mu_{ks}$ — then the eigenvalues take the form

Zero and negative roots are no good: for $\mu$ $= 0$ the eigenfunction $J_{k + 1/2}(0 \cdot r)$ $\equiv 0$ — identically zero, which is not an eigenfunction (the index $k$ $+ 1/2$ $> 0$), and negative roots give no new solutions. Therefore we take only the positive roots: $\mu_{ks}, s$ $\in (1..\infty)$.

Thus, the eigenvalues and eigenfunctions have the form

where $A_{kms}$ is an arbitrary constant factor: the eigenfunction is defined up to it. Let us set $A_{kms}$ $= 1$ — this factor is absorbed into the expansion coefficients anyway, as in the general solution.

To expand functions into a Fourier series in $\Psi_{1kms}(r, \theta, \phi)$ and $\Psi_{2kms}(r, \theta, \phi)$ we need to compute the norms $\|\Psi_{1kms}(r, \theta, \phi)\|^2$ and $\|\Psi_{2kms}(r, \theta, \phi)\|^2$; the weight for spherical coordinates is $\rho(r, \theta, \phi)$ $= r^2 \cdot \sin(\theta)$.