## 4 Computation of PDs of Surface-Wave Group Velocities
### Fully Analytical Method
 \tag{12}\]  
Therefore, we differentiate the last expression above with respect to \(c\) and thus obtain  
\[\partial^{2}F/\partial c^{2}=\frac{\partial\overline{\textbf{Y}}_{n+1}}{ \partial c}\overline{\textbf{V}}_{n+1}+2\overline{\textbf{Y}}_{n+1}\frac{ \partial\overline{\textbf{V}}_{n+1}^{*}}{\partial c}+\overline{\textbf{X}}_{ n+1}\frac{\partial^{2}\overline{\textbf{V}}_{n+1}^{*}}{\partial c^{2}}. \tag{13}\]  
Then, in order to obtain \(\partial\overline{\textbf{Y}}_{n+1}/\partial c\), we shall resort to the second expression in (12):  
\[\partial\overline{\textbf{Y}}_{m+1}/\partial c=\frac{\partial\overline{\textbf {Y}}_{m}}{\partial c}\overline{\textbf{T}}_{m}^{*}+2\overline{\textbf{Y}}_{ m}\frac{\partial\overline{\textbf{T}}_{m}^{*}}{\partial c}+\overline{\textbf{X}}_{ m}^{*}\frac{\partial^{2}\overline{\textbf{T}}_{m}^{*}}{\partial c^{2}},m=1,2, \cdots,n. \tag{14}\]  
Similarly, if we introduce a new vector \(\overline{\textbf{Z}}_{m}^{*}=\partial\overline{\textbf{Y}}_{m}/\partial c\), then (14) can be reformulated as follows:  
\[\overline{\textbf{Z}}_{m+1}^{*}=\overline{\textbf{Z}}_{m}^{*}\overline{ \textbf{T}}_{m}^{*}+2\overline{\textbf{Y}}_{m}^{*}\frac{\partial\overline{ \textbf{T}}_{m}^{*}}{\partial c}+\overline{\textbf{X}}_{m}^{*}\frac{\partial^{2 }\overline{\textbf{T}}_{m}}{\partial c^{2}},m=1,2,\cdots,n. \tag{15}\]  
Obviously, \(\overline{\textbf{Z}}_{1}^{*}\) is a zero vector as well. Thus, by taking advantage of (12), (15), and (13), we can finally implement the analytical computation of \(\partial^{2}F/\partial c^{2}\) (Table 3). See Appendix C for the elements of \(\partial^{2}\overline{\textbf{T}}_{m}^{*}/\partial c^{2}\) and of \(\partial^{2}\overline{\textbf{V}}_{n+1}^{*}/\partial c^{2}\) for Love waves. \begin{table}
\begin{tabular}{l} \hline
**Pseudo-code 3** The computation of \(\partial^{2}F/\partial c^{2}\) \\ \hline
**Input**: angular frequency \(\omega\), phase velocity \(c\), earth model \(\textbf{Q}_{\text{R}}\) or \(\textbf{Q}_{\text{L}}\) \\ \hline
1: Initialize six vectors: \(\textbf{xold}=\overline{\textbf{X}}_{1}\), \(\textbf{yold}=\overline{\textbf{Y}}_{1}\), \(\textbf{zold}=\overline{\textbf{Z}}_{m}^{*}\), \(\textbf{xnew}=\textbf{0}\), \(\textbf{ynew}=\textbf{0}\), and \(\textbf{znew}=\textbf{0}\) \\
2: **for**\(m=1\)-\(n\)**do** \\
3: calculate \(\overline{\textbf{T}}_{m}^{*}\), \(\partial\overline{\textbf{T}}_{m}^{*}/\partial c\) and \(\partial^{2}\overline{\textbf{T}}_{m}^{*}/\partial c^{2}\) \\
4: update **xnew**: \(\textbf{xnew}=\textbf{xold}\times\overline{\textbf{T}}_{m}\) \\
5: update **ynew**: \(\textbf{ynew}=\textbf{yold}\times\overline{\textbf{T}}_{m}+\textbf{xold}\times \partial\overline{\textbf{T}}_{m}/\partial c\) \\
6: update **znew**: \(\textbf{znew}=\textbf{zold}\times\overline{\textbf{T}}_{m}+2\times\textbf{ yold}\times\partial\overline{\textbf{T}}_{m}/\partial c+\textbf{xold}\times \partial^{2}\overline{\textbf{T}}_{m}/\partial c^{2}\) \\
7: update **xold**: \(\textbf{xold}=\textbf{xnew}\) \\
8: update **yold**: \(\textbf{yold}=\textbf{ynew}\) \\
9: update **zold**: \(\textbf{zold}=\textbf{znew}\) \\
10: **end for** \\
11: compute \(\overline{\textbf{V}}_{n+1}\), \(\partial\overline{\textbf{V}}_{n+1}^{*}/\partial c\) and \(\partial^{2}\overline{\textbf{V}}_{n+1}/\partial c^{2}\) \\
**Output**: \(\partial^{2}F/\partial c^{2}=\textbf{znew}\times\overline{\textbf{V}}_{n+1}+2 \times\textbf{ynew}\times\partial\overline{\textbf{V}}_{n+1}/\partial c+\textbf{ xnew}\times\partial^{2}\overline{\textbf{V}}_{n+1}/\partial c^{2