B.4. Solving integral and differential equations¶
We encounter a variety of integral and differential equations in the study of galaxies. While these can never be solved analytically for fully realistic galaxy models, they are often amenable to analytical treatment under assumptions that are quite reasonable and, therefore, analytical solutions play an important role in gaining insights about the astrophysics of galaxies. Here, we discuss some techniques for solving these equations that occur in a variety of contexts in this book.
B.4.1. Abel integral equation¶
We will start with an integral equation that comes up a surprising amount in the study of stellar systems: \begin{equation}\label{eq-math-intdiffeq-abel-original} f(x) = \int_0^x\,\mathrm{d}y\,{g(y)\over \sqrt{x-y}}\,. \end{equation} This is an Abel integral equation, because it appears in the solution of the Abel problem, which is the problem of finding the shape of a wire (\(\propto g[y]\)) such that a bead slides down for a given vertical height \(x\) in a given time (\(\propto f[x]\)). This integral can be inverted to \begin{equation}\label{eq-math-intdiffeq-abel-inverted} g(y) = {1\over \pi}\,{\mathrm{d} \over \mathrm{d} y} \int_0^y\,\mathrm{d}x\,{f(x)\over \sqrt{y-x}}\,. \end{equation} To show this, we plug Equation (B.96) into Equation (B.97) and switch the order of integration, giving \begin{equation} g(y) = {1\over \pi}\,{\mathrm{d} \over \mathrm{d} y} \int_0^y\,\mathrm{d}t\,g(t)\,\int_t^y\mathrm{d}x\,{1\over \sqrt{(y-x)(x-t)}}\,. \end{equation} The inner integral can be written as \(\int_0^1\mathrm{d} z/\sqrt{z\,(1-z)} = \pi\) and \(g(y) = g(y)\) then immediately follows. A more general inversion also holds: \begin{align}\label{eq-math-intdiffeq-abel-original-general} f(x) & = \int_0^x\,\mathrm{d}y\,{g(y)\over (x-y)^\alpha}\,,\quad g(y) = {\sin(\pi\alpha)\over \pi}\,{\mathrm{d} \over \mathrm{d} y} \int_0^y\,\mathrm{d}x\,{f(x)\over (y-x)^{1-\alpha}}\,, \end{align} for \(0 < \alpha < 1\). This can be demonstrated in the same way as above, except that you now need that \(\int_0^1\mathrm{d} z/[z^\alpha\,(1-z)^{1-\alpha}] = \pi/\sin(\pi \alpha)\). Similar inversions exist when the integrals go from a lower value to \(\infty\), except that the inversion includes a minus sign then, that is \begin{align}\label{eq-math-intdiffeq-abel-original-general-xtoinf} f(x) & = \int_x^\infty\,\mathrm{d}y\,{g(y)\over (y-x)^\alpha}\,,\quad g(y) = -{\sin(\pi\alpha)\over \pi}\,{\mathrm{d} \over \mathrm{d} y} \int_y^\infty\,\mathrm{d}x\,{f(x)\over (x-y)^{1-\alpha}}\,, \end{align} again for \(0 < \alpha < 1\). All of these inversions are referred to as Abel inversions.
B.4.2. Ordinary differential equations¶
We encounter many ordinary differential equations in this book—differential equations that only involved derivatives with respect to a single variable. Some of these can be solved analytically using the techniques discussed in this section. Most of the remaining ones can be solved numerically without great difficulty with standard tools; we discuss some of these in Chapter 4.5.1 in the context of numerical orbit integration.
Ordinary differential equations are called linear when the unknown function and its derivatives only appear linearly. They are called homogeneous when the equation does not include any terms that do not involve the unknown function. Linear homogeneous differential equations with constant coefficients have the form \begin{equation} a_0\,y+a_1\,y'+a_2\,y''+\ldots + a_n\,y^{(n)} = 0\,, \end{equation} where the \(a_i\) are constants and the prime(s) and \((n)\) superscript denote differentiation with respect to the dependent variable \(x\) (\(n\) times in the case of the \((n)\) superscript). These equations can be solved by obtaining the characteristic equation \begin{equation} a_0+a_1\,t+a_2\,t^2+\ldots + a_n\,t^n = 0\,. \end{equation} If the \(n\) roots of this equation are \(t_i\), then the solutions of the differential equation are \(y(x) = e^{t_i\,x}\) (note that when a root appears multiple times, it corresponds to a set of solutions \(y(x) = x^k\,e^{t_i\,x}\), where \(k\) ranges from zero to the root’s multiplicity minus one). For example, the harmonic-oscillator equation \begin{equation} \omega^2\,y+y'' = 0\,, \end{equation} has the characteristic equation \begin{equation} \omega^2+t^2= 0\,, \end{equation} with roots \(\pm i\,\omega\) and, thus, solutions \(y(x) = e^{\pm i \omega\,x} = \sin(\omega x)\) and \(\cos(\omega x)\). Because the differential equation is linear in the unknown function, we can multiply each solution with an arbitrary number and still obtain a solution. The \(n\) coefficients of the solutions are, therefore, determined by the initial or boundary conditions of the problem. For example, the general solution of the harmonic-oscillator equation is \(y(x) = A\,\cos(\omega x) + B\,\sin(\omega x)\); if we specify that the oscillator is maximally extended up to \(y=C\) at positive \(y\) at \(x=0\), then the solution becomes \(y(x) = C\,\cos(\omega x)\).
Inhomogeneous linear differential equations contain one or more terms that do not involve the unknown function (but that can involve the dependent variable): \begin{equation} a_0\,y+a_1\,y'+a_2\,y''+\ldots a_n\,y^{(n)} = f(x)\,. \end{equation} Such differential equations can be solved by combining a particular solution with the solutions of the homogeneous form of the equation (that with \(f = 0\)). A particular solution is any solution that does not contain arbitrary constants and satisfies the inhomogeneous equation. For example, adding a constant-forcing term to the harmonic-oscillator equation \begin{equation} \omega^2\,y+y'' = F\,, \end{equation} a particular solution is \(y(x) = F/\omega^2\) and the full solution of the equation is \(y(x) = A\,\cos(\omega x) + B\,\sin(\omega x)+F/\omega^2\). The value of the arbitrary constants \(A\) and \(B\) is again set by the initial or boundary conditions of the problem. We often encounter versions of this equation where the forcing term is sinusoidal as well \begin{equation} \omega^2\,y+y'' = F\,\cos\left(\Omega x\right)\,. \end{equation} A particular solution of this equation is \(y(x) = F\,\cos\left(\Omega x\right)/[\omega^2-\Omega^2]\) and the full solution is \(y(x) = A\,\cos(\omega x) + B\,\sin(\omega x)+F\,\cos\left(\Omega x\right)/[\omega^2-\Omega^2]\), with constants \(A\) and \(B\) set by the initial/boundary conditions. The \(\omega^2-\Omega^2\) in the denominator gives rise to resonances at \(\omega \approx \Omega\) where the sinusoidal forcing with frequency \(\Omega\) interferes constructively with the system’s natural frequency \(\omega\) to create a very large response. When the inhomogeneous contribution includes multiple terms representing constant forcing and sinusoidal forcing or multiple different sinusoidal components, we can simply add up the particular solutions of each forcing term separately.
Ordinary differential equations with non-constant coefficients are much harder to solve analytically and generally require a numerical treatment. However, some important special cases can be solved analytically. For example, if we can re-write a differential equation as \begin{equation}\label{eq-math-intdiffeq-diffeq-linear-separate} f(y)\,y' = g(x)\,, \end{equation} where \(f(\cdot)\) and \(g(\cdot)\) are arbitrary functions, then we can solve this equation by writing it as \(f(y)\,\mathrm{d}y = g(x)\,\mathrm{d}x\), and integrating each side \begin{equation} \int \mathrm{d}y\,f(y) = \int \mathrm{d}x\,g(x) + \mathrm{constant}\,. \end{equation} If these integrals can be performed analytically, then we obtain an analytical solution. We can leverage this result to solve any first-order (involving only \(y'\) and \(y\)) linear differential equation. First, we reduce the equation to the form \begin{equation}\label{eq-math-intdiffeq-diffeq-linear-general} y' + p(x)\,y = g(x)\,, \end{equation} where \(p(x)\) and \(g(x)\) are arbitrary functions. Then we define the integrating factor \(\gamma(x)\) through \begin{equation}\label{eq-math-intdiffeq-diffeq-linear-integrating-factor} \gamma(x) = \exp\left[\int\mathrm{d} x'\,p(x')\right]\,, \end{equation} which is such that \(\gamma'(x) = \gamma(x)\,p(x)\). Multiplying Equation (B.110) by \(\gamma(x)\), we find that \begin{equation}\label{eq-math-intdiffeq-diffeq-linear-general-deriv1} \gamma(x)\,y' + \gamma(x)\,p(x)\,y = \gamma(x)\,y' + \gamma'(x)\,y =\gamma(x)\,g(x)\,, \end{equation} such that we have \begin{equation}\label{eq-math-intdiffeq-diffeq-linear-general-deriv2} \tilde{y}' = \gamma(x)\,g(x)\,, \end{equation} where \(\tilde{y} \equiv \gamma(x)\,y\). This equation is of the form of Equation (B.108) and can thus be solved as \begin{equation}\label{eq-math-intdiffeq-diffeq-linear-general-result} y(x) = {\int\mathrm{d}x\,\gamma(x)\,g(x)+ \mathrm{constant} \over \gamma(x)}\,, \end{equation} with \(\gamma(x)\) given by Equation (B.111). The integration constant involved in the definition of \(\gamma(x)\) can be absorbed into the general constant in Equation (B.114).
.
.
.