Preface

The Resolvent method for linear partial differential equations was developed in mid 1980s by Vladimir Dobrushkin.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VI of the course APMA0340
Introduction to Linear Algebra with Mathematica

Glossary

The Resolvent method for heat equation

Previously, we discussed the separation of variables method that can be applied to certain classes of linear partial differential equations (PDEs for short). However, this method has serious limitations because not every linear PDE is acceptable for separation of variables. Also it requires the given domain, PDE, and boundary conditions to be separable as well. Furthermore, it expresses the solution as either an integral or a series, neither of which are uniformly convergent on the boundary of the domain (for appropriate boundary conditions), which renders such expressions unsuitable for numerical computations.

A further extension of separation of variables is the integral transform method. However, it suffers the same drawback as the separation of variables---integrals or series (through which solutions are expressed) do not converge uniformly at the boundary. The integral transform method is successful only when its inverse transform is known, which heavily depends on the spectrum of the corresponding operator. Therefore, many practical boundary value problems are not suitable for separation of variables and integral transform methods.

The resolvent method reduces the boundary value problem to the integral equation of the second kind with compact operator. Such integral equations always have a unique solution. Then the solution of the boundary value problem is expressed explicitly through the solution of the integral equation. Till some extend, the resolvent method is an extension of the boundary element method.

We illustrate the resolvent method by considering the heat equation on the half line with the boundary conditions of the third kind. Similar problem was treated by the Fokas method (or unified transform) in the previous section. So our main concern is the following initial boundary value problem for positive half-line:

\begin{align*} &\mbox{Heat equation:} \qquad &u_t &= u_{xx} \qquad\mbox{for } x > 0 \mbox{ and } t > 0, \\ &\mbox{Initial condition:} \qquad &u(x,0) &= f(x) , \qquad x > 0, \\ &\mbox{Boundary condition:} \qquad &a\,u(0,t) - b\,u_x (0,t) &= g(t) , \end{align*}
with nonzero coefficients 𝑎 and b. It is convenient to introduce the projection operators:
\begin{align*} \left[ B_1 u \right] (t) &= u(0,t) = \lim_{x\to 0} \, u(x,t) , \\ \left[ B_2 u \right] (t) &= u_x (0,t) = \lim_{x\to 0} \, \frac{\partial u}{\partial x} . \end{align*}
Then the above boundary condition at x = 0 can be written in compact form:
\[ a\,B_1 u - b\,B_2 u = g . \]

We separate the given problem into two auxiliary boundary problems subject to the Dirichlet and Neumann conditions. Since the initial conditions play no role in the boundary conditions, we can assume they homogeneous without any loss of generality. Therefore, we consider the following two problems.

Auxiliary Problem 1:

\[ \begin{split} u_t &= \alpha\, u_{xx} , \qquad x \in (0 ,\infty ), \quad t > 0, \\ u(x,0) &= 0 , \\ u(0,t) &= g(t) \qquad \mbox{or} \qquad B_1 u = g; \end{split} \]
and

Auxiliary Problem 2:

\[ \begin{split} u_t &= \alpha\, u_{xx} , \qquad x \in (0 ,\infty ), \quad t > 0, \\ u(x,0) &= 0 , \\ u_x (0,t) &= g(t) \qquad \mbox{or} \qquad B_2 u = g. \end{split} \]
Solutions of these problems were found in section of this tutorial. They are
\begin{align*} S_1 g (x,t) &= \frac{x}{2\sqrt{\pi\alpha}} \, \int_0^t g(\tau )\left( t- \tau \right)^{-3/2} e^{-x^2/4\alpha (t-\tau )}\,{\text d}\tau = \frac{2}{\sqrt{\pi}} \,\int_{v(x,t)}^{\infty} g \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) e^{-\eta^2} \,{\text d}\eta , \\ S_2 g (x,t) &= -\frac{\sqrt{\alpha}}{\sqrt{\pi\,t}} \, \int_0^t g(\tau )\left( t- \tau \right)^{-1/2} e^{-x^2/4\alpha (t-\tau )}\,{\text d}\tau = -\frac{x}{\sqrt{\pi t}} \, \int_{v(x,t)}^{\infty} g \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) e^{-\eta^2} \,\eta^{-2} \,{\text d}\eta , \end{align*}
where \( \displaystyle v(x,t) = \frac{x}{2\sqrt{\alpha\,t}} . \)

With the projection operators

\[ \left[ B_1 u(x,t) \right] (t) = \lim_{x\to 0} \, u(x,t) \qquad \mbox{and} \qquad \left[ B_2 u(x,t) \right] (t) = \lim_{x\to 0} \, \frac{\partial u}{\partial x} , \]
we have obviously the following relations:
\[ B_1 S_1 = I \quad\mbox{or} \quad \left[ B_1 S_1 \right] g(t) = g(t) \qquad\mbox{and} \qquad B_2 S_2 = I \quad\mbox{or} \quad \left[ B_2 S_2 \right] g(t) = g(t) , \]
where I is the identity operator.

According to the resolvent method, we seek a solution to the given initial mixed boundary value problem in the form

\[ u(x,t) = \left( S_1 + S_2 \right) \varphi , \]
where unknown function φ(t) will be chosen shortly. Recall that the heat equation and the initial conditions are homogeneous, therefore, the above function u(x,t) will also satisfy them. Hence, we need to choose φ so that the mixed boundary conditions will be valid. In other words, we need to satisfy the boundary condition
\[ \left( a\,B_1 - b\,B_2 \right) \left( S_1 + S_2 \right) \varphi (t) = g(t) \]
Using established previously identities \( B_1 S_1 = I \quad\mbox{and} \quad B_2 S_2 = I , \) we conclude that the required function φ should satisfy the equation
\[ \left( a\,B_1 S_2 - b\,B_2 S_1 \right) \varphi (t) + \left( a-b \right) \varphi (t) = g(t) . \]
Therefore, we need to determine two limits:
\begin{align*} \left( B_1 S_2 \right) \varphi (t) &= -\lim_{x\to 0} \, \frac{\sqrt{\alpha}}{\sqrt{\pi\,t}} \, \int_0^t \varphi (\tau )\left( t- \tau \right)^{-1/2} e^{-x^2/4\alpha (t-\tau )}\,{\text d}\tau = -\sqrt{\frac{\alpha}{\pi\, t}} \, \int_0^t \varphi (\tau )\left( t- \tau \right)^{-1/2} \,{\text d}\tau , \\ \left( B_2 S_1 \right) \varphi (t) &= \lim_{x\to 0} \, \frac{\partial}{\partial x} \frac{x}{2\sqrt{\pi\alpha}} \, \int_0^t \varphi (\tau )\left( t- \tau \right)^{-3/2} e^{-x^2/4\alpha (t-\tau )}\,{\text d}\tau = \lim_{x\to 0} \, \frac{\partial}{\partial x}\,\frac{2}{\sqrt{\pi}} \int_{v(x,t)}^{\infty} \varphi \left( t - \frac{x^2}{4\alpha \eta^2} \right) e^{-\eta^2} \,{\text d}\eta . \end{align*}
Assuming[{t > 0 && aa > 0 && x > 0}, 1/2/Sqrt[Pi*aa]* Integrate[tau^(-3/2)*Exp[-x^2/4/aa/(tau)], {tau, 0, t}]]
Erfc[x/(2 Sqrt[aa t])]/x
Assuming[{t > 0 && aa > 0 && x > 0}, aa/Sqrt[Pi*aa]* Integrate[tau^(-1/2)*Exp[-x^2/4/aa/(tau)], {tau, 0, t}]]
(Sqrt[aa] (2 E^(-(x^2/(4 aa t))) Sqrt[t] - ( Sqrt[\[Pi]] x Erfc[x/(2 Sqrt[aa t])])/Sqrt[aa]))/Sqrt[\[Pi]]
The first operator B1S2 is an integral operator with weak singularity, so it converges uniformly with respect to x and we can interchange the order of integration and limit. Indeed, the first integral can be broken in two parts
\[ \sqrt{\frac{\alpha}{\pi}}\int_0^t g(t-\tau )\,\tau^{-1/2} e^{-x^2 /4\alpha \tau} \,{\text d}\tau = \sqrt{\frac{\alpha}{\pi}}\int_0^{\varepsilon} g(t-\tau )\, \tau^{-1/2} e^{-x^2 /4\alpha \tau} \,{\text d}\tau + \sqrt{\frac{\alpha}{\pi}}\int_{\varepsilon}^{\infty} g(t-\tau )\, \tau^{-1/2} e^{-x^2 /4\alpha \tau} \,{\text d}\tau , \]
where ε is a small number. In the interval [0,ε], the function g(t - τ) is approximately g(t), so it can be moved out of integral:
\[ \sqrt{\frac{\alpha}{\pi}}\int_0^{\varepsilon} g(t-\tau )\, \tau^{-1/2} e^{-x^2 /4\alpha \tau} \,{\text d}\tau \approx g(t)\, \sqrt{\frac{\alpha}{\pi}}\int_0^{\varepsilon} \tau^{-1/2} e^{-x^2 /4\alpha \tau} \,{\text d}\tau . \]
The latter is a small number because the integration is performed along a small interval. So we conclude that the main contribution comes from the integral over semi-infinite interval [ε,∞), which converges uniformly with respect to x and we can perform finding the limit as x → 0.
FindMaximum[1/Sqrt[t] Exp[-(1/t)], {t, 0.01`}]
{3.72008*10^-43, {t -> 0.01}}

We know that

\[ \lim_{x\to 0} \frac{x}{2\sqrt{\pi\alpha}} \, \tau^{-3/2} e^{-x^2/4\alpha \tau} = \delta (\tau ) ,\qquad\mbox{the delta function of Dirac}, \]
so this limit is not uniform with respect to x. First, we find the derivative assuming that g(0) = 0:
\[ \frac{\partial}{\partial x} \,S_1 g(t) = \frac{\partial}{\partial x} \, \frac{2}{\sqrt{\pi}} \,\int_{v(x,t)}^{\infty} g \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) e^{-\eta^2} \, {\text d}\eta = \frac{2}{\sqrt{\pi}} \,\int_{v(x,t)}^{\infty} g' \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) \left( -\frac{x}{2\alpha\,\eta^2} \right) e^{-\eta^2} \, {\text d}\eta , \]
where \( \displaystyle v(x,t) = \frac{x}{2\sqrt{\alpha\,t}} . \) Next, we integrate by parts using the identity
\[ \frac{\partial}{\partial \eta}\, g \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) = g' \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) \left( \frac{x^2}{2\alpha\,\eta^3} \right) ; \]
this yields
\begin{align*} \frac{\partial}{\partial x} \,S_1 g(t) &= \frac{2}{x\sqrt{\pi}} \, \int_{v(x,t)}^{\infty} \left[ \frac{\partial}{\partial\eta}\,g \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) \right] \frac{2\alpha \eta^3}{x^2} \left( -\frac{x}{2\alpha\,\eta^2} \right) e^{-\eta^2} \,{\text d}\eta \\ &= \frac{2}{x\sqrt{\pi}} \, \int_{v(x,t)}^{\infty} g \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) \left( 1 - 2\eta^2 \right) e^{-\eta^2} \,{\text d}\eta . \end{align*}
So
\[ \frac{\partial}{\partial x} \,S_1 g(t) = \frac{2}{x\sqrt{\pi}} \, \int_{v(x,t)}^{\infty} g \left( t - \frac{x^2}{4\alpha\,\eta^2} \right) \frac{\text d}{{\text d}\eta} \left[ \eta\,e^{-\eta^2} \right] \,{\text d}\eta . \]
Assuming that the function g(t) is analytic and g(0) = 0,
\[ g(t) = \sum_{k\ge 1} a_k t^k , \]
we find the limit
\begin{align*} \lim_{x\to 0} \, \frac{\partial}{\partial x} \,S_1 g(t) &= \lim_{t\to 0} \, \frac{2}{x\sqrt{\pi}} \,\int_{v(x,t)}^{\infty} \sum_{k\ge 1} a_k \left( t - \frac{x^2}{4\alpha\,\eta^2} \right)^k \left( \frac{\text d}{{\text d}\eta} \, \eta\, e^{-\eta^2} \right) {\text d}\eta \\ &= \frac{2}{\sqrt{\pi}} \,\sum_{k\ge 1} a_k \, \lim_{x\to 0} \, \frac{1}{x}\, \int_{v(x,t)}^{\infty} \left( t - \frac{x^2}{4\alpha\,\eta^2} \right)^k \left( \frac{\text d}{{\text d}\eta} \, \eta\, e^{-\eta^2} \right) {\text d}\eta . \end{align*}
Using the binomial theorem,
\[ \left( t -a \right)^k = t^k - t^{k-1} a + \sum_{i=2}^k \binom{k}{i} t^{k-i} \left( -a \right)^i , \]
we see that every term will be
\begin{align*} &\lim_{x\to 0} \, \frac{1}{x}\, \int_{v(x,t)}^{\infty} \left( t - \frac{x^2}{4\alpha\,\eta^2} \right)^k \left( \frac{\text d}{{\text d}\eta} \, \eta\, e^{-\eta^2} \right) {\text d}\eta \\ &= \lim_{x\to 0} \, \frac{1}{x}\,t^k \int_{v(x,t)}^{\infty} \left( \frac{\text d}{{\text d}\eta} \, \eta\, e^{-\eta^2} \right) {\text d}\eta - \lim_{x\to 0} \, x\, \int_{v(x,t)}^{\infty} \left( \frac{1}{4\alpha\,\eta^2}\right)\left( \frac{\text d}{{\text d}\eta} \, e^{-\eta^2} \right) + \cdots . \end{align*}
So each limit will be zero as having multiple x except the first one. Hence,
\begin{align*} \lim_{x\to 0} \, \frac{\partial}{\partial x} \,S_1 g(t) &= \frac{2}{\sqrt{\pi}} \,\sum_{k\ge 1} a_k \, \lim_{x\to 0} \, \frac{1}{x}\, \int_{v(x,t)}^{\infty} t^k \left( \frac{\text d}{{\text d}\eta} \, \eta\, e^{-\eta^2} \right) {\text d}\eta \\ &= \frac{2}{\sqrt{\pi}} \,g(t)\, \lim_{x\to 0} \, \frac{1}{x}\, \int_{v(x,t)}^{\infty} \left( \frac{\text d}{{\text d}\eta} \, \eta\, e^{-\eta^2} \right) {\text d}\eta \\ &= - \frac{2}{\sqrt{\pi}} \,g(t)\, \lim_{x\to 0} \, \frac{1}{x}\, v(x,t)\, e^{-v^2} = - \frac{1}{\sqrt{\pi\alpha\,t}}\, g(t) . \end{align*}
Finally, we obtain for φ(t) the integral equation of second kind:
\[ -a\, \sqrt{\frac{\alpha}{\pi\,t}} \, \int_0^t \varphi (\tau )\left( t- \tau \right)^{-1/2} \,{\text d}\tau + \frac{b}{\sqrt{\pi\alpha\,t}}\, \varphi (t) + \left( a-b \right) \varphi (t) = g(t) . \]
Note that the above Volterra integral equation has a weak singularity and therefore the above equation has a unique solution. Therefore, this integral equation (which has a unique solution) can be solved numerically.

When a = b, the above integral equation becomes

\[ -\sqrt{\frac{\alpha}{\pi}} \, \int_0^t \varphi (\tau )\left( t- \tau \right)^{-1/2} \,{\text d}\tau + \frac{1}{\sqrt{\pi\alpha}}\, \varphi (t) = g(t) \,\sqrt{t} . \]
It can be solved explicitly using the Laplace transform:
\[ -\sqrt{\frac{\alpha}{\pi}} \,\varphi^L (\lambda )\, \sqrt{\frac{\pi}{\lambda}} + \frac{1}{\sqrt{\pi\alpha}}\, \varphi^L (\lambda ) = {\cal L} \left[ g(t) \,\sqrt{t} \right] , \]
where
\[ \varphi^L (\lambda ) = {\cal L} \left[ \varphi \right] = \int_0^{\infty} \varphi (t)\, e^{-\lambda\,t} \,{\text d}t \]
is the Laplace transform of the function φ. Solving the above algebraic equation and applying the inverse Laplace transform, we arrive
\[ \varphi (t) = {\cal L}^{-1} \left[ \frac{1}{\frac{1}{\sqrt{\pi\alpha}} - \sqrt{\frac{\alpha}{\lambda}}} \, {\cal L} \left[ g(t) \,\sqrt{t} \right] \right] = \sqrt{\pi\alpha} \,{\cal L}^{-1} \left[ \frac{\sqrt{\lambda}}{\sqrt{\lambda} - \alpha\sqrt{\pi}}\, {\cal L} \left[ g(t) \,\sqrt{t} \right] \right] . \]

Example: We consider the initial boundary value problem

\[ \begin{split} u_t = \alpha\, u_{xx} , \qquad x,t \in (0,\infty ), \\ u(0,t) - u_x (0,t) = g(t) = \sin (\omega\, t) , \\ u(x,0) = f(x) = 0 . \end{split} \]
where ω is a positive constant.

According to the resolvent method, the solution of the above problems is given by

\[ u(x,t) = \left( S_1 + S_2 \right) \varphi , \]
where unknown function φ(t) is a solution of the integral equation
\[ -\sqrt{\frac{\alpha}{\pi}} \, \int_0^t \varphi (\tau )\left( t- \tau \right)^{-1/2} \,{\text d}\tau + \frac{1}{\sqrt{\pi\alpha}}\, \varphi (t) = \sqrt{t}\,\sin (\omega\,t) . \]
Application of the Laplace transformation to the above integral equation of the second kind gives
\[ \varphi (t) = \sqrt{\pi\alpha} \,{\cal L}^{-1} \left[ \frac{\sqrt{\lambda}}{\sqrt{\lambda} - \alpha\sqrt{\pi}}\, \frac{\sqrt{\pi} \,\sin \left[ \frac{3}{2}\,\arctan \left( \frac{\omega}{\lambda} \right) \right]}{2 \left( \omega^2 + \lambda^2 \right)^{3/4}} \right] . \]
because Mathematica knows how to find the Laplace transform:
LaplaceTransform[Sqrt[t]*Sin[omega*t], t, s]
(omega Sqrt[\[Pi]] Sin[3/2 ArcTan[Sqrt[omega^2]/s]])/(2 Sqrt[omega^2] (omega^2 + s^2)^( 3/4))
This formula has only theoretical value because finding the inverse Laplace transform is a changeable task because the integrand has also an issential singularity at λ = 0 besides the poles, and a branch point.    ■

 

  1. Cargniard, L., Reflection and Refraction of Progressive Seismic Waves, New York, McGraw-Hill Book Company, 1962.
  2. Dobrushkin, V.A., Boundary Value problems for elastodynamics in wedge shaped domains, Minsk, Science and Engineering, 1988.

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions