A variety of numerical methods may be used to solve (dispersionrelation), including for example Crout's reduction method (Crout, 1941). However, since the system is relatively small () and since our purposes include gaining some physical insight into the processes involved, it will prove instructive to do some more analysis on the problem prior to the ultimate numerical calculations.

First, note that our analysis will be considerably simplified by an obvious rearrangement of the determinant (dispersionrelation) so that

(- v^2) = 0, where is the identity matrix and

= ^-1^-1^-1,
having complex matrix elements *g*_{ij}, for *i*,*j* = 1,2,3.
The matrix was given previously in (variablechange)
and is clearly invertible for finite values of the two porosities.
Matrix is also always invertible
for realistic choices of material parameters.

Using the properties of determinants, it is not hard to show that the determinant (rearrangement) can be reduced to a determinant in either of two convenient forms:

&&
(g_11-v^2)(g_22-v^2) - g_12g_21 & -1
g_12g_23g_31 + g_13g_32g_21
- (g_22-v^2)g_13g_31 - (g_11-v^2)g_23g_32 & g_33 - v^2
&=&
(g_11-v^2)(g_33-v^2) - g_13g_31 & -1
g_12g_23g_31 + g_13g_32g_21
- (g_33-v^2)g_12g_21 - (g_11-v^2)g_23g_32 & g_22 - v^2
&=& 0.
It is useful to write the determinant in each of these two ways
in order to make connections with single-porosity models.
The first version in (newdet) has as its upper left-hand
element an expression that corresponds precisely to the determinant
for the system when the only porosity present is the matrix (or storage) porosity.
Similarly, the second version in (newdet) has as its element in the upper
left-hand position a term that corresponds to the determinant for the
system when the only porosity present is the
fracture/crack porosity. In each case the remaining terms determine
the effect of coupling to the second type of porosity. Depending on
other system parameters such as permeabilities and porosities, either of
these limits may be useful to consider, and may also provide good
starting points for the iterative (Newton-Raphson) procedure that we will
use to solve the determinant equation for the complex velocity *v*.

If we express the determinant (rearrangement)
as a polynomial *D*(*x*) where *x* = *v ^{2}*, then

D(x) = -x^3 &+& (g_11+g_22+g_33)x^2 &-& (g_22g_33+g_33g_11+g_11g_22 &-& g_23g_32 - g_13g_31 - g_12g_21)x &+& (g_11g_22g_33 - g_23g_32g_11 - g_13g_31g_22 &-& g_12g_21g_33 + g_12g_23g_31 + g_13g_32g_21) = 0.
The Newton-Raphson method (Hildebrand, 1956)
for solving this equation for *x* is an
iteration process starting from some initial choice *x ^{(0)}* and computing

x^(i)=x^(i-1)-D(x^(i-1))/D'(x^(i-1)), for i = 1,2,
..., N_c
(where *D*' is the first derivative of *D* with respect to *x*)
until some convergence criterion has been met at *i* = *N*_{c}.
The fact that the coefficients and the solution of the problem are
complex adds no special complication to this procedure. However,
since the polynomial is complex, it is important that good
starting values *x ^{(0)}* be obtained for at least two of the three roots of
(polynomial), as a search procedure in the complex plane
would be considerably more difficult to implement than is pure
Newton-Raphson iteration.

From (newdet), we see that two choices of starting values for the
complex parameter *x* are given by (for example)

x_1,2^(0) = 12(g_11+g_22
(g_11-g_22)^2 + 4g_12g_21).
Once the Newton-Raphson iteration has converged from both of these
starting values to their final values of *x _{1}* and

D(x) _n=0^3 D_n(-x)^n = -(x-x_1)(x-x_2)(x-x_3).
Using the linear independence of the terms in powers of *x* in
(factoring), we have three equations showing that

x_3 = D_2 - x_1 - x_2 = D_1-x_1 x_2x_1+x_2 =
D_0x_1 x_2.
Any one of these three formulas may be used to calculate *x _{3}* directly, or
combinations of them (such as the geometric mean of any two)
may be used to reduce the errors that might be introduced
by premature termination of the Newton-Raphson iteration process
for

x_2,3 = 12[(D_2 - x_1) (D_2 - x_1)^2 -
4D_0/x_1].
In this approach it is clearly preferable to solve for the fast wave
solution as *x _{1}*, and then use the formulas in (p2p3) to obtain the
two slow wave solutions since the fast wave solution will virtually
always be well-behaved.

For each value of *x* that solves the dispersion equation,
we can then compute the wave velocity and inverse of the quality factor
*Q* using the definition

1x 1v()
1v_r(1 + i2Q),
where the choice of square root is determined by
the requirement of physical realizability as discussed in the
following section.
This definition of 1/*Q* is accurate when attenuation losses are low,
but should be carefully interpreted for high loss situations
(*e.g.,* low frequency diffusive modes).
In particular, it is often stated that, when losses are high,
(vandQ) needs to be modified so that if

1v()
1v_r(1 + iv_r),
where is the attenuation coefficient (having units of inverse length),
then [as Hamilton (1972) and Bourbié *et al.* (1987) show]

1Q 2v_r- ^2 v_r^2/. However, for diffusive modes the imaginary and real parts of the velocity are of comparable size, and therefore it is possible that the formula (highlossQ) will become singular at low frequencies for such modes. This complication can and does happen in practice. To avoid this complication, we use the definition (vandQ) in all cases, and then interpret those situations in which as an indication that the mode under consideration is actually diffusive rather than propagatory.

4/20/1999