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

1D Heat Equation: Dirichlet BC


We consider the initial boundary value problem for the heat equation subject to the Dirichlet 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} u(+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 \\ u(0,0) &= u_0 (+0) = g(+0) , \qquad \mbox{corner 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 heat equation on the half-line (0, ∞) with homogeneous boundary conditions:
\begin{align*} u_t &= \alpha u_{xx} , & x > 0, \quad 0 < t < \infty , \\ u(+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}), \\ u(0,0) &= u_0 (+0) = 0 \qquad (\mbox{corner condition}). \end{align*}
If the solution to the above mixed initial/boundary value problem exists, then we know that it must be unique from an application of the maximum principle. In terms of the heat conduction, one can think of u in (1.1) as the temperature in an infinite rod, one end of which is kept at a constant zero temperature. The initial temperature of the rod is then given by u₀(x). Our goal is to solve the IBVP (1.1), and derive a solution formula, much like what we did for the heat IVP on the whole line. But instead of constructing the solution from scratch, it makes sense to try to reduce this problem to the IVP on the whole line, for which we already have a solution formula. This is achieved by extending the initial data u₀(x) to the whole line. We have a choice of how exactly to extend the data to the negative half-line, and one should try to do this in such a fashion that the boundary condition of (1.1) is automatically satisfied by the solution to the IVP on the whole line that arises from the extended data. This is the case, if one chooses the odd extension of φ(x), which we describe next. By definition, a function ψ(x) is odd, if ψ(−x) = −ψ(x). But then plugging in x = 0 into this definition, one gets ψ(0) = 0 for any odd function. Recall also that the solution u(x, t) to the heat IVP with odd initial data is itself odd in the x variable. This follows from the fact that the sum [u(x, t) + u(−x, t)] solves the heat equation and has zero initial data, hence, it is the identically zero function by the uniqueness of solutions. Then, by our above observation for odd functions, we would have that u(0, t) = 0 for any t > 0, which is exactly the boundary condition of (1.1)

This shows that if one extends u₀(x) to an odd function on the whole line, then the solution with the extended initial data automatically satisfies the boundary condition of (1.1). Let us then define \[ \mathfrak{\phi_{\mathrm{odd}}}(x) = \begin{cases} u_0 (x), &\quad\mbox{for }\ x> 0 , \\ - u_0 (-x), &\quad\mbox{for }\ x< 0 , \\ 0, &\quad\mbox{for }\ x = 0. \end{cases} \] It is clear that φodd is an odd function, since we defined it for negative x by reflecting the u₀(x) with respect to the vertical axis, and then with respect to the horizontal axis. This procedure produces a function whose graph is symmetric with respect to the origin, and thus it is odd. One can also verify this directly from the definition of odd functions. Now, let v(x, t) be the solution of the following IVP on the whole line \[ \begin{cases} v_t &= \alpha\,v_{xx} , \qquad -\infty < x < \infty , \quad 0 < t < \infty , \\ v(x,0) &= \mathfrak{\phi_{\mathrm{odd}}}(x) , \qquad -\infty < x < \infty . \end{cases} \tag{1.2} \] From previous section we know that the solution to (1.2) is given by the formul \[ v(x,t) = \int_{-\infty}^{+\infty} S(x-\xi , t)\,\mathfrak{\phi_{\mathrm{odd}}}(\xi )\,{\text d}\xi , \quad t > 0, \] where the kernel is \[ S(x,t) = \frac{1}{\sqrt{4\pi\alpha t}}\, e^{- x^2/(4\alpha t)} \] Restricting the x variable to only the positive half-line produces the function \[ u(x,t) = v(x,t) \big\vert_{x\ge 0} . \] We claim that this u(x, t) is the unique solution of IBVP (1.1). Indeed, u(x, t) solves the heat equation on the positive half-line, since so does v(x, t). Furthermore, \begin{align*} v(x,t) &= \int_0^{\infty} S(x-\xi , t) \,\mathfrak{\phi_{\mathrm{odd}}}(\xi )\,{\text d}\xi + \int_{-\infty}^0 S(x-\xi , t) \,\mathfrak{\phi_{\mathrm{odd}}}(\xi )\,{\text d}\xi \\ &= \int_0^{\infty} S(x-\xi , t) \,u_0 (\xi )\,{\text d}\xi - \int_{-\infty}^0 S(x-\xi , t) \,u_0 (-\xi )\,{\text d}\xi . \end{align*} Making the change of variables ξ ↦ −ξ in the second integral on the right hand-side, and flipping the integration limits gives \[ v(x,t) = \int_0^{\infty} S(x-\xi , t) \,u_0 (\xi )\,{\text d}\xi - \int_{-\infty}^0 S(x+\xi , t) \,u_0 (\xi )\,{\text d}\xi \] Using definition of ϕodd and the above expression for v(x, t), as well as the expression of the heat kernel S(x, t), we can write the solution formula for the IBVP (1) as \[ u(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 . \tag{1.3} \] The method used to arrive at this solution formula is called the method of odd extensions or the reflection method. We can make a physical sense of formula (1.3) by interpreting the integrand as the contribution from the point ξ minus the heat loss from this point due to the constant zero temperature at the endpoint.

   ■

End of Example 1

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

\[ \frac{{\text d} u^s}{{\text d}t} = \alpha \int_0^{\infty} \frac{\partial^2 u}{\partial x^2} \,\sin (kx) \,{\text d}x - \gamma \,u^s (k,t) + f^s (k,t) , \]
and the initial condition
\begin{equation} \label{eqHeat.4} u^{s} (k, +0) = u^s_0 (k) = \int_0^{\infty} u_0 (x)\,\sin (kx)\, {\text d} x \end{equation}
for the sine Fourier transform of the unknown function
\[ u^{s} (k, t) = ℱ^s_{x \to k} \left[ u(x,t) \right] = \int_0^{\infty} u(x,t)\,\sin (kx)\,{\text d} x . \]
The sine Fourier transforms of the given function in Eq. \eqref{eqHeat.1} is
\[ f^{s} (k, t) = ℱ^s_{x \to k} \left[ f(x,t) \right] = \int_0^{\infty} f(x,t)\,\sin (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} \,\sin (kx)\,{\text d}x &= \left. \frac{\partial u}{\partial x}\,\sin (kx) \right\vert_{x=0}^{x\to +\infty} - \int_0^{\infty} \frac{\partial u}{\partial x} \,k\, \cos (kx)\,{\text d}x \\ &= -u(x,t)\,k\,\cos (kx) \Large\vert_{x=0}^{x\to +\infty} - \int_0^{\infty} u(x,t)\, k^2 \sin (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]

We integrate by parts each integral:

\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 )\,\frac{\text d}{{\text d}\xi} \left( -\frac{\cos (k\xi )}{k} \right) \\ &= \frac{2}{\pi} \int_0^{\infty} \sin (kx)\,e^{- (\alpha k^2 + \gamma )t} \,{\text d} k \left[ -\left. u_0 (\xi )\, \frac{\cos (k\xi )}{k} \right\vert_{\xi =0}^{+\infty} + \int_0^{+\infty} \frac{\cos (k\xi )}{k} \, \frac{{\text d} u_0 }{{\text d}\xi}\,{\text d}\xi \right] \\ &= \frac{2}{\pi}\, u_0 (+0) \,\int_0^{\infty} \frac{\sin (kx)}{k}\,e^{- (\alpha k^2 + \gamma )t} \,{\text d} k + \frac{2}{\pi} \int_0^{\infty} \sin (kx)\,e^{- (\alpha k^2 + \gamma )t} \,{\text d} k \,\int_0^{+\infty} \frac{\cos (k\xi )}{k} \, \frac{{\text d} u_0 }{{\text d}\xi}\,{\text d}\xi \\ &= u_0 (+0)\,e^{-\gamma t}\,\mbox{erf} \left( \frac{x}{2\sqrt{\alpha t}} \right) + \frac{1}{2} \, \int_0^{+\infty} {\text d}\xi \left[ \mbox{erf} \left( \frac{x - \xi}{2\sqrt{\alpha t}} \right) + \mbox{erf} \left( \frac{x + \xi}{2\sqrt{\alpha t}} \right) \right] \frac{{\text d} u_0 }{{\text d}\xi} , \end{align*}
Integrate[Sin[k*x]*Exp[-a*k^2]/k, {k, 0, Infinity}]
ConditionalExpression[1/2 \[Pi] Erf[x/(2 Sqrt[a])], Re[a] > 0]
Integrate[Sin[k*x]*Cos[k*xi]*Exp[-a*k^2]/k, {k, 0, Infinity}]
ConditionalExpression[ 1/4 \[Pi] (Erf[(x - xi)/(2 Sqrt[a])] + Erf[(x + xi)/(2 Sqrt[a])]), Re[a] > 0]
where
\[ \mbox{erf} (t) = \frac{2}{\sqrt{\pi}}\,\int_0^t e^{-\eta^2} \,{\text d}\eta \]
is the error function. Since x²/t has no limit as x,t approach zero, the integral Iu0(x, t) has not limit value at the origin unless u0(0) = 0.

Next we make a substitution in integral Ig

\[ \tau = t - \frac{x^2}{4\alpha T} \qquad \Longrightarrow \qquad {\text d}\tau =\frac{x^2}{4\alpha T^2} \,{\text d}T . \]
This yields
\begin{align*} I_{g} (x, t) &= \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 \\ &= \frac{x}{2\sqrt{\pi\alpha}} \,\int_{x^2/(4\alpha t)}^{+\infty} e^{-\gamma x^2 /(4\alpha T)} \, e^{-T} \left( \frac{4 \alpha T}{x^2} \right)^{3/2} \,\frac{x^2}{4\alpha T^2} \,g \left( t - \frac{x^2}{4\alpha T} \right) {\text d}T \\ &= \frac{1}{\sqrt{\pi}} \, \int_{x^2/(4\alpha t)}^{+\infty} e^{-\gamma x^2 /(4\alpha T)} \, e^{-T} \, g \left( t - \frac{x^2}{4\alpha T} \right) \frac{{\text d}T}{\sqrt{T}} . \end{align*}
Upon adding and subtructing g(t), we transfer the integral to the following form:
\[ I_g (x,t) = \frac{g(t)}{\sqrt{\pi}} \, \int_{x^2/(4\alpha t)}^{+\infty} e^{-\gamma x^2 /(4\alpha T)} \, e^{-T} \,\frac{{\text d}T}{\sqrt{T}} + \frac{1}{\sqrt{\pi}} \, \int_{x^2/(4\alpha t)}^{+\infty} e^{-\gamma x^2 /(4\alpha T)} \, e^{-T} \left[ g \left( t - \frac{x^2}{4\alpha T} \right) - g(t) \right] \frac{{\text d}T}{\sqrt{T}} . \]
Using Mathematica, we simplify tis expression
Integrate[ Exp[-gamma*x^2 /(4*a*t)]*Exp[-T]/Sqrt[T], {T, x^2 /(4*a*t), Infinity}]/Sqrt[Pi]
E^(-((gamma x^2)/(4 a t))) Erfc[1/2 Sqrt[x^2/(a t)]]
\[ I_g (x,t) = g(t)\, e^{- \gamma \,x^2 /(4\alpha t)} \,\mbox{erf}\left( \frac{x}{2\sqrt{\alpha t}} \right) + \frac{1}{\sqrt{\pi}} \, \int_{x^2/(4\alpha t)}^{+\infty} e^{-\gamma x^2 /(4\alpha T)} \, e^{-T} \left[ g \left( t - \frac{x^2}{4\alpha T} \right) - g(t) \right] \frac{{\text d}T}{\sqrt{T}} . \]
In integrals Iu0 and If we make substitution span class="math">\( \displaystyle \quad \xi \pm x = X\.\sqrt{4\alpha t} ; \quad \) this yields
\begin{align*} I_{u0} (x,t) &= \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] \\ &= e^{-\gamma t} \,\int_{-x/\sqrt{4\alpha t}}^{+\infty} \, e^{-X} \, u_0 \left( x + X\,\sqrt{4\alpha t} \right) {\text d}X - e^{-\gamma t} \,\int_{x/\sqrt{4\alpha t}}^{+\infty} \, e^{-X} \, u_0 \left( -x + X\,\sqrt{4\alpha t} \right) {\text d}X ; \end{align*}
and
\begin{align*} I_{f} (x,t) &= \frac{1}{2\sqrt{\pi\alpha}} \int_0^{\infty} f(\xi, \tau ) \,{\text d}\xi \int_{-x/\sqrt{4\alpha t}}^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 )} . \\ &= \int_0^t {\text d}\tau \,e^{-\gamma (t-\tau )} \left[ \int_{-x/\sqrt{4\alpha t}}^{\infty} {\text d}X \, f \left( x + X\,\sqrt{4\alpha t} \right) e^{-X} - \int_{x/\sqrt{4\alpha t}}^{\infty} {\text d}X \, f \left( -x + X\,\sqrt{4\alpha t} \right) e^{-X} e^{-X} \right] . \end{align*}

Numerical Experiment


We calculate upon evaluation integrals with pecewise input constant functions. Let us start with u₀.
\[ u_0 (x) = \begin{cases} 1 , & \quad\mbox{if} \quad 0 \le x \le 1 , \\ 0, & \quad \mbox{otherwise} . \end{cases} \]
Integrate[ Exp[-(x - xi)^2 /4/t] - Exp[-(x + xi)^2 /4/t], {xi, 0, 1}]/(2* Sqrt[Pi*t])
1/2 (-Erf[(-1 + x)/(2 Sqrt[t])] + 2 Erf[x/(2 Sqrt[t])] - Erf[(1 + x)/(2 Sqrt[t])])
\[ I_0 (x,t) = \mbox{erf} \left( \frac{x}{2\sqrt{t}} \right) - \frac{1}{2}\,\mbox{erf} \left( \frac{x-1}{2\sqrt{t}} \right) - \frac{1}{2}\,\mbox{erf} \left( \frac{1+x}{2 \sqrt{t}} \right) , \]
wjere erf is the error function. We plot the function
Plot3D[1/ 2 (-Erf[(-1 + x)/(2 Sqrt[t])] + 2 Erf[x/(2 Sqrt[t])] - Erf[(1 + x)/(2 Sqrt[t])]), {x, 0, 5}, {t, 0, 1}, PlotStyle -> Orange, Mesh -> None, PlotRange -> {0, 1}]
Function u₀.

As it is seen from the graph, the corresponding integral I0(x, t) is unbounded at x = 0. To check the problem, we consider another initial function that vanishes at the origin:

\[ u_0 (x) = \begin{cases} 1 , & \quad\mbox{if} \quad 0.5 \le x \le 1.5 , \\ 0, & \quad \mbox{otherwise} . \end{cases} \]
Integrate[ Exp[-(x - xi)^2 /4/t] - Exp[-(x + xi)^2 /4/t], {xi, 0, 1}]/(2* Sqrt[Pi*t])
1/2 (-Erf[(-3 + 2 x)/(4 Sqrt[t])] + Erf[(-1 + 2 x)/(4 Sqrt[t])] + Erf[(1 + 2 x)/(4 Sqrt[t])] - Erf[(3 + 2 x)/(4 Sqrt[t])])
\[ I_0 (x,t) = \frac{1}{2}\,\mbox{erf} \left( \frac{2x-1}{4\sqrt{t}} \right) - \frac{1}{2}\,\mbox{erf} \left( \frac{2x-3}{4\sqrt{t}} \right) + \frac{1}{2}\,\mbox{erf} \left( \frac{2x+1}{4\sqrt{t}} \right) - \frac{1}{2}\,\mbox{erf} \left( \frac{2x+3}{4 \sqrt{t}} \right) , \]
wjere erf is the error function. We plot the function
Plot3D[1/ 2 (-Erf[(-3 + 2 x)/(4 Sqrt[t])] + Erf[(-1 + 2 x)/(4 Sqrt[t])] + Erf[(1 + 2 x)/(4 Sqrt[t])] - Erf[(3 + 2 x)/(4 Sqrt[t])]), {x, 0, 5}, {t, 0, 1}, PlotStyle -> Orange, Mesh -> None, PlotRange -> {0, 1}]
Function u₀ vanishes at the origin

Now we consider function ug, where g is again fiecewise constant.

(x/2/Sqrt[Pi])* Integrate[ Exp[-(t - tau)]*Exp[-x^2/4/(t - tau)]/(t - tau)^(3/2), {tau, 0, Max[1, t]}]
1/(2 Sqrt[x^2])E^-Sqrt[ x^2] x (-Erfc[-((-2 t + Sqrt[x^2])/(2 Sqrt[t]))] + Erfc[(2 t - Sqrt[x^2] - 2 Max[1, t])/(2 Sqrt[t - Max[1, t]])] + E^(2 Sqrt[ x^2]) (Erfc[(2 t + Sqrt[x^2])/(2 Sqrt[t])] - Erfc[(2 t + Sqrt[x^2] - 2 Max[1, t])/(2 Sqrt[t - Max[1, t]])])
Since this expression is cumbersome, we use numerical integration. Note that the integrand has a weak singularity as 1/√X as X → 0.
Plot3D[(x/2/Sqrt[Pi])* NIntegrate[ Exp[-(t - tau)]*Exp[-x^2/4/(t - tau)]/(t - tau)^(3/2), {tau, 0, Min[1, t]}], {x, 0, 5}, {t, 0, 3}, PlotStyle -> Purple, Mesh -> None, PlotRange -> {0, 2}]
Function ug.
Finally, we plot function uf for
\[ f (x,t) = \begin{cases} f(t) , & \quad \mbox{if} \quad 0 \le x \le 1 , \\ 0 , & \quad \mbox{otherwise}. \end{cases} \]
We substitute the intermedian integral (with α = 1 for simplicity)
(1/2/Sqrt[Pi])* Integrate[ Exp[(x - xi)^2 /4/(t - tau)] - Exp[(x + xi)^2 /4/(t - tau)] , {xi, 0, 1}]
-1/2 Sqrt[ t - tau] (Erfi[(-1 + x)/(2 Sqrt[t - tau])] - 2 Erfi[x/(2 Sqrt[t - tau])] + Erfi[(1 + x)/(2 Sqrt[t - tau])])
into the general formula,
\[ I_f (x, t) = \int_0^t f(\tau )\,{\text d}\tau \left\{ \frac{1}{\sqrt{t - \tau}} \left[ \mbox{erf} \left( \frac{x}{2\sqrt{t - \tau}} \right) - \frac{1}{2} \,\mbox{erf} \left( \frac{x-1}{2\sqrt{t - \tau}} \right) - \frac{1}{2} \,\mbox{erf} \left( \frac{x+1}{2\sqrt{t - \tau}} \right) \right] \right\} . \tag{f.1} \]
For numerical integration, we intrduce ε = 0.1 amd integrate with respect to τ from 0 to t − ε because the integrand approach 0 exponentially near τ = t; so
\[ I_f (x, t) = \int_0^{t-\epsilon} f(\tau )\,{\text d}\tau \left\{ \frac{1}{\sqrt{t - \tau}} \left[ \mbox{erf} \left( \frac{x}{2\sqrt{t - \tau}} \right) - \frac{1}{2} \,\mbox{erf} \left( \frac{x-1}{2\sqrt{t - \tau}} \right) - \frac{1}{2} \,\mbox{erf} \left( \frac{x+1}{2\sqrt{t - \tau}} \right) \right] \right\} . \tag{f.2} \]
Upon choosing f(t) = χ[0, 2] (t), we obtain
Plot3D[ NIntegrate[ -(1/2) Sqrt[ t - tau] (Erfi[(-1 + x)/(2 Sqrt[t - tau])] - 2 Erfi[x/(2 Sqrt[t - tau])] + Erfi[(1 + x)/(2 Sqrt[t - tau])]), {tau, 0, Min[2, t - 0.1]}], {x, 0, 5}, {t, 0.1, 3}, PlotStyle -> Brown, Mesh -> None, PlotRange -> {0, 1}]
Function uf.

Boundary operator


 

  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