4 separation theory

Separation Theory PhL 8.17.11 1. Introduction 1 2. Notation, Cast of Players and References 2 3. Initial Processing of the Helmholtz Equation. 5 Comments on (3.5) and example of toroidal coordinates 6 Resume Flow: the traditional separation does not work at this point 7 Resume Flow and define Q 8 Compute Q for the toroidal system using (3.9) 9 Finish Flow 10 4. Starting from the other end: The Stackel Matrix. 11 Introduction of the Stackel Matrix. 11 5. Simple separation 12 Solving the Cofactor/Robertson Conditions: Non-uniqueness of the Stackel Matrix Φ 14 Theorem of Equivalent Stackel Matrices: 16 A Viability Condition for Simple Separability 16 How do we find the first column of the Stackel matrix Φ ? 17 Toroidal: 17 Spherical: 18 Summary of Simple Separation 18 6. R-separation 19 1. Introduction The subject here is separation of the Helmholtz equation (and its special case the Laplace equation) in various 3D curvilinear coordinate systems. The term "separation" means that you start with the 3D Helmholtz equation (2 + k12)ψ=0, which is of course a PDE, and you try to find solutions of the form ψ(ξ1,ξ2,ξ3) = X1(ξ1)X2(ξ2)X3(ξ3)/ R(ξ1,ξ2,ξ3) where you can produce an ODE for each of the functions Xn(ξn). Thus, you can search for solutions "separately" for each Xn. If Xn is any solution to its ODEn for n = 1,2,3, then the ψ shown above is a solution of the Helmholtz equation. This is useful as a way to solve the problem in the first place, and secondly, because boundary conditions are often specified separately for the three coordinates, this separation turns each ODE into its own little 1D boundary value problem. In fact, the solution ODE's always contain self-adjoint differential operators and are therefore amenable to normal 1D "Sturm-Liouville theory". If the separation can be done as outlined above for a certain orthogonal curvilinear coordinate system, then we say that the system is R-separable. If this can be done with R=1, it is simple separable. It turns out that only the 11 "classical" 3D Euclidean orthogonal coordinate systems are simple separable for the Helmholtz equation (and therefore also for the Laplace equation). These systems are discussed in the first chapter (called Section I) of Moon & Spencer. It also turns out that with R ≠ constant, the Helmholtz PDE is not separable unless k12 = 0 in which case it is the Laplace PDE. In other words, no Helmholtz equation with k12≠0 is R-separable unless R = constant. In this case one might as well assume R=1 since the rest of a constant could be absorbed into the Xi functions. The upshot is that there are systems like toroidal coordinates (not one of the classical 11 systems) in which the Laplace equation is R-separable but not simple separable. The Helmholtz equation is extremely significant because it arises very naturally in problems involving the heat conduction equation and the wave equation, where the time derivative term in the PDE is replaced by a constant parameter by applying a Laplace or Fourier time transform to the PDE. A huge swath of mathematical physics is dominated by these two PDE equation types, one always needs to solve the Helmholtz equation that results, and that then involves the notion of "separation". 2. Notation, Cast of Players and References The entire analysis is a study of functional forms and this makes it a bit slippery. By functional form we mean the way a function of multiple variables depends on those variables in terms of possible factorization of the form. For example, some u(x,y,z) might be expressible as v(x)w(y,z) and this would then be the functional form of u. (Most functions of course cannot be factored in this manner.) Another sense of functional form is simply what variables a function is a function of. We might write v(x) simply as v if our expressions get complicated, but we must remember than that its functional form is v(x). The symbols ξ1, ξ2, ξ3 (these are "xi", pronounced "zeye", different from ζ = zeta and χ = chi ) are often used for ellipsoidal coordinates since those coordinates are so closely related to each other. Here we shall use ξ1, ξ2, ξ3 to represent an arbitrary triplet of orthogonal curvilinear coordinates. We shall use the notation ∂n for a partial derivative, = ∂f/∂ξn = ∂nf n = 1,2,3 and the following other symbols Σn ≡ Σn=13 LHS = left hand side of some equation RHS = right hand side PDE = partial differential equation ODE = ordinary diff eq There will be many function symbols used below, and each symbol has an implied functional form. At first one needs to show the forms in full detail, but eventually one learns to work without doing this. Here is a triplet of function symbols where we show four different notations for each, g1(ξ2, ξ3) = g1(23) = g1(≠1) = g1 g2(ξ3, ξ1) = g2(31) = g2(≠2) = g2 g3(ξ1, ξ2) = g3(31) = g3(≠3) = g3 Each row is the forward cyclic permutation of the previous row. Only the notations of the last two columns are amenable to the generic notation gn(≠n) = gn n = 1,2,3 Notice these facts: ∂1g1 = 0 ∂1(g1F) = g1(∂1F) ∂ngn = 0 ∂n(gnF) = gn(∂nF) These follow, for example, since g1 is a function only of ξ2 and ξ3. "Facts" like these will be crucial in the analysis below. The functions gn are "helper functions" which have no particular significance. Only one other function symbol set will have this same functional form Mn(≠n) = Mn n = 1,2,3 The letter M stands for "minor" as in the minor of a 3x3 matrix, but in fact M is really a cofactor. This confusion seems to go back to Morse and Feshbach who use the word "minor" to mean what we usually call "cofactor". In current terminology cofactorpq = (-1)p+q minorpq where p and q label the rows and columns of a matrix, as in Apq. To be a little more specific, the Mn will be the cofactors of the elements of the first column of a certain 3x3 matrix called the Stackel matrix Φ which we will discuss soon below. Certain function symbols imply a very simple functional form as follows fn(ξn) = fn(n) = fn n = 1,2,3 Xn(ξn) = Xn(n) = Xn n = 1,2,3 We already know what the Xn are from the Introduction. The fn functions are just more helper functions that will be used in our analysis and which appear in the separated ODE equations. There are several function symbols that are assumed to have in general no factored functional form, and here they are R(ξ1,ξ2,ξ3) = R(123) = R // the function appearing in (1) above ("modulation factor") Q(ξ1,ξ2,ξ3) = Q(123) = Q // a helper function (called u by Morse & Feshbach) S(ξ1,ξ2,ξ3) = S(123) = S // the determinant of Φ (coming soon) hn(ξ1,ξ2,ξ3) = hn(123) = hn // the curvilinear system scale factors hn2 = gnn (metric tensor) H(ξ1,ξ2,ξ3) = H(123) = H ≡ h1h2h3 The symbol H is "non standard" but I found it convenient to use it as the product of the three curvilinear scale factors. We might say that all the quantities list here have a generic (123) functional form. This does not mean they cannot also have some kind of factored form. The symbol ψ, our Helmholtz equation solution function, has the special functional form noted in the Introduction (that is, we seek solutions ψ of this functional form) ψ(ξ1,ξ2,ξ3) = X1(ξ1)X2(ξ2)X3(ξ3)/ R(ξ1,ξ2,ξ3) ψ(123) = X1(1)X2(2)X3(3)/ R(123) ψ = X1X2X3/R Historically it seems that R was defined "in the denominator of ψ". Certain constants shall appear below: k12 = the parameter appearing in the Helmholtz equation (2 + k12)ψ=0 k22, k32 = two generic constants which will be called "separation constants" Although these three kn2 are written "squared", they are really meant to be arbitrary real numbers, so if some ki2 < 0 then we imagine that the corresponding ki is imaginary. The squared notation arises from the way the Helmholtz equation looks when it arises from a transformed wave equation. We now come finally to the Stackel matrix Φ (and its determinant S) Φ = S = det(Φ) The subscripts on Φnm are "standard" where the first index n is the row index and the second m the column index, and the first value of each index is 1. A key fact to recognize is that the three functions of row n are functions only of ξn ! Sometimes we will indicate this fact by using this notation Φnm(n) = Φnm(ξn) so the first index matches the argument index. The cofactors of the elements of the first column of the Φ matrix are these, Mn ≡ cof(Φn1) = (-1)n+1 minor(Φn1) For example, M2 = (-1)2+1 minor(Φn1) = – = – [Φ12(ξ1) Φ33(ξ3) – Φ13(ξ1) Φ32(ξ3)] Notice that M2 has the functional form M2(≠2) = M2(ξ1,ξ3), and in general Mn = Mn(≠n), as it was presented earlier in this section. Sometimes we might use the following notations Mn(Φ) S(Φ) to stress that these quantities are functions of the Stackel matrix elements. This Φ matrix is named after German mathematician Paul Stäckel (shtay'kle), 1962-1919. Our symbols match exactly those of Morse and Feshbach except for Q which they call u, and except for our added symbol H for h1h2h3. Our symbols match Moon & Spencer except they use Un(un) in place of our Xn(ξn) Morse and Feshbach address this 3D separation subject in two places in their Volume 1: simple separation is treated pp 508-511, while R-separation is treated pp 518-519. They do not use the term R-separation and call the R function a "modulation factor". Moon and Spencer discuss simple separation on pp 5-7, and R-separation on page 96. They provide more detail in some of their other books. 3. Initial Processing of the Helmholtz Equation. We are going to treat R-separation and simple separation at the same time, and branch off later into the two cases. Remember that simple separation just means R = 1. In this section, we shall set the Helmholtz parameter to K12 instead of k12 for a reason which will be seen later. In orthogonal curvilinear coordinates, using our notation as defined above, the Helmholtz equation can be written (2+K12)ψ = 0 : H-1{ ∂1[(H/h12)(∂1ψ)] + cyclic}ψ + K12ψ = 0 H-1 Σn ∂n[(H/hn2)(∂nψ)] + K12ψ = 0 (3.1) where sometimes the cyclic form is more useful for observing functional forms. We are seeking a solution of this functional form ψ = X1X2X3/R (3.2) where the four functions Xn and R are as yet unknown. It is they we seek to find! Inserting this form into (3.1) gives H-1{ ∂1[(H/h12)(∂1{ X1X2X3/R })] + cyclic}ψ + K12ψ = 0 H-1{ X2X3∂1[(H/h12) ∂1{ X1/R }] + cyclic}ψ + K12ψ = 0 (3.3) Sometimes we shall "show every line" so the reader can clearly see the changes made. Here we have just used a kind of "fact" noted earlier and you see how the implicit functional forms come into play. We next insert into (3.3) with ∂1(X1/R) = ∂1(R-1X1) = R-1∂1X1 - R-2X1∂1R = R-2{ (∂1X1)R - X1(∂1R) } to get H-1{ X2X3∂1[ (H/[R2h12]){ (∂1X1)R - X1(∂1R) } ] + cyclic}ψ + k12ψ = 0 (3.4) In order to move toward a "separated form" wherein the terms above are less coordinate-entangled, we shall attempt to select function R so the following equation is satisfied in terms of functional form (comments below): (H/[R2hn2]) = fn(n)gn(≠n) n = 1,2,3 (3.5) (H/[R2h12]) = f1(1)g1(23) // for example because, if we insert this into (3.4), we can pull g1(23) to the left through ∂1 to get H-1{ X2X3∂1[ f1(1)g1(23) { (∂1X1)R - X1(∂1R) } ] + cyclic}ψ + K12ψ = 0 H-1{ X2X3 g1∂1[ f1 { (∂1X1)R - X1(∂1R) } ] + cyclic}ψ + K12ψ = 0 and, recalling that ψ = X1X2X3/R from (3.2), we divide by ψ to get (R/H){(1/X1) g1∂1[f1 { (∂1X1)R - X1(∂1R) }] + cyclic} + K12 = 0 (R/H) Σn[(gn/Xn) ∂n[fn { R(∂nXn) - Xn(∂nR) }] + K12 = 0 If we now use (3.5) to replace gn in favor of fn, we get (1/R) Σn[ (1/[hn2Xn]) (1/fn)∂n[fn { R(∂nXn) - Xn(∂nR) }] + K12 = 0 (3.6) This is at least starting to look like a separated form. Comments on (3.5) and example of toroidal coordinates We now comment as promised on our assumed conditions (3.5) above (1/R2) (H/hn2) = fn(1)gn(23) n = 1,2,3 which requires that the LHS factor in a certain specific manner. We can certainly find an R(123) that works for the first equation with n=1: we could just select two arbitrary functions f1(1) and g1(23) and then define R by [R(123) ]-2 ≡ f1(1)g1(23) (H/h12) But if we do this, it is unlikely that (3.5) will be viable for n=2 and 3! So we have to assume that we can find a set of 7 functions (three fn, three gn and one R) which make (3.5) valid. Any curvilinear system for which we cannot find a happy set of 7 functions is therefore not R-separable. So we view (3.5) as a restriction or condition on our curvilinear system which must be met to obtain any separation. We do not address the question of whether there might be some other possible separation solution where (3.5) is not assumed. Since this is perhaps a confusing concept, let's pause to examine the toroidal coordinate system as an example. We will use the Moon & Spencer page 112 notation where η labels toroids, θ labels bowls, ψ labels azimuthal half planes, a is the radius of the limiting toroid, and we assume (ξ1,ξ2,ξ3) = (η,θ,ψ). We find that (azimuthal ψ does not even appear in the scale factors) R-2 ≡ [ch(ξ1)-cos(ξ2)] h1 = h2 = a/[ch(ξ1)-cos(ξ2)] => h1 = h2 = aR2 h3 = a sh(ξ1)/[ch(ξ1)-cos(ξ2)] => h3 = a sh(ξ1)R2 where R is our candidate for a viable R function. We then find that (H/h12) = h2h3/h1 = h3 = a sh(ξ1)/[ch(ξ1)-cos(ξ2)] (H/h22) = h3h1/h2 = h3 = a sh(ξ1)/[ch(ξ1)-cos(ξ2)] (H/h32) = h1h2/h3 = a2/[ch(ξ1)-cos(ξ2)]2 * [ch(ξ1)-cos(ξ2)] /a sh(ξ1) = a/{[ch(ξ1)-cos(ξ2)] sh(ξ1)} We now attempt to find a set of 6 functions fn and gn which satisfy (3.5), fn(1)gn(23) = (1/R2) (H/hn2) = [ch(ξ1)-cos(ξ2)] (H/hn2) Here we go: f1(1)g1(23) = [ch(ξ1)-cos(ξ2)] a sh(ξ1)/[ch(ξ1)-cos(ξ2)] = a sh(ξ1) = [sh(ξ1) ] [a ] f2(2)g2(31) = [ch(ξ1)-cos(ξ2)] a sh(ξ1)/[ch(ξ1)-cos(ξ2)] = a sh(ξ1) = [ 1 ] [a sh(ξ1) ] f3(3)g3(12) = [ch(ξ1)-cos(ξ2)] a/{[ch(ξ1)-cos(ξ2)] sh(ξ1)} = a /sh(ξ1) = [ a ] [1/ sh(ξ1) ] Thus with this candidate R we have found our 6 functions. Here is a solution set f1(1) = sh(ξ1) g1(23) = a R = [ch(ξ1)-cos(ξ2)]-1/2 f2(2) = 1 g2(31) = a sh(ξ1) f3(3) = a g3(12) = 1/ sh(ξ1) Therefore, the toroidal system at least has a chance of being R-separable. It is in fact R-separable, but we state it this way since we might find other restrictions that will need checking. If in the above toroidal discussion we had tried R = 1, our problem would have been to find 6 functions fn and gn which satisfy (3.5), which says for the case n = 1, f1(1)g1(23) = (H/h12) = a sh(ξ1)/[ch(ξ1)-cos(ξ2)] Our only chance is to select f1(1) = a sh(ξ1) but then we are stuck with g1(23) = 1/[ch(ξ1)-cos(ξ2)] which involves coordinate ξ1 which violates the functional form g1(23). Therefore the toroidal system is not simple-separable! From now on the gn helper functions will no longer appear in our equations, but we must remember that they must be identified at some point so that (3.5) may be verified. Resume Flow: the traditional separation does not work at this point We now resume our processing of the Helmholtz equation where we left off, which was here (1/R) Σn[ (1/[hn2Xn]) (1/fn)∂n[fn { R(∂nXn) - Xn(∂nR) }] + K12 = 0 (3.6a) where we have assumed (3.5). Now ψ is "gone" and we can see that we are slowly moving toward some form that might be separable, but is not yet so since variables other than n appear in the nth summand. These appear in R(123) and hn(123). One approach would be to somehow manipulate (3.6) to get it into this form L1(X1)/X1 + L2(X2)/X2 + L3(X3)/X3 = constant (3.7) where Ln is a differential operator in just ξn. As just noted, (3.6a) is most definitely NOT in this form -- the coordinates are still "entangled". But if we could somehow obtain this form (3.7), we could conclude by the usual argument that each term on the LHS is a constant ci2 so we would then have Ln(Xn)/Xn = cn and c1 + c2 + c3 = constant and we would identify two of the three ci2 as "separation constants" in the traditional sense. The "usual argument" referred to above is as follows. If in (3.7) we fix ξ2 and ξ3 and allow ξ1 to vary wildly in its coordinate range (whatever that might be), we can conclude that L1(X1)/X1 must be a constant. If it varied, that would be a contradiction since the other 3 terms in (3.7) are constant in this little scenario. Thus Ln(Xn)/Xn = cn . We are not going to take this search path right now and we return to the Stackel approach. Resume Flow and define Q To further process (3.6a), we compute the ∂n derivatives, ∂n[R fn(∂nXn)] = ∂n[R{fn(∂nXn)}] = R∂n{fn(∂nXn)} + (∂nR) {fn(∂nXn)} ∂n[fnXn(∂nR)] = ∂n[Xn{fn(∂nR)}] = Xn∂n{fn(∂nR)} + (∂nXn) {fn(∂nR)} Therefore we find for the quantity appearing in (3.6), ∂n[fn { R(∂nXn) - Xn(∂nR) })] = R∂n{fn(∂nXn)} – Xn∂n{fn(∂nR)} since the two second terms cancel. Then (3.6a) becomes (1/R) Σn (1/[hn2Xnfn]) [R∂n{fn(∂nXn)} – Xn∂n{fn(∂nR)}] + K12 = 0 (1/R) Σn[ (1/[hn2Xnfn]) [R∂n{fn(∂nXn)}] – (1/R) Σn[ (1/[hn2Xnfn]) [Xn∂n{fn(∂nR)}] + K12 = 0 Σn[ (1/[hn2Xnfn]) [∂n{fn(∂nXn)}] – Σn[ (1/[hn2Rfn]) [∂n{fn(∂nR)}] + K12 = 0 (3.8) Notice that we have managed to segregate the Xn part from the R part. Since we have already determined our viable set of 4 functions fn and R by the method illustrated above in the toroidal example, we can compute the R sum shown in (3.8). It is going to be some general function of 123 which we are going to write in a fairly strange manner: Σn (1/[hn2fnR]) [∂n{fn(∂nR)}] = - k12/Q(123) (3.9) Here, it is the RHS which is fully determined by the LHS and we choose to partition the RHS into two factors, a constant -k12 (M&S call this -α1) and 1/Q(123). Remember that our Helmholtz parameter is K12 and has nothing to do with this new constant k12 we just introduced. In theory, k12 can be any constant, but one usually chooses it to make the resulting Q(123) have some simple form. When we consider the simple separation case with R = 1, we will choose Q(123) = 1 and k12 = 0 in (3.9), which of course makes it still valid since (∂nR) = 0. Compute Q for the toroidal system using (3.9) We start with the numerators. First, (∂1R) = ∂1[ch(ξ1)-cos(ξ2)]-1/2 = (-1/2)sh(ξ1) [ch(ξ1)-cos(ξ2)]-3/2 = (-1/2)sh(ξ1) R3 (∂2R) = ∂2[ch(ξ1)-cos(ξ2)]-1/2 = (-1/2)sin(ξ2) [ch(ξ1)-cos(ξ2)]-3/2 = (-1/2)sin(ξ2) R3 (∂3R) = 0 => entire third term in the sum is 0, so ignore it from now on Then we do this very tedious algebra which we show in detail for this example: [∂1{f1(∂1R)}] = [∂1{ sh(ξ1) * (-1/2)sh(ξ1) R3}] = (-1/2) ∂1[sh2(ξ1)R3] = (-1/2) [ 2sh(ξ1)ch(ξ1) R3 + sh2(ξ1)3R2(∂1R) ] = (-1/2) [ 2sh(ξ1)ch(ξ1) R3 + sh2(ξ1)3R2(-1/2)sh(ξ1) R3 ] = (-1/2) [ 2sh(ξ1)ch(ξ1) R3 – sh3(ξ1)(3/2)R5 ] = (-1/4)sh(ξ1)R5 [ 4ch(ξ1) R-2 – 3sh2(ξ1) ] = (-1/4)sh(ξ1)R5 [ 4ch(ξ1) [ch(ξ1)-cos(ξ2)] – 3sh2(ξ1) ] = (-1/4)sh(ξ1)R5 [ 4ch2(ξ1) - 4 ch(ξ1) cos(ξ2)] – 3sh2(ξ1) ] = (-1/4)sh(ξ1)R5 [ 3 + ch2(ξ1) - 4 ch(ξ1) cos(ξ2)] [∂2{f2(∂2R)}] =[∂2{1 *(∂2R)}] = ∂2[(-1/2) sin(ξ2) R3 = (-1/2)[ cos(ξ2)R3 + sin(ξ2)3R2 (∂2R) ] = (-1/2)[ cos(ξ2)R3 + sin(ξ2)3R2 (-1/2)sin(ξ2) R3 ] = (-1/2)[ cos(ξ2)R3 – sin2(ξ2)(3/2)R5 ] = (-1/2) R5 [ cos(ξ2)R-2 – sin2(ξ2)(3/2)] = (-1/4) R5 [ 2cos(ξ2) [ch(ξ1)-cos(ξ2)] – 3sin2(ξ2)] = (-1/4) R5 [ 2cos(ξ2) ch(ξ1)-2cos2(ξ2) – 3sin2(ξ2)] = (-1/4) R5 [ 2cos(ξ2) ch(ξ1) - 2 - sin2(ξ2)] which we summarize as [∂1{f1(∂1R)}] = (-1/4) R5 [ 3 + ch2(ξ1) - 4 ch(ξ1) cos(ξ2)] sh(ξ1) [∂2{f2(∂2R)}] = (-1/4) R5 [ 2cos(ξ2) ch(ξ1) - 2 - sin2(ξ2)] Next we shall need the denominator factors h12f1R = a2R4 sh(ξ1) R = a2 sh(ξ1) R5 h22f2R = a2R4 1 R = a2 R5 Finally we can evaluate the LHS of (3.9) to get Σn (1/[hn2fnR]) [∂n{fn(∂nR)}] = [a2 sh(ξ1) R5]-1 * (-1/4)sh(ξ1)R5 [ 3 + ch2(ξ1) - 4 ch(ξ1) cos(ξ2)] + [a2 R5]-1 * (-1/4) R5 [ 2cos(ξ2) ch(ξ1) - 2 - sin2(ξ2)] = (-1/4a2) [ 3 + ch2(ξ1) - 4 ch(ξ1) cos(ξ2)] + (-1/4a2) [ 2cos(ξ2) ch(ξ1) - 2 - sin2(ξ2)] = (-1/4a2) [3 + ch2(ξ1) - 4 ch(ξ1) cos(ξ2)+ 2cos(ξ2) ch(ξ1) - 2 - sin2(ξ2) ] = (-1/4a2) [1 + ch2(ξ1) - 4 ch(ξ1) cos(ξ2)+ 2cos(ξ2) ch(ξ1) - sin2(ξ2) ] = (-1/4a2) [cos2(ξ2) + ch2(ξ1) - 2 ch(ξ1) cos(ξ2)] = (-1/4a2) [ch(ξ1)-cos(ξ2)]2 = (-1/4a2) R-4 Therefore we have shown that -k12/Q(123) = Σn (1/[hn2fnR]) [∂n{fn(∂nR)}] = (-1/4a2) R-4 = (-1/4) / {a2R4} so we can make this simple partitioning k12 = α1 = (1/4) Q = a2R4 = a2[ch(ξ1)-cos(ξ2)]2 which we are happy to see agrees with Moon & Spencer page 97. We have done a brute force computation of Q(123) in the toroidal system and when lots of algebraic dust settles out we find that it is a very simple function. We can summarize our toroidal information to this point in case we want to further develop this example: h1 = h2 = a/[ch(ξ1)-cos(ξ2)] => h1 = h2 = aR2 H = sh(ξ1)a3R6 h3 = a sh(ξ1)/[ch(ξ1)-cos(ξ2)] => h3 = a sh(ξ1)R2 f1(1) = sh(ξ1) g1(23) = a R = [ch(ξ1)-cos(ξ2)]-1/2 f2(2) = 1 g2(31) = a sh(ξ1) Q = a2[ch(ξ1)-cos(ξ2)]2 = a2R4 f3(3) = a g3(12) = 1/ sh(ξ1) k12 = (1/4) Finish Flow So here is where we were prior to the above example: Σn(1/[hn2fn Xn]) [∂n{fn(∂nXn)}] –Σn 1/[hn2fnR] [∂n{fn(∂nR)}] + K12 = 0 (3.8) Σn (1/[hn2fnR]) [∂n{fn(∂nR)}] = - k12/Q (3.9) which we now combine to get Σn(1/[hn2Xn]) (1/fn)∂n[fn(∂nXn)] + k12/Q + K12 = 0 (3.10) Notice in (3.10) the position of the two copies of function fn. Equation (3.10) is still not separated because hn= hn(123) and Q = Q(123). At least in the simple separation case the Q term goes away, since as noted above we select k12 = 0 and Q = 1. This concludes our "initial processing of the Helmholtz equation". Equation (3.10) essentially is the Helmholtz equation (2+K12)ψ=0 with the various assumptions we have made: the form ψ = X1X2X3/R, equation (3.5), and the definition of Q and k12 from (3.9). At this point, then, for a given curvilinear coordinate system, we know all these items in full detail: K12 hn(123) H(123) R(123) fn(n) gn(≠n) Q(123) k12 4. Starting from the other end: The Stackel Matrix. Our goal is to find differential operators Ln in terms of which we can separate (3.10). The most general second order linear Ln can be written this way LnXn = (1/pn)∂n[pn(∂nXn)] + rn (∂nXn) + qnXn = 0 (4.1) where the first derivative term is here mixed into two terms. If we want this Ln operator to be formally self-adjoint, which means (u,Lv) = (Lu,v) where we ignore the "parts" in the parts integrations that make this be true, where (u,v) indicates an appropriate scalar product in a function Hilbert space, then we must at once set rn = 0 since L = ∂nXn and L* = - ∂nXn and so this term cannot be self-adjoint. We are then reduced to this candidate form for our separated equations: LnXn = (1/pn)∂n[pn(∂nXn)] + qnXn = 0 Looking now at (3.10) we are highly motivated to set pn = fn and then our candidate form is LnXn = (1/fn)∂n[fn(∂nXn)] + qnXn = 0 (4.2) which says (1/Xn)(1/fn)∂n[fn(∂nXn)] = - qn (4.3) If we insert (4.3) into (3.10), we get - Σn(1/hn2)qn + k12/Q + K12 = 0 (4.4) This is the form of our processed Helmholtz equation if we make the assumption that Xn satisfies the equation (4.2) Introduction of the Stackel Matrix. Now comes a fairly inspired and unexpected step in the development. Suppose we choose to write the qn(n) function in this manner, which is to say, the unknown function qn(n) is simply written as a linear combination of three other unknown functions, -qn(n) = [ κ12Φn1(n) + k22Φn2(n) + k32Φn3(n)] n = 1,2,3 (4.5) where κ12, k22 and k32 are at the moment three arbitrary real constants, and the Φnm(n) are the 9 functions of what is called the Stackel matrix mentioned in the Notation section above, Φ = S(Φ) ≡ det(Φ) (4.6) Mn(Φ) ≡ cof(Φn1) = (-1)n+1 minor(Φn1) (4.7) With (4.5) used in (4.2), our assumed form for Ln becomes LnXn = (1/fn)∂n[fn(∂nXn)] + [ κ12Φn1(n) + k22Φn2(n) + k32Φn3(n)]Xn = 0 (4.8) and (4.4) becomes Σn(1/hn2) [ κ12Φn1(n) + k22Φn2(n) + k32Φn3(n)] + k12/Q + K12 = 0 (4.9) As a reminder, this is the form our processed Helmholtz equation takes if we expand the qn of (4.4) as shown in (4.5). Our goal is to find a set of Φnm which makes (4.9) be true. It is not obvious that such a set of Φnm exists. The sum is still entangled due to hn(123). At this point we finally need to branch into our two paths: simple separation and R-separation.. 5. Simple separation In this case we have R = 1, of course, and from (3.9) we then make the choice Q = 1 and k12= 0 which is consistent with LHS (3.9) since (∂nR) = 0. Then the processed Helmholtz equation (4.9) becomes Σn(1/hn2) [ κ12Φn1(n) + k22Φn2(n) + k32Φn3(n)] + K12 = 0 For simple separation we choose κ12 = K12 (the Helmholtz parameter) to get Σn(1/hn2) [ K12Φn1(n) + k22Φn2(n) + k32Φn3(n)] + K12 = 0 (5.1) and then our separated equation (4.8) becomes LnXn = (1/fn)∂n[fn(∂nXn)] + [ K12Φn1(n) + k22Φn2(n) + k32Φn3(n)]Xn = 0 (5.2) If we can find the 9 element functions Φnm of the Stackel matrix, then (5.2) gives us the separated ODE for function Xn and our simple separation of the Helmholtz equation will have succeeded! Moreover, our ODE will contain two free real parameters k22 and k32 which we shall call "separation constants". Since all three constants are free parameters, in that we expect (5.1) to be satisfied for arbitrary real values of all three, we must set the coefficients of the constants to 0 to get Σn (1/hn2) Φn1(n) = 1 should be -1 ? Σn (1/hn2) Φn2(n) = 0 Σn(1/hn2) Φn3(n) = 0 To put this into standard matrix equation, we need to use [ΦT]mn(n) ≡ Φnm (n)so we have Σn ΦT1n(n) (1/hn2) = 1 Σn ΦT2n(n) (1/hn2) = 0 Σn ΦT3n(n) (1/hn2) = 0 and our matrix equation is then ΦT V = W where (5.3) V = W= Φ = S ≡ det(Φ) As shown in Appendix A, from linear algebra we know that, as long as S ≠ 0, our matrix equation can be inverted to give V = [ΦT]-1 W = [1/det(Φ)] C W = (1/S) C W where C is the "cofactor matrix" Cpq ≡ cof(Φpq) (5.4) The elements of matrix C are the cofactors of elements of Φ. In components we then have Vn = (1/S) ΣmCnmWm In our case we have from (4.12) (5.3) Vn = (1/hn2) and Wm = δm,1 so our solution is then (1/hn2) = (1/S) Cn1 Following tradition, we define Mn = Cn1 = cof(Φn1) = the cofactor of the first element in the nth row of Φ (5.5) and finally we arrive at the Cofactor Condition for simple separation, Mn/S = (1/hn2) (5.6) Alternatively, we can write (4.13) as Mn(Φ)/S(Φ) = (1/hn2) " Cofactor condition" (5.6a) and this emphasizes that we have 3 conditions that the Φ matrix must satisfy for a given coordinate system with scale factors hn. Notice that this condition arises from our assumed form for Ln in (5.2), it was not "imposed". Our next step is to in fact impose an extra condition on the Φ matrix known as the Robertson Condition. The requirement is that S = det(Φ) take this form, S(Φ) = H / (f1f2f3) "Robertson condition" (5.7) It is not obvious that, for a given coordinate system, we can find a matrix Φ which satisfies the Cofactor and Robertson conditions, but let's try to construct a solution Φ. Combining (5.6) and (5.7) we get Mn = (S/hn2) = (H/hn2) /(f1f2f3) (5.8) We see that the Robertson condition lets us write the solution Mn directly in terms of known quantities. If we can continue on to find a full solution matrix Φ, then this condition is retroactively justified. Using (3.5) we can write this last result in an alternate form Mn = gnfn/(f1f2f3) eg M1 = g1/(f2f3) (5.8a) Solving the Cofactor/Robertson Conditions: Non-uniqueness of the Stackel Matrix Φ To reduce our symbol writing load, let's write (beware overloaded symbols) Φ = S = det(Φ) so that we can read off the cofactors of the first column elements, M1(23) = e(2)i(3)-f(2)h(3) -M2(31) = b(1)i(3)-c(1)h(3) M3(12) = b(1)f(2)-c(1)e(2) which we combine with (4.17) to get (H/h12) /(f1f2f3) = e(2)i(3)-f(2)h(3) - (H/h22) /(f1f2f3) = b(1)i(3)-c(1)h(3) (H/h32) /(f1f2f3) = b(1)f(2)-c(1)e(2) Can we find 6 functions {b,c,e,f,h,i} that satisfy these three equations? The answer hinges on the functional form of the left hand sides. In order for there to be a solution, the LHS's must have these functional forms: (H/h12) /(f1f2f3) = E(2)I(3)-F(2)H(3) - (H/h22) /(f1f2f3) = B(1)I(3)-C(1)H(3) (H/h32) /(f1f2f3) = B(1)F(2)-C(1)E(2) (5.9) which using 3.5 with R=1 can also be written (H/h12) /(f1f2f3) = E(2)I(3)-F(2)H(3) - (H/h22) /(f1f2f3) = B(1)I(3)-C(1)H(3) (H/h32) /(f1f2f3) = B(1)F(2)-C(1)E(2) (5.9) If (4.18) is valid, then our problem is to solve these equations for {b,c,e,f,h,i} e(2)i(3)-f(2)h(3) = E(2)I(3)-F(2)H(3) b(1)i(3)-c(1)h(3) = B(1)I(3)-C(1)H(3) b(1)f(2)-c(1)e(2) = B(1)F(2)-C(1)E(2) (5.10) We consider this generic form of the above problem (using temporary new symbols) a(1)b(2)-c(1)d(2) = A(1)B(2)-C(1)D(2) If we didn't have the functional forms to worry about, we would say that the solution to this problem was that the matrix could be any matrix which has the same determinant as , which would certainly allow many solutions. But the functional form requirement severely limits the number of solutions to the following, which are pretty self-evident, solution #1: a(1) = α A(1) b(2) = B(2)/α c(1) = β C(1) d(2) = D(2)/β solution #2: a(1) = -α C(1) b(2) = D(2)/α c(1) = -β A(1) d(2) = B(2)/β where in either solution α and β are arbitrary constants. If we apply the conclusion of this generic example to the three lines of (5.10), we conclude that, given a solution set of functions {b,c,e,f,h,i}, we can generate new solution sets by using either or both of the following operations: (1) swap the last two columns of Φ and then negate either of these columns. (2) multiply one of the last two columns by a constant α and the other by 1/α. These operations each maintain the values of the three Mn and of S. We shall now add another operation to the above list. Suppose we have found a set of solutions { b,c,e,f,h,i } for the rightmost two columns of Φ Φ = We can generate a new set of solutions { b',c',e',f',h',i' } by adding a multiple of either of the last two columns to the other. We know that such an operation does not change the determinant S, and we can show that the cofactors of the last column are also unchanged. For example, let's make a new third column by adding a multiple of the second column to it, c'(1) = c(1) + α b(1) f'(2) = f(2) + α e(2) i'(3) = i(3) + α h(3) Then for example M1'(23) = e(2)i'(3)-f'(2)h(3) = e(2)[ i(3) + α h(3)]-[ f(2) + α e(2)]h(3) = e(2)i(3)-f(2)h(3) + α [e(2) h(3) - e(2) h(3)] = e(2)i(3)-f(2)h(3) = M1(23) and this will be similarly true for the other two cofactor Mn. This gives a new solution set { b,c',e,f',h,i' }. Let's now summarize our "lack of uniqueness" rules for the last two columns of the Φ matrix: (1) swap the last two columns of Φ and then negate either of these columns. (2) multiply one of the last two columns by a constant α and the other by 1/α. (3) add any constant multiple of one of the last two columns to the other. All three of these operations maintain both the Mn and S. We could also add a multiple of one of the last two columns to the first column: this gives a new solution Φ matrix, but does not change { b,c,e,f,h,i }. We cannot add a constant multiple of the first column to one of the last two columns, however. So let's state our results one more time: Theorem of Equivalent Stackel Matrices: If we find a Stackel matrix Φ which satisfies our Cofactor and Robertson conditions for S and the Mn, then we can generate alternative equivalent Stackel matrices using the following operations: (1) swap the last two columns of Φ and then negate either of these columns. (2) multiply one of the last two columns by a constant α and the other by 1/α. (3) add any constant multiple of one of the last two columns to a different column. These rules are useful if you compute a Stackel matrix and it does not appear to agree with a Stackel matrix appearing in the literature. A Viability Condition for Simple Separability The functions shown on the left of these equations (which are determined from the choice of curvilinear system as outlined above) must have the functional form shown on the right. (H/h12) /(f1f2f3) = E(2)I(3)-F(2)H(3) - (H/h22) /(f1f2f3) = B(1)I(3)-C(1)H(3) (H/h32) /(f1f2f3) = B(1)F(2)-C(1)E(2) (5.9) Not only does each line declare a required functional form, but there is correlation between the three lines since each function appears on two different lines. Since the 11 classical curvilinear systems are simple-separable for the Helmholtz equation, we expect that each of them will meet the above viability condition. How do we find the first column of the Stackel matrix Φ ? In terms of our simplified notation matrix above Φ = S = det(Φ) We can write S = a(1) M1(23) + d(2) M2(31) + g(3) M3(12) and we can insert the Mn values from (5.8) to get S = a(1) (H/h12) /(f1f2f3) + d(2) (H/h22) /(f1f2f3) + g(3) (H/h32) /(f1f2f3) But the Robertson condition (4.16) says H / (f1f2f3) = S, so the above becomes 1 = a(1) (1/h1)2 + d(2) (1/h2)2 + g(3) (1/h3)2 (5.11) which we have to regard as a non-trivial functional form requirement. So this is another condition like (5.9) on the viability of a curvilinear coordinate system as a candidate for simple separation. We can examine now two systems: Toroidal: We already learned that toroidal coordinates are not R-separable , but they provide a good illustration of the notion of trying to solve a "functional form equation". We have from earlier h1 = h2 = aR2 h3 = sh(ξ1) a R2 R-2 ≡ [ch(ξ1)-cos(ξ2)] so (4.20) becomes 1 = a(1) a-2R-4 + d(2) a-2R-4 + g(3) a-2R-4 /sh2(ξ1) a2 = R-4 [{a(1) + d(2)} + g(3) /sh2(ξ1) ] a2/ [ch(ξ1)-cos(ξ2)]2 = [{a(1) + d(2)} + g(3) /sh2(ξ1) ] Since 1/(A-B)2 has no partial fraction expansion, the RHS can only have a single term and it must be a function of both 1 and 2 (that term must be, in fact, the left hand side). But none of the RHS terms shown even allows a function of 1 and 2, so Toroidal coordinates fails this condition. This gives us another reason why this system is not R-separable. Spherical: 1,2,3 = r,θ,φ h12 = 1 h22= ξ12 h32 = ξ12sin2ξ2 1 = a(1) (1/h1)2 + d(2) (1/h2)2 + g(3) (1/h3)2 (5.11) 1 = a(1) + d(2) /ξ12 + g(3)/ (ξ12sin2ξ2) An obvious solution is this a(1) = 1 d(2) = 0 g(3) = 0 which agrees with the first column of the Stackel matrix shown in M&S page 24. In general, we can see that if any of the three scale factors hi is a constant, we can get a simple first column solution. For example, if h1 = K, then we set a(1) = 1/K and the other two functions to 0, just as we did for Sphericals above. All cylindrical systems have g3 = 1 so they are all OK on (5.11). Not all 11 classical systems have a constant hi, however, so they would have to be examined. Summary of Simple Separation We attempt to solve the Helmholtz equation (2 + K12)ψ = 0 assuming ψ(123) = X1(1)X2(2)X3(3). We find that the separated ODE's for the Xn must have the following form, where k22 and k32 are arbitrary real parameters (separation constants) LnXn = (1/fn)∂n[fn(∂nXn)] + [ K12Φn1(n) + k22Φn2(n) + k32Φn3(n)]Xn = 0 (5.2) We have shown simple separation if we can find the 9 functions Φnm(n) of a Stackel matrix Φ as well as certain helper functions fn and gn such that the following five requirements are satisfied. Once the first requirement below is satisfied, we can regard it as determining the gn and we can then exclude gn from the list of "data items" we need to enumerate for a simple separating system. The first requirement is that the three equations (3.5) must be satisfied with R=1: (H/hn2) = fn(n)gn(≠n) n = 1,2,3 (3.5) with R=1 where H = h1h2h3. This represents three functional-form requirements directly on the set of scale factors. As we showed above equation (3.6), this requirement is not met in toroidal coordinates, for example. Equation (3.5) was assumed "in order to have a chance" of doing separation at all. The Cofactor Condition (5.6) says Mn(Φ)/S(Φ) = (1/hn2) and is a condition on the Stackel matrix Φ. It arises from our assumed form for the separated differential operators Ln shown above, and from the processed form of the original Helmholtz equation. The Robertson Condition (5.7), which says S(Φ) = H / (f1f2f3), is an imposed condition on Φ (not a condition on the curvilinear system parameters). It is fair to ask whether, if some other condition were imposed, could some formerly excluded coordinate system be made to join the list of simple-separable systems? We can regard this as part of a larger question which is: could new simple-separable coordinate systems be found by some different method of processing than that used in this write-up? We leave these questions to the reader. No doubt they have been addressed pretty thoroughly in the last 150 years. A fourth condition is that stated in (5.9), and it is another functional form condition on the functions fn and hn which we deduce from the coordinate system. A fifth and final condition is the functional-form condition (5.11) which is given directly in terms of the scale factors and relates to finding the first column of the Stackel matrix Φ. Here is a statement of all 5 of these conditions: (H/hn2) = fn(n)gn(≠n) is solvable for the fn and gn " fg conditions" (3.5) with R=1 Mn(Φ)/S(Φ) = (1/hn2) " Cofactor conditions" on Φ (5.6) S(Φ) = H / (f1f2f3) "Robertson condition" on Φ (5.7) (H/h12) /(f1f2f3) = E(2)I(3)-F(2)H(3) "Column 23 Conditions" - (H/h22) /(f1f2f3) = B(1)I(3)-C(1)H(3) are solvable for all 6 cap functions shown (H/h32) /(f1f2f3) = B(1)F(2)-C(1)E(2) (5.9) 1 = Φ11(1) (1/h1)2 + Φ21(2) (1/h2)2 + Φ31(3) (1/h3)2 is solvable for the 3 Φn1(n) (5.11) "Column 1 Condition" The Cofactor and Roberson conditions do not impose any restrictions on the curvilinear coordinate scale factors hi, whereas the other three conditions do. They must all be checked. 6. R-separation We shall use the exact same form for operator Ln shown in (4.8), and we are facing our same processed Helmholtz equation (4.9). We rewrite both those equations here LnXn = (1/fn)∂n[fn(∂nXn)] + [ κ12Φn1(n) + k22Φn2(n) + k32Φn3(n)]Xn = 0 (6.1) Σn(Q/hn2) [ κ12Φn1(n) + k22Φn2(n) + k32Φn3(n)] + k12 + QK12 = 0 (6.2) For simple separation we had k12 = 0, cancelled the Q's, and proceeded from there setting κ12 = K12. For R-separation we want to do the same matrix inversion we did before, but this time with Q as part of the summand, and we set κ12 = k12. This method cannot support K12 ≠ 0, so we are forced at this point to give up on the Helmholtz equation and try to be successful with the Laplace equation (K12 = 0). We then have from (6.2), Σn(Q/hn2) [ k12Φn1(n) + k22Φn2(n) + k32Φn3(n)] + k12 = 0 (6.3) We now proceed through the same steps as in simple separation but with Q sitting in the sum. We set the coefficients of the kn2 to 0 in (6.3) to get three equations Σn (Q /hn2) Φn1(n) = 1 Σn (Q /hn2) Φn2(n) = 0 Σn(Q /hn2) Φn3(n) = 0 which we write as a matrix equation exactly as in (5.3) but with Q sitting over the hn2. We invert the matrix equation in the same manner and write out the inverted equation to get a new version of our Cofactor Condition, analogous to (5.6) but with this extra Q, Mn/S = (Q/hn2) " Cofactor condition" (6.4) For our Robertson Condition we make a slightly different ansatz this time which is S(Φ) = { H / (f1f2f3)} [1/(QR2)] "Robertson condition" (6.5) where the simple-separation Robertson condition was just the {...} part. As noted in the simple separation section, if we can find a solution for Φ with this assumed condition, then it is justified de facto. Remember that the Robertson condition is just a condition on matrix Φ, not on any of the underlying curvilinear scale factors. What now are the expressions for Mn analogous to (5.8) ? From two equations above we get Mn = SQ/hn2 = [ H / (f1f2f3)] [1/(QR2)] (Q/hn2) = { (H/hn2) /(f1f2f3) } (1/R2) (6.6) where the object in {...} was the simple separation result for Mn. Since we know all the objects on the RHS of (6.6), we know the three Mn cofactors. Using (3.5) which says (H/hn2) = R2 fn(n)gn we get this alternate form Mn = { fn(n)gn /(f1f2f3) } (6.6a) which we might better express as M1 = g1 /(f2f3) M2 = g2 /(f3f1) M3 = g3 /(f1f2) (6.6b) Next, to get a new version of our "first column equation", we start with S = a(1) M1(23) + d(2) M2(31) + g(3) M3(12) and we can insert the Mn values from (6.6) to get R2S = a(1) { (H/h12) /(f1f2f3) } + d(2) { (H/h22) /(f1f2f3) } + g(3) { (H/h32) /(f1f2f3) } But the Robertson condition (6.5) says QR2S = { H / (f1f2f3)}, so the above becomes R2S = a(1) (1/h12) QR2S + d(2) (1/h22)QR2S + g(3) (1/h32)QR2S 1/Q = a(1) (1/h12) + d(2) (1/h22) + g(3) (1/h32) (6.7) Question: For R separation, what are the corresponding "Column 23 Conditions" ? They were given in (5.9) for simple separation. We still have this underlying problem to solve in R separation M1(23) = E(2)I(3)-F(2)H(3) -M2(31) = B(1)I(3)-C(1)H(3) M3(12) = B(1)F(2)-C(1)E(2) which is a functional form requirement on the Mn . Installing from (6.6b) above we get g1 /(f2f3) = E(2)I(3)-F(2)H(3) - g2 /(f3f1) = B(1)I(3)-C(1)H(3) g3 /(f1f2) = B(1)F(2)-C(1)E(2) (6.8) So this is then our "Column 23 Condition". Calculation of Φ for Toroidal Coordinates We want of course to try this out for our ongoing Toroidal coordinates example, where we had h1 = h2 = a/[ch(ξ1)-cos(ξ2)] => h1 = h2 = aR2 H = sh(ξ1)a3R6 h3 = a sh(ξ1)/[ch(ξ1)-cos(ξ2)] => h3 = a sh(ξ1)R2 f1(1) = sh(ξ1) g1(23) = a R = [ch(ξ1)-cos(ξ2)]-1/2 f2(2) = 1 g2(31) = a sh(ξ1) Q = a2[ch(ξ1)-cos(ξ2)]2 = a2R4 f3(3) = a g3(12) = 1/ sh(ξ1) k12 = (1/4) We first compute S from the new Roberson condition (6.5) S(Φ) = { H / (f1f2f3)} [1/(QR2)] = { sh(ξ1)a3R6/ (ashξ1) }[1/(QR2)] = a2R4/Q = 1 Next the Mn from the Cofactor condition (6.4) M1 = S(Q/h12) = [a2R4/Q] (Q/h12) = a2R4/h12 = a2R4/(a2R4) = 1 M2 = S(Q/h22) = [a2R4/Q] (Q/h22) = a2R4/h22 = a2R4/(a2R4) =1 M3 = S(Q/h32) = [a2R4/Q] (Q/h32) = a2R4/h32 = a2R4/(a2R4sh2(ξ1)) = 1/(sh2(ξ1) We now try to obtain the last two columns of Φ Φ = S = det(Φ) M1(23) = e(2)i(3)-f(2)h(3) = 1 -M2(31) = b(1)i(3)-c(1)h(3) = -1 M3(12) = b(1)f(2)-c(1)e(2) = 1/(sh2(ξ1) Start with b(1) = 1/[ sh2(ξ1)] and f(2) = 1 in line 3 , then need i(3) = 0 in line 2, then need e(2) = 0 in line 3 giving at this juncture M1(23) = -h(3) = 1 -M2(31) = -c(1)h(3) = -1 M3(12) = b(1) = 1/(sh2(ξ1) The first line says h(3) = -1, so then the second line tells us c(1) = -1. Results for our last two Stackel matrix columns are: b(1) = 1/(sh2(ξ1) c(1) = -1 e(2) = 0 f(2) = 1 h(3) = - 1 i(3) = 0 Next we attempt to solve (6.7) for the first column: S = a(1) M1(23) + d(2) M2(31) + g(3) M3(12) 1 = a(1) 1 + d(2) (-1) + g(3) 1/(sh2(ξ1) and one solution to this is g(3) = 0, d(2)=0, a(1) = 0, so we now have a(1) = 1 b(1) = 1/(sh2(ξ1) c(1) = -1 d(2) = 0 e(2) = 0 f(2) = 1 g(3) = 0 h(3) = - 1 i(3) = 0 Φ = = Now to get an equivalent Φ, we swap the last two columns and then negate the third column Φeq = = and this is the Stackel matrix which appears in Moon & Spencer page 112. Aside: we could to the first row stuff using (6.7) 1/Q = a(1) (1/h12) + d(2) (1/h22) + g(3) (1/h32) (6.7) h12/Q = a(1) + d(2) + g(3) (h12/h32) a2R4/ (a2R4) = a(1) + d(2) + g(3) 1/(sh2(ξ1) 1 = a(1) + d(2) + g(3) 1/(sh2(ξ1) and I must have (6.7) still wrong after 4 repairs to it!!! Appendix A: A small digression on inverting a matrix From linear algebra, if we have a matrix equation in En X = AY we know that if detA ≠0 we can invert to get Y = A-1Y where [A-1]pq = cof(Aqp)/det(A) Notice that the subscripts are switched on the two sides, and cof(Aqp) = (-1)p+q minor(Aqp) where minor(Aqp) is the determinant of the (n-1)x(n-1) matrix obtained by crossing out row q and column p of the matrix A. Now suppose we are given this problem instead X = BTY where T means transpose. If det(BT) ≠ 0, we can invert to get Y = (BT)-1X where, noting that det(BT) = det(B), (BT)-1pq = cof(BTqp)/det(BT) or (BT)-1pq = cof(Bpq)/det(B) where now the indices are in the same order on both sides. To clarify notation, we might then define a matrix C where Cpq ≡ cof(Bpq) " the cofactor matrix" so we then have (BT)-1pq = [1/det(B)]Cpq so that our inversion is given by Y = [1/det(B)] C X inversion of X = BTY