Preface

This section provides an introduction to Cagnard's method.

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

Cagniard Method

Louis Cagniard

Louis Cagniard (1900-1971) was a French geophysicist who was noted for his seminal works in divergent branches of mathematical geophysics. His classic text Reflection and refraction of progressive seismic waves (1938) introduced a clever transformation in the Fourier-Laplace domain (later modified by A. de Hoop) that allows exact solutions of wave equations for multilayer media to be obtained analytically through a process of contour integration.

Adrianus Hoop

Adrianus T. “Adrian” de Hoop received his M.Sc. degree in electrical engineering with the distinction cum laude from Delft University of Technology, the Netherlands, in 1950. In 1956–1957, he enjoyed a one-year research assistantship at the Institute of Geophysics of the University of California–Los Angeles on the invitation of its director, Louis B. Slichter. De Hoop carried out research on elastodynamic wave propagation and scattering in the Seismic Scattering Project, funded by a consortium of U. S. oil companies under the supervision of Leon Knopoff. Adrianus T. de Hoop is the founder of the Laboratory of Electromagnetic Research at Delft University. He pioneered a modification of the Cagniard technique for calculating impulsive wave propagation in layered media, now known as the Cagniard–de Hoop technique, a standard in the industry to analyze time-domain wave propagation.

We demonstrate the Cagniard--de Hoop method on the following example of the initial value problem for the two-dimensional wave equation:

\[ \square_c u =0 , \qquad (x,y) \in \mathbb{R}^2 , \qquad u(x,y,0) = f(x,y), \quad \dot{u}(x,y,0) = g(x,y) , \]
where □ is the d'Alembert operator
\[ \square_c u = \frac{\partial^2 u}{\partial t^2} - c^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) , \]
and overdot denotes the derivative with respect to time variable (Newton's notation). Application of the Laplace transform yields
\[ \lambda^2 u^L - g - \lambda f = c^2 \left( \frac{\partial^2 u^L}{\partial x^2} + \frac{\partial^2 u^L}{\partial y^2} \right) , \]
where
\[ u^L = {\cal L} \left[ u \right] (\lambda ) = \int_0^{\infty} u(x,y,t)\,e^{-\lambda\,t} \,{\text d} t \]
is the Laplace transform of unknown function u. Next, we apply the Fourier transform with respect to variable x. This leads to the following problem
\[ \left( \lambda^2 + c^2 k^2 \right) v = c^2 \,\frac{{\text d}^2 v}{{\text d} y^2} + g^F + \lambda\, f^F , \qquad y\in \mathbb{R} , \]
where
\[ v(y) = \int_{-\infty}^{\infty} u^L (x,y,\lambda )\, e^{{\bf j} x\,k}\,{\text d} x \]
is the Fourier transform of uL, and
\[ f^F (y) = 𝔉ℱ\left[ f \right] = \int_{-\infty}^{\infty} f (x,y)\, e^{{\bf j} x\,k}\,{\text d} x , \qquad g^F (y) = \int_{-\infty}^{\infty} g (x,y)\, e^{{\bf j} x\,k}\,{\text d} x . \]
We solve the nonhomogeneous differential equation for v using variation of parameters. For sake of compactness, we denote the inhomogeneous term as
\[ \varphi (y) = g^F (y) + \lambda \, f^L (y) . \]
If
\[ \delta = \sqrt{c^{-2} \lambda^2 + k^2} \]
denotes the positive branch of the square root, which is defined for positive entries, then its general solution becomes
\[ v(y) = A(y)\, e^{-\delta \,y} + B(y)\, e^{\delta\, y} . \]
For derivatives of functions A and B, we get the system of equations:
\[ \begin{cases} A' (y)\,e^{-\delta \,y} + B'(y)\, e^{\delta\, y} &= 0 , \\ -\delta\,A' (y)\,e^{-\delta \,y} + \delta\, B'(y)\, e^{\delta\, y} &= c^{-2} \varphi (y) . \end{cases} \]
Upon adding these two equations, we obtain
\[ 2\delta\, B'(y)\, e^{\delta\, y} = c^{-2} \varphi (y) , \qquad -2\delta\, A'(y)\, e^{-\delta\, y} = c^{-2} \varphi (y) . \]
Integration yields
\[ B(y) = - \frac{1}{2\,\delta}\, \int_y^{\infty} {\text d}\xi \,c^{-2} \varphi (\xi ) \,e^{-\delta\, \xi} , \qquad A(y) = - \frac{1}{2\,\delta}\, \int_{-\infty}^{y} {\text d}\xi \,c^{-2} \varphi (\xi ) \,e^{\delta\, \xi} . \]
So the solution for v becomes
\[ v(y) = - \frac{1}{2c^2\delta}\, e^{-\delta\, y} \int_{-\infty}^{y} {\text d}\xi \, \varphi (\xi ) \,e^{\delta\, \xi} - \frac{1}{2c^2 \delta}\, e^{\delta\, y} \int_y^{\infty} {\text d}\xi \, \varphi (\xi ) \,e^{-\delta\, \xi} . \]
Application of the inverse Fourier transform gives
\[ u^L (\lambda ) = - \frac{1}{2c^2}\,\frac{1}{2\pi} \, \int_{-\infty}^{\infty} {\text d}k\, e^{-{\bf j} xk} \,\frac{1}{\delta} \, \int_{-\infty}^{\infty} {\text d}\xi \, \varphi (\xi ) \,e^{-\delta\,|y- \xi |} . \]
The main idea of Cagniard--de-Hoop method is to transform the integral expression for uL in the form of the Laplace transform. Then the integrand will be the original required function u(x,y,t). To achieve it, we change the variable of integration
\[ k = \lambda \alpha \qquad \delta = \lambda\,\theta , \qquad {\text d}k = \lambda\,{\text d}\alpha , \]
where \( \displaystyle \theta = \sqrt{c^{-2} + \alpha^2} . \) If we set φ to be the Diral delta function δ. we get an auxiliary integral that will allow us to recover the full expression for u(x,y,t):
\[ w^L = \frac{1}{2\pi} \, \int_{-\infty}^{\infty} \frac{{\text d}k}{\delta}\, e^{-{\bf j} xk} \, e^{-\delta\, A} , \qquad A = |y-\xi | \ge 0. \]
Then
\[ w^L = \frac{1}{2\pi} \, \int_{-\infty /\lambda}^{\infty /\lambda} \frac{{\text d}\alpha}{\theta}\, e^{-{\bf j} x\,\lambda\,\alpha} \, e^{-\lambda \theta\, A} = \frac{1}{2\pi} \, \int_{-\infty /\lambda}^{\infty /\lambda} \frac{{\text d}\alpha}{\theta}\, e^{- \lambda \left( {\bf j} x\,\alpha + \theta\,A \right)} , \]
where
\[ \theta = \sqrt{c^{-2} + \alpha^2} . \]
Changing the variable α to the time t, we put
\[ t= \theta\,A + {\bf j} x\,\alpha \qquad \Longrightarrow \qquad \alpha = \frac{-c\,{\bf j}\,t \pm \sqrt{\left( A^2 + x^2 \right) \left( c^2 t^2 - A^2 \right) - c^2 t^2}}{c^2 \left( A^2 + x^2 \right)} . \]
Solve[t^2 - x^2 *a^2 - 2*t*I*a == (1/c^2 + a^2)*A , a]
a -> -((I (2 t + I Sqrt[-4 t^2 - 4 (-(A/c^2) + t^2) (-A - x^2)]))/( 2 (A + x^2)))

We deform the contour of integration, which in our case is a straight line \( \displaystyle \left( -\frac{1}{\lambda}\,\infty , \frac{1}{\lambda}\,\infty \right) , \) into a such contour for which the expression for t will have positive (real) values. It is defined by equation \( \displaystyle t= \theta\,A + {\bf j} x\,\alpha , \) where t is a parameter. Such transformation to a new contour, called Cagniard's contour, is very simple in our case because no poles or branch points are crossed.

Let us denote by L contour defined by the equation \( \displaystyle t= \theta\,A + {\bf j} x\,\alpha \) in complex plane ℂ. Upon setting in this equation α = 0, we obtain t = θ A. It is possible only when the radical in the formula for α is a pure imaginary number, so \( c^2 t^2 < A^2 . \) This condition a;;pws us to define contour L analytically:

\[ \alpha (x,A,t) = \begin{cases} \frac{-c\,{\bf j}\,t \pm \sqrt{\left( A^2 + x^2 \right) \left( c^2 t^2 - A^2 \right) - c^2 t^2}}{c^2 \left( A^2 + x^2 \right)} , & \quad \mbox{for} \quad ct > A, \\ \frac{-c\,{\bf j}\,t - {\bf j}\sqrt{\left( A^2 + x^2 \right) \left( A^2 - c^2 t^2 \right) + c^2 t^2}}{c^2 \left( A^2 + x^2 \right)} , & \quad \mbox{for} \quad ct < A, \end{cases} \]
Note that in the above formula, all roots take arithmetic values.
half = ParametricPlot[{0.2*Cos[t], 0.2*Sin[t]}, {t, Pi, 2*Pi}, PlotStyle -> {Thickness[0.02], Black}, Axes -> False];
left = Plot[1*t^2 + 1, {t, 0.2, 1}, PlotStyle -> {Thickness[0.02], Black}, Axes -> False, PlotRange -> {{0, 2}, {0, 2}}];
right = Plot[1*t^2 + 1, {t, -0.2, -1}, PlotStyle -> {Thickness[0.02], Black}, Axes -> False];
line1 = Graphics[{Thickness[0.02], Line[{{0.2, 0}, {0.2, 1}}]}];
line2 = Graphics[{Thickness[0.02], Line[{{-0.2, 0}, {-0.2, 1}}]}];
arrow1 = Graphics[{Arrowheads[0.07], Arrow[{{-1.5, 0}, {1.5, 0}}]}];
arrow2 = Graphics[{Arrowheads[0.07], Arrow[{{0, -0.3}, {0, 2.2}}]}];
txt1 = Graphics[ Text[StyleForm["Re \[Alpha]", FontSize -> 16, FontWeight -> "Bold"], {1.3, 0.2}]];
txt2 = Graphics[ Text[StyleForm["Im \[Alpha]", FontSize -> 16, FontWeight -> "Bold"], {0.3, 2.1}]];
pt = Graphics[{PointSize[Large], Point[{0, 1}]}];
tx = Graphics[ Text[StyleForm["M", FontSize -> 16, FontWeight -> "Bold"], {0.1, 1.15}]]; Show[line1, line2, left, right, half, arrow1, arrow2, txt1, txt2,pt,tx]
   Mathematica codes,     Cagniard's contour.

 

Therefore, the Cagniard's contour L is symmetric with respect to imaginary axis. The corresponding formula for α contains sign minus before radical for ℜα<0, and sign plus for ℜα>0. To fix a square root in θ, we make a cut along vertical axis connecting branch points -j/c and j/c through infinity.

Using Jordan lemma, we deform our line of integration into Cagniard's contour. With such new contour, we integrate twice along interval OM in opposite directions (compensate each other). As a result, we get

\begin{align*} w^L &= \frac{1}{2\pi} \int_{L_1} \frac{1}{\theta_1}\,\frac{{\text d}\alpha_1}{{\text d}t}\, e^{-\lambda\, t} {\text d}t + \frac{1}{2\pi} \int_{L_2} \frac{1}{\theta_2}\,\frac{{\text d}\alpha_2}{{\text d}t}\, e^{-\lambda\, t} {\text d}t \\ &= \frac{1}{2\pi} \int_{-\infty}^{-\sqrt{c^{-2} + \alpha^2}} \frac{1}{\theta_1}\,\frac{{\text d}\alpha_1}{{\text d}t}\, e^{-\lambda\, t} {\text d}t + \frac{1}{2\pi} \int_{\sqrt{c^{-2} + \alpha^2}}^{\infty}\frac{1}{\theta_2}\,\frac{{\text d}\alpha_2}{{\text d}t}\, e^{-\lambda\, t} {\text d}t , \end{align*}
where
\begin{align*} \alpha_1 &= \frac{-c\,{\bf j}\,t - \sqrt{\left( A^2 + x^2 \right) \left( c^2 t^2 - A^2 \right) - c^2 t^2}}{c^2 \left( A^2 + x^2 \right)} , \\ \alpha_2 &= \frac{-c\,{\bf j}\,t + \sqrt{\left( A^2 + x^2 \right) \left( c^2 t^2 - A^2 \right) - c^2 t^2}}{c^2 \left( A^2 + x^2 \right)} , \end{align*}
and the derivatives are
\[ \frac{1}{\theta_{1,2}} \,\frac{{\text d} \alpha_{1,2}}{{\text d}t} = \pm \frac{t}{\sqrt{t^2 - c^{-2} \left( x^2 + A^2 \right)}} . \]
Please do not ask me how I derived the above formulas---I just did not have time for this job, and I dedicated it to a computer algebra system. However, some genius people, such as Louis Cagniard and Adrianus Hoop did all calculations by hand.

The key point of the method is to modify above integrals into one of the form

\[ w^L = \int_{t_0 \ge 0}^{\infty} F(t)\,e^{-\lambda\, t}\,{\text d}t . \]
Then the inverse Laplace transform of wL is exactly F. The final answer is given below.
\[ u(x,y,t) = \frac{1}{2\pi} \,\frac{\partial}{\partial t} \left( t \int_{\xi_1^2 + \xi_2^2 \le 1} \frac{f \left( x+ ct\,\xi_1 , y+ ct\,\xi_2 \right)}{\sqrt{1 - \xi_1^2 - \xi_2^2}}\,{\text d}\xi_1 {\text d} \xi_2 \right) + \frac{t}{2\pi} \int_{\xi_1^2 + \xi_2^2 \le 1} \frac{g \left( x+ ct\,\xi_1 , y+ ct\,\xi_2 \right)}{\sqrt{1 - \xi_1^2 - \xi_2^2}}\,{\text d}\xi_1 {\text d} \xi_2 . \]
Upon introducing polar coordinates (r,φ), we obtain so called the Poisson formula
\begin{align*} u(x,y,t) &= \frac{1}{2\pi\, c}\,\frac{\partial}{\partial t} \int_0^{ct} r\,{\text d}r \int_0^{2\pi} {\text d}\varphi \, \frac{f \left( x+r\,\cos\varphi , y + r\,\sin\varphi \right)}{\sqrt{c^2 t^2 - r^2}} \\ & \quad + \frac{1}{2\pi\, c}\, \int_0^{ct} r\,{\text d}r \int_0^{2\pi} {\text d}\varphi \, \frac{g \left( x+r\,\cos\varphi , y + r\,\sin\varphi \right)}{\sqrt{c^2 t^2 - r^2}} . \end{align*}

 

Initial value problem for 3D wave equation

Consider the three-dimensional wave equation subject to the initial conditions: two-dimensional wave equation:
\[ \square_c u =0 , \qquad (x,y,z) \in \mathbb{R}^3 , \qquad u(x,y,z,0) = f(x,y,z), \quad \dot{u}(x,y,z,0) = g(x,y,z) , \]
where □ is the d'Alembert operator
\[ \square_c u = \frac{\partial^2 u}{\partial t^2} - c^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \right) . \]
Application of the Laplace transform with respect to time variable t yields
\[ \frac{\partial^2 u^L}{\partial x^2} + \frac{\partial^2 u^L}{\partial y^2} + \frac{\partial^2 u^L}{\partial z^2} - c^{-2} \lambda^2 u^L = c^{-2} \varphi (x,y,z) , \]
where
\[ u^L = {\cal L} \left[ u(x,y,z,t) \right] = \int_0^{\infty} u(x,y,z,t)\,e^{-\lambda\, t}\,{\text d} t \]
is the Laplace transform of the unknown function u(x,y,z,t) and
\[ \varphi (x,y,z) = \lambda\,f(x,y,z) + g(x,y,z) . \]
Next, we apply the double Fourier transform with respect to variables x,y:
\[ \frac{{\text d}^2 v}{{\text d}z} - \left( c^{-2} \lambda^2 + k^2 + q^2 \right) v = c^{-2} \varphi^{FF} , \]
where
\[ v(z) = 𝔉ℱ\left[ u^L \right] = \int_{\mathbb{R}^2} u^L (x,y,z,\lambda )\,e^{{\bf j} \left( k\,x + q\,y \right)} \,{\text d} x\,{\text d} y , \]
and
\[ \varphi^{FF} (z) = \int_{-\infty}^{\infty} {\text d} x \int_{-\infty}^{\infty} {\text d} y\, \varphi (x,y,z)\, e^{{\bf j} \left( k\,x + q\,y \right)} . \]
Solution of the above second order differential equation that satisfies the regularity condition at infinity is exactly the same as it was done previously:
\[ v(z) = - \frac{1}{2c^2 \delta} \,\int_{-\infty}^{\infty} {\text d}\zeta \, \varphi^{FF} (\zeta )\, e^{-\delta |z- \zeta |} , \]
where
\[ \delta (k,q) = \sqrt{c^{-2} + k^2 + q^2} . \]
Now we apply the inverse Fourier transform and obtain
\[ u^L (x,y,z,\lambda ) = \frac{1}{4\pi^2} \int_{\mathbb{R}^2} {\text d}k\,{\text d}q \, v(z) \, e^{-{\bf j} \left( x\,k + q\,y \right)} . \]
To simplify our job, we set φ to be the Diral delta function to obtain an auxiliary integral:
\[ w^L (x,y,z,\lambda ) = \frac{1}{4\pi^2} \int_{\mathbb{R}^2} {\text d}k\,{\text d}q \, \frac{1}{\delta}\,e^{-\delta\,z} \, e^{-{\bf j} \left( x\,k + q\,y \right)} . \]
In the above integral, we make change of variables
\[ k=\lambda\,\alpha , \qquad q = \lambda\, \beta ,\qquad \delta = \lambda\,\sqrt{c^{-2} + \alpha^2 + \beta^2} . \]
This allows us to rewrite the integral in compact form
\[ w^L (x,y,z,\lambda ) = \frac{\lambda}{4\pi^2} \int_{\mathbb{R}^2} {\text d}\alpha\,{\text d}\beta \, \frac{1}{\theta}\,e^{-\lambda \left( \theta\,z + {\bf j}\,\alpha \,x + {\bf j}\,\beta \,y \right)} . \]
Since multiplication by λ corresponds differentiation of the original, we consider another auxiliary function
\[ w^L_0 (x,y,z,\lambda ) = w^L (x,y,z,\lambda ) /\lambda \qquad \Longrightarrow \qquad w(x,y,z,t) = \frac{\partial}{\partial t}\, w_0 (x,y,z,t) . \]
In double inetgral for w0, we make change the variables
\[ \alpha = \xi\,\frac{x}{r} - \eta\,\frac{y}{r}, \qquad \beta = \xi\,\frac{y}{r} + \eta\,\frac{x}{r} \qquad\mbox{if } y≥0, \]
and
\[ \alpha = \xi\,\frac{x}{r} - \eta\,\frac{|y|}{r}, \qquad \beta = -\xi\,\frac{|y|}{r} + \eta\,\frac{x}{r} \qquad\mbox{if } y≤0, \]
where
\[ r = \sqrt{x^2 + y^2} . \]
The Jacobian of such transformation is 1. Then
\[ w^L_0 (x,y,z,\lambda ) = \frac{1}{4\pi^2} \int_{\mathbb{R}^2} {\text d}\xi\,{\text d}\eta \, \frac{1}{\theta}\,e^{-\lambda \left( \theta\,z + {\bf j}\,\zi\,r \right)} \]
because \( \displaystyle x\,\alpha + y\,\beta = r\,\xi . \) Now we deform the line of integration over variable ξ into Cagniard's contour:
\[ w^L_0 (x,y,z,\lambda ) = \frac{1}{4\pi^2} \int_L e^{-\lambda\,t} \,{\text d}\xi \,\int_{-\infty}^{\infty} {\text d}\eta \,\frac{1}{\theta}\,{\text d}\eta , \]
where
\[ t = z \sqrt{c^{-2} + \xi^2 + \eta^2} + {\bf j}\, r\,\xi \qquad\mbox{and} \qquad \theta = \sqrt{c^{-2} + \xi^2 + \eta^2} . \]
Then ξ is expressed through t as
\[ \xi = \frac{-{\bf j}\,r\,t \pm \sqrt{t^2 - R^2 \left(\eta^2 + c^{-2} \right)}}{R^2} , \qquad R^2 = z^2 + r^2 = x^2 + y^2 + z^2 . \]
Then we proceed exactly in the same way as in two dimensional case. The only difference is that the branch point for θ depends on the variable of integration. Using the formula that the derivative of the Heaviside function H(t - a) with respect to (t) is the delta function δ(t - a), we obtain

 

In spherical coordinates (θ,φ), we get so called the Kirkhhoff formula
\begin{align*} u(x,y,z,t) &= \frac{1}{4\pi} \,\frac{\partial}{\partial t}\,t \int_0^{2\pi} {\text d}\varphi \,\int_0^{\pi} \sin\theta\,{\text d}\theta\,f\left( x+ct\,\sin \theta\,\cos\varphi , y+ ct\,\sin\theta\,\sin\varphi , z + ct\,\sin \theta\, \sin\varphi \right) \\ &\quad + \frac{1}{4\pi} \,t \int_0^{2\pi} {\text d}\varphi \,\int_0^{\pi} \sin\theta\,{\text d}\theta\,g\left( x+ct\,\sin \theta\,\cos\varphi , y+ ct\,\sin\theta\,\sin\varphi , z + ct\,\sin \theta\, \sin\varphi \right) . \end{align*}

 

  1. Achenbach, J.D., Wave Propagation in Elastic Solids, North Holland Publishing Company, 1973.
  2. Cagniard, L. (1962) Reflection and refraction of progressive seismic waves McGraw-Hill. Translation of Reflexion et refraction des ondes sismiques progressives" translation by Edward A. Flinn and C. Hewitt Dix.
  3. Berenger, J.-P., Three-dimensional perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics, 1996, Volume 127, Issue 2, September 1996, Pages 363-379; https://doi.org/10.1006/jcph.1996.0181
  4. Cagniard, L. Basic theory of the magneto-telluric method of geophysical prospecting, Geophysics, 1953, Vol. 18, Issue 3, pp. 605--635.
  5. De Hoop, A.T., List of publications A. T. De Hoop - Papers in Refereed Journals
  6. De Hoop, A.T., A modification of Cagniard's method for solving seismic pulse problems, Applied Scientific Research, B8 (1960) pp. 349--356.
  7. De Hoop, A.T., Frankena, H. J., Radiation of pulses generated by a vertical electric dipole above a plane, non-conducting, earth, Applied Scientific Research, B8 (1960) pp. 369--377.
  8. Diaz, J., Joly, P., An Analysis of Higher Order Boundary Conditions for the Wave Equation. [Research Report] RR-4962, INRIA. 2003. ￿inria-00071617￿
  9. Dobrushkin, V.A., The boundary value problems in the dynamic theory of elasticity for the wedge-shaped domains, 1988, Minsk: Nauka i Technika (in Russian).
  10. Jingu, T., Matsumoto, H., Nezu, K., The transient stress in an elastic half-space excited by impulsive loading over one quarter of its surface, Bulletin of JSME, 1986, Volume 29, Issue 247, pp. 44--51.
  11. Parnell, W.J., Nguyen, V-H., Assier, R., Naili, S., and Abrahams, I.D., Transient thermal mixed boundary value problems in the half-space, 1918, https://arxiv.org/pdf/1507.06291.pdf
  12. Štumpf, M. and De Hoop, A.T., An application of the Cagbiard--de Hoop technique for solving initial-boundary value problems in bounded regions, The Quarterly Journal of Mechanics and Applied Mathematics, 2013, Volume 66, Issue 2, pp. 185--197. https://doi.org/10.1093/qjmam/hbs025

 

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