"Science, Faculty of"@en .
"Earth, Ocean and Atmospheric Sciences, Department of"@en .
"DSpace"@en .
"UBCV"@en .
"Johnson, Ian Mayhew"@en .
"2011-03-21T19:25:57Z"@en .
"1972"@en .
"Doctor of Philosophy - PhD"@en .
"University of British Columbia"@en .
"The inversion problem of the geomagnetic secular variation, that is the recovery of the electrical conductivity distribution of the earth's lower mantle from the attenuation of the magnetic field across the lower mantle, is shown to have a unique solution. This implies that the electrical conductivity is uniquely determined by the spatial and time history of the magnetic field at the base of the mantle. The derivation of the inversion method is based on the inverse Sturm-Liouville theory of I. M. Gel'fand and B. M. Levitan.\r\nNumerical experiments on the theory with synthetic electrical conductivity profiles indicate a computationally stable inversion method. The resolving power of the inversion is shown to have a predictable dependence on the extent of the attenuation information. The attenuation properties required for the inversion are shown to be recoverable from the spatial deconvolution of two autocorrelation functions which are expressions of the statistical characteristics of the magnetic field components at the boundary of the core and the mantle."@en .
"https://circle.library.ubc.ca/rest/handle/2429/32653?expand=metadata"@en .
"INVERSION OF THE GEOMAGNETIC SECULAR VARIATION UNIQUENESS AND FEASIBILITY by Ian Mayhew Johnson B.Sc, University of Western Ontario, 1968 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Department of GEOPHYSICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July, 1972 In p r e s e n t i n g t h i s t h e s i s in p a r t i a l f u l f i l m e n t o f the requirements f o r an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia , I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e fo r reference and s tudy . I f u r t h e r agree t h a t p e r m i s s i o n fo r e x t e n s i v e copy ing o f t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the Head o f my Department o r by h i s r e p r e s e n t a t i v e s . I t i s understood that copying or p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l ga in s h a l l not be a l lowed wi thout my w r i t t e n p e r m i s s i o n . Department of GEOPHYSICS AND ASTRONOMY The U n i v e r s i t y of B r i t i s h Columbia Vancouver 8, Canada Date AUGUST 14, 1972 (i) ABSTRACT The inversion problem of the geomagnetic secular variation, that i s the recovery of the e l e c t r i c a l conductivity distribution of the earth's lower mantle from the attenuation of the magnetic f i e l d across the lower mantle, is shown to have a unique solution. This implies that the e l e c t r i c a l conductivity is uniquely determined by the spatial and time history of the magnetic f i e l d at the base of the mantle. The derivation of the inversion method is based on the inverse Sturm-Liouville theory of I. M. Gel'fand and B. M. Levitan. Numerical experiments on the theory with synthetic e l e c t r i c a l conductivity profiles indicate a computationally stable inversion method. The resolving power of the inversion i s shown to have a predictable dependence on the extent of the attenuation information. The attenuation properties required for the inversion are shown to be recoverable from the spatial deconvolution of two autocorrelation functions which are expressions of the s t a t i s t i c a l characteristics of the magnetic f i e l d components at the boundary of the core and the mantle. ( i i ) TABLE OF CONTENTS ABSTRACT (i) LIST OF FIGURES ( i i i ) ACKNOWLEDGEMENTS (iv) CHAPTER I INTRODUCTION 1.1 The Problem 1 1.2 Review of Related Work . 3 1.3 Scope of this Thesis 6 CHAPTER II THEORY 2.1 Theoretical Background 8 2.2 Attenuation Properties 12 2.3 Review of Inverse Sturm-Liouville Theory 21 2.4 The Inversion 23 2.5 Generalized Transfer Functions and the Inversion 36 CHAPTER III NUMERICAL TESTS OF THE INVERSION 3.1 General Aims 40 3.2 Transfer Functions to Spectral Function 41 3.3 Spectral Function to Conductivity Distribution 48 CHAPTER IV SUGGESTIONS FOR THE APPLICATION 82 CHAPTER V CONCLUSIONS ' 86 APPENDIX A-l SOLUTION OF THE INTEGRAL EQUATION 88 APPENDIX A-2 SPATIAL CONVOLUTION OF P0L0IDAL FIELD COMPONENTS 92 APPENDIX A-3 PROGRAM LISTING: GENERATION OF A SPECTRAL FUNCTION 97 APPENDIX A-4 PROGRAM LISTING: THE INVERSION ,103 REFERENCES 118 ( i i i ) LIST OF FIGURES Figure Page Power transfer function curves of different harmonic degrees for 1 43 different lower mantle ele c t r i c a l conductivity models. Results of the numerical inversions 2 51 for different lower mantle 3 53 conductivity models, different 4 54 harmonic degrees, and different 5 55 spectral function types. 6 58 7 59 8 61 9 63 10 65 11 66 12 67 13 69 14 71 15 72 16 73 17 74 18 75 19 76 Average profile from the numerical inversion with the two spectral functions for the same conductivity model. (iv) ACKNOWLEDGEMENTS By f a r the greatest part of my appreciation goes to Dr. D. E. Smylie. For him, I reserve a s p e c i a l thankyou f o r h i s guidance and encouragement. I would l i k e to thank Dr. G. W. Bluman of the Department of Mathematics f o r c r i t i c i s m s and suggestions. I am also g r a t e f u l to the Department of Geophysics under Dr. R. D. Russell and to the University of B r i t i s h Columbia for support during t h i s research. A s p e c i a l note of thanks goes to Ms. Rosanne Rumley f o r her very able help i n the typing of the thesis. CHAPTER I INTRODUCTION 1.1 The Problem The magnetic f i e l d at the surface of the earth i s observed to be continually varying with time and with location, and is known to be the result of the superposition of the effects of electromagnetic f i e l d sources which are both above and below the earth's surface. That part of the geomagnetic f i e l d which is of internal origin exhibits a variation with time which is slow, that i s , changes commonly have periods of years. This behavior of the internal f i e l d is called the geomagnetic secular variation. Contemporary geomagnetic theory suggests that the secular variation is generated within the earth's liquid core, the radius of the core being just over half that of the earth. The region between the core and the surface is divided into mantle and crust, the crust being a rela-tively thin shell near the surface. The mantle therefore represents the largest part of the earth by volume; a part about which l i t t l e is known and much is speculated concerning the composition and the el e c t r i c a l properties. One of the few methods available to increase our knowledge of the ele c t r i c a l conductivity of the mantle i s based on the study of how the internal and external parts of the geomagnetic f i e l d are each affected by the mantle's conductivity structure. Of particular importance i n this approach is the role of inversion theorems. In the case of the secular variation, a time and space varying -2-geomagnetic f i e l d generated in the earth's core passes through the mantle and crust to the surface. The attenuation i t suffers in passing through the mantle and the crust is determined by simple geometrical spreading and the so-called physical attenuation determined entirely by the elec-t r i c a l conductivity structure. The inversion problem of the geomagnetic secular variation may be stated as follows: is the el e c t r i c a l conductiv-it y distribution of the mantle uniquely defined by the attenuation of the magnetic f i e l d in diffusing through the mantle and crust? And i f so, how is i t recovered from the attenuation properties? The lack of direct information on the magnetic f i e l d at the boundary of the core and the mantle, and therefore on the attenuation properties of the earth from the core to the surface, is an obvious constraint on the usefulness of such a theory. There i s , however, a large body of theoretical understanding of the fl u i d motions in the core which would make the solution to the inversion problem less academic than might be expected. Fundamental to any solution of this problem i s i t s directness. Given the attenuation properties, an inverse theory i s desired which shows that the conductivity profile can be recovered uniquely and also, hopefully, an algorithm which shows promise of direct numerical inversion poss i b i l i t i e s with real data. Models and iterative approaches to non-unique solutions are to be avoided. This i s a particularly important distinction to make as others maintain a different definition for geomag-netic inversion. References to their work, which does not overlap the work presented here, may be found i n the next section. The problem then is to build a workable inverse theory for the -3-secular variation and to see how i t may be applied to improve our under-standing of the ele c t r i c a l properties of the earth's mantle. 1.2 Review of Related Work A review of the important contributions to geomagnetic research should start with the \"General Theory of Terrestrial Magnetism\" by Carl Friedrich Gauss (1841). With potential theory and spherical harmonic analysis, he gave to geomagnetism the firm mathematical footing which i s the foundation of any attack on our ignorance of the e l e c t r i c a l conductivity structure of the mantle. The study of geomagnetic variations, with the purpose of determining the ele c t r i c a l conductivity of the mantle, has traditionally taken one of two possible approaches. The transmission problem of the diffusion of the secular variation through the mantle is the subject to which this thesis contributes. There i s also the reflection problem from the external component of the magnetic f i e l d , penetrating the upper part of the mantle, and the earth's response measured at the surface. The relations between input and response fields determines indirectly the e l e c t r i c a l conductivity of the upper part of the mantle. The reflection problem is also an inversion problem in the st r i c t sense. It has only recently been solved by Bailey (1970). His results show that the inversion i s unique: that i s , given the external f i e l d variations and the earth's response, the ele c t r i c a l conductivity distribution of the mantle, assumed to be spherically symmetric, may be recovered uniquely. The applicability of the theory appears, however, to be jeopardized by an overly sensitive dependence- on the completeness -4-of the information about the geomagnetic f i e l d through space and time at the surface. The most complete analysis to date of the external f i e l d v a r i a t i o n s f o r the determination of the e l e c t r i c a l c onductivity has been done by Banks (1969). He employed the t r a d i t i o n a l method, developed by S. Chapman and used by L a h i r i and P r i c e (1939), of suggesting a paramet-r i z e d c o n d u c t i v i t y model f o r the mantle of a s p h e r i c a l l y symmetric earth and comparing the t h e o r e t i c a l l y predicted response of the earth with that a c t u a l l y observed. Adjusting the model such that the observed and pred i c t e d values agree, gives the parameters of the model. His r e s u l t s confirmed a jump i n the e l e c t r i c a l c onductivity of approximately two orders of magnitude, to a value of 1 (ohm m.)\" 1 , or mhos per meter, over a depth of 400 to 600 kms. The depth at which the conductivity e x h i b i t s such a large increase i s conventionally used as the boundary between the upper and the lower mantle. The usefulness of sounding the conductivity s t r u c t u r e with the external component of the geomagnetic f i e l d i s l i m i t e d by the very poor r e s o l u t i o n of the e l e c t r i c a l c onductivity with the method i n the lower mantle. Most of the estimates of the conductivity i n t h i s region have been a r r i v e d at by applying the modelling approach to the d i f f u s i o n of the secular v a r i a t i o n . The r e s u l t s s u f f e r from the use of unsophisticated models of the magnetic f i e l d at the core-mantle boundary. A comprehensive i n v e r s i o n theory was unavailable, making d e t a i l e d t h e o r e t i c a l estimates unwarranted. E l s a s s e r (1950), Runcorn (1955), McDonald (1957), Smylie (1965), and Currie (1968) used simple conductivity models to make e s t i -mates of the e l e c t r i c a l c onductivity over the lower mantle. The most -5-extensive analysis is that of McDonald, who used a two-parameter conduc-t i v i t y model. A l l of these results suffer from large error bounds, although an average el e c t r i c a l conductivity for the lower mantle of something like 100 mhos per meter seems to have found general agreement. An inverse theory for the secular variation would eliminate some of the restrictions and ambiguities inherent i n the modelling approach, an approach which has been responsible for most of our ideas of lower mantle conductivity structure to date. Any review of inverse theories and work done to determine the el e c t r i c a l conductivity of the mantle would not be complete without some mention of the so-called inversion method of Backus and Gilbert (1970) and the employment of that method by Parker (1970). Despite the t i t l e , \"Uniqueness and Inversion of Inaccurate Gross Earth Data\" of the work by Backus and Gilbert, the result seems to be the addition of some logic and order into general geophysical modelling procedures but not, in fact, a direct inversion technique. Parker has reworked Banks' data of the external f i e l d variations for some interesting results on the accuracy-resolution trade-off of the recovered conductivity pr o f i l e , following Backus and Gilbert, but the result remains essentially unchanged from that of Banks. The most recent work on this 'inversion' method and that which gives a complete l i s t of references for this expanding area of geophysical interest i s by Wiggins (1972). Hence, that area of geophysics, which deals with the determina-tion of the el e c t r i c a l conductivity of the earth's lower mantle, is much in need of a useable inverse theory of the geomagnetic secular variation. The need i s both for completeness of our theoretical understanding of -6-the problem and f o r improved estimates of the lower mantle conductivity-s t r u c t u r e . 1.3 Scope of This Thesis The problem to be in v e s t i g a t e d i s , again: The magnetic i n d u c t i o n , B , generated i n the earth's core, passes through the mantle and appears at the surface of the earth i n some attenuated form. Given a d e s c r i p t i o n of the attenuation of the i n t e r n a l magnetic f i e l d , i s i t p o s s i b l e to recover the e l e c t r i c a l conduc-t i v i t y of the lower mantle uniquely? And i f so, how i s the i n v e r s i o n c a r r i e d out and how much information i s required? The b a s i s of the s o l u t i o n to t h i s problem i s an inverse theory which guarantees that the conductivity s t r u c t u r e , which i s assumed to be s p h e r i c a l l y symmetric, i s uniquely defined by the attenuation properties of the mantle. In f a c t , the r a t i o of two t r a n s f e r functions of any harmonic degree i s a l l that i s required f o r a unique recovery. The two t r a n s f e r functions describe the attenuation i n the harmonic c o e f f i c i e n t s of the r a d i a l and transverse parts of the p o l o i d a l f i e l d across the earth's lower mantle. Extensive computer t e s t i n g with model conductivity d i s t r i b u t i o n s has been c a r r i e d out with the i n v e r s i o n theory. The r e s u l t s , presented i n chapter I I I , i n d i c a t e an algorithm with enough numerical s t a b i l i t y to suggest that the method w i l l give u s e f u l c o n d u c t i v i t y information with the eventual a p p l i c a t i o n to r e a l data. A v a i l a b l e data would of course s u f f e r from a lack of complete s p a t i a l and time information at the surface and estimates of the magnetic -7-f i e l d at the core-mantle boundary are nec e s s a r i l y speculative. With a view to the a p p l i c a t i o n of the method, suggestions are included on some of the mathematical development and numerical problems which may a r i s e when the in v e r s i o n of the geomagnetic secular v a r i a t i o n i s c a r r i e d out with r e a l data. CHAPTER' II THEORY 2.1 Theoretical Background The mathematical expression of some basic geomagnetic theory i s reviewed here and the notation used In the thesis i s introduced. The s p i r i t and style of the theoretical concepts are those of Smylie (1965), from which only those parts necessary to a theoretical background to the inversion are extracted. The basic assumption, common to most geomagnetic theory of the secular variation, that the el e c t r i c a l conductivity of the mantle i s a real positive function of radius only, i s also made here. The magnetic induction, B , in the earth's mantle and crust, where relative motions of material are assumed to have negligible electro-magnetic affects, may be expressed as B = S\" + T (1) where S is called the poloidal part and T the toroidal part of the induction B. They are defined as S = V x (VxrS) T = VxrT (2) where r is a radius vector and S and T are scalars, determined to within an additive constant by the induction B. The relatively low e l e c t r i c a l conductivity^ of the upper mantle assures that the toroidal f i e l d , which has no radial component, forms a negligible part of the f i e l d at the surface. The secular variation observed -9-at the surface ts therefore almost entirely poloidal and the attenuation properties w i l l describe the attenuation of the poloidal part of the total f i e l d B. The spherical symmetry of the mantle suggests the spherical harmonic expansion S = I I s^(r,t)P* 0 , n 0 n n+m \u00E2\u0080\u0094 2 n! du and for a l l m by p - m = m i n ^ l l pm n (n+m) i n The condition that the induction, B , i s real implies that the complex coef f ic ients , s \u00E2\u0084\u00A2 ( r , t ) , obey the relat ion * \u00E2\u0080\u00A2 ' -m m (n+m)! m*> . s n (r . t ) - (-1) s n (r ,t) where the asterisk denotes the complex conjugate. When the spherical harmonic expansion for the scalar S (equation (3)) i s substituted into the vector d i f fe rent ia l form for the poloidal part of the induction (equations (2)), the poloidal vector has components, each of which i s a double sum over i t s respective vector elements. The polo idal vector elements -10-are gtyen by n r n n CS\"1) - i f (rs111) % e n 6 r 3r n -TT\u00E2\u0080\u0094 do dP m im /rnu im 3 , mN _,m im** - n(n+l) s m - . 3 m = y a(r) \u00E2\u0080\u0094 s o 3t n ( 5 ) where P q i s the permeability of free space and a(r) i s the e l e c t r i c a l conductivity. The corresponding dif f e r e n t i a l equation for the scalar function S i s V2S = u 0(r) |- S 0 d t The l a s t point to be made i n this section concerns an assumption about the e l e c t r i c a l conductivity of the upper mantle, and the boundary m, conditions on s (r,t) which result from this assumption. F i r s t some n - l i -no tat ion: let b represent the radius of the core-mantle boundary, c the radius of the boundary between the lower and upper mantle, and d the radius of the earth. The dist inct ion between, the upper and lower mantle is made because the region of the earth between the radius r=c and the surface i s thought to have an average e lec t r i ca l conductivity which is very much smaller than that of the lower mantle. The div is ion is not c r i t i c a l for the inverse theory. It does provide a convenient and physically reasonable point at which to apply the condition that the e lec t r i ca l conductivity becomes zero. Between r=c and r=d , under the assumption of zero e l e c t r i c a l conductivity, the scalar S sat is f ies Laplace's equation. That is, the magnetic f i e l d becomes harmonic, the poloidal f i e l d coeff icients sat isfy ing the equation 2 1 3 , nu n(n+l) m . 7 ~ 2 - ~ T - S n = \u00C2\u00B0 ' r \u00C2\u00B1 - C -3r r Subject to the condition that the internal f i e l d should weaken with increasing radius, the solution in the upper mantle i s m / /d.n+1 m., . s n ( r , t ) = (-) s n (d , t ) The continuity of the magnetic induction across the boundary between the lower and upper mantle, that i s across the radius c, implies that m / . \ /d,n+l m, ,. N s n ( c , t ) = (-) s n (d , t ) 3 . m. ,dvn+l m,, . \u00E2\u0080\u0094 (rs ) = -n (\u00E2\u0080\u0094) s (d,t) 3r n r=c c n (7) -12-The theory of the diffusion of the poloidal f i e l d from the core-mantle boundary to the surface has been mathematically cast into the form of a part i a l d i f f e r e n t i a l equation describing the diffusion of the poloidal harmonic coefficients,(equation (5)),and the boundary conditions on these coefficients applied at r=c (equations (7)). The next task is to define what is meant by attenuation properties and to relate the above theoretical basis to functions which describe the attenuation of the poloidal f i e l d over the lower mantle. 2.2 Attenuation Properties In introducing the inversion problem, terminology was used, (input, output, transfer function) which i s particular to a way of looking at many similar problems. The diffusion of a given poloidal f i e l d coefficient of a particular spherical harmonic can be thought of as a distributed linear system (Schwarz and Friedland, 1965). Distributed means a system for which the spatial configuration i s important and the dynamics are mathematically represented by a parti a l differential equation. The diffusion of the radial component of a given harmonic degree n and order m, for example, across the lower mantle, the 'system' for the inversion problem, i s represented by the pa r t i a l differential equation (5). For this system, the change i n , or the attenuation of any radial poloidal coefficient, s\u00E2\u0084\u00A2(r,t) , from the core-mantle boundary to the surface i s completely described by either an impulse response function in the time domain, or by a transfer function in the frequency domain. These two functions make up a Fourier transform pair. In the time domain, the radial ouput coefficients of the poloidal f i e l d at the surface, may each be expressed as a convolution of the input radial coefficients, s m(b,t) , with a -13-corresponding impulse response function. In the frequency domain, the system transfer function i s the ratio of the Fourier transforms of the output and the input coefficients of the same harmonic degree and order. Attenuation properties may similarly be assigned to the transverse part of the poloidal f i e l d . It would seem that transfer functions, or impulse response functions, for a l l possible combinations of the harmonic degree and order are needed to describe completely the attenuation of any component, r, 8, or of the poloidal f i e l d . The theory w i l l reveal, however, that the attenuation properties, either transfer functions or impulse response functions, are independent of the harmonic order m. The physical attenuation, as opposed to the simple geometrical attenuation, of a l l f i e l d components is determined by the physical parameters of the system, in this case the e l e c t r i c a l conductivity structure of the lower mantle. The poloidal f i e l d , as seen in the vector element forms (equations (4)), i s composed of two independent parts with respect to the harmonic coefficients. The radial part of the f i e l d i s determined by the coefficients s m ( r , t ) . The transverse part of the f i e l d i s determined by the coefficients 3/8r(rs\u00E2\u0084\u00A2). It i s necessary, therefore, to specify two transfer functions for each harmonic degree to describe the attenuation of a l l vector components of the poloidal f i e l d The necessity of two transfer functions, or two impulse response functions for each harmonic degree n, implies the presence of two systems. The attenuation, however, of a l l coefficients of a given harmonic degree of a l l the vector components of the poloidal f i e l d i s completely determined by the e l e c t r i c a l conductivity structure of the mantle. The confusion i s a result of the fact that the two parts, radial and transverse, of the -14-poloidal f i e l d , and therefore the two transfer functions, are related through the differential equation (5). The need for two system transfer functions for one system i s a r t i f i c i a l as both transfer functions are determined by the e l e c t r i c a l conductivity. The inverse theory w i l l show that, in fact, the ratio of the two transfer functions alone i s required to recover the conductivity profile. That i s , the ratio of the transfer functions describing the attenuation of the radial and transverse poloidal f i e l d coefficients of some harmonic degree, i s required to define the e l e c t r i c a l conductivity uniquely. To get to that point, i t i s necessary to relate the poloidal f i e l d coefficients and the attenuation properties. The derivation of the relationship between the coefficients s\u00E2\u0084\u00A2(r,t) and 3/3r(rs\u00E2\u0084\u00A2) , and the corresponding transfer functions, begins with the transform into the frequency domain, of the differential equation (5). The choice of the type of integral transform, either Fourier or Laplace, to be used, i s guided by the understanding that the i n i t i a l f i e l d distribution in the mantle should be unimportant today. The Fourier transform of the function s m ( r , t ) is defined here as n Z n(r,iu) - ^ ^ ' ^ e d t \u00E2\u0080\u0094oo This form i s used, contrary to the conventions of time-series analysis, as i t simplifies the notation of the inverse theory. When the Fourier transform i s applied to the parti a l differential equation (5) and the associated boundary conditions (equations (7)), and the resulting equations are divided by the normalization, c\u00C2\u00A3 m(c,iw), n the ordinary differential equation and boundary conditions -15-2 \u00E2\u0080\u0094 . M (r,Iu) +\u00E2\u0080\u00A2 { u aCr ) iw - n C n + 1 ) }M C r ' i u ) = 0 , ^ n o ^ n dr r M nCc,iw) = 1 ' (8) \u00E2\u0080\u0094 - M (r,iw) = -n/c dr n r=c r e s u l t . The functions M^Cr.iu) are r e l a t e d to the Fourier transforms of the p o l o i d a l f i e l d c o e f f i c i e n t s by \u00C2\u00A3\"Cr,iw) = ( | ) n + 1 ( ^ ) M n(r,io))Z^(d,i U) . (9) The trans f e r function, which describes the attenuation of the r a d i a l p o l o i d a l harmonic c o e f f i c i e n t s of degree n and a r b i t r a r y order m, across the lower mantle, H n ( b , c o ) , i s given by H C b . u O , O f f i . . ( 1 0 ) n The t r a n s f e r function, G (b, oi) , which describes the attenuation of the n c o e f f i c i e n t s of the transverse p o l o i d a l f i e l d c o e f f i c i e n t s of degree n, i s d , \u00E2\u0080\u009Em. \u00E2\u0080\u0094 ( r Z ) dr n r=c . . G n ( b \u00C2\u00BB u ) = d m \" H \u00E2\u0080\u00A2 <1\u00C2\u00AB f - ( r Z m ) 4 M C r . i i o ) v d r v n r=b dr n r=b -16-t t t s reasonable to expect that the transfer functions w i l l be known for real frequencies only when real data is used. In anticipation of the requirements of the inverse theory, therefore, i t i s necessary to establish the relationship between the transfer functions, known along the real u> axis, and their representation i n the rest of the complex u plane. This i s done by expanding the transfer functions, which are meramorphic functions, as i n f i n i t e series. The representation of the transfer functions i n terms of their poles, requires some consideration of the nature of the zeros of their denominators. This leads to Sturm-Ltouville theory for which Ince (1956) i s used as the main reference. The substitution of X for ia> i n the ordinary d i f f e r e n t i a l system (equations (8)), and the addition of an arbitrary condition on the solutions at the core-mantle boundary, gives d 2 M (r,X) + {y a(r)X - n ( n t 1 ) } M (r,X) = 0 \u00E2\u0080\u0094 T n o I n dr M (c,X) = 1 n d M (r,X) n r=c dr -n (12) d M (r,X) 3'M (b,A) + 3 - r - 5 - ^ = 0, n dr a boundary value problem of Sturm-Liouville type where 3' and 3 are arbitrary constants. As the boundary condition at r=b i s homogeneous, 3' and 3 represent one arbitrary real constant. The number of boundary conditions on the system (12) dictates that only particular values of the parameter X are allowed -17-tn the aolutton, M QO:,X). The particular values of X form an i n f i n i t e sequence of real numbers, X^ ,X^ . Each member of this sequence i s called a characteristic frequency or characteristic value or eigenvalue. To the numbers, X , correspond the solutions M (r,X ), called the characteristic functions or eigenfunctions of the Sturm-Liouvllle problem (12) . For this particular problem, the zeros of f(X) = S'M^CbjX) + 3 4-M (r,X) , , are distinct. That i s , the poles i n the complex X plane dr n r=b of l/f(A) are a l l simple poles (see Ince, pg. 241). The characteristic frequencies are a l l positive i f B'/B < 0 (see Ince, pg. 235). The condition that the characteristic frequencies of the system are a l l greater than zero i s physically important as any values of these frequencies which are equal to or less than zero, make the impulse response function a constant or an increasing function of time. (The forms of the impulse response functions are given i n equations (15).) The general results on the distribution of the characteristic values of the Sturm-Liouville problem (see Ince, pg. 233) may also be used to relax the condition that B'/B < 0 . Positive characteristic values are insured when B'/B < n/b. From the accepted results of Sturm-Liouville theory, therefore, the zeros X , i=l,2,\"- , of B'M (b,X) + B 4~M (r,X) , , are simple n. J n dr n r=b J and are i n an ordered sequence up the positive real X axis (negative imaginary co axis), where B'/B i s an arbitrary constant, less than n/b. The cases which are most applicable to the inversion problem are 3 = 0, specifying the location of the poles of the transfer function H (b,co), n and B' =0, specifying the location of the poles of the transfer function G (b ,co). n These results for the nature of the characteristic values of -18-a Sturm-Liouville problem, permit the expansion of the transfer functions in terms of their poles, using the Mittag-Leffler theorem (Whittaker and Watson, 1927, pg. 134). If 1 , j=l,2,-*' , i s the sequence of poles of j H (b,w), that Is from the boundary condition M (b,X) = 0, and u ,.j=l,2,' i s the sequence of poles of Gn(b,w), that i s from the boundary condition \"f r^ n( r>^) r =b = 0, the system transfer functions may be expressed as H (b,w) = H (b,o) + I (rr- . n ' n ' . L, X -ico R R n. n. \u00E2\u0080\u00A21 1 j-1 n \u00C2\u00B0\u00C2\u00B0 \"n. G (b,u) = G (b,o) + I ( ?-n n . -, y -it J = 1 n. ) n. n. n. 3 where R are the residues of H (b,u) at ico=X , and Q are the n. n n. n. J 2 3 residues of G (b,a)) at iw=y . The residues have the form n n. 2 R = \u00E2\u0080\u0094 { d Mn C b' X )X=X T 1 , \" j C dx\" n j n _ -n , d2M (r,X) ,-1 Q = \u00E2\u0080\u0094 { n r=b > j drdX X=u n. J The differential equation in Mn(r,X) implies that the transfer functions should go to zero as \u00C2\u00B1a approaches i n f i n i t y along the real oi (imaginary X) axis. Hence the transfer functions take on the form H (b,to) n Gn(b,co) R y - S C X - iu>) J=l n. j - I n, J-1 ( yn.\" \u00C2\u00B1 ( 0 ) J (13) -19-Prom the system of di f f e r e n t i a l equations (12) , i t may be seen that the transfer functions have no zeros i n the f i n i t e part of the co plane. This fact, combined with the analytically known expressions for the zero frequency responses, HQ(b,o) and G n(b,o), gives the alternate expression for the transfer functions Hn(b,u>) = G (b,aj) = (b_) n + 1 c ( 1 - i ^ - ) ( 1 - i ^ - ) c (1-i-^-) (1-i-^-) n l n2 (14) The attenuation properties of the lower mantle, for a given harmonic degree, are therefore completely determined by the locations of the poles, i n the complex frequency plane, of the two appropriate transfer functions. The bias towards describing the attenuation of the poloidal f i e l d coefficients i n terms of frequency i s intended as most of the inversion theory rests on this approach. Although the time domain characteristics of the attenuation are not used directly, they are included here. The significance of the poles and residues of the transfer functions i n describing the attenuation of the poloidal f i e l d , i s perhaps better understood when the attenuation properties are i n the form of impulse response functions. For a given harmonic degree, the responses observed at the radius r=c of the earth, from an impulse **~~ for t > o, respectively. These expressions are derived from the inverse Fourier transforms of the respective transfer functions, H (b,to) and The ratio of the transfer functions i s a general function of frequency, dependent on the harmonic degree n, and given, from equations (14),by This ratio i s characterized, i n the complex frequency plane, by an alternating series of simple zeros and poles down the negative imaginary to axis. In anticipation of the requirements of the inverse theory, i t i s important to be able to determine the pole and zero locations of this ratio, knowing i t for real frequency only. Analytic continuation of the function H (b,to)/G (b,co) from the real co axis insures that the two n n characteristic value sequences, A , j=l,2,*\" , and u , j=l,2,**** , n. n. J j are, i n theory, uniquely recoverable. The inversion problem of the secular variation i s thus transformed into an inverse Sturm-Liouville problem. The relationship n G (b,co) . n (16) between the two transfer functions, radial and transverse, and the two corresponding characteristic value sequences of two Sturm-Liouville problems, is established. The connection i s made as inverse Sturm-Liouville theory requires two distinct characteristic value sequences i n order to recover the e l e c t r i c a l conductivity uniquely. For simplicity, the subscript n w i l l now be assumed and notation such as M (r,A ) w i l l become M(r,A.). n n. J J 2.3 Review of Inverse Sturm-Liouville Theory From the behavior of the eigenfunctions at r=b, in terms of the parameter A, that is from the transfer functions, the e l e c t r i c a l conductivity profile>is sought. This operation is the 'inverting' of the boundary value problem (equations (12)). The technique for inverting linear ordinary homogeneous second order boundary value problems used here, is from work by I. M. Gel'fand and B. M. Levitan (1955). They were instrumental in developing an inverse theory for the Sturm-Liouville problem. Modifications by Levitan (1968) and others, completed the inverse theory, showing that a unique inversion requires two characteristic value sequences, or two spectra, of the boundary value problem. Other geophysical applications of this type of inverse theory are found in seismology (Ware and Aki, 1969) and in geomagnetic depth sounding (Weidelt, 1970). To outline the results of the inversion theory of Gel'fand and Levitan, i t is convenient to use the modified form of the Sturm-Liouville problem -22-d^ 0(u,\u00C2\u00A3) +' U2 - q ( u ) } 0 - ( u , i O = 0 , du 2 (17) A .+ H 0 ( 1 , J I ) =0 , 0 < u < 1. du u=I - -They- established what i s necessary and how to recover the function q(u) uniquely. Given the normalization of the solutions, 0(o,&) for example, this form of the Sturm-Liouville problem (equations (17)), i s seen later to be equivalent to the differential system (12), after a simple variable change. The constants h and H are, as before, real. One of the intriguing properties of the eigenfunctions, 0(u,\u00C2\u00A3_.), is their orthogonality. They form a complete orthogonal set (see Ince, pg. 273). Hence the orthogonality property f 0(u,fc )0(u,A )du = p 6jj where 6^ i s the Kronecker delta, may be used in the expansions of arbitrary functions, which are at least piecewise continuous, i n a uniformly convergent series. This property of the eigenfunctions i s the same for cosines and sines of Fourier series, a special case, (q(u) = 0), of the Sturm-Liouville problem (17). The original work by Gel'fand and Levitan established that, given the normalizations or weights, p^ . , j=l,2,*\"' , the characteristic values \u00C2\u00A3j , j=l,2,\"** , and their asymptotic values, of the boundary value problem (17), the function q(u) may be recovered uniquely. The two sequences of numbers, p., j=l,2,\u00C2\u00AB--, and I., j=l,2,---, make up the -23-s p e c t r a l function of the system. The Gel'fand-Levitan i n v e r s i o n method, based on eigenfunction expansions, i s used to solve the e l e c t r i c a l conductivity i n v e r s i o n problem. Levitan (1968) and others showed that .the weights, p^ . , j=l,2,*-' , are determined by any two d i s t i n c t c h a r a c t e r i s t i c value sequences. The r e l a t i o n of inverse Sturm-Liouville theory to the d i f f u s i o n problem of the secular v a r i a t i o n i s now c l e a r . The two c h a r a c t e r i s t i c value sequences are known from the r a t i o of the frequency attenuation of the r a d i a l and transverse p o l o i d a l f i e l d c o e f f i c i e n t s of a given harmonic degree. The i n v e r s i o n method, applied to the recovery of the e l e c t r i c a l conductivity of the lower mantle from two c h a r a c t e r i s t i c value sequences, should complete the inverse theory f o r the geomagnetic secular v a r i a t i o n . 2.4 The Inversion The mantle attenuation of the p o l o i d a l f i e l d determines two sequences of c h a r a c t e r i s t i c values of two Sturm-Liouville problems. The eigenfunctions of e i t h e r boundary value problem form a complete orthogonal set of functions, with the e l e c t r i c a l conductivity as a weighting function. The orthogonality, i n the case of the boundary condition which generates the poles of a r a d i a l t r a n s f e r function, i s The eigenfunctions M(r,y_.), where y , j=l,2,--\u00C2\u00AB , are the poles of the t r a n s f e r function of the transverse p o l o i d a l f i e l d c o e f f i c i e n t s , describe another set of orthogonal functions. b (18) c 1 9 6 5 ) by -24-The normalizations or weights of equation (18) are given (Smylie, P ^ i / ^ V W C ^ ) \u00E2\u0080\u00A2 . (19) The f i r s t term of this expression i s derived from the transfer function G(b,u>) (equations 0-4)) of the transverse f i e l d coefficients. The second term i n equation (19) i s obtained by differentiating the i n f i n i t e product expansion form of H(b,u)) (equations (14)). The weights then take on the form A. A. j. n -o. zn-t-i - >-P \u00E2\u0080\u0094 T ( r > 2 n + 1 \u00C2\u00AB i - ^ > C i - ^ > . . . > j Vj C b . Vl V2 (20) A A A. \u00E2\u0080\u00A2 { ( l - ' r 1 ) \u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2 a - ' i ^ - x i - r 3 - ) \u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2>\u00E2\u0080\u00A2 \u00E2\u0080\u00A2 i A j - i A j + i This expression for the weights shows how the characteristic value sequences of two boundary value problems, determine the weights associated with either problem. The form of the weights may seem to indicate some preference for the Sturm-Liouville problem related to the radial transfer functions. The choice i s arbitrary. The use of the sets of weights and characteristic values related to the transfer functions of the transverse f i e l d for the inversion, involves only minor modifications i n the inverse theory. The result i s essentially the same. The inversion procedure i s developed with the spectral function of the 'radial' problem. Changes in the theory with alternate spectral functions, such as the one for the 'transverse' problem, are outlined later. -25-The expansion, in the eigenf unctions M(r,A..), of a piecewise continuous function f(r) may be written as b { 00 r \u00E2\u0080\u00A2' \u00E2\u0080\u00A2 . /a(r) M(r,A.) f ( r ) = I { f ( r)v^ ( r T M C r,A . ) d r } , (21) c for which both Parseval's formula b 00 ( f 2 ( r ) d r = I k I fCr)^CrT M(r,A )dr} 2 , (22) and the analogous relation b j=l \"3 a b f f(r)g(r)dr = . J { f f ( r ) ^ ) M(r,A )dr} / 1-1 P . i J J \u00E2\u0080\u00A2 { j g(r)/a(r) MCr,Xj#dr} , (23) c follow. As suggested earlier, the asymptotic values of the weights and characteristic values are required for the inversion. If the mean root squared e l e c t r i c a l conductivity i s defined as f b , /a(r) dr , a =\u00E2\u0080\u00A2 { r } , (24) mrs b-c then, for large j , the characteristic values A^ approach 2 2 ( 2 j - l ) V \u00E2\u0080\u00A2 , s \u00E2\u0080\u00A2 v - v \u00E2\u0080\u00A2[ > . (25) J 4u (b-c) a 0 mrs - 2 6 -The eigenf unctions are asymptotic (see E r d e l y i , 1956) to the sequence of functions r , . 1/4 r N C r . v ^ = C f ^ y ) cos{^Tvj\" J /cJ(rTdr} ,\u00E2\u0080\u00A2 j - 1 , 2,--* , which obey the orthogonality r e l a t i o n aO:)NCr,v )NCr,v k)dr = ~ \" / a C c ) a m r s 6 k = a 8k ' ( 2 6 ) The functions N(r,v^) also form a complete orthogonal s e t , which w i l l be used f o r expansions s i m i l a r i n form to equation (21) and i t s attendant r e l a t i o n s (equations (22) and (23)). The s t a r t i n g point of the Gel'fand-Levitan i n v e r s i o n , and b a s i c to i t , i s the expression f o r the eigenfunctions i n terms of t h e i r asymptotic values, and an unknown k e r n e l L ( r , s ) , i n the i n t e g r a l equation r \u00E2\u0080\u00A2o(r) M(r,A^) = /a(r) N(r,A..) + j L ( r , s ) / a ( s ) N(s,A..)ds . (27) c From t h i s p o i n t , the development of the i n v e r s i o n method takes two paths. The function L(r,s) i s determined from the s p e c t r a l function and the asymptotic values of i t s terms. The e l e c t r i c a l c onductivity i s r e l a t e d to the function L ( r , s ) . Following the theory of Gel'fand and Levitan, the d e r i v a t i o n of an equation f o r L(r , s ) i s a seemingly d i r e c t i o n l e s s development based on expansions i n the eigenfunctions and i n t h e i r asymptotic forms. The r e s u l t , a f t e r many manipulations, i s an i n t e g r a l equation, (equation (34)), i n L ( r , s ) , to which the reader may skip i f so moved. Required f o r the d e r i v a t i o n , Is the representation f o r the N(r,A^) as -27-/ c r ( r ) N(r,Aj) = VaXrT MCr.Aj) - j L ^ r ^ y a l s T M(s,A..)ds c The proof of the uniqueness of th i s representation, as with that of equation (27) , waits f o r the uniqueness of L ( r , s ) . To derive the i n t e g r a l equation, the following three r e l a t i o n s are required x y 00 ( f 1 \u00E2\u0080\u0094 { /JCTT M(r,A.)dr}-{ /a(r) N(r,A.)dr} = 0 , (28) x \u00E2\u0080\u00A2 I ^ { | /a(r) N(r,v. . )dr}-{ j /a(r) N ( r ,v . . )dr } = 0 , a e co a x I \ ^ T r T N ( r , v ) [ L ( s , r ) d s dr x (29) (30) + | / 0 ( r ) N(r,v.) | L(s,r)ds dr}\u00C2\u00AB{ j /a(r) N(r,v^)dr} .y x L ( s , r ) d r , a l l f o r c > e > y > x > a > b . (31) e a To obtain (29), equation (28) i s integrated to give -28-)dr - | /crCc) M(>,Aj)dr y r ** | L 1 ( r , s ) / a C s ) MCs,Xj)-ds dr e c y e y f r ' i^(rTM(r.A )dr - j /a(r) M(r,A ) L^(s,r)ds dr y / - |\" V a C r J M(r,Aj) | L 1 ( s , r ) d s dr r r = I ^aCrTMCr.Aj) \u00E2\u0080\u00A2 {1 - j L 1(s,r)ds}dr - j v\u00C2\u00A3 (r) M(r,A^) I L ^ s . ^ d s dr Therefore, the RHS of equation (29) may be written as b b I \u00C2\u00B1-{ [ f(r)\u00C2\u00BB^(rT M(r,A.)dr}-{ [ g(r)vo\"(rT M(r,A )dr} j = l j J J J J -29-vfaere a n d g C r ) = 0 1 0 c > r > x x > r > a a > r L 1 C s , r ) d s 1 - I L ^ C s . r J d s , c > r > e e > r > y , y > r Using the relation (23), this sum i s seen to vanish. Thus, (29) holds. Similarly, the LHS of (30) can be written co D b I \Z i [ f(r)yoTrT N(r,v.)dr}-{ ( g(r) N(r,v .)dr} j=l a ) J ) J (32) w h e r e and f(r) = 0 , c > r > x g(r) = 1 0 0 + 1 0 , x > r > a a > r , c > r > e , e > r > y y > r (33) Using the relation similar to (23) for the asymptotic approximations to the eigenfunctions, the sum (32) i s found to be zero and (30) i s valid. -30-The LHS of (31) can be giyen the form (32) where, in place of C33), fCr) = x L(s ,r) ds , c > r > x X X L(s,r)ds + L(s,r)ds , x > r > a , a > r and gCr) = 0 1 0 c > r > e e > r > y y > r Again using the relation similar to (23) for the asymptotic approximations to the eigenfunctions, the RHS of (31) is obtained. Integration of (27) yields x x | /a(r) M(r,A..)dr = j /a(r) N(r,A^)dr x + J VaTrY N(r,Aj) c a x x f / f L(s,r)ds dr L(s,r)ds dr S u b s t i t u t i n g the RHS of this equation i n t o (29), subtracting the LHS of (30) and the LHS of (31), and adding the RHS of (31) gives x y CO . f f I \u00E2\u0080\u0094 { yalrT N(r,A.)dr}-{ /a(r7 N(r,X .) dr} oo 3 X X X + I J-l [ /cr(r) N ( r , ^ ) j L(s,r)ds dr + j M r ) N(r,X^.) j L(s,r)dsdr} 3 J c a a r { f /a(r) NU . X^dr} x - \u00C2\u00A3 - { f /a(r) N(r,v.)dr}-{ f /a(r) N(r,v.)dr} j = l a . . J J J J J a e ob a x x x ^ ^ { | /a(r) N(r,v..) j L(s,r)ds dr + | /a(r) N(r,v..) j L ( s , r ) d s dr} y x \u00E2\u0080\u00A2{ | /a(r) N(r,v^)dr} + j j L(s,r)ds dr - 0 . e e a On d i f f e r e n t i a t i n g with respect to both x and y, the i n t e g r a l equation f o r L(x,y) r e s u l t s , x f F(x,y) + L(x,y) + F(s,y) L(x,s)ds = 0 , (34) for x < y , -32-where ^ \u00C2\u00AB NCx,V)N(y,A.) - N(5C,v.)NCy,v.) FCx,y) - A(x)a(y) J { ^- ^ L } P J (35) Equation (34) i s a Volterra equation of the second kind, with a symmetric kernel and a unique solution (Lovltt, 1950). The uniqueness of L(r,s) insures that the representations (equations (27) and (28)), assumed earlier, are also unique. From the above representation for F(x,y), the solution of the integral equation (34) is seen to depend only on the weights and characteristic values, that i s the spectral function, and their asymptotic approximations. The relationship between L(r,s) and the e l e c t r i c a l conductivity is derived from the substitution of the integral representation of the eigenfunctions (equation (27)), into the differential equation (see equations (12)). The notation i s simpler and the parallel with. Gel'fand-Levitan theory i s more apparent, i f the system of equations (12) is subjected to the variable transformations u = r \ \u00E2\u0080\u0094 f dr , (b-c)vk J mrs c Z = (c-b) y/a o mrs ( \ 1 / 4 0Cu,A) = (^-) M(r,A) (36) The differential equation for M(r,A) then becomes -33-2 : \u00E2\u0080\u0094 , 0 G l \u00C2\u00BB A ) +'U 2 -qGil) - 0 , C37) da ' where , v 1 d , , .1/4-, n(n+l) mrs .2 / 0 0 . < C u ) \" -77174 T T { a ( r ) } + \u00E2\u0080\u0094 T \" 7M CB\"C) \u00E2\u0080\u00A2 ( 3 8 ) a(r) du r The eigenfunctions are transformed i n t o u 0Cu,Jt^) = cosAjU + j K(u,v) cosAjV dv , (39) 0 with where KCu,v) = Cb-c)y^7 U r ' s ) . m r s ( a ( r ) o ( s ) ) 1 / 4 s v = [ /a(s) ds (b-c) v\u00C2\u00A3 J mrs c The i n t e g r a l equation f o r L ( r , s ) now becomes an i n t e g r a l equation f o r K(u,v), u f(u,v) + K(u,v) + j f(w,v)K(u,w)dw = 0 , v < u , (40) o f o r 0 < u < 1, where f(u,v) i s G O . f(u,v) = I { -^os^.ucosA.v - 2cos(2j-l)^ucos(2j-l)|v}. (41) j = l P j 2 . 3 .At t h i s point, the in v e r s i o n equations are i n the same form as those treated hy Gel'fand and Levttan. Arguments used here are borrowed l a r g e l y from t h e i r work. To r e l a t e the e l e c t r i c a l c onductivity, or at l e a s t the function q(u) (equation C 3 8 ) ) , to the s o l u t i o n K(u,v) of the i n t e g r a l equation (40), the i n t e g r a l representation f o r the 0(u , iL^) i s substituted i n t o the d i f f e r e n t i a l equation (37). From Gel'fand-Levitan inverse theory, t h i s s u b s t i t u t i o n uniquely implies the p a r t i a l d i f f e r e n t i a l equation 9 2K(u,v) . , s v , ^ S 2K(u,v) - q(u)K(u,v) = 2 ? \u00E2\u0080\u00A2 , 3u 3v i n the function K(u,v), and the conditions on K(u,v), m 0 3v v=o i f ^ L = fq(u) . (42) The d i f f e r e n t i a t i o n of 0(u,l/) with respect to u, i n equation (39), shows that K(0,0) = 3/3u(0(u,\u00C2\u00A3.) _ , which i s known from the boundary J \"CI 0 conditions. Equation (42) may then be integrated to 3 0 ( 1 1 , * , ) ; U + \u00C2\u00B1 q(u)du u=o ' 0 The function K(u,v), therefore, and i t s normal d e r i v a t i v e are s p e c i f i e d along the boundary u=v. The normal d e r i v a t i v e of K(u,v) , alone, i s given along the boundary v=0. The general theory of hyperbolic d i f f e r e n t i a l equations insures that these boundary conditions are- s u f f i c i e n t to uniquely -35-define KCu,y). The combination of equations (38) and (42) results i n the formidable non-linear d i f f e r e n t i a l equation , 2 z, x J / 4 , s mrs du r (u) a (r) du (43) for the recovery of the e l e c t r i c a l conductivity distribution. If the variable substitutions u y, (u) = r(u) = c + (b-c)Ja I d u * l v ' ' v mrs J r~7~\ 1 /o(u) y 0(u) = a 1 / 4 ( r ) are made, the di f f e r e n t i a l equation (43) becomes the system of three linked d i f f e r e n t i a l equations dy. (b-c) j/q mrs ( y 2 V 2 - y. du dy du\" d y 3 _ n(n+l) (b-c) 2 du \" ~ . .2 , .3 (y x) ( y 2 ) . _ dK(u,u) a +2 \u00E2\u0080\u0094 ^ ? y_ mrs du 2 (44) subject to the boundary conditions -36-y^Co) \u00C2\u00AB c , y 1 C D \u00E2\u0080\u00A2 b , C45) where the e l e c t r i c a l conductivity at r=c i s assumed to be known. The boundary value problem (equations (44) and (45)) may he integrated numerically as an i n i t i a l value problem. The value of y^CO) is. i t e r a t i v e l y adjusted to satisfy the condition that y^(l) = b. Hence, the proof that the e l e c t r i c a l conductivity distribution i s uniquely determined by two characteristic value sequences of two given boundary value problems, i s complete. The application of the Gel'fand-Levttan inverse theory has also outlined an inversion procedure which w i l l hopefully be numerically stabile enough to warrant its, application to real data. The numerical behavior of the inversion method i s the subject of the next Chapter. 2.5 Generalized Transfer Functions and the Inversion The uniqueness of the inversion of the geomagnetic secular variation i s established. The inversion has been proven with the characteristic values which correspond to the location of the poles of a radial transfer function. The theory of Gel'fand and Levitan allows that any two distinct characteristic value sequences are sufficient for the inversion. The characteristic values of the boundary value problem (12) are determined by the condition on the solution at r=b. This boundary condition, in i t s most general form, generates the poles of some arbitrary combination of the radial and transverse transfer functions, H(b,co) and G(b,to). To -37-generalize the inversion theory, the attenuation of two linearly independent combinations of the transformed coefficients, E m(r,iw) and n d/dr{E m(r,iu)} , over the lower mantle, determines the e l e c t r i c a l n conductivity structure uniquely. This is not of great theoretical significance as any new transfer functions are linear combinations of those already assigned. The uniqueness has been established and the method i s ready for numerical testing: the addition of a more general approach would seem to be of theoretical interest only. The motivation is based on a suspected improvement in the numerical behavior of the inversion method. This outline of some of the ideas and notation involved, i s included, partly as a prelude to the numerical experiments of the next Chapter. The relation between a general linear combination of the transforms of the poloidal f i e l d coefficients, at the core-mantle boundary r=b, and, at r=c , may be given, from equation (9), by be'Sm(b,IuO + 6 ^-{rZ mCr,ia3)} , n dr n r=b = {&'M(b,iw) + 3 4-MCr,iaj) }.cZm(c,iw) dr r=b n The f i r s t term of the RHS of this equation i s a combination of the two established transfer functions, oi-vr^c - \ _i_ o d \u00E2\u0080\u00A2 \ 3'b/c 3n/c 3 MCbjiu)) + 3 -r-M(r,:u)) - \u00E2\u0080\u00A2 -dr v ' 'r=h HCb,u) G(b,u) The l e f t hand side of this equation, is the general expression for the boundary condition at r=b on the Sturm-Liouvilleproblem (12). The -38-arguments used to relate the transfer functions to the characteristic values, of an associated boundary value problem, may- also be used here. The result i s an expression, F(b,w) say, describing the attenuation of an arbitrary combination of input coefficients. F(b,to) is related to H(b,a>) and \u00E2\u0080\u00A2 G(b,u) by \u00E2\u0080\u00A2 1 = g'b/c gn/c F(h,co) E(b,co) \" G(b,u>) (46) Previous reasoning, based on Sturm-Liouville theory, concerning the representation of the transfer functions as i n f i n i t e sums (equations (13)), or as i n f i n i t e products (equations (14)), is also applicable to F(b ,co). The formulae for the residues and the zero-frequency response, F(b,0), only, are different. The ratio of F(b,w) to the attenuation of any other linearly independent combination of transformed input coefficients, of the same harmonic degree, would complete the information necessary for the inversion. The ratio F(b,co)/E(b ,0)) for example, determines two characteristic value sequences and therefore, the e l e c t r i c a l conductivity distribution. The two possible approaches are equivalent, the use of the transfer functions, for the radial and transverse poloidal fields being a particular case of an array of possible inverse Sturm-Liouville problems. There are a few changes in the formulae used to develop the inversion method. Representing the pole locations of the transfer function F(b,co) as ico = rij , j=l,2,\"* , the orthogonality of the appropriate eigenfunctions of the boundary value problem (12) is given by -39-where the normalizations, p^, are, i n this case, _ 1 \u00E2\u0080\u00A2\u00E2\u0080\u00A2rdMCr,*)1 . d r M f , 8 dMCr,A)i , , p j \" v {dr }r=b d T i M C b ' A ) + 8 V s The function f(u,v) i n the integral equation (41), becomes f(u,v) = ~\u00E2\u0080\u0094\"COS&.. U C O S ^ i . v - 1 P j 1 1 + \ {\u00E2\u0080\u0094cosA .ucosA .v - 2cos(j-l)T7ucos(j-l)irv} . (48) j=2 p j - 3 J The inverse theory i s mathematically complete, i f physically confused, i n this more general treatement. CHAPTER III NUMERICAL TESTS OF THE INVERSION 3.1 General Alms The inversion theory guarantees that the e l e c t r i c a l conductiv-i t y of the lower mantle may be recovered uniquely from the ratio of two transfer functions. These two transfer functions describe the attenua-tion of the poloidal f i e l d of a given harmonic degree. Hidden i n such statements of uniqueness, however, are assumptions of arbitrary accuracy. Many questions of practical importance to the usefulness of the method remain to be answered. The analytical continuation of the ratio of the transfer functions from the real frequency axis may pose numerical problems. How sensitive i s the whole inversion procedure to inaccurate and limited information of the attenuation properties? What kind of resolution i n the conductivity structure i s recoverable with how much attenuation information? These and other questions are investigated with synthetic conductivity profiles. An understanding of the numerical limitations of the inversion method i s essential before any attempt at the inversion of the secular variation, using real data, should be made. The starting point of a l l the numerical inversions i s the assumption of a lower mantle conductivity p r o f i l e . E l e c t r i c a l conductiv-i t y profiles consisting of an exponential increase with depth and super-imposed sinusoids are used. The height and wavelength of the sinusoids, or 'bumps' i n the model, the mean conductivity over the lower mantle, -41-and the harmonic degree n are varied independently in anticipation of possible variations in the numerical behavior of the method. Different values of 3 ' and 3 , which determine the boundary condition on the eigenfunctions at r=b , are also used to generate different spectral functions for the same conductivity model. Banks' (1969) value of 1 mho per meter for the e l e c t r i c a l conductivity at r = c and a lower mantle thickness of 2400 kilometers are assumed in a l l conductivity models. The general numerical procedure for the synthetic inversions is as follows. For a given conductivity model and a given harmonic degree, the boundary value problem (12) is integrated numerically to generate the sequences of characteristic values and associated weights, to some arbitrary number. At this point, the spectral function of the Sturm-Liouville problem i s defined, the conductivity model may be forgotten and the inversion is carried out numerically. Knowing the characteristic values and weights is not, however, in practice, equivalent to knowing the ratio of two transfer functions for real frequencies only. It is f i r s t necessary to look at the transfer functions and the ratio of transfer functions, for requirements on their accuracy and extent in the recovery of the two characteristic value sequences and, therefore, the spectral function. 3.2 Transfer Functions to Spectral Function As a prelude to the recovery of the spectral function, i t is interesting to look at the general nature of the attenuation properties of the lower mantle. Much can be deduced from the mathematical forms of -42-the transfer functions - equations (14). Whether representing the attenuation of the radial part of the poloidal f i e l d of a certain harmonic degree, or the transverse part, or a combination of the two, the transfer function reflects the fact that the earth's mantle, as any conductor, i s , electromagnetically, a low pass f i l t e r . The amplitude attenuation of either poloidal coefficient, sn(*,t) or (rs m) , and therefore of any poloidal f i e l d vector component, across the lower mantle increases with increasing frequency. The transfer functions should also show a growing independence of the conductivity of the mantle, at a given frequency, as the harmonic degree increases. The power transfer function curves should 'flatten' as n increases. These and other properties are illustrated i n Figure 1 which shows the dependence of different normalized power transfer functions, F(b,0))-F*(b,a))/F(b,0) 2 , or H(b ,u) -H*(b ,w)/H(b,0)2 or G(b,u)-G*(b,u>)/G(b,0)2 for the combined, or radial, or transverse f i e l d coefficients respectively, on frequency. As the mean conductivity across the lower mantle increases, the frequency, beyond which the surface expression i s a negligible frac-tion of an input f i e l d harmonic component, a cutoff frequency, decreases. From Figure 1, curve (a) shows that with a conductivity p r o f i l e with a mean value of about 50 mhos per meter, the amplitude response at the surface is near zero at a frequency of about one cycle per year. A crude comparison of the spectrum of the secular variation at the surface with an input f i e l d , which i s assumed to contain equal power at a l l frequencies for any harmonic, favours this type of transfer function. Therefore something l i k e this value for the mean e l e c t r i c a l conductivity Figure 1. Normalized power transfer functions, (a) a -60 mhos per mrs meter, characteristic values of spectral function generated by B',g^o , a 4-=? 5ir, n=3.(b) Power ratio curve with (a) and one for 3=o. (c) a - 6 mhos per meter, 3 ' , 3 4 o, a^ = 9 rr, n=3. (d) same as (c) but a^ = 57r.s(e) a - 1.6 mhos per meter, 3 = o, ai+ = 5TT, n=l. (f) same as (e) but n= 3. ak specifies the frequency of the super-imposed sinusoids (see page 48). -44-over the lower mantle i s suspected. One of the anticipated benefits from an inversion theory i s the resolution of structural detail, beyond the simple models used to date, i n the conductivity p r o f i l e of the mantle. C r i t i c a l to the resolv-ing power of the method i s the effect with which structural details are registered in the transfer functions. The curves (c) and (d) i n Figure 1 are those for two conductivity models of approximately the same mean conductivity (curve (c) i s for the larger \u00C2\u00B0 m r s ) \u00C2\u00BB a n <* for the same harmonic degree. The length scale of the structural 'bumps', the super-imposed sinusoids, in the model used to generate the curve (c) i s half that used for the curve (d). The difference between these two power transfer function curves w i l l seem small in view of the. complex numerical operations which follow. As predicted, the power transfer function curves flatten out as the degree of the harmonic component of the f i e l d studied, increases. Curves (e) and (f) in Figure 1 are the attenuation curves for a conduc-t i v i t y model with a a -1.5 mhos per meter and n = 1 and n = 3 J mrs * respectively. The amplitude attenuation at any frequency, however, increases with increasing n , as the power transfer functions are pro-portional to ( b / c ) n + ^ . The inversion requires the ratio of two transfer functions. Curve (b) is the ratio of the power transfer functions, H(b ,OJ) \u00C2\u00ABH*(b,OJ) /(F(b ,co) *F*(b ,a>)) , where the normalized power curve of F ( b , u ) , shown i n curve (a), i s used. The sensitivity of the inversion method towards conductivity structures depends on the accurate recovery of the locations of the zeros and poles of the ratio of two transfer functions. The zeros and poles, -45-that i s the two characteristic value sequences, are also recoverable from the ratio of the power transfer functions. If , j - l , 2 , * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 , are the zeros and X^ ,j=l,2,\u00C2\u00AB \u00E2\u0080\u00A2 \u00E2\u0080\u00A2, are the poles of the ratio H(b ,co)/G(b ,co) of two transfer functions, the ratio of the power transfer functions, herein called the power ratio, i s \u00C2\u00AB b . , . ) H * ( b . \u00C2\u00AB ) . g + - 2 / - ? ( i + - 2 / f l - t ( w G(b,o))G*(b,aj) (1+ < u 2 / A 2Ml+ co2/*2) \u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2 for a given harmonic degree. Analytic continuation insures that the power ratio's poles and zeros, which are now mirrored by the real u> axis, are uniquely determined by the power ratio on the real axis. In theory, each pole and each zero of the power ratio w i l l have some effect on the power ratio along the real axis, no matter how far away. If the power ratio i s known to arbitrary accuracy for some f i n i t e length along the real frequency axis, the zero and pole locations should, theoretically, be recoverable to arbitrary accuracy and number. In practice, i t i s a d i f f i c u l t numerical problem: i n s t a b i l i t i e s and slow convergence are not unexpected. The accuracy and extent of the power ratio spectrum, needed to recover a certain number of zeros and poles to predictable accuracy, are d i f f i c u l t to estimate. For real frequencies, the effects on the power rat i o , of the zeros and poles further from the real axis, are relatively more at greater frequencies. To recover the pole and zero locations which are further away, therefore, given certain accuracy limitations, the power ratio must be known to proportionately greater real frequencies. For high frequencies, the power ratio should be known to frequencies -46-greater in magnitude than the furthest pole or zero location desired, in order to determine the pole and zero locations up to this maximum. This implies, for example, that ten zeros and ten poles, that is two characteristic value sequences of ten terms each, require the knowledge of the power ratio of curve (b) (\u00C2\u00B0\" m r s - 50 mhos per meter) , for up to more than 60 cycles per year. The magnitude of the problem of the recovery of the two charac-t e r i s t i c value sequences from the power ratio, known for real frequency only, can only be fully appreciated with real data. Speculation on the eventual application of the inversion method would predict this step to be a hazardous one. There are associated numerical problems before the inverse Sturm-Liouville theory is tested. The other half of the spectral func-tion, the sequence of weights, is determined from an i n f i n i t e product expansion (equation (20)), with the two characteristic value sequences. The problem of the errors, arising out of the analytic continuation, transferred into the weights i s augmented by the need for long character-i s t i c value sequences for convergence of this i n f i n i t e product. This i s alleviated by a decreased demand on the accuracy of the larger character-i s t i c values. The asymptotic form of the terms, and p.., i n the spectral function (equations (24), or (47), and (25)), required in the inversion, depend upon the mean root squared value of the e l e c t r i c a l conductivity over the lower mantle. This quantity may be estimated directly from the analytically determined form of M(b,iu) for high frequency, compared through equation (9), with the high frequency behavior of the input and output f i e l d coefficents. -47-In short, the recovery of the spectral function of the d i f f e r -ential system demands much of the appropriate power ratio. Spectral functions known to only ten terms have been used as representative of the type needed to recover the el e c t r i c a l conductivity. It remains to be shown that this number of terms, that is this amount of information, i s , in fact, enough to achieve reasonable results for the inversion. -48-3.3 Spectral Function to Conductivity Distribution The recovery of the spectral function, that i s the sequence of characteristic values and weights, from the ratio of two power transfer functions, i s a demanding numerical operation. Unavoidable errors in the terms, X_. and p_. , of the spectral function increase as the order, j , increases. The sensitivity of the inversion method, from the spectral function to the conductivity distribution, to such errors is investigated later i n this section. For the moment, the terms of the spectral function are as accurate as the numerical integration of the boundary value problem (equations (12)) allows, thus isolating any numerical problems in this part of the inversion. Many different conductivity models are used in testing the numerical behavior of the inversion. A l l models are of the general form -2(~P) crCr) n a\u00C2\u00B1 + a 2e + a 3 s i n C a 4 -g^ ) The constants a^, a 2, and a 3 determine the relative amplitudes of the exponential and sinusoidal parts of the model, and a^ specifies the wavelength of the sinusoids superimposed on the exponential. As the power transfer function curves change with the harmonic degree, the dependence of the numerical inversion on n is also investigated. The spectral function is generated by the numerical integration of the boundary value problem (12). (The program to calculate the spectral function for a given conductivity model, degree n , and boundary conditions, given by 8' and 3 , is liste d in the third Appendix.) Given the characteristic value and weight sequences, X. and -49-, j=l,2,''\u00C2\u00ABN, for a given harmonic degree, the inversion is carried out numerically as follows. The integral equation (40) is solved for dK(u,u)/du by truncating the series expansion form for f(u,v) (equations (41) or (48)). (The method used for the numerical calculation of dK(u,u)/du is outlined in Appendix A-l.) The values of dK(u,u)/du are substituted into the non-linear system of equations (44) . The boundary value problem (equations (44) and (45)) is numerically integrated to recover the e l e c t r i c a l conductivity. (The computer program for the inversion from spectral function to conductivity distribution, i s l i s t e d in the fourth Appendix.) Comparison of the recovered conductivity profile with the assumed model measures the success of the inversion. The accuracy of the recovered profile echoes the effects of, among other things, the truncated series f(u,v) . The validity of the truncation depends upon how close the characteristic values and weights come to their asymptotic values. As more terms are included in the series, the solution of the integral equation is more accurate and the recovered conductivity distribution converges to i t s true value. The whole inversion method may require very long spectral functions for even moderate s t a b i l i t y . This would seriously limit the usefulness of the inversion method. The relationship between the length of the spectral function and the accuracy of the recovered conductivity profile i s fundamental to an understanding of the limitations of the inversion method. There are many factors which determine the rate at which the terms of the spectral function converge to their asymptotic values. For attenuation properties of higher harmonic degree, this convergence is -50-slower. This conclusion is based on the flattening of the power transfer function curves as n increases. Changes in the height and width of the conductivity structures of the model are also expected to influence the length of the spectral function required to achieve a reasonable recovery. Depending upon which set of characteristic values is used, the spectral function i s different f o r the same conductivity model and harmonic degree. The numerical experiments are designed to establish what amount of accuracy and resolution i s possible with how much of the spectral function and how this dependence varies with different conduc-t i v i t y models, harmonic degrees, and types of spectral functions. The f i r s t test of the inversion from spectral function to ele c t r i c a l conductivity is for a model conductivity with a mean value of approximately 1 mho per meter. The superimposed sinusoids are long i n wavelength, 1000 kilometers, and short in height, 0.15 mhos per meter. The result of the inversion, for n=l , with the mean root squared conductivity, a =1.6 mhos per meter and with three bumps in the mrs model (there are actually only two and one half complete sinusoids), i s shown in figure 2. This recovery results when two terms of the spectral function are used, that is when the series f(u,v) is truncated at two (N=2) terms. The characteristic values in the spectral function are, in this case, the poles of the transfer function for the radial c o e f f i -cients of the poloidal f i e l d . The recovered conductivity distribution compares poorly with the original model, especially, near the base of the mantle. The inab i l i t y of the numerical inversion to resolve the three bumps in the model might have been expected. The characteristic values had .units of (time) ^ -51-0.0 50.0 100.0 150.0 , 200.0 250.0 DEPTH FROM C (KM.) CX103 ) Figure 2. Inversion result (solid line) and model (dashed lin e ) . N = n = 1, characteristic values of spectral function from 8 = 2, o. -52-and were called the characteristic frequencies of the system. The role of the transformed characteristic values (see equations (36)), Is as spatial characteristics, as they have units of 1/u where u 'is measure of distance. As with Fourier series, i t may be possible to recover only those components of the conductivity model which are 'normal' to the eigenfunctions of the f i r s t two characteristic values. The analogy i s not exact as the eigenfunctions i n this problem are only asymptotic to cosines. The differences between the behavior of Fourier expansions and the inversion are d i f f i c u l t to predict. The type of conductivity model used may be thought of as a perturbation of the constant conductivity case, for which the eigenfunctions, with n=0 , are cosines. Discontin-uities in the conductivity model would therefore make high spatial fre-quencies of the spectral function important to the inversion; a behavior similar to Gibb's phenomenon. On the basis of this explanation for the significance of the length, N , of the spectral function, the inversion for the same conductivity model should recover, for N=3 , the three bumps in the original model. The result of this inversion test are shown in figure 3. The sinusoids are recovered, i f displaced and exaggerated i n height. The recovered values of the e l e c t r i c a l conductivity near the base of the mantle are s t i l l poor. Figures 4 and 5 show the results of the inversion for the same conductivity model as in figure 3, but with the series in the integral equation truncated at N=5 and N=10 terms respectively. The recovery improves as more terms of the spectral function are included, -53-50.0 100.0 150.0 , 200.0 DEPTH FROM C (KM.) (X101 ) 250.0 Figure 3. Inversion result (solid line) and model (dashed li n e ) . N n = 1, characteristic values of spectral function from 8 (Flat section indicates truncation for graphing.) = 3, = o. -54-Figure 4. Inversion result (solid line) and model (dashed l i n e ) . N = 5, n = 1, characteristic values of spectral function from 6 = o. Figure 5. Inversion result (solid line) and model (dashed lin e ) . N = 10, n = 1, characteristic values of spectral function from 6 = o . -56-the agreement being good for N=10 at a l l r a d i i except near r=b . Thus, i t appears that the inversion method requires a spectral function of up to at least N terms to resolve conductivity structures which have length scales of (b-c)/N, or longer. The effect of increasing the mean conductivity of the model, that is with the e l e c t r i c a l conductiv-i t y distribution less like a constant over the lower mantle, may alter this conclusion. The small errors in the terms of.the spectral function are determined by the numerical integration of the boundary value problem (12). These characteristic values and weights are accurate to at least 0.1%. As i t is unrealistic to imagine such accuracy to be available with real data, the inversion is tested under more adverse circumstances. To imitate terms of the spectral function which have an unpre-dictable error, a uniform random number generator is used. To the true value of the terms in the spectral function, is added some fraction of a random number. The average value of the random number, 0.5, and the fraction, dictate the average percentage error introduced into the weights and characteristic values. To simulate what is expected with a frequency limited power transfer function, this percentage error is increased linearly from zero to some maximum, as the order of the terms in the spectral function increases. A 1% random error, for example, would mean no error in the f i r s t characteristic value and weight but the percent error is increased to an average value of \u00C2\u00B1 1% for the last, t i l the N , term of the spectral function. The result of the inversion with the e l e c t r i c a l conductivity, harmonic degree,and boundary conditions of figure 5, where the spectral function terms have an average random error of less than \u00C2\u00B1 1% is shown in -57-figure 6. The comparison is favourable. A tolerance of \u00C2\u00B1 1% in the higher terms of the spectral function does not, in this case, have any serious effects on the recovered profile. Figure 7 is the same inversion problem as that of Figure 6, except the percentage error reaches a maximum mean value of \u00C2\u00B1 5% in the characteristic values and weights. This situation, which allows for as much as \u00C2\u00B1 10% in and , results i n a high spatial frequency error in the recovered conductivity distribution. The recovery i n figure 7 shows that the effect of the \u00C2\u00B1 5% random error, is to a r t i f i c i a l l y enhance a high frequency mode of the conductivity to where i t almost swamps the longer wavelength information. As the higher terms are the most inaccurate, a general correlation persists between the resolving power of the inversion method and the length of the spectral function used. The flattening of the power transfer function curves, as the harmonic degree increases, i s a symptom of the migration of the poles of the transfer functions away from the real frequency axis. As the asymptotic approximations to the characteristic values are independent of n , the convergence of the terms in the spectral function to their asymptotic values i s expected to get worse as n increases. There are practical reasons for studying the numerical effects on the inversion with the attenuation properties of higher harmonics of the f i e l d . The case n=l corresponds to the dipole f i e l d , which is known to undergo only slow changes with time.. The frequency information about the secular variation is therefore limited, thus jeopardizing the recovery of higher order terms of the spectral function. ... -58-Figure 6. Inversion result (solid line) and model (dashed lin e ) . N = 10, n = 1, terms of spectral function have a maximum average \u00C2\u00B1 1% random error and are from 8 = o. -59-50.0 DEPTH 100.0 FROM C 150.0 , 200.0 (KM.) (X101 ) 250.0 Figure 7. Inversion result (solid line) and model (dashed lin e ) . N = 10. n = 1, terms of spectral function have a maximum average \u00C2\u00B15% random error and are from 8 = o. (Flat section indicates truncation for graphing.) . . . - 6 0 -Th e recovered profile for the harmonic degree n=3 and N=10 , \u00E2\u0080\u00A2with the same conductivity model and spectral function type as in figures 2 through 7, is shown in figure 8. This recovery should be compared with figure 5, the inversion result for the dipole component of the poloidal f i e l d . The recovery is poorer for larger n , the difference being most pronounced near the core-mantle boundary. The error in many of the recovered conductivity profiles has the appearance of an exponentially modulated sinusoid. The spatial frequency of this error roughly corresponds to the value of N . The dependence of the sinusoid error in the recovered conductivity on the truncation point, is found to be consistent. This behavior is similar to Fourier series where the function to be expanded contains higher modes than those used in the truncated series. The effects, on the inversion, of truncating the spectral function, are inconsistent with Fourier analysis, as figure 5.shows, i n that long wavelength structures require higher modes of the eigenfunctions for resolution. An i n s t a b i l i t y of the inversion in the region of the core-mantle boundary, is a general feature of a l l the synthetic inversions so far. The only way to improve the recovery near r=b , i t appears, is by including more terms of the spectral function, that is by using higher spatial frequency information in the inversion. The behavior of the error near r=b is similar to that of the eigenfunction expansion of a discontinuous function. When the eigenfunctions, corresponding to the boundary condition M(b,X^ .) = 0 , are used, this i n s t a b i l i t y i s a result of expanding a non-zero function in a series of orthogonal functions, which are zero at one boundary. Many terms of the eigenfunction expan--61-- i i i i I r -> 0.0 50.0 100.0 150.0 200.0 250.0 DEPTH FROM C (KM.) (X10] ) Figure 8. Inversion result (solid line) and model (dashed l i n e ) . N = 10, n = 3, characteristic values of spectral function from 8 = o. -62-sion are required, i n the area of the discontinuity, as the series is asked tb reproduce an i n f i n i t e slope. To improve the recovery near r=b , the eigenfunctions which are non-zero at r=b are used. That i s , the poles of the transfer functions of the transverse f i e l d or of the combined fields are used as the characteristic values of the spectral function. The result of the inversion with the conductivity model and harmonic degree of figure 8, with the spectral function of the 'transverse' problem, that is with characteristic values for the boundary condition dM(r,X)/r=b = 0 , Is dr shown i n figure 9. The error of the recovered conductivity distribution is of the same magnitude as in figure 8, but the error sinusoid i s shifted by a phase of \u00E2\u0080\u00A2-ir near the core-mantle boundary. The error has simply changed sign from the 'radial' inversion. The average of the two profiles (figures 8 and 9) would remove much of the error i n the recovered conductivity distribution near the base of the mantle. The possibility of superimposing the two recoveries with the possible spectral functions, from the same attenuation properties, to remove the Gibb's phenomenon behavior of the inversion, is investiga-ted at the end of this section. The use of the spectral function from the 'transverse' problem does not correct this convergence problem but causes a similar error effect as i n the 'radial' inversion. The apparent dependence of the error in the recovered conductivity profile, on the nature of the charac-t e r i s t i c values, suggests a combined boundary condition to generate the characteristic value sequence of the spectral function. The result of the inversion for the boundary conditions, given by 8'/B= -n/b , (see equation (12)), which should give equal weight to Figure 9. Inversion result (solid line) and model (dashed l i n e ) . N = 10, n = 3, characteristic values of spectral function from B'= o. -64-the radial and transverse transfer functions, is shown in figure 10. The recovery is almost identical to that in figure 9. This result is reasonable as the asymptotic approximations to the characteristic values of a l l boundary value problems, -00<8'/3 > e> figures 17, 18, and 19, respectively. The performance of the whole Inversion method with limited information is the most important for this type of e l e c t r i c a l conductivity model, as just this type of mean conduc-t i v i t y i s favoured for the earth's lower mantle. The results show that the recovery of long wavelength structures improves as N increases. Any direct relationship between N and the a b i l i t y of the inversion to recover structures of frequency up to (N-l) cycles over the lower mantle, no longer applies. The information necessary for the resolution of longer wave-length structures is contained i n many terms of the spectral function. The correlation between Fourier modes of the conductivity and the trunca-tion of the spectral function was based on the results for conductivity models which were much closer to the Fourier case r.that i s , closer to a -71-Figure 14. Inversion result (solid line) and model (dashed lin e ) . N = 20, n = 3, characteristic values of spectral function from 8 ' , 8 ^ o . Figure 15i Inversion result (solid line) and model (dashed l i n e ) . N = 1 0 , n = 3 , characteristic values of spectral function from B',8^o. -73-Figure 16. Inversion result (solid line) and model (dashed l i n e ) . N = 20, n = 3, characteristic values of spectral function from 8 ' , 8 ^ o . -74-Figure 17. Inversion result (solid line) and model (dashed l i n e ) . N = 10, n = 3, characteristic values of spectral function from. B',8^o. (Flat sections indicate truncation for graphing.) , -75-Figure 18. Inversion result (solid line) and model (dashed l i n e ) . N = 20, n = 3, characteristic values of spectral function from 8 ' , B ^ o . -76-C3 a 0.0 50.0 100.0 150.0 200.0 250.0 DEPTH FROM C (KM.) (X101 ) Figure 19. Inversion result (solid line) and model (dashed l i n e ) . N = 30, n = 3, characteristic values of spectral function from $',6^o. -77-constant value. When the el e c t r i c a l conductivity at r=c is tied to 1 mho per meter and the model has relatively large rates of change in the conductivity with depth, the analogy with Fourier analysis breaks down. Using the rough guide for analytic continuation, the use of 10, 20, and 30 terms of the spectral function implies the ratio of the power transfer functions is known to up to more than 60, 250 and 600 cycles per year, respectively. The recovered conductivity distribution is comparable, for N=20, to the result for N=3 with a =1.6 mhos e ' * mrs per meter (figure 3). The frequency requirements on the ratio of the power transfer functions is of the same order in both cases. The recovered values of the e l e c t r i c a l conductivity near the base of the mantle are as inaccurate i n this series of synthetic inver-sions as before. Earlier results showed that i f one set of the charac-t e r i s t i c values, of the two required for the inversion, i s that from the poles of the radial transfer function, the two possible inversion results are alike except for a change in sign of the error near r=b . The recovered profiles from the spectral functions could be averaged to dampen a high spatial frequency error. The result of such a superposition is shown in figure 20. This result has used the recovered profile of figure 18 as half of its average. The error is attenuated enough to make this extra work worth-while. There is s t i l l a large error in a(r) at the base of the mantle. Given the p o s s i b i l i t i e s for numerical error in the inversion, the i n s t a b i l i t i e s in the superposition of the two profiles and no guarantee of exact error cancellation, the improvement in the. recovered conductiv-i t y i s sufficient to justify general use of this average conductivity -78-Figure 20. Inversion result for the average of two recoveries (solid line) and model (dashed l i n e ) . N = 20, n = 3, Figure 18 and another from the spectral function with characteristic values from 8 = o. -79-p r o f i l e as the f i n a l resul t . It should be noted that the ins tab i l i t y in the inversion method near r=b i s part icular ly regrettable in view of geophysical \"interest in the physical properties of the lower mantle in this region. The resolution of the e l e c t r i c a l conductivity at the base of the mantle is important to considerations of core-mantle electromagnetic interactions. A large number of different conductivity models have been used to generate spectral functions for testing the inversion method. The results show that the accuracy of the inversion for a given amount of attenuation information i s independent of the type of character ist ic values used as those of the spectral function. The recovery i s dependent on the harmonic degree and roughly independent of a m r s ' It i s possible to extract a rough rule relat ing the resolving power of the inversion method and the frequency length of the rat io of the power transfer functions. A conductivity with a mean of 1 mho per meter, may be used as a base for which a one to one relationship ex is ts , through the inversion, between the Fourier modes of the conductivity d is t r ibut ion, and the order of the terms of the spectral function. The minimum number of terms of the spectral function, required to accurately recover sinusoids in the conductivity p r o f i l e , which have a spat ia l frequency of M/(c-b) or l ess , i s roughly N \u00C2\u00AB (M - h v \u00C2\u00A3 ~ ~ + | I mrs 2 Hence a conductivity p rof i le with a mean root squared value of 60 mhos per meter, such as used for figures 17 through to 20, requires more than twenty terms in the spectral function to resolve sinusoidal structures -80-with a wavelength of 800 kilometers. This is a consequence of the assumption, supported by the numerical tests, that the ratio of the power transfer functions, known up to some frequency, determines any elec t r i c a l conductivity distribution to the same relative.accuracy. Applying the rule for the analytic continuation, the wave-length, L. , of conductivity structures, recoverable from power ratios known to up to cycles per year, i s , from the asymptotic form of the characteristic.values (equations (25) or (41)), (y (b-c)2w ) 1 / 2 L - ( c - b ) - { \u00E2\u0080\u00A2 \u00E2\u0080\u0094 2 \u00E2\u0080\u0094 j ^ + I}\" 1 If the values of the constants are used, L -(1+ 0.15*^) where u> is in cycles per year. This means that a power ratio, known for frequencies up to 50 cycles per year, is enough information to recover two sinusoids in the el e c t r i c a l conductivity profile. These estimates of the a b i l i t i e s of the inversion method are necessarily crude, valid for large M , and subject to decay as the harmonic degree increases. It is established that the types of transfer functions in the required ratio used do not affect the inversion. There is a particular advantage to including the 'radial' transfer function as one of the two, as the two possible inversion results may be superimposed to reduce the error in the recovered electrical conductivity near the base of the mantle. . -81-Hence a rough understanding of the numerical demands of the inversion i s available. It' is now possible to use what is known, or suspected, about the attenuation of the poloidal f i e l d across the earth's lower mantle, with some assurance as to the accuracy of the outcome of the inversion. CHAPTER IV SUGGESTIONS TOR THE APPLICATION The. inverse theory of the diffusion of the secular variation has established that the e l e c t r i c a l conductivity distribution of the lower mantle Is uniquely determined by the ratio of two power transfer functions. Numerical testing with many conductivity models has shown that the inversion algorithm is numerically stable enough to permit the recovery of a meaningful conductivity profile from a limited knowledge of the attenuation properties of the lower mantle. Short of carrying out the inversion with real data, i t i s important to consider sources of information, actual and theoretical, available about the attenuation. What i s the simplest and most accurate way of estimating the power ratio and what would be the limitations of such a method? In order to complete an evaluation of the f e a s i b i l i t y of the inversion of the secular variation, these and other questions must be investigated. The ratio of the transfer functions, radial and transverse, i s given in terms of the Fourier transforms of the harmonic coefficients of the poloidal f i e l d , from equations CIO) and (11), as \u00E2\u0080\u0094 - f r m v H (b,cj) = _ JL dr n r=b G (b,u) U Zm(b,iu>) n ' n ' Hence, the e l e c t r i c a l conductivity of the lower mantle is uniquely determined by the differences between the radial and transverse parts of the poloidal f i e l d at the base of the mantle. -83-'The assumption of zero e l e c t r i c a l conductivity at c > b, i s responsible f o r what i s , i n f a c t , the a l t e r a t i o n of the i n v e r s i o n problem from a study of the transmission through the mantle, to the r e f l e c t i o n at the base of the mantle, of the p o l o i d a l f i e l d . An a p p l i c a t i o n of the i n v e r s i o n method could r e l y , therefore, e n t i r e l y on t h e o r e t i c a l estimates of the p o l o i d a l f i e l d at the core-mantle boundary. The secular v a r i a t i o n measured at the surface of the earth would not be involved. A l o g i c a l approach f o r estimating the r a t i o of the power tr a n s f e r functions i s through a s t a t i s t i c a l d e s c r i p t i o n of the f l u i d motions i n the core. For high harmonic degrees, the f l u i d motions i n the core might be considered as turbulent flow, and s t a t i s t i c a l considerations would apply. A convenient expression f o r the s t a t i s t i c a l p roperties of the p o l o i d a l f i e l d components, i s through s p a t i a l autocorrelations. The a u t o c o r r e l a t i o n of a f u n c t i o n , f (b ,0 ,<(> ,ico) , over the the surface of the core-mantle boundary, i s given by $ f f(b,e,,6 f , + ' ,i\u00C2\u00AB) = E {f (b , 9 , ,ito) \u00E2\u0080\u00A2 f *(b , 9 1 ,' ,ico) } where E i n d i c a t e s the expected value. Then the geometical convolution, (see Appendix A-2), - 8 4 -$ f f Q j , 8,$,6' ,<|. r , iu) TT TT \u00E2\u0080\u00A2if 2rr d(j)\" 5>p C B ^ \" , ^ ' ' ^ ' , ,iu)) i s related to the Fourier transforms of the transverse components of the poloidal f i e l d at r=b, by 3 2 PAP\u00E2\u0080\u009E . 3036' f f 2 P.P.,. sin0sin0' 343d)' f f 9 9 . The power ratio of the radial and transverse transfer functions of a given harmonic degree n, is then recoverable from the deconvolution of equation ( 5 0 ) and then through H Cb,u>)H*(b,a)) 2 / = (n+1) . E(0,co)P (cos0)sin0d0 G (b,w)G (b,oj) . . J n n n o These equations are just one way of relating the poloidal f i e l d components through a convolution integral. The number of possible convolution expressions for the poloidal f i e l d , and therefore, possible deconvolutions' to recover the attenuation properties, i s almost endless. -85-The overdeterminacy of the inversion method, due to the dependence of the attenuation properties, and the ratio thereof, on the harmonic degree, remains unexploited. The relating of the stati s t i c s of the poloidal f i e l d components, through a convolution Integral, i s then one suggestion for the recovery of the power ratio of the transfer functions. Its application to the secular variation' and the determination of the s t a t i s t i c a l error limits on the recovered conductivity distribution, w i l l have to await further study. CHAPTER V CONCLUSIONS The poloidal part of the magnetic induction E , generated in the earth's core, diffuses through the mantle and i s measured on the earth's surface. Within the limits of the assumptions of spherical symmetry and of a non-conducting upper mantle, the attenuation, through the mantle, of the radial and transverse coefficients of any harmonic of the poloidal f i e l d , uniquely determines the e l e c t r i c a l conductivity distribution of the lower mantle. A unique recovery of the conductivity i s , i n fact possible, with no more information than the ratio of the amplitude spectra of two transfer functions. Two transfer functions describe the attenuation with frequency of the poloidal f i e l d of a given harmonic degree. The power ratio of amplitude spectra i s determined solely by the poloidal f i e l d components at the core-mantle boundary. The derivation of the inversion theory is based on work by I. M. Gel'fand and B. M. Levitan with the inverse Sturm-Liouville problem. The marriage of the frequency characteristics of the attenuation of the poloidal f i e l d , and general Sturm-Liouville theory, i s basic to the theoretical development of the inversion method. The proof of uniqueness also outlines a numerical method which may be applied to the inversion of the secular variation from real data. Numerical experiments with model e l e c t r i c a l conductivity distributions, indicate an algorithm with enough s t a b i l i t y to give good results for the - 8 7 -inversion with limited attenuation information. Given the power spectrum of the ratio of two transfer functions, up to a certain frequency, the resolution and accuracy of the inversion method, are found to be predictable. A power ratio, known for up to f i f t y cycles per year, for example, would insure the detection of 1000 kilometer wavelength structures i n the ele c t r i c a l conductivity distribution of the lower mantle. The implementation of this inverse theory rests on what i s known or suspected about the attenuation properties of the lower mantle. There are many possible methods for estimating the attenuation of the poloidal f i e l d , the most useful involving s t a t i s t i c a l considerations of the magnetic fields in the core. The information required to recover the el e c t r i c a l conductivity distribution, the amplitude ratio of two transfer functions, i s recoverable from the spatial deconvolution of two autocorrelation functions. These functions of space and frequency are expressions of the statistics of the radial and transverse poloidal fields over the core-mantle boundary. The numerical complexity of the whole inversion, the indeterminacies in the power ratio, the approximations involved, a l l suggest that the application of'the inversion method w i l l be a lengthy computational effort. It i s , nevertheless, an application which should significantly increase our understanding of the el e c t r i c a l conductivity distribution of the earth's lower mantle. APPENDIX' A-l SOLUTION OF THE INTEGRAL EQUATION \u00E2\u0080\u00A2 The Integral equation (40) i s u f(u,v) + K(u,v) + | f (w,v)K(u,w)dw = 0 , v < u, (A-l-1) o where f(u,v) (equations (41) or (47)) may be written f(u,v) = \u00C2\u00A3 {2r.cos*.ucosft.v - 2s.cosy.ucosy.v} , (A-1-2) js=]_ ^ ^ J J 3 J where Yj = ( 2 j - l ) 2 L or Yj \" CJ-1)T , j=l,2,\u00C2\u00AB-' , with the condition that s^ = i f = 0. As the difference between the characteristic values and weights and their asymptotic values can be made as small as desired, the i n f i n i t e series for f(u,v) is truncated, differentiated and integrated with j u s t i f i e d abandon. The result of substituting the truncated form of f(u,v) (equation (A-1-2)) into the integral equation (A-l-1), i s -89-N , fCu,v) + K.(u,v) + \u00E2\u0080\u00A2\u00C2\u00A3 \u00C2\u00A32r^cosAHy cos& wK(u,w)dw J o -2s^cosYjV j cosYjWK(u,w)dw} = 0. If the symbols u Cj(u) = i\ I cosJljWKCu,w)dw o u dj(u) = -s_j J cosYjWK(u,w)dw , j=l,2,**'N , are used, K(u,v) Is given by N K(u,v) = -f(u,v) - I {2cos\u00C2\u00A3.vc.(u) + 2cosY.vd.(u)} . (A-1-3) j=l 1 3 1 J To determine the c^ . (u), equation (A-l-3) is multiplied by r^cos\u00C2\u00A3/vdv and integrated from o to u. The result i s the set of linear equations c.(u) N - 4 \u00E2\u0080\u0094 + I {c.(u)a,.(u) + d.(u)b..(u)} i j=l cosYjVf (u,v)dv , i=l,2,\"\u00C2\u00ABN . (A-1-4) The multiplication of equation (A-l-3) by s^cosy/vdv and integration from o to u, gives the other N equations -90-d Cu) N - T + ' I {c , (u )b , ,Cu) + d . (u ) e f , (u ) } i 1=1 3 . . u cosy^fCu,v)dv , i=l,2,-\u00C2\u00AB\u00C2\u00ABN , (A-1-5) where u a. . (u) = 2 cosft.vcosfl, .vdv i l J i J o u b^j Cu) = 2j cosJl/vcosy.. vdv , o u V - 1 cosy^vcosYjVdv Hence, the c^(u) and d^ Cu) , for a given value of u, are found through the solution of 2N simultaneous linear equations. Substitution into equation (A-1-3) , gives K(u,u) which may, through interpolation be numerically differentiated for dK(u,u)/du. The alternative, which is used here, (see the computer l i s t i n g for the inversion in Appendix A-4), is to calculate dK(u,u)/du by differentiating equation (A-1-3). The result i s dK(u,u) = _ df(u,u) du du N + I {2i!..sinJ!-.uc. (u) .+ 2y.siny.ud.(u)} J 3 3 \u00E2\u0080\u00A2 J I I 3 3 N 11 ' - I {2cosA.ue,(u) + 2cosy.uf.(u)} , (A-1-6) j=l J J. J J -91 -where An operation s i m i l a r to that used to determine the c,. Cu) and the d,. (u), Is used to generate the 2N simultaneous equations e. Cu) ' N ~ \u00E2\u0080\u0094 + T {a (u)e.Cu) + b. .(u)f .(u)} r ^ j = l u = cosJl .uf Cu.u) - I cos&.v ^^ U ? V^dv x J x 3u o N da db I {c.Cu) -r^1 + d.Cu) - r ^ } j\u00E2\u0080\u009E! J du 3 du f. Cu) N \u00E2\u0080\u00A2+ \u00C2\u00A3 {e,Cu)b4.(u) + f .Cu)B, .Cu)}; u \" \ ( 9f(u,v) , - cosy^uf Cu,u) - | C O S Y / V \u00E2\u0080\u0094 ^ ' ' dv o N db dB . 3=1 J J The recovered values f o r e.Cu) and f . (u), along with c.Cu) and 3 J , 3 d.Cu) f o r j = l , 2 , * \" N, are used i n equation (A-1-6) to determine J .i dK(u,u)/du f o r a given value of u, o < u < 1. APPENDIX A-2 SPATIAL CONVOLUTION OF POLOIDAL FIELD COMPONENTS An a r b i t r a r y f unction of p o s i t i o n may be expanded i n a s e r i e s of s p h e r i c a l harmonics as oo n S C r ^ e , * , t ) = I I s mCr t)P mCcose)e n=l m=-n im (A-2-1) where the T^tcosd) are the associated Legendre functions as defined on n the second page of Chapter I I . The c o e f f i c i e n t s s m ( r , t ) are given by IT ir s m C r 1 t t ) = C-D\"(n+l/2) ( ( g ( r e,<(),t)p- mCcose)e- : L m < f >sinede. (A-2-2) n \u00C2\u00B1 Zir j j l n -ir o m. The c o e f f i c i e n t s s ( r . t ) are assumed to be the product of two other n c o e f f i c i e n t s m, * , , . m, . s ( r ,t) = h (r , t ) g (r ,t) n 1 n z n 2 (A-2-3) (n+1/2) m, Both h ( r - , t ) and g (r\u00E2\u0080\u009E,t) are considered as the c o e f f i c i e n t s of two n Z n z harmonic expansions of the form H(r 6,t) = I h (r t)P (cos6) z _ n Z n n=l oo n G(r 9,6,(|),t) = I I g m ( r 2 , t ) P m ( c o s e ) e lm<$> n=l m=-n (A-2-4) - 9 3 -where the coefficients h (r\u00E2\u0080\u009E,t) are related to H(r o,0,t) by n 2 z Tf 1 f h n ( r 2 , t ) = (n+ -) H(r 2 >8,t)P n(cose)sin0de . (A-2-5) The relationship between the three harmonic coefficients, s\u00E2\u0084\u00A2 , n h^ , and g\u00E2\u0084\u00A2 , (equation (A-2-3)), is now used to relate the function 8(^,6,9,0 (equation (A-2-1)), to H(r 2,8,t) and G(r ,9,9,t) (equations (A-2-4)). Substituting the RHS of equation (A-2-5), and a similar one to (A-2-2) for g\u00E2\u0084\u00A2 , into equation (A-2-1), changing the order of integration and summation, gives TT Tf S ( r l f e , 9 , t ) - j-ri | H(r2,8',t)d8'}-{ | j G(r 2,8\",9\",t)sin8\"de\"d9\"} 0 - Tf 0 *{ I (n+-i)sin9'P (cos8')P (cos8 )} , (A-2-6) , 2 . n n o n=l where the addition theorem of the Associated Legendre functions, I: P (cos8 ) = I (-l) mP m(cose)P T O(cos8\")e i m (*- < t ,\" ) n o \" n n m=-n where C O S 0 Q = cos8cos8\" + sin6sin6\"cos(9-9\") , has been used td reduce the equations to the form of (A-2-6). The sum, over the harmonic degree n, in equation (A-2-6), is the expansion, in Legendre functions, of the Dirac delta function, 6(6-0').\u00E2\u0080\u00A2' Therefore, equation (A-2-6) may be further reduced to the desired convolution -94-Tf IT S(r1>0,ct),t) - \u00E2\u0080\u0094 | | GO:2,O,,>(l>\",t)H(.r2,6o>t)sinO,,d0,,d(j),, , (A-2-7) TT 0 for two arbitrary functions , S(r^, 0, cj), t) and G(r 2,G, <$>, t) , defined, at different r a d i i , over a spherical surface. The number of possible applications of this convolution to the vector components of the poloidal f i e l d , or their autocorrelations, at the earth's surface and at the core-mantle boundary, are almost endless. A particular case, that of the convolution of the autocorrelations of the transverse and radial poloidal fields at the core-mantle boundary, r=b, is given here. As the power ratio of two transfer functions i s the aim of the convolution, the convolution in time, a linear operator, i s omitted. As the Fourier transforms of the transverse poloidal f i e l d components, Pg or P^ , may be expressed as PQ(r,e,0,io>) - f g d I i| 7{rE^(r,ico)}P^(cos8)e i^} , n=l m=-n 3 = -^Q- f(r,e,,ico) V ^ 9 ' ^ - ^ - ^ f ( r > 6 > * > i u ) the convolution of the autocorrelation of f (b, 0, ,ico) and P (b , 0 , ,ico) , is more convenient. The spatial autocorrelation of the function f (b , 6,,0'\u00C2\u00BB'.i\") = E{f(b,0,,iu)f*(b,e' ,<(>' ,iw)} -95-where E indicates the expected value. The assumption that there i s zero probability of any spatial correlation between different harmonic coefficients of a given f i e l d component, allows the autocorrelation of the Fourier transform of the poloidal radial f i e l d , for example, to be expressed as $ p p (b,e,9,e',9',iu>) r r 0 0 n 2 2 r r n (n+1) c/_m. . m* . . -.^m. Q. -m. QiN im( \u00E2\u0080\u009E E{Z (b,io>)\u00C2\u00A3 (b,iw)}P (cos0)P (cos\u00C2\u00A9 )e Y Y n n n n n=l m=-n b In the.scheme of equation (A-2-2), the two autocorrelations, and \u00C2\u00A7 p , are related through the coefficients of their respective r r spherical harmonic expansions. If m 1 , 9 , \u00E2\u0080\u009Em, . . , 3 , vm. , .-, * lT.-m. ... -ima' s = \u00E2\u0080\u0094\u00E2\u0080\u009EE{ \u00E2\u0080\u0094 { r Z (r,io))} ,-r-{rZ (r,ia>)}*}P (cos0')e Y , n , 2 3r n r=b3r n r=b n b 2 2 * m n (n+1) . , . \u00E2\u0080\u009Em ., . . , -m, A l N -imd> g = \ E{Z (b,iu)Z (b,iw)}P (cos6')e y n , z n n n b m 1 * s (n+ f) H (b,oi)H (b.aj) , n , ., l . Z n n h = \u00E2\u0080\u0094 (n+ \u00E2\u0080\u0094) = : , g\" (n+l) Z G (b)tl.)G (b.o)) ; n n n where Hn(b,oj) and Gn(b,w) are the radial and transverse transfer functions, the convolution of the autocorrelations is given by $f(b,e,9,6,,' ,iu) rr TT 1 _ f f 2TT -Tf O <5p (b,0\",9,l,e' ,9* ,iu)H(6 ,u)sine\"d0\"d9\" , r r where - 9 6 -H ( b , c o ) H * ( b , u ) . / \u00E2\u0080\u0094 = (n+1) H(6 ,OJ)P Ccose)sin6d0 G n ( b , u ) G ( b , u ) . J n n , n o APPENDIX A-3 PROGRAM LISTING: GENERATION OF A SPECTRAL FUNCTION The input to this program i s a specified e l e c t r i c a l conductivity model (given by the constants SGC, SGB, CC, DC), the harmonic degree (NORD), and the boundary conditions, given by 8' and 8 (the constants BP and BO respectively), on the eigenfunctions. A fourth order Runge-Kutta integration scheme i s used to solve the boundary value problem (equations (12)) for a given value of the parameter A (EIGEN). The characteristic values, the zeros of 8'M(b,A) + 8 \u00E2\u0080\u00A2^^(rjA) k > a r e f Q u n d through iterative adjustment of the parameter A , using the method of false position. The Runge-Kutta algorithm used, in SUBROUTINE DRK, is l i s t e d in Appendix A -4 . - 9 8 -C M A I N PROGRAM TO I N T E G R A T E A G I V E N S T U R M - L I 0 0 V I L L E P R O B L E M C TO G E N E R A T E T H E CHAR A C T E R I S T I C V A L U E S AND W E I G H T S , I . E . T H E C S P E C T R A L F U N C T I O N . N E C E S S A R Y I N T E G E R P A R A M E T E R S A R E : C NORD = HARMONIC D E G R E E C NS = ORDER OF F I R S T TERM TO B E D E T E R M I N E D C N-NS+1 = TOTAL NUMBER OF S P E C T R A L F U N C T I O N TERMS G E N E R A T E D C M = NUMBER OF S T E P S I N RUNGE-KUTTA I N T E G R A T I O N C P H Y S I C A L P A R A M E T E R S N E C E S S A R Y A R E : C B = R A D I U S OF C O R E - M A N T L E BOUNDARY C C = R A D I U S OF UPPER-LOWER MANTLE BOUNDARY C SGMRS = MEAN ROOT SQUARED C O N D U C T I V I T Y ( O B T A I N E D C BEFOREHAND FROM I N T E G R A T I O N OF MODEL) C S G C , S G B , C C , D C = P A R A M E T E R S OF C O N D U C T I V I T Y MODEL C BP,BO = BOUNDARY C O N D I T I O N C O N S T A N T S , B E T A ( P R I M E ) AND C B E T A , ON THE E I G E N F U N C T I O N S AT R = B C N.B. THE C O N D U C T I V I T Y MODEL I S G E N E R A T E D I N T H E M A I N PROGRAM C AND STORED I N COMMON FOR AUXRK. A L L A R I T H M E T I C I S I U C DOUBLE P R E C I S I O N C I M P L I C I T R E A L * 8 { A - H , 0 - Z ) R E A L * 4 F L O A T COMMON S I G ( 4 0 0 1) , F A C T (4 0 0 1 ) , E I G E N, H, D3 , FUN , Z DH S T , C A N T , 1 N MR K NORD=3 NS=1 N=10 M=1000 H=0.001DO READ ( 6 , 7 7 5 ) S G M R S , I S G C , I S G B , I C C , I D C 7 7 5 FORMAT ( D 2 6 . 1 6 , 4 1 3 ) P I = 3 . 1 4 1 5 9 2 6 5 3 5 8 9 7 4 D 0 B C = D B L E ( F L O A T ( I S G B - I S G C ) ) / ( 1 . DO-DEXP (-2.D0) ) C C = D B L E ( F L O A T ( I C C ) ) / 1 . D2 D C = P I * D B L E ( F L O A T ( I D C ) ) BP=1.DO BO=0.DO U M B E R = 3 . 1 5 5 6 9 2 5 9 7 4 7 D 7 B = 3 4 8 8 . D 0 C = 5 8 7 1 . D 0 FUN=B-C C O N S T = 0 . 4 D 0 * P I * F U B * F U N / U M B E R C A N T = P I * 1 . D 4 * U M B E R / ( 1 6 . D0*FIJN*FUN* S 3 M R S ) FT= D B L E ( F L O A T ( W O R D ) ) AC=1.DO+BC D3=-FT*FUN/C S A S = F T * ( F T + 1 . D 0 ) D 1 = 2 . D 0 * C / F U N D2 = C * C / ( F U N * F U N ) WRITE ( 6 , 6 6 6 ) NS, N , NOR D, H ,BO -99-666 FORMAT {//5X, 'EIGENVALUES AND WEIGHTS FROM I=',I2, 1 1X,',TO I = \u00E2\u0080\u00A2 , 12,2 X,* FOR N=\u00E2\u0080\u00A2,12,1X,\u00E2\u0080\u00A2ST1PSIZE (RK) = ', 2 F8.5,2X, 'FH=',F8.5,/, 'SIGMA (R) = AC-BC*EXP (-2X } , 3 \u00E2\u0080\u00A2+CC*SIN (DC*X) ') WRITE (6,667) AC,BC,CC,DC 667 FOR MAT (1 OX, 1 AC =\u00E2\u0080\u00A2,1PD10.2,3X,\u00E2\u0080\u00A2BC =\u00E2\u0080\u00A2,1PD10.2,3X, 1 'CC =',1PD10.2,3X,'DC =',1PD10.2) J=1 X1 = 0.D0 Z=AC-BC*DEXP(-2.D0*X1)+CC*DSIN (DC*X1) SIG (J)=Z*CONST FACT (J) = SAS/(X1*X1+D1*X1 + D2) 739 X1=X1+0.5D0*H Z= AC-BC*DEXP(-2.D0*X1)+CC*DSIN (DC*X1) SIG (J + 1)=Z*CONST SIG (J + 2 ) = S I G ( J + 1) FACT (J+1) =SAS/ (X1*X1 + D1*X1 + D2) FACT (J + 2) =FACT (J + 1) X1=X1+0.5D0*H Z=AC-BC*DEXP (~2.D0*X1) +CC*DSIN (DC*X1) SIG (J+3)=Z*CONST SIG (J + 4)=SIG(J+3) FACT (J+3)=SAS/(X1*X1 + D1*X1+D2) FACT (J + 4 ) = FACT (J + 3) J=J + 4 IMK=4*M~3 IF ( J . GT. IMK) GO TO 740 GO TO 739 740 CONTINUE CALL EIG (NS,N,M,BP,BO) STOP END SUBROUTINE EIG (NS,N,M,BP,BO) C C SUBROUTINE TO GENERATE THE SPECTRAL FUNCTION OF A GIVEN C STURM-LIOUVILLE PROBLEM. RUNGE-KUTTA ALGORITHM INTEGRATES C THE I N I T I A L VALUE PROBLEM FOR A GIVEN VALUE OF THE PARAMETER C LAMBDA. FALSE POSITION ALGORITHM ADJUST LAMBDA (EIGEN) , TO C SATISFY BOUNDARY CONDITIONS AT R=B. WEIGHTS ARE CALCULATED C WITH TRAPEZOIDAL RULE. FIRST ESTIMATES FOR CHARACTERISTIC C VALUES COME FROM ASYMPTOTIC APPROXIMATIONS (ASEIG). C IMPLICIT REAL*8(A-H,O-Z) REAL*4 FLOAT -100-DIMENSION X(50,4) ,Y (3) ,F(3) ,Q (3) ,FNCSQ (1001) , 1 STORE (1001) ,SUM (10) ,VALUE (5) COMMON SIG(4001) , FACT (4001) , EIGEN, H , D3 , FUN , CONST, CANT , 1 NMRK MP1=M+1 IFORM=1 IF (BO.LT. 1.D-6) IFORM=2 FORM=DBLE (FLOAT (IFORM) ) DO 99 INDEX=NS,N BAG= DBLE (FLOAT(IFORM*INDEX-1) ) /FORM ASEIG=4.D0*BAG*BAG *CANT FTR= 0.DO IF (IFORM.EQ. LAND-INDEX. EQ.1) FTR = 0. 1 DO IF (INDEX.GE. 4) GO TO 98 X (INDEX, 1) =0.9 D0*ASEIG X(INDEX,2)=1.1 DO* AS EIG+ FTR GO TO 99 98 X (INDEX, 1) =ASEIG-0. 2 DO X (INDEX,2) =ASEIG + 0.2D0 99 CONTINUE DO 3 I=NS,N DO 3 K= 1,2 K1=K+2 EIGEN=X (I,K) Y(1) = 0.DO Y (2) = 1, DO Y (3) =D3 NMRK=0 CALL DRK(Y,F,Q,H,3,M) X (I,K1)=BP*Y (2) + B0*Y (3) /FUN 3 CONTINUE DO 11 I=NS,N WRITE (6,613) 613 FORMAT (//18X,'EIGENVALUE 1,10X, \u00E2\u0080\u00A2HM+DM/DR\u00E2\u0080\u00A2,5X, 1 'DERIVATIVE WRT LAMB DA', 5X , 'DERIVAIIV E W RT X') WRITE (6,5) I 5 FORMAT(5X,\u00E2\u0080\u00A21=1,12) L=1 7 CHI= (X (1,3) *X (I,2)-X (1,1) *X (1,4) ) / (X (1,3) -X (1,4) ) Y ( 1 ) = 0 . D 0 Y (2) = 1 . D 0 Y (3) =D3 EIGEN=CHI NMRK=0 DO 6 J=1,M CALL DRK (Y,F,Q,H,3, 1) FNCSQ (J+1) =Y (2) 6 STORE (J+ 1) =Y (1) FUNC=BP*Y(2) +B0*Y (3)/FUN - 101-DERIV=Y (3) CHI=EIGEN IF (FUNC) 13,77,12 13 SLOPE= (X (1,3) -FUNC) / (X (1,1) -CHI) GO TO 77 12 SLOPE=(FUNC-X ( 1 , 1 ))/(CHI-X ( 1,2) ) GO TO 77 77 CONTINUE WRITE (6,4) CHI, FU NC, SLOPE, DER IV 4 FORMAT(10X,1P4E20.9) IF (DABS (FUNC) . LE. 1 0. D- 10) GO TO 15 IF (L.GE.10) GO TO 14 IF (FUNC) 25,11,26 C C PREVIOUS THREE STATEMENTS ARE KICK-OUTS FROM FALSE C POSITION - RK - CHARACTERISTIC VALUE GENERATOR. C 25 X(I,2)=CHI X(I,4)=FUNC L=L+1 GO TO 7 26 X(I,1)=CHI X (I,3)=FUNC L=L+1 GO TO 7 14 WRITE(6,83) 83 FORMAT (//'TOO SLOW CONVERGENCE'/) 15 FNCSQ (1) =1 .DO STORE (1) =0.D0 DO 16 J=1,MP1 IL=4*J-3 16 FNCSQ (J) =FNCSQ (J) *FNCSQ (J) *1.D3*SIG (IL) /CONST NO=M/8 DO 97 IP=1,4 JN=2** (4-IP) SUM (JN) =0.D0 NIP=JN+1 DO 96 INDEX=NIP,M,JN 96 SUM (JN)=SUM (JN) +FNCS Q (INDEX) VALUE (IP) =FUN* (FNCSQ (1) /2.D0+FNCSQ (MP1) /2. DO +SUM (JN) ) / 1 DBLE (FLOAT (NO) ) 97 NO=2*NO WRITE (6,81) FNCSQ (MP 1) 81 FORMAT(/5X,'COMPOSITE T RULE, SIGMA*M*M=',1PE15.8/) WRITE (6,82) (VALUE (IG) ,IG= 1, 4) 82 FORMAT (10X,1PE20.9) 11 CONTINUE RETURN END -102-SUBROOTINE AUXRK (Y, F) C C SUBROUTINE REQUIRED BY SUBROUTINE DRK (A PUBLIC PROGRAM C AT THE UNIVERSITY OF BRITISH COLUMBIA COMPUTING CENTER), C A DOUBLE PRECISION FOURTH ORDER RUNGE-KUTTA ALGORITHM. C IMPLICIT REAL*8{A-H, O-Z) DIMENSION Y (3) ,F (3) COMMON SIG(4001) , FACT (4001) ,E IG EN,H,D3,FUN,CONST,CANT, 1 NMRK F (2)=Y (3) NMRK=NMRK+1 F (3) = (FACT (NMRK) -SIG (NMRK) *EIGEN) *Y (2) RETURN END APPENDIX A-4 PROGRAM LISTING: THE INVERSION The input to t h i s program i s the s p e c t r a l function, that i s the c h a r a c t e r i s t i c value sequence, X_. , j = l , 2 , N, and the sequence of weights, Py j = l , 2 , N . The harmonic degree, the mean root squared e l e c t r i c a l conductivity, a m r s \u00C2\u00BB and the type of asymptotic approximation to the c h a r a c t e r i s t i c values, whether ( 2 j - l ) - j or (J-I )T T, are also required f o r the in v e r s i o n . The program i s designed to carry out the in v e r s i o n , i f desired, f o r two s p e c t r a l functions and to average the two recovered d i s t r i b u t i o n s . For e i t h e r recovery, dK(u,u)/du i s evaluated, i n SUBROUTINE DKDU, at selected values of u, o < u < 1, chosen by SUBROUTINE SPLINE. SPLINE interpolates dK(u,u)/du to other values of u, values c a l l e d f o r by the main program, f o r use i n the Runge-Kutta i n t e g r a t i o n of the boundary value problem i n the e l e c t r i c a l conductivity. The value of ^'^ r^ | i s i t e r a t i v e l y adjusted, through the f a l s e p o s i t i o n algorithm, to f i n d the zero of (b - rl u_^)\u00C2\u00AB The systems of simultaneous l i n e a r equations are solved by Gaussian el i m i n a t i o n and double back s u b s t i t u t i o n . The subroutines LRD (double p r e c i s i o n , t r i a n g u l a r decomposition of the c o e f f i c i e n t matrix i n DKDU), and DBS (double p r e c i s i o n , double back s u b s t i t u t i o n i n DKDU) are not l i s t e d . They are p u b l i c programs supplied by the un i v e r s i t y computing center. -104-C THIS IS THE MAIN PROGRAM FOR THE INVERSION FROM SPECTRAL C FUNCTIONS TO CONDUCTIVITY DISTRIBUTION. I T IS DESIGNED TO C PRODUCE AN AVERAGE PROFILE FROM TWO RECOVERIES (TWO SPECTRAL C FUNCTIONS). INTEGER PARAMETERS NEEDED ARE: C N = NUMBER OF TERMS IN SPECTRAL FUNCTION C NORD = HARMONIC DEGREE C NIX = NUMBER OF STEPS IN RUNGE-KUTTA ALGORITHM C NSPL = NUMBER OF EVLUATION POINTS OF DK(U,U)/DU C IFOR M = PARAMETER SPECIFYING THE TYPE OF ASYMPTOTIC C CHARACTERISTIC VALUES, = 1, OR = 2. C PHYSICAL PARAMETERS HEEDED ARE: C B1 = RADIUS OF CORE-MANTLE BOUNDARY C C1 = RADIUS OF UPPER-LOWER MANTLE BOUNDARY C SGMRS = MEAN ROOT SQUARED CONDUCTIVITY C SGC = ELECTRICAL CONDUCTIVITY AT R = C1. C N.B. MOST ARITHMETIC IS IN DOUBLE PRECISION. C IMPLICIT REAL*8(A-H,O-Z) REAL*4 FLOAT DIMENSION Y (4) ,F (4) ,Q (4) ,EIG (50) , RATE (50) , FACT (50) , BIG (50) DIMENSION Y1P (201) , Y2P (201) ,X1P (201) ,X3P (20 1) ,U3 (201) DIMENSION A A (40,40) , BA (4 0,40) ,CA (4 0, 1) , DA (40,1) ,EA (2 0,20) , 1 FA (2 0,20) ,GA (20,2 0) , AB (20) , BB (2 0) ,CB(20) , DB (20) ,EB (20) , 2 GB (2 0) ,HB (20) ,OB (20) ,BIT (50) ,U4 (201) , FB (20) , IB (80) DIMENSION ZA (201) ,ZB (201) ,X (4) COMMON EIG,FACT,BIG,BIT,BIN,CAKE,IFORM CALL PLOTS DATA NORD,B1,C1/3,3488.00,5871.DO/ M=2 0 NIX=2 01 MD=2*M MQ=4*M NSPL=50 FUN=B1-C1 PI=3. 14 159265358974D0 UMBER=3.15569259747D7 C C UMBER IS THE NUMBER OF SECONDS PER YEAR. ISGC, ISGB, C ICC, IDC, IPOW, ARE INTEGER FORMS OF THE MODEL CONDUCTIVITY C PARAMETERS C READ (6,775) SGMRS,ISGC,ISGB,ICC,IDC,IPOW 775 FORMAT (D26.16,513) SGB= DBLE (FLOAT (ISGB) ) SGC=DBLE (FLOAT (ISGC) ) DIG=SGB-SGC BC=DIG/ (1. DO-DEXP (-2. DO) ) AC=SGC+BC -105-CC=DBLE (FLOAT (ICC) ) / (10. **IPOW) DG=PI*DBLE (FLOAT (IDC) ) SGC=SGC*1.D3 SGB=SGB*1.D3 SQMRS=DSQRT(SGMRS) SQSGB=DSQRT(SGB) SQSGC=DSQRT(SGC) QUSGC=DSQRT(SQSGC) FAD=-FUN*SQMRS*DSQRT (4.D0*PI/ (1.D4*0MB ER)) GAD=1.D3*SQMRS/(4.D0*SGC*QUSGC) ALPHA=FUN*SQSGC*SQMRS/2.DO BIN= DBLE (FLOAT (NOR D) ) BIN=BIN*(BIN +1.DO)*FUN*FUN*SGMRS CAKE=FUN*SQMRS WRITE (6,12) NORD, M 12 FORMAT(/5X,'N= ',I2,5X,'NO. OF TERMS= *,I3,5X,'SIGMA (R) 1 \u00E2\u0080\u00A2 = ' , 1X, \u00C2\u00BBAC-BC*EXP (-2.*X) +CC*DSIN (DC*X) . \u00E2\u0080\u00A2/) WRITE(6,667) AC, BC, CC, DC ,SGMRS 667 FORMAT (10X,'AC =',1PD 10,2,3X, \u00E2\u0080\u00A2BC =', 1PD10.2,3X, 1 'CC = ' , 1PD10 . 2, 3 X, 1 DC =\u00C2\u00AB,1PD10.2,3X,'SGMRS =',1PD10.2) DO 917 J1=1,2 C C THE ALLOWED VALUE OF J1 DETERMINE THE NUMBER OF C RECOVERIES TO BE MADE. J1 = 1, OR J1 = 2. J1 = 2 MEANS C THE RESULT WILL BE AN AVERAGE PROFILE. C READ (6,916) IFORM 916 FORMAT (12) FORM=DBLE (FLOAT (IFORM) ) READ (5,11) (EIG (I) , FACT (I) ,1=1 ,M) C C EIG (I) AND FACT(I) REPRESENT THE SPECTRAL FUNCTION. C EIG {I) CORRESPONDS TO LAMBDA (I) AND FACr(I) CORRESPONDS C TO RHO(I). (SEE APPENDIX A-1). C 11 FORMAT (2D26.16) DO 13 1=1,M EIG (I) =FAD*DSQRT (EIG (I) ) FACT (I) =ALPHA/FACT (I) BIG (I) =DBLE (FLOAT (IFORM*I-1) ) *PI/FORM RATE (I) =EIG(I) -BIG (I) BIT(I)=1.DO 13 CONTINUE BIT (1)=FORM/2.D0 C C EIG(I) NOW CONTAINS THE TRANSFORMED CHARACTERISTIC VALUES, C AND FACT (I) CONTAINS THE ALPHA/RHO (I) . BIT (I) CONTAINS C THE S ( I ) . (SEE APPENDIX A-1). -106-IF (IFORM.EQ. 1) BO = 1.D0 IF (IFORM.EQ. 2) BO=0. DO WRITE(6,777) BO 777 FORMAT(/20X,\u00E2\u0080\u00A2M PRIME IS MULTIPLIED BY \u00C2\u00BB,F8.3) WRITE (6,78) 78 FORMAT (//14X, ' L ( I ) - (2*1-1) *PI/2',6X, * A LP HA/RHO (I) '/) WRITE (6,79) (RATE (I) , FACT (I) ,1=1 , M) 79 FORMAT (7X,2F20.8) EPSLN=1.D-7 CALL SPLINE(M,MD,MQ,EPSLN,NSPL,AA,BA,CA,DA,EA,FA,GA, 1 A B,B B,CB, D B , E B, F B , G B , HB,OB,IB) WRITE (6,45) 45 FORMAT (//) H=0.005DO DERIVC=2,D0*BC+CC*DC C C NOW ENTERING A FALSE POSITION ALGORITHM TO ITERATE THE C BOUNDARY VALUE PROBLEM. D (SIGMA (R) )/DR IS ADJUSTED TO C MEET THE CONDITION Y (R) AT U=1 IS EQUAL TO B1. NOTE C THAT THE NOTATION DIFFERS FROM THE APPENDIX A-1. C IF(J1.EQ.1) DERIVC=11.D0 IF(J1.EQ.2) DERIVC=10.5D0 F1 1 = 0.01DO END=140.DO X (1) =DERIVC X (2) =DERIVC* (1 . D0-F1 1) C C F11 IS THE RELATIVE SISE OF THE FIRST STEP IN THE C DERIVATIVE VS. (B1-RADIUS) ITERATIONS. END IS C A MAXIMUM CUTOFF FOR PLOTTING PURPOSES. C NIP=NIX-1 DO 4 1 1= 1,2 DERIVC=X(I) Y (1)=0.D0 Y (2)=C1 Y (3) =QUSGC Y (4) =GAD*DERIVC CALL DRK(Y,F,Q,H,4,NIP) X (1 + 2)=B1-Y (2) IG=I-2 WRITE(6,44) IG,DERIVC,X (1 + 2) 41 CONTINUE L=1 42 DERIVC= (X (3) *X (2) - X (1) *X (4) ) / (X (3)-X (4) ) Y(1)=0.D0 Y (2) =C1 Y (3) =QUSGC -107-Y (4) =GAD*DERIVC U4 (1)=Y (1) Y1P (1)=SGC/1.D3 X1P(1)=Y{2) U3 (1)=DERIVC/FUN DO 43 1=2,NIX CALL DRK (Y(, F, Q, H, 4,1) Y1P (I) = (Y (3) **4) /1. D3 X1 P (I)=Y (2) U3 (I) =4.D0*Y1P (I) *Y (3) *Y (4) / (FUN*SQMRS) 04 (I)=Y (1) 43 CONTINUE DIFER=B1-Y (2) WRITE (6,44) L, DERIVC, DIFER 44 FORMAT (10X, \u00C2\u00AB ITERATION=\u00C2\u00BB,12,5X, 'DERIVC= 1PE1 5, 5, 1 5X,'DIFER=',1PE15.5) IF (DABS (DIFER) . LE. 1. DO) GO TO 59 IF (L. GE. 10) GO TO 58 IF (DIFER) 46,59,48 C C THE PREVIOUS THREE STATEMENTS ARE KiCK-OuTS FROM THE C BOUNDARY VALUE PROBLEM IN CONDUCTIVITY. C 46 X(2)=DERIVC X (4) =DIFER L-L+1 GO TO 42 48 X(1)=DERIVC X (3) =DIFER L=L+1 GO TO 42 58 WRITE (6,50) 50 FORMAT(//'TOO SLOW CONVERGENCE'/) 59 CONTINUE DO 51 1=1,NIX X3P (I) =-X1P (I) +C1 FAM=-X3P (I) /FUN Y2P (I) =AC-BC*DEXP (-2. D0*FAM) +CC*DSIN (DC*FAM) 51 CONTINUE WRITE (6,4) 4 FORMAT (1H1) WRITE (6,1) 1 FORMAT(//20X,\u00E2\u0080\u00A2U',15X,'R',10X, 'SIGMA (R) \u00E2\u0080\u00A2 ,5X, 1 \u00E2\u0080\u00A2D/DR(SIGMA(R))',3X,'SIGMA EXACT'/) WRITE (6,2) ( (U4 (I) ,X1P(I) ,Y1P (I) ,03 (I) ,Y2P (I) ), 1=1,NIX, 5) WRITE (6,4) CALL SUMFUN(Y1P,X3P,Y2P,NIX,END,J1) IF (J1.EQ.2) GO TO 714 DO 713 J=1,NIX -108-ZA (J) =C1-X1P (J) 713 Z B ( J ) = Y 1 P ( J ) 917 CONTINUE IF (J1.EQ. 1) GO TO 888 714 CONTINUE CALL SPINE(NIX,EPSLN,ZA,ZB) C C ONE OF THE TWO POSSIBLE ( I F . J 1 . E Q . 2 ) , RECOVERED C PROFILES, IS FED INTO A SPLINE INTERPOLATION C SUBROUTINE, FOR SUPERPOSITION. C 2 FORMAT (10X,5F15.4) DO 712 J=1,NIX RAD=C1-X1P (J) CALL SPINV(RAD,FAT) Y1P (J) = (Y1P (J) +FAT) /2. DO 712 CONTINUE WRITE (6,1) WRITE (6,2) ( (U4 (I) ,X1P (I) ,Y1P (I) , U3 (I) ,Y2P (I) ) , 1=1, NIK, CALL SUMFUN(Y1P,X3P,Y2P,NIX,END,J1) C C SUMFUN CONTAINS THE PLOTTING SUBROUTINES. C WRITE (6,4) 888 CONTINUE CALL PLOIND STOP END SUBROUTINE AUXRK (Y, F) C C AUXRK IS THE SUBROUTINE REQUIRED BY DRK (A PUBLIC C F I L E SUPPLIED BY THE UNIVERSITY OF BRITISH COLUMBIA C COMPUTING CENTER), A FOURTH ORDER VERSION, WITH GILLS C COEFFICIENTS, OF THE RUNGE-KUTTA ALGORITHM. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y (4) ,F (4) , EIG (50) , FACT (50) ,BIG (50) ,BIT (50) COMMON EIG,FACT,BIG,BIT,BIN,CAKE,IFORM F (2) =CAKE/(Y (3) *Y (3) ) F(3) = Y ( 4 ) U=Y (1) CALL SPLINV(U,DK) F (4) =-BIN/(Y (3) *Y (3) *Y (3) *Y (2) *Y (2) ) +2. D0*DK*Y (3) RETURN END -109-SUBROUTINE SPLINE (MI,HD,MQ,EPSLN,N4,AA,BA,CA,DA,EA,FA, C C SPLINE INTERPOLATOR FOR THE FUNCTION DK(U,U)/DU. C SERVES TO CALL THE SIMULTANEOUS LINEAR EQUATION C SOLVER, SUBROUTINE DKDU, AT N4 EQUALLY SPACED C POINTS IN THE VARIABLE U. NOTE THAT DK(U,U)/DU IS C NOT EVALUATED AT U= 1 , BUT INTERPOLATED BY SPLINE. C 1 GA,AB,BB,CB,DB,EB,FB,GB,HB,OB,IB) IMPLICIT REAL*8 (A-H,0-Z) REAL*4 FLOAT DIMENSION EIG(50),FACT (50) ,BIG (50),BIT (50) DIMENSION T(1) ,SS (1) ,SS1 (1) ,SS2 (1) DIMENSION AA (MD, MD) , BA (MD,MD) ,CA (MD, 1) , DA (MD, 1) , 1 FA (MI,MI) ,GA (MI, MI) ,AB (HI) ,BB (MI) ,CB(MI) , DB (MI) , 2 GB (MI) , HB (MI) ,OB (MI) , EA (MI, MI) , EB (MI) , FB (MI) , IB (MQ) DIMENSION X(51),Y (51) , DELI (51) ,DELSQY(51) ,S2 (51) , 1 S3 (5 1) ,C (51) , B (51) ,H (51) , H2 (51) COMMON EIG,FACT,EIG,BIT,BIN,CAKE,IFORM WRITE (6 , 105) 105 FORMAT(//5X,'INPUT VALUES TO SPLINE'/) M=0 N=N4 N1 = N- 1 STEP=1.D0/DBLE (FLOAT (N) ) DO 1 1=1,N X(I)= DBLE(FLOAT (I))*STEP-STEP U=X (I) CALL DKDU(U,DK,MI,MD,MQ,AA,BA,CA,DA,EA,FA,GA,AB,BB, 1 CB,DB,EB,FB,GB,HB,OB,IB) 1 Y(I)=DK WRITE (6 , 104) (Y(I),I=1,N) WRITE(7,104) (Y(I),I=1,N) 104 FORMAT (1P5E20.9) C C VALUES OF DK(U,U)/DU ARE WRITTEN IN FILE C DEVICE 7, AS A HEDGE AGAINST FAILURE TO ACHIEVE AN C ACCURATE ENOUGH SOLUTION TO THE BOUNDARY VALUE C PROBLEM IN THE CONDUCTIVITY. C H(1)=X(2) -X(1) DELY (1) = (Y (2) -Y (1) )/H (1) 4 DO 52 1=2,N1 H (I) =X (1+1) -X (I) -110-H2 (I) =H (1-1) +H (I) B (I) =0.5D0*H (I-1)/H2 (I) DELY (I)= (Y (1 + 1)-Y (I) )/H(I) DELSQY (I) = (DELY (I)-DELY (1-1) ) /H2 (I) S2 (I) -2. DQ*DELSQY (I) 52 C(I)=3.DO*DELSQY (I) S2 (1)=S2 (2) S2 (N) =S2 (N-1) OMEGA=1.071797D0 IG=0 5 ETA=0. IG=IG+1 6 DO 10 1=2,S1 7 W= (C (I)-B (I) *S2 (1-1)- (0.5D0-B (I) ) *S2 (1 + 1) -S2 (I) ) *OMEGA 8 IF (DABS (W) , LE. ETA) GOTO 10 9 ETA=DABS (W) 10 S2 (I)=S2 (I) +W IF (IG. GE.N) GO TO 101 13 IF (ETA. GE. EPSLN) GO TO 5 GO TO 103 101 WRITE(6,102) 102 FORMAT(/5X,'QUESTIONABLE CONVERGENCE'/) 103 CONTINUE 14 DO 53 1=1,N1 53 S3 (I)= (S2 (1 + 1)-S2 (I))/H(I) AR=DELY (N1) +H (N1) *S2 (N1) /6.D0 AL= DELY (1)-H(1)*S2 (2)/6.D0 IF(M.LE.O) RETURN GO TO 15 ENTRY SPLINV(ARG,SP) M=1 T (1) =ARG 15 DO 6 1 J=1,M 16 1=1 54 IF (T (J)-X (1) ) 58,17,55 55 IF (T (J) -X (N) ) 57,59,580 56 IF (T (J) -X (I) ) 60,17,57 57 1=1+1 GO TO 56 58 SS (J)=AL*(T(J)-X (1) )+Y(1) SS 1 (J) =AL SS2 (J) =0.D0 GO TO 61 580 SS (J) =AR* (T(J)-X (N) )+Y (N) SS 1 (J)=AR SS2 (3) =0 . DO GO TO 61 59 I=N 60 1=1-1 - 1 1 1 -17 HT1 = T (J) -X (I) HT2=T (J) -X (1 + 1) PR0D=HT1*HT2 SS2 (J)=S2 (I) +HT1*S3 (I) DELSQS= (S2 (I) +S2 (1+1) +SS2 (J) ) /6. DO SS (J) = Y (I) + HT 1 *DELY (I) +PROD*DELSQS SS1 (J) = DEL Y(I) + (HT1+HT2)*DELSQS+PROD*S3(I)/6.D0 61 CONTINUE SP=SS (1) SP1=SS 1(1) SP2 = SS2 (1) RETURN END SUBROUTINE DKDU(U,DK,M,N,MQ,A,T,B,X,AH,BH,BP,SL,CL, C C THE BASIC SUBROUTINE TO CALCULATE DK(U,U)/DU FOR A C GIVEN U. SUBPROGRAMS SUPPLIED BY THE COMPUTING CENTER C ARE: LRD - TRIANGULAR DECOMPOSITION OF COEFFICIENT C MATRIX, STORED IN THE ARRAY A (N,N) , AND DBS - DOUBLE C BACK SUBSTITUTION. C 1 SJ,CJ,TIFF,CLF,GIF,SLI,SJB,IPERM) IMPLICIT REAL*8(A-H,O-Z) COMMON EIG,FACT,BIG,BIT,BIN,CAKE,IFORM DIMENSION EIG (50) , FACT (50) , B IG (50 ) , BIT (50) DIMENSION A (N,N) ,T (N,N) ,B (N, 1) ,X (N, 1) , AH (M,M) , BH (M,M) , 1 SL (M) ,CL (M) , SJ (M) , CJ (M) ,TIFF (M) , CLF (M) , GLF {M) , SLI (M) , 2 BP (M,M) ,SJB (M) , IPERM (MQ) MEND=M DO 996 1=1,MEND SL (I) =DSIN (EIG (I) *U) CL (I) =DCOS (EIG (I) *U) SJ (I) =DSIN (BIG (I) *U) CJ (I) =DCOS (BIG (I) *U) CLF (I) =CL (I) * FACT (I) GLF (I) =CJ (I) *BIT (I) SLI (I) =SL (I) *EIG (I) SJB (I) =SJ (I) *BIG (I) JUG=U* (BIG (I) -EIG (I) ) 996 TIFF (I)=JUG*JUG DO 3 J=1,M DO 3 1=1,M ID=I-J IF (ID) 4,5,4 -112-4 AH (I, J) = (SL {I) *CL (J)-CL (I) *SL (J) ) / (EIG (I) -EIG (J) ) 1 (SL (I) *CL (J) +CL (I) *SL (J) ) / (EIG (I) +EIG (J) ) BH (I,J) = (SL (I) *CJ (J)-CL (I) *SJ (J) )/ (EIG (I)-BIG (J) ) 1 (SL (I) *CJ (J) +CL (I) *S J (J) ) / (EIG (I) +BIG (J) ) BP (I,J) = (SJ(I) *CJ (J)-CJ (I) *SJ (J) ) / (BIG (I)-BIG (J) ) 1 (SJ (I) *CJ (J) +CJ (I) *SJ (J) ) / (BIG (I) +BIG (J) ) GO TO 3 5 CONTINUE 35 AH (I,J)=U+SL(I) *CL (I)/EIG (I) BH (I, J) = (SL (I) *CJ (J)-CL (I) *SJ (J) ) / (EIG (I) -BIG (J) ) 1 (SL (I) *CJ (J) +CL (I) *SJ (J) ) / (EIG (I) +BIG (J) ) IF (IFORM.EQ. 2) GO TO 995 IF (I.EQ.1) GO TO 34 995 BP (I,J)=U + S J ( I ) * C J (I)/BIG(I) GO TO 3 34 BP (I, J) =2. D0*U 3 CONTINUE DO 6 J=1,M DO 6 1=1,M IP=I-J IF (IP) 7,8,7 7 A (I, J)=AH (I, J) GO TO 6 8 A(I,J)=AH(I, J) +1. DO/FACT (J) 6 CONTINUE MP1 = M+1 DO 9 J=MP1,N DO 9 1=1,M K=J-M 9 A (I, J ) =BH (I, K) DO 10 J=1,M DO 10 I=MP1,N K=I-M 10 A (I, J)=BH (J,K) DO 11 J=MP1,N L=J-M DO 11 I=MPi1,N K=I-M IN=I-J IF(IN) 12,13,12 12 A (I, J) =BP(K,L) GO TO 11 13 A(I,J)=BP(K,L) -1.DO/BIT (K) 11 CONTINUE DO 14 1=1 , M B (1,1) =0. DO DO 14 K=1,M 14 B (I, 1)=B (I, 1) +GLF (K) *BH (I,K) -CLF (K)*AH (I,K) DO 16 I=MP1,N - 1 1 3 -L=I-M B(I,1)=0.D0 DO 16 K=1,M 16 B(I,1)=B (1,1) +GLF (K) *BP(L,K)-CLF (K) *BH (K,L) CALL LBD (N,N,N,A,IPERM,N,T) DO 47 1=1,N DO 4 7 J=1,N 47 A (I, J) =T (I, J) CALL DBS(N,1 , N, B, X, I PERM , N , T) SUM13=0.D0 SDM12=0.0D0 SUM5=0.0D0 DO 26 K=1,M KPM=K+M SUM12=SUM12+CL (K) *CLF (K) -CJ(K) *GLF (K) SUM13=SUM13+SLI (K) *X (K,1) + SJB (K) *X (KPM ,1) 26 SUM5=SUM5+X (K, 1) *CL (K) +X (KPM,1)*CJ (K) GOAT=2.D0*(SUM12+SUM5) DO 27 1=1,M SUM7=0.D0 SOM8=O.D0 DO 28 K=1,M SUM7=SDM7 + SLI (,K) * AH (I,K) *FACT (K) -SJB (K ) *B H (I, K) *BIT (K) 28 SOM8=SUM8 + SLI (K) *BH (K,I) *FACT (K) -SJB (K) *BP (I,K) *BIT (K) B (I, 1) =-CL (I) *GOAT+SUM7 27 B (I + M, 1) =-CJ (I) *GOAT + S0M8 CALL DBS (N,1,N,B,X,IPERM,N,A) SUM9=0.D0 SUM10=0.D0 BO 77 J=1,M JPM=J+M SUM9=SUM9+SLI (J) *CLF (J) -SJB (J) *CJ (J) *BIT (J) 77 SUM 1 0=SUM 10-CL (J) *X (J , 1) -C J (J) *X (JPM,1) DK=4.D0*SUM9+2.D0* (SUM 10+SUM1 3) RETURN END SUBROUTINE SPINE (N,EPSLN,X,Y) C C SUBROUTINE SPINE, SIMILAR TO SPLINE, INTERPOLATES C ONE OF THE TWO CONDUCTIVITY RECOVERIES FOR SUPERPOSITION. C IMPLICIT REAL*8 (A-H,0-Z) REAL*4 FLOAT DIMENSION EIG (50) , FACT (50) ,BIG (50) , BIT (50) -114-EIMENSION T(1) ,SS (1) rSS1 (1) ,SS2 (1) DIMENSION X (201) , I (201) , DELY (201 ) , DELSQY (201 ) , 1 S2 (201) , S3 (201) , C(201) ,B(201) ,H (201) ,H2(201) COMMON EIG,FACT,BIG,BIT,BIN,CAKE M=0 N1 = N- 1 H(1)=X(2)-X(1) DELY (1) = (Y (2) -Y (1) )/H (1) 4 DO 52 1=2,N1 H (I) =X (1+1) -X (I) H2 (I) =H (1-1) +H (I) B (I) =0.5D0*H(1-1) /H2 (I) DELY (I) = (Y (I+1)-Y (I) )/H (I) DELSQY (I) = (DELY (I) - DELY (1-1) ) /H2 (I) S2 (I) =2.D0*DELSQY (I) 52 C (I) =3.D0*DELSQY (I) S2 (1) =S2 (2) S2 (N)=S2 (N-1) OMEGA=1.071797D0 IG=0 5 ETA=0. IG=IG+1 6 DO 10 1=2,N1 7 W= (C (I)-B (I) *S2 (1-1)- (0. 5D0-B (I) ) *S2 (1+1) -S2 (I) ) *OMEGA 8 IF (DABS (W) .LE.ETA) GO TO 10 9 ETA= DABS (W) 10 S2 (I) =S2 (I) +W IF (IG.GE.N) GO TO 101 13 IF (ETA.GE.EPSLN) GO TO 5 GO TO 103 101 WRITE (6,102) 102 FORMAT(/5X,'QUESTIONABLE CONVERGENCE'/) 103 CONTINUE 14 DO 53 1=1,N1 53 S3 (I) = (S2 (1 + 1)-S2 (I) )/H (I) AR=DELY (N1) +H (N1) *S2 (N1) /6.D0 AL=DELY (1) -H (1) *S2 (2) /6. DO IF(M.LE.O) RETURN GO TO 15 ENTRY SPINV(ARG,SP) M=1 T (1) =ARG 15 DO 61 J= 1 , M 16 1=1 54 IF (T (J) \u00E2\u0080\u0094 X (1) ) 58,17,55 55 IF (T (J)-X (N) ) 57,59,580 56 IF (T (J)-X(I) ) 60,17,57 57 1=1+1 GO TO 56 - 1 1 5 -58 SS (J)=AL*(T(J)-X (1) )+Y (1) SS 1 (J)=AL SS2 (J)=O.DO GO TO 61 580 SS (J)=AR* (T(J) -X "Thesis/Dissertation"@en .
"10.14288/1.0052853"@en .
"eng"@en .
"Geophysics"@en .
"Vancouver : University of British Columbia Library"@en .
"University of British Columbia"@en .
"For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en .
"Graduate"@en .
"Inversion of the geomagnetic secular variation: uniqueness and feasibility"@en .
"Text"@en .
"http://hdl.handle.net/2429/32653"@en .
~~