Consider the second order ODE \[ \mathcal L y = \frac{d}{dx}\left(p(x) \frac{dy}{dx} \right) + s(x) y = f(x) \] subject to boundary conditions at \(x=a\) and \(x=b\) such that the operator \(\mathcal L\) be Hermitian. Note that the operator has the same form as the Sturm-Liouville problem, but the ODE is inhomogeneous and doesn't define an eigenvalue problem.
The preceding definitions mean that we are only interested at the solution \(y(x)\) for \(a\le x \le b\). Assume that all functions take in reals as inputs and complex numbers as outputs, and that we have homogeneous boundary conditons (if we multiply the solution by a constant we still satisfy the boundary conditions).
Let's consider now the related ODE \[ \mathcal L G(x,\xi) = \delta(x-\xi), \] that is, the inhomogenous part is the Dirac delta. The advantage of doing this is that, if we solve this ODE instead of the original problem, we can find the solution to the original problem as \[ y(x) = \int_a^b G(x,\xi) f(\xi) d\xi \] because if we apply opon this equation the operator \(\mathcal L\) we have \[ \mathcal L y(x) = \int_a^b \mathcal L G(x,\xi) f(\xi) d\xi = \int_a^b \delta(x-\xi) d\xi = f(x). \] The function \(G(x,\xi)\) is called the Green function. If we solve for the Green function of an inhomogeneous ODE, we can then build the solution \(y(x)\) for any inhomogeneous function \(f(x)\) without repeating the calculation as long as \(\mathcal L\) and the boundary conditions remain the same.
Before we can study how to solve a Green function problem, we have to explore the properties of such a function around \(x=\xi\), where the delta goes to infinity. First, we take the Green ODE and integrate in a small interval around \(x=\xi\) (with \(\epsilon \to 0\)): \[ \int_{\xi-\epsilon}^{\xi+\epsilon} \frac{d}{dx} \left[ p(x) \frac{dG(x,\xi)}{dx} \right]dx + \int_{\xi-\epsilon}^{\xi+\epsilon} s(x) G(x,\xi) dx = \int_{\xi-\epsilon}^{\xi+\epsilon} \delta(\xi-x)dx \]
The first term on the l.h.s. can be integrated through the fundamental theorem of calculus, and the r.h.s. is 1 because of the definition of the delta: \[ \left. p(x) \frac{dG(x,\xi)}{dx} \right|_{\xi-\epsilon}^{\xi+\epsilon} + \int_{\xi-\epsilon}^{\xi+\epsilon} s(x) G(x,\xi) dx = 1 \] The remaining integral is zero in the limit \(\epsilon\to 0\) (the area under the curve is zero from \(\xi\) to \(\xi\)). If, in general, the function \(p(x)\) is continuous at \(\xi\) (it converges to the same value when evaluated from the left and the right) then we must have \[ \lim_{\epsilon \to 0} \left[ \left.\frac{dG(x,\xi)}{dx}\right|_{x=\xi+\epsilon} - \left.\frac{dG(x,\xi)}{dx}\right|_{x=\xi-\epsilon} \right] = \frac{1}{p(\xi)}, \] which means that the derivative of \(G(x,\xi)\) is discontinuous at \(x=\xi\). The derivative is discontinuous but the \(G\) is continuous. To picture the situation, imagine the graph of the function \(|x|\) (the derivative or slope has two constant values, -1 and 1, for negative and positive values of x, respectively, but both lines meet at \(x=0\) so the function is continuous while its derivative is discontinuous).
Another property of the Green function can be found if we expand it as a series of base eigenfunctions \(\varphi_n(x)\) of the Sturm-Liouville operator \(\mathcal L\). The expansion of \(G(x,\xi)\) is \[ G(x,\xi) = \sum_{mn} g_{nm}\varphi_n(x) \bar\varphi_m(\xi) \] and the expansion of the Dirac delta is (see PDF "Exercises of Dirac delta") \[ \delta(x-\xi) = \sum_m \varphi_m(x) \bar\varphi_m(\xi). \] We now apply the Sturm-Liouville operator \(\mathcal L \varphi_n = \lambda_n \varphi_n\) and find for the Green ODE \[ \sum_{mn}\lambda_n g_{mn} \varphi_n(x) \bar\varphi_m(\xi) = \sum_{m}\varphi_m(x) \bar\varphi_m(\xi) \] \[ \implies g_{nm} = \delta_{nm}/\lambda_n \] \[ \implies G(x,\xi) = \sum_{n} \frac{\bar\varphi_n(\xi) \varphi_n(x)}{\lambda_n} \] Finally, we see the symmetry property of the Green function \[ G(x,\xi) = \bar G(\xi,x) \]
For solving the inhomogeneous Green equation \(\mathcal L G(x,\xi) = \delta(x-\xi)\) we use the method of variation of parameters, the symmetry of the Green function and the discontinuity of the derivative to find the following recipe:
A 1D boundary-value problem \(\mathcal L y = f\) can be solved with via the Green function method with the following procedure:
Suppose we have an initial value problem \[ \mathcal L y \equiv a(t)y''(t) + b(t)y'(t) + c(t)y(t) = f(t) \] subject to the initial conditions \(y(0) = y_0\) and \(y'(0)=v_0\). We can find the Green function \(G(t,\tau)\) for the initial value problem by the method of variation of parameters. Let's assume we know the homogenous solution \[ y_h(t) = c_1y_1(t)+c_2y_2(t) \] and we want the particular solution \(y_p(t)\). We vary the parameters \[ y_p(t) = c_1(t) y_1(t) + c_2(t) y_2(t) \] and substitute. In the PDF "Exercises of Green function", there is the demonstration of how this method produces the Green function for this case. However, we can also directly apply the two properties of the Green function (symmetry and discontinuity of the derivative) to find the same answer. Let's call this method finding the Green function by the definition.
We have the two base functions \(y_1(t)\) and \(y_2(t)\), found by solving the homogeneous ODE. In an initial value problem, we only have one boundary value, at \(t=0\). This can be used to find the Green function \(G(t,\tau)\) for the region \(t<\tau\). For the region \(t>\tau\), we are free to vary the parameters as \[ G(t,\tau) = c_1(\tau) y_1(t) + c_2(\tau) y_2(t) \] Then the conditions to be imposed are:
The conditions yield a system of two equations with two unknowns, \(c_1\) and \(c_2\). The notation \(X_<\) and \(X_>\) is very common in Green functions and means the criterion for \(t<\tau\) and \(t>\tau\) respectively, for a variable \(X\).
In 3D problems, the following techniques have to be applied: