Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to computing page for the fourth course APMA0360
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to Mathematica tutorial for the fourth course APMA0360
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340
Return to the main page for the course APMA0360
Introduction to Linear Algebra with Mathematica

Glossary

1D Heat Equation: Neumann BC

We consider the initial boundary value problem for the heat equation subject to the Neumann boundary conditions on half line:

\begin{align} \label{eqHeat.1} u_t &= \alpha u_{xx} - \gamma\, u + f(x,t) , & x > 0, \quad 0 < t < \infty , \\ \label{eqHeat.2} \frac{\partial u}{\partial x} (+0,t) &= -g(t) , & 0< t < \infty , \\ \label{eqHeat.3} u(x,+0) &= u_0 (x) , & x> 0 , \\ u &\in 𝔏^2 (\mathbb{R}_{+} \times \mathbb{R}_{+}) , \qquad \mbox{regularity conditions}, \notag \\ \int_0^{\infty} g(t)\,{\text d}t &= 0 , \qquad \mbox{compatibility condition}. \notag \end{align}
where subscripts indicate partial derivatives. For example, \( \displaystyle \quad u_t = \frac{\partial u}{\partial t} . \quad \) The heat equation contains (phisical) real parameters, so α > 0 and γ ≥ 0. All functions in the heat equation \eqref{eqHeat.1}, boundary condition \eqref{eqHeat.2}, and the initial condition \eqref{eqHeat.3} are assumed to be in the Hilbert space ℌ². Also all derivatives (ut, uxx) are assumed to be continuous in the open domain 0 < x < ∞, 0 < t < ∞

Example 1: Let us consider the Neumann problem on the half-line with homogeneous boundary conditions:
\begin{align*} u_t &= \alpha u_{xx} , & x > 0, \quad 0 < t < \infty , \\ \frac{\partial u}{\partial x}(+0,t) &= 0 , & 0< t < \infty , \tag{1.1} \\ u(x,+0) &= u_0(x) , & x> 0 , \\ u &\in 𝔏^2 (\mathbb{R}_{+} \times \mathbb{R}_{+}) \qquad (\mbox{regularity conditions}) . \end{align*}
We reduce the given IBVP to an IVP on the whole line by extending the initial data u₀(x) to the negative half-axis in such a fashion that the boundary condition is automatically satisfied. Notice that if ψ(s) is an even function, i.e., ψ(−x) = ψ(,i>x), then its derivative function will be odd. Indeed, differentiating in the definition of the even function, we get −ψ′(−x) = ψ′(x), which is the same as ψ′(−x) = −ψ′(x). Hence, for an arbitrary even function ψ(x), ψ′(0) = 0. It is now clear that extending the initial data so that the resulting function is even will produce solutions to the IVP on the whole line that automatically satisfy the Neumann condition of (9).

We define the even extension of u₀(x), \[ \mathfrak{\phi_{\mathrm{even}}}(x) = \begin{cases} u_0 (x), &\quad\mbox{for }\ x> 0 , \\ u_0 (-x), &\quad\mbox{for }\ x \le 0 , \end{cases} \] and consider the following IVP on the whole line \[ \begin{split} v_t &= \alpha\,v_{xx} , \qquad -\infty < x < \infty , \quad 0 < t < \infty , \\ v_x (x, 0) = \mathfrak{\phi_{\mathrm{even}}}(x) . \end{split} \tag{1.2} \] It is clear that the solution v(x, t) of the IVP (1.1) will be even in x, since the difference [v(−x, tv(x, t)] solves the heat equation and has zero initial data. We then use the solution formula for the IVP on the whole line to write \[ v(x,t) = \int_{-\infty}^{\infty} S(x-\xi , t)\,\mathfrak{\phi_{\mathrm{even}}}(\xi )\,{\text d}\xi . \] Then function \[ w(x,t) = v(x,t) \big\vert_{x\ge 0} . \] solves the IBVP (1.1), and use the expression for the heat kernel, as well as the definition (1.2), to write the solution formula as follows: \[ w(x,t) = \frac{1}{\sqrt{4\pi\alpha t}} \int_0^{\infty} \left[ e^{-(x-\xi )^2 /(4\alpha t)} + e^{-(x+\xi )^2 /(4\alpha t)} \right] u_0 (\xi )\,{\text d}\xi . \] In terms of heat conduction, the Neumann condition in (1.1) means that there is no heat exchange between the rod and the environment (recall that the heat flux is proportional to the spatial derivative of the temperature). The physical interpretation of formula for w is that the integrand is the contribution of u₀(ξ) plus an additional contribution, which comes from the lack of heat transfer to the points of the rod with negative coordinates.

As a numerical example, we take u₀(x) = χ[0, R], the characteristic function of interval [0, R]. The the solution of corresponding IBVP (1.1) becomes \begin{align*} u(x,t) &= \frac{1}{\sqrt{4\pi\alpha t}} \int_0^R \left[ e^{-(x-\xi )^2 /(4\alpha t)} + e^{-(x+\xi )^2 /(4\alpha t)} \right] {\text d}\xi \\ &= \frac{1}{\sqrt{\pi}} \int_{-x/\sqrt{4\alpha t}}^{(R-x)/\sqrt{4\alpha t}} e^{-p^2} {\text d}p + \frac{1}{\sqrt{\pi}} \int_{x/\sqrt{4\alpha t}}^{(R+x)/\sqrt{4\alpha t}} e^{-q^2} {\text d}q , \end{align*} where we made the changes of variables \[ p = \frac{\xi -x}{\sqrt{4\alpha t}} , \qquad q = \frac{\xi +x}{\sqrt{4\alpha t}} . \] Using the error function, we simplify \[ w(x,t) = \mbox{erf}\left( \frac{R-x}{\sqrt{4\alpha t}} \right) + \mbox{erf}\left( \frac{R+x}{\sqrt{4\alpha t}} \right) + C , \] where C is an arbitrary constant.

   ■
End of Example 1

To solve the given initial boundary value problem, we use the cosine Fourier transformation. So we multiply the differential equation \eqref{eqHeat.1} and the initial condition \eqref{eqHeat.3} by cosine (kx) and integrate with respect to x from 0 to ∞. This results in the equation

\[ \frac{{\text d} u^c}{{\text d}t} = \alpha \int_0^{\infty} \frac{\partial^2 u}{\partial x^2} \,\cos (kx) \,{\text d}x - \gamma \,u^c (k,t) + f^c (k,t) , \]
and the initial condition
\begin{equation} \label{eqHeat.4} u^{c} (k, +0) = u^c_0 (k) = \int_0^{\infty} u_0 (x)\,\cos (kx)\, {\text d} x \end{equation}
for the cosine Fourier transform of the unknown function
\[ u^{c} (k, t) = ℱ^s_{x \to k} \left[ u(x,t) \right] = \int_0^{\infty} u(x,t)\,\cos (kx)\,{\text d} x . \]
The cosine Fourier transforms of the given function in Eq. \eqref{eqHeat.1} is
\[ f^{c} (k, t) = ℱ^s_{x \to k} \left[ f(x,t) \right] = \int_0^{\infty} f(x,t)\,\cos (kx)\,{\text d} x . \]
Since the differential equation contains second derivative under integration, we integrate by parts twice and obtain
\begin{align*} \int_0^{\infty} \frac{\partial^2 u}{\partial x^2} \,\cos (kx)\,{\text d}x &= \left. \frac{\partial u}{\partial x}\,\cos (kx) \right\vert_{x=0}^{x\to +\infty} - \int_0^{\infty} \frac{\partial u}{\partial x} \,k\, \sin (kx)\,{\text d}x \\ &= -u(x,t)\,k\,\sin (kx) \Large\vert_{x=0}^{x\to +\infty} - \int_0^{\infty} u(x,t)\, k^2 \cos (kx)\,{\text d}x \\ &= k\, u(+0,t) - k^2 u^s (k,t) . \end{align*}
This leads to the the ordinary differential equation for the sine Fourier transform of the unknown function:
\begin{equation} \label{eqHeat.5} \frac{{\text d} u^s}{{\text d}t} = - \left( \alpha\,k^2 + \gamma \right) u^s (k,t) + k\alpha\, g (t) + f^s (k, t) . \end{equation}
We ask Mathematica to solve the initial value problem \eqref{eqHeat.5}, \eqref{eqHeat.4}.
DSolve[{u'[t] + (a*k^2 + gamma)*u[t] == F[t], u[0] == u0}, u[t], t]
{{u[t] -> E^(-gamma t - a k^2 t) (u0 + Inactive[Integrate][ E^(gamma K[1] + a k^2 K[1]) F[K[1]], {K[1], 0, t}])}}
\[ u^s (k, t) = e^{- (\alpha k^2 + \gamma )t} \left[ u_0^s (k) + \int_0^t e^{ (\alpha k^2 + \gamma ) \tau} \left( k\alpha\, g (\tau ) + f^s (k, \tau ) \right) {\text d} \tau \right] . \]
Application of the inverse sine Fourier transform to the latter gives the required solution
\[ u(x,t) = \frac{2}{\pi} \int_0^{\infty} u^s (k, t)\,\sin (kx)\,{\text d} k = I_{u0} + I_g + I_f , \]
where
\begin{align*} I_{u0} &= \frac{2}{\pi} \int_0^{\infty} \sin (kx)\,e^{- (\alpha k^2 + \gamma )t} \,{\text d} k \,\int_0^{\infty} {\text d}\xi u_0 (\xi )\,\sin (k\xi ) , \\ &= \frac{1}{2\sqrt{\pi\alpha t}}\, e^{-\gamma t} \,\int_0^{+\infty} {\text d}\xi u_0 (\xi ) \left[ e^{- (x-\xi )^2 /(4\alpha t)} - e^{- (x+\xi )^2 /(4\alpha t)} \right] \\ I_{g} &= \frac{2}{\pi} \int_0^{\infty} \sin (kx)\,{\text d} k \,k\alpha \, \int_0^t e^{-\gamma (t- \tau )} \,e^{- \alpha k^2 (t-\tau )} \,g(\tau ) \,{\text d} \tau \\ &= \frac{x}{2\sqrt{\pi\alpha} } \,\int_0^{t} e^{-\gamma (t- \tau )} \, e^{-x^2 /(4\alpha (t-\tau ))} \left( t- \tau \right)^{-3/2} g(\tau )\,{\text d}\tau , \\ I_{f} &= \frac{2}{\pi} \int_0^{+\infty} \sin (kx)\,e^{- (\alpha k^2 + \gamma )(t- \tau )} \,{\text d} k \,\int_0^t f^s (k, \tau )\,{\text d}\tau \\ &= \frac{2}{\pi} \int_0^{\infty} \sin (kx)\,e^{- (\alpha k^2 + \gamma )(t- \tau )} \,{\text d} k \,\int_0^t {\text d}\tau \int_0^{\infty} f(\xi , \tau )\,\sin (k\xi )\,{\text d}\xi \\ &= \frac{1}{2\sqrt{\pi\alpha}} \int_0^{\infty} f(\xi, \tau ) \,{\text d}\xi \int_0^t {\text d}\tau \left[ e^{(x-\xi )^2 /(4\alpha (t-\tau ))} - e^{(x+\xi )^2 /(4\alpha (t-\tau ))} \right] e^{-\gamma (t-\tau )} . \end{align*}
Integrate[Sin[k*x]*Sin[k*xi]*Exp[-a*k^2], {k, 0, Infinity}]
ConditionalExpression[((E^(-((x - xi)^2/(4 a))) - E^(-((x + xi)^2/(4 a)))) Sqrt[\[Pi]])/(4 Sqrt[a]), Re[a] > 0]
Integrate[Exp[-a*k^2]*k*Sin[k*x], {k, 0, Infinity}]
ConditionalExpression[(E^(-(x^2/(4 a))) Sqrt[\[Pi]] x)/(4 a^(3/2)), Re[a] > 0]
Integrate[Exp[-a*k^2]*Sin[k*x]*Sin[k*xi], {k, 0, Infinity}]
ConditionalExpression[((E^(-((x - xi)^2/(4 a))) - E^(-((x + xi)^2/(4 a)))) Sqrt[\[Pi]])/(4 Sqrt[a]), Re[a] > 0]

 

  1. Boyd, J.P. and Flyer, N., Compatibility conditions for time-dependent partial differential equations and the rate of convergence of Chebyshev and Fourier spectral methods, Computer Methods in Applied Mechanics and Engineering, 1999, Volume 175, Issues 3–4, Pages 281--309. https://doi.org/10.1016/S0045-7825(98)00358-2
  2. Chatziafratis, A., Boundary behaviour of the solution of the heat equation on the half line via the Fokas unified transform method, 2024, https://doi.org/10.48550/arXiv.2401.08331
  3. Chatziafratis, A., Fokas, A., Aifantis, E.C., Variations of heat equation on the half-line via the Fokas method, 2024, First published: 08 September 2024 https://doi.org/10.1002/mma.10303
  4. Chatziafratis, A., Mantzavinos, D., Boundary behavior for the heat equation on the half-line, 2022, https://doi.org/10.1002/mma.8245

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Basic Concepts
Return to the Part 2 Fourier Series
Return to the Part 3 Integral Transformations
Return to the Part 4 Parabolic PDEs
Return to the Part 5 Hyperbolic PDEs
Return to the Part 6 Elliptic PDEs
Return to the Part 6P Potential Theory
Return to the Part 7 Numerical Methods