. Nassau; ~r§=rrv-w :4»? 2.‘{..a.-.\\:\ m—u :x: «nu-gm.» 1‘? . 5-3.. .3".th $.13!!! "\.4 x a: :23: R. 11.2; ‘3'}; ' 4.5;“ 1"S'?’H‘* X; ’: 2.13 |- {5: "UV 213‘. k: -';.‘. I 4 a‘- - K. . ’ " 4" . LUIS w? W4: “VJ: “R ’ " ‘.’=\ 2. ..tK‘f; \‘4 1%.? 2%.! 55. k." . .2. w " $3 '3 1214st ’36 1:19 ”leg?” 3: 3... .. . 13;!" ‘ 7" f .1 $13,551: 91-. 5‘ .3':_2.fl‘t‘- ,9 I‘d. fink}.\r’ ki‘aig V'-\‘:!~'1'v ‘3.-X .o. P " r, m, gr) #2: t; 0 “ :3 ‘ (92!. h- , '.“‘fgfl:. o, 4-. “wad‘dn @3111“: a This is to certify that the thesis entitled A NUMERICAL METHDD FOR DETERMINING CHARACTERISTIC VALUES OF THE WAVE EQUATION presented by RICHARD CALVIN HAVEN-S has been accepted towards fulfillment of the requirements for I'IASTER OF SCIENCE degree in ELECTRICAL ENGINEmINC I I / .o / - / .‘ . 1» pl" _ ( ( ‘4. I .. 1. A.~‘-( -\ . 'v Major professor Date August 19, 1959 0-169 LIBRARY Michigan State University A NUMERICAL METHOD FOR DETERMINING CHARACTERISTIC VALUES OF THE WAVE EQUATION By RICHARD CALVIN HAVENS AN ABSTRACT Submitted to the College of Engineering Michigan State University of Agriculture and Applied Science in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Electrical Engineering . {322 Approved: 1 9 64/- MM) ear/A. ABSTRACT Methods of solving the two dimensional wave equation assuming wave guide boundary conditions are discussed. In particular, numerical methods are mainly stressed. Difference equations, relaxation methods, and iteration methods receive the most emphasis. The objective of the thesis is to find a method for deter- mining characteristic values. The result is an iteration program which will determine characteristic values for a uniform cross section wave guide. The method approximates the wave equation by difference equations. The eigenvalues are then computed by an iteration process. Alternative methods of solving the wave equation and methods of improving the existing program are given. The results of two applications are given as a check on the method. Finally, more difficult problems to which the same method might be applied are discussed. A NUMERICAL METHOD FOR DETERMINING CHARACTERISTIC VALUES OF THE WAVE EQUATION By RICHARD CALVIN HAVENS A THESIS Submitted to the College of Engineering Michigan State University of Agriculture and Applied Science in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Electrical Engineering 1959 ACKNOWLEDGEMENTS The author would like to sincerely thank Dr. L. W. Von Tersch for his assistance and would also like to extend the most genuine thanks to Mr. Richard J. Reid for his untiring help and guidance. ii TABLE OF CONTENTS ACKNO‘IILEDGEI‘ENTS . O O O O O O O O O O 0 Chapter I. INTRODUCTION . . . . . . . . . . A. B. Guided Waves . . . . . . . The Wave Equation . . . . II. DIFFERENCE EQUATIONS . . . . . . A. B. C. III. METHODS FOR SOLVING THE DIFFERENCE EQUATIONS A. B. C. First Order Difference Equations . Higher Order Difference Equations Special Type Difference Equations Matrix Methods . . . . . . Relaxation . . . . . . . . Iteration and Eigenvalue Methods IV. DIGITAL COMPUTER PROGRAMS . . . A. B. C. D. CONCLUSION . . Appendix An Iteration Program . . . Application of the Program Program Improvements . . . Results . . . . . . . . . I. THE ITERATION PROGRAM . . . . . II. PACKING ROUTINE . . . . . . . . III. INTERCONNECTION PROGRAM . . . . REFERENCES . . iii Page ii 45 53 56 61 CHAPTER I INTRODUCTION A. Guided Waves As is obvious from the name itself, the purpose of a wave guide is to transfer electromagnetic radiation or energy from one place to another by means of some sort of guiding structure. The design of the wave guide may be as varied as the applications to which it is applied. In many cases the wave guide is a length of hollow pipe which has as its physical characteristics the critical dimensions of its cross section, the conductivity of the metal structure, and the magnetic permeability, dielectric constant, and conductivity of the dielectric material in its interior. This type of wave guide is not the only type of interest, but it is this general type of wave guide which will be considered in this thesis. In other words, the wave guides considered in this thesis are wave guides for which: (1) the cross section is uniform, (2) the dielectric material in the wave guide has constant dielectric constant e, constant magnetic permeability/u, and zero conductivity, (3) the metal of which it is made is assumed to be a perfect conductor. This last stipulation is of no great practical concern since most wave guides are constructed out of metal with high enough conductivity so that the metal can be assumed to be a perfect conductor for all practical purposes. A wave guide has many modes of operation and, in fact, infinitely many. However, it is usually only the lower order modes which are of interest. Propagation of electromagnetic energy in the wave guide is governed by a homogeneous partial differential equation called the wave equation. This equation has solutions only for certain discrete values of one of its parameters. These discrete values are the characteristic values or eigenvalues of the wave equation. Each characteristic value corresponds to its own mode as was mentioned in the preceeding paragraph. Once the characteristic values are known, it is then a somewhat simpler matter to solve for the eigen- vectors of eigenfunctions which represent the field in the wave guide. As it turns out for a uniform cross section wave guide, it is only necessary to solve a two dimensional wave equation in order to determine the characteristic values. The dimensions involved are of course the X and Y dimensions of the cross section. Solution of the two dimensional wave equation then becomes a matter of 3 finding the characteristic values for the modes of interest. The wave equation is not unique in form but is found in problems involving many different branches of science. It so happens, that the wave equation is identically the same as the equation governing the vibrations of thin elastic membranes of arbitrary cross section, which are rigidly clamped on some or all of their edges. The vibrating membrane is brought up here because it has certain analogies which help understand the wave guide. The analogies will be explained after a few brief remarks on the wave guide. From the characteristic value of a mode of prOpagation in a wave guide, it is possible to calculate from classical formulas (l) the associated character- istics of the wave guide such as the cut-off frequency, the attenuation of the guide, the power handling capacity, and the mode separation. The latter is naturally dependent on more than one eigenvalue. Now for the analogies. The characteristic values associated with the vibrating membrane are related to the extreme positions in which it is possible for the membrane to vibrate. The lowest possible frequency at which the membrane will vibrate is analogous to the lowest cut-off frequency or cut-off frequency of the dominant mode of the wave guide. 4 The lowest frequency at which the membrane will vibrate in a more complex nature corresponds, clearly, to the cut-off frequencies of the higher order modes. The cut—off frequency of the wave guide, is the frequency below which prOpagation of energy in the wave guide is not possible for that mode. The rate at which the oscillations of the membrane are damped is analogous to the attenuation of the fields in the wave guide. The power handling capacity, which is determined by the maximum fields allowable before the dielectric material breaks down, is analogous to the extreme amplitude in which the membrane can vibrate without damage. Obvi- ously, mode separation is its own analogy. It is these properties which are of final interest, however, it is not the intent of this thesis to discuss these properties of the wave guide but rather the characteristic values or eigenvalues, from which these properties can be calculated. In order to analyze a wave guide of particular shape, it is first necessary to determine the eigenvalues. This is a very difficult job in all but the very simplest cases. The purpose of this thesis can now be stated. The object of this thesis is to develop a general program for a digital computer, in particular the MISTIC, which will, when given the necessary details about the boundary conditions, compute the desired characteristic values. B. The Wave Equation In the introduction the wave equation was mentioned as a homogeneous partial differential equation which governs the electromagnetic fields in a wave guide. The wave equation in general form is VP:- /«- g3}: or Viv-7% :1 (1-1) These equations are easily derived from Maxwell's equations. They are based on the assumption that the dielectric material has zero conductivity, is homo- geneous, linear, and isotrOpic. They are identical for the electric or magnetic field and must be satisfied by all components of the electric and magnetic field. If it is now assumed that the fields vary with hat time as eJ , the wave equations become VIE-v —-w:ae F/ 72/“: 'wfaefl (1-2) If it is further assumed that the propagation in the Z direction (the direction of prOpagation) varies as eF’2, where J is known as the prepagation constant, the equations then become: 1' 1. 1. "" 7’ 1- 1.. va E: ’(k {OJ/.Lé)£‘/ ng H: ’[’+w/‘¢)H (1-3) 6 In these equations ViyE is now the two dimensional Laplacian in the transverse plane; the third partial derivative having been set equal to r23 because of the (Z assumption that Ez varies as e- . Finally, it is convenient to replace 62 +tdgfltrby K3 and the equations are then, simply I- VHF—14:? [7,; //.—.- —/«‘// (1-4) ) c where K: =f2 +W2/‘6’ These equations have solutions only for discrete values of Ki. It is the Kc's with which the thesis is 2‘ 2" 2 4. =K=K = c1 02 c3 . . . . etc., where K0 is the eigenvalue of the dominant 1 mode and K0 is the characteristic value of the next 2 highest order mode, etc. Once the Kc is known the cut- concerned. These numbers are such that K off frequency fc, for example, is easily found by letting X'= 0 so that: K : ._ (5) .{i 2?W'W7¢e- It is common to classify solutions to the wave equation inside a wave guide into three general types. These modes of operation are called transverse electro- magnetic, transverse electric, and transverse magnetic; TEM, TE, and TM respectively. The TEM mode has no E or H component in the Z direction (i.e. direction of prOpa- gation) and may in many cases not even be a possible mode 7 of Operation. Its prepagation usually requires at least two separate conductors. The TE and TM modes are the usual modes of propagation in the wave guide. The TE modes are called transverse electric modes because they have no E component in the axial Z direction, and likewise, the TM has no H component in the Z direction. As was mentioned before, the wave equation must be satisfied by all of the individual components of the electric and magnetic field present in the problem being analyzed. For this reason, it is not necessary to solve the wave equation for the entire electric or magnetic field if the main interest is merely the characteristic values. One is at liberty to pick the component of E or H field which is most convenient in applying the boundary conditions. What can be done in practice is to decide what mode is to be analyzed and pick for solution the component of the electric or magnetic field, for example, Ex’ Ey’ Ez’ Hx’ Hy’ or Hz, which allows the easiest application of boundary conditions. Thus far the wave equation has been considered only in a rectangular coordinate system. This is because the detailed method given in the following chapters uses the rectangular coordinate system. However, it is just as logical to write the wave equation in a cylindrical or spherical coordinate system or any other coordinate system. A common method of solution of the two dimensional or the three dimensional wave equation in any coordinate system is known as separation of variables (1, p. 145) or as the product solution method. The solution is assumed to be a product of functions, each of which is a function of only one variable of the coordinate system. The assumed product solution is substituted into the differential equation, and it is then possible to separate it into ordinary differential equations which can be solved individually. Of course, there still remains the application of the boundary conditions. When this is done the individual solutions are multiplied together giving the final product solution. This is easily done when the geometry of the wave guide matches nicely the coordinate system. For instance, the solution to the wave equation for a retangular wave guide can be found quite easily by writing the wave equation in the rectangular coordinate system. The solution turns out to be trigonometric functions which are easily matched to the rectangular geometry of the wave guide. In the case of the cylindrical wave guide, Bessel functions result for the radial variation in the solution, having written the wave equation in cylindrical coordinates. Again, the Bessel functions are easily matched to the circular boundary. In a like manner it 9 is even possible to solve a parabolic wave guide (2) using a parabolic coordinate system. However, when the cross section of the wave guide does not fit the coordinate system, applying the boundary conditions becomes quite difficult. As a result, solution of the wave equation is very difficult except for these few special cases above. A great deal of work has been done during the past decade or so on analyzing the ridge wave guide. S. B. Cohn has shown (3) the ridge wave guide to have a lower cut-off frequency and higher mode separation than the corresponding rectangular wave guide. This partic- ular feature has lead to its commercial use in airline weather penetration radar. The only difference between the rectangular wave guide and the ridge wave guide is the superposition of a rectangular indentation on one of the sides of the rectangular cross section. Yet, this simple change is enough to make the problem extremely difficult to solve. More recently, work has been done on a semi- circular ridge wave guide. The solution of this problem has shown (4) the semicircular ridge wave guide to be superior in power handling capicity as compared with a rectangular ridge guide of the same mode separation factor. There are other means of solving the wave equation lO analytically besides a product solution, but each seems to require a great deal of ingenuity and specialization to a particular problem. It would be valuable to have a general method of solution to the wave equation. A convenient subterfuge is a numerical solution. Many types of numerical solutions are possible. Some of these will be discussed broadly so as to give a wider outlook on the problem, and then the particular method used will be discussed in detail. One numerical method is to substitute difference equations for the differential equation. It is then possible to write a difference equation for each node of a grid system imposed on the cross section of the wave guide. The next step is to effect a simultaneous solution of all the difference equations. Since the number of equations necessarily must be great in order to obtain accuracy, the simultaneous solution of these finite difference equations is not always easy. The following chapters will discuss the difference equations and methods of solving the difference equations. Finally, Chapter IV will discuss a program for MISTIC, embodying these ideas, which will solve the wave equation given any arbitrary cross section of wave guide. CHAPTER II DIFFERENCE EQUATIONS A. First Order Difference Equations Writing difference equations for a differential equation involves two procedures. First, it is necessary to impose a grid system on the domain of the problem so that a difference equation can be written for each node. The second is to find a suitable approximation or difference equation for the original differential equation. Thus, instead of having to solve one differ- ential equation, the problem reduces to a matter of finding a simultaneous solution for the numerous difference equations. The wave equation for which one would like to find a difference equation is 1'_ ‘27 L" 1 Actually, once one finds a difference expression for 72E the entire equation (2-1) is easily written in difference equation form. Therefore, the main objective of this chapter will be to discuss difference expressions for VQE. To begin with one imposes on the cross section of 11 12 the wave guide a grid system. A portion of this grid system is shown in Figure 2.1. «r—ln—a- u 71 .1 r I 5" it I l‘h 4 8 0 ,6 z i 54 i . . 7 r. r 3 -———a.x I Figure 2.1 The solid lines show the actual grid lines; while the dotted lines are only for the purpose of deriving the difference equation. A derivation of the first order difference equation in rectangular coordinates (5, p. 51) will now be given. Since E is a function of two variables, the definition of jL—-is (BX . "' x+.‘.‘. E x ~19. ‘5) i, __ Ltwt t( o 31 Ho) — ( o 1. j o (2-2) ex 0 " W” k as . . q; . where -——- means the partial derivative at the p01nt ‘AX.O .3x 0, etc. If h is taken small enough, .95., z 5.21:5. (2-3) ex 0 k and similarly, E7 _. ’ - 53—71 % ___F2 F0 anJ 5.9.5) A, E; F“ (2_4) ‘9 G k .9)( g h Now, since 3’? = 3;... at? 9x1 9x “532 it follows that 3E] ”35' (32; ex A ”7 a §1~ A; (2.5) ° )1 By substituting equations (2-4) into equation (2-5) one obtains 915' A, EJE‘I "259 ‘— 1, ~ (2-6) 29x c, ha. 1 This is the finite difference equation for 3%... I. In a similar manner §}§wbecomes: 92? (V F‘I‘F - 25: .——-—1I N I 3 (2_7) 9;; o ht. Combining equations (2-6) and (2-7) one obtains £3'*£;.+EEI+E;{—-VEEI ha. In equation (2-8) the expression E1+E2+35+E4-4EO is the “‘5 V7151 (2-8) finite difference approximation for the Laplacian, VZE. As a check on the error (5, p. 52) of this approximation, assume the function E(x,y) can be expanded in a Taylors series. First, check the error in the l4 approximation for 6X49. Assume E(x,y) to be a function of X with Y constant; one may write F09) ‘-'" 4.. +4:(X.'Y43+ 01094;)" 1 (2-9) +a3(x.-X;) + aq(g-x,-)+ ... Where i = 2,4,0 0 o o . With X = X0 as the origin, one then obtains a 4055:” 959‘1/ “42:15:- 3|a=9£ I ext 0 ' ' 3 OX3) Now letting i = 2 and 4, one obtains: t 4 E1 ”‘0 +4.14 “'1" +43A’1—q,é + " ' (2-10) L E, =4o—Q,A+az‘1 ~43A3+aqlty+ . ' ' (2-11) Adding equations (2-10) and (2-11), I A“ ‘/ ‘ 1+5]: 2410 +944 +2441: {- .. ° (2-12) Now solving equation (2-12) for a2, one obtains 2 E: f- : «9)! 0 vhl )2, git Therefore, the error in the finite difference equation (2- 6) is given by the algebraic sum of the terms ”k:L $E'ooo 0 etc. H N" By an analogous procedure, one can find an analogous error associated with‘;Y4; 'Hence, the total error in approximatingV2 is just twice that assumed x 15 by %%)0 alone . B. Higher Order Difference Equations Thus far only first order difference equations have been mentioned. A great deal more accuracy can be accomplished by using higher order difference equations. , For example, in an auxiliary method of solving the wave equation for a rectangular wave guide with a grid system of ten pivotal points, the error in determining the lowest characteristic value was 2.26% for the first order difference equation. However, applying the second order difference equation in the same method produced an error of only 0.08%; while a third order difference equation produced a low error of 0.0lOS%. The second order difference equations can easily be found by assuming a power series expansion. Given the points of a grid system as shown in Figure 2.2, first, compute %§%i . «L 5' .. i ‘8 ,fi 0 .2 ‘6 T I T V 3 Jr» 3 L X ‘P' 7 Figure 2.2 Assume that the E field existing at these points is given 16 as a power series expanded about the center point. i.e. E(X’J) 2: ”0+ ‘4,(Xo*&)+a1(xo —x’" )14‘ a3 (xg"X,(' )3 L (2—14) 4- (MUCH) + ' ' ' Where i = 2, 4, 6, and 8. It then follows that F = E, +a,h + a1h1+43h3+ at, h 4 (2-15) E9 :: E;*Q'I’l+allt1-a3 L3+4qhq (2-16) ’2 A; =- 5, +2q,h #14111 + 843134144qu (2-17) x H E8 :- Eo—JQ,A~ 4:40.114 '843112't-IGQ‘II" (2-18) Adding equations (2-15) and (2—16) and adding equations (2-17) and (2-18) it follows that .— 7. 5+5, = 2601524,). +,?a.,l.‘/ (2-19) F ~ 9+8 lit-r324 ‘1" 54E; - 3 o a, ‘I (2-20) By multiplying equation (2-19) by sixteen and then subtracting equation (2-20) from equation (2-19), a4 is eliminated, and the resulting equation can be solved for 2a2. Finally one has 20 __ at} = /6(E;+€v)-[E;+Eg) ”305. (2 21) 2 9X‘ 0 lent Equation (2—21) is the second order difference equation ‘5 for %}%u,. Similarly the second order partial derivative in -+>- I T 1 F1 1» l (a) rst Order -<>- -l -0- H. -l 16 lb -1 1 i '60 i 1 T T T T -0- lb - b Second Order 1? t ( ) A».z -"~ -27 -0- 270 Z ‘27 270 780 210 -37 2 i 1 ' 1 1 i T v if r + 1 4» 270 -27 Z (c) Third Order Figure 2.3 18 the other direction is 922:] -.- /6(E,'+E3)—-(Erf€7) “-7050 <9?g ° 131." Combining equation (2-22) and (2-21) one obtains ,2 715: “(an-9&5.) -— .(E{+E;+§,+E,.) —éoEo h1- The V72 Operator in difference equation form has (2-22) (2-23) been derived for the first and second order cases. V 2 operators which consider higher orders, points on the diagional, or even three dimensions can be derived in a perfectly analogous manner. Some of these Operators (6, p. 170) are shown in Figure 2.3. C. Special Type Difference Equations Besides difference equations of higher orders there are many other ways in which a difference expression can be written for the Laplacian. For one thing, it is not absolutely necessary to use a grid system with constant h. One might use a rectangular grid with h as the distance between nodes in the X direction and k as the distance between nodes in the y direction. It would be nice if h could be entirely variable. This would allow fitting a very close grid in regions of high irregularity and allow a very coarse grid in regions where there is little variation of the fields. Many types of coordinate systems may be used, 19 for example, a skew rectangular coordinate system or a triangular coordinate system. One other system which does hold a great deal of interest is the cylindrical coordinate system. The two dimensional Laplacian, to which one is restricted here, in the cylindrical coordinate system is 2.— 1” at: ’"E’ Vt 39—52:? Jgfs'F‘l‘JFI 953—1. (2-24) The grid system for this equation is shown in Figure 2.4. The difference expressions for the derivatives of equation (2-24) can be derived in a similar manner to those in the rectangular coordinate system (6, p. 224) 54¢, a“ o ’1 + 2 f & ‘fl 4‘ 43’ / {p c. Figure 2.4 20 and are ’37; v- at" q E -- 5’ 5‘11; :: E; +Ey “'95:: a¢t Y‘ where h =Ar, =4¢ Substituting equation (2-25) into equation (2-24) one Obtains in cylindrical coordinates. v72” =(t+—’—‘-— )2’. (J:— )1:- +(I——‘—*'- )E 0 are I to! 1 rat 3 * tale ”04:13)? E where rO is the distance from the origin to E0. Equation (2-26) is the first order difference equation (2—26) for the Laplacian.V72EIin a cylindrical coordinate system. It is entirely possible to write higher order difference equations in cylindrical coordinates or with a variable h or X'. However, the equations become very complex with these increasing variations. CHAPTER III METHODS FOR SOLVING THE DIFFERENCE EQUATIONS A. Matrix Methods Many methods of solving difference equations are possible, but just a few will be discussed here. The first to be discussed are two matrix methods which were investigated by the author. They were not found too helpful because of the lack of digital computer programs to handle large enough matrices or determinants. But they will be discussed since the only missing link is a computer with sufficient memory. By using a computer with larger memory, solution by these methods may be feasible. The wave equation, equation (2-1), can be written in matrix form as: , FE.” ’ 1 “PE." 2 z 2 t . 2 — 0 (3-1) v : Kc ’H : i. i , E31 _ J _ Egd In equation (3-1) {V1} represents the V2 Operator matrix. 21 22 Its typical entries can be filled in directly from the difference equation (2-8) and the grid system of the particular cross section of wave guide for which the wave equation is to be solved. Each row of the V72 matrix corresponds to a difference equation for a single node, and the columns correspond to the coefficients of the E values. The [u] is, of course, the unit matrix. The column matrices are n dimensional space vectors_ representing the E field in the wave guide. By multiplying Kg times the unit matrix and subtracting the right hand side of the matrix equation from both sides of equation (3-1); one Obtains 1p 1 T? r 'fi q'I’Kc qua qua ' " am E’ all sz'k: . .1 . :r C? an ' ' . ' (3-2) - ' . 2 F . H This matrix equation has non zero solutions only when the determinant of the coefficient matrix is zero. Obtaining an approximate solution to the wave equation is, then, just a matter of finding the eigenvalues and eigenvectors of the coefficient matrix. Programs for the MISTIC are available to compute eigenvalues and eigenvectors of symmetric matrices up to forty by forty. Once the wave equation is written in this form its 23 solution follows simply by using these programs. This method was used by the author to determine the characteristic values of a rectangular wave guide, Of aspect ratio two to one, using a grid mesh of ten points. Very good accuracy for the dominant mode was obtained by using second and third order difference equations, as was mentioned in Chapter II. However, there is an inherent problem. A grid system of more than just a few points, in this case ten, can not be used because the boundary conditions make the matrix asymmetrical even if the boundary of the guide is symmetrical. This is because boundary conditions enter the difference equations for nodes on the boundary of the mesh but not for those equations in the interior of the mesh. A program to compute the determinant of up to forty by forty asymmetrical matrices is available for MISTIC. This routine could be used by first making an approximation to the value of an eigenvalue. Using this value one could compute the determinant of the matrix. If the approximation were correct, which is highly unlikely, the determinant would be zero. Most likely the determinant would not be zero, so another approximation would be made. Again the determinant would be computed. If this approximation is not correct either, some method would be used to pick a better value 24 on the basis of the previous two. This process could be continued until a sufficiently accurate answer is obtained. A program of this nature is conceivable. However, one would still be limited to a grid system of forty points. This is certainly better than ten points, but it is desirable to have as many nodes as possible so as to be able to handle irregular shapes with more accuracy. B. Relaxation Since 1938 Sir Richard Southwell and a small group of his followers have been busy developing their numerical method known as relaxation. Relaxation (7) is a method for solving a number of simultaneous equations by successive approximations. It is a very convenient way Of solving the difference equations representing the wave equation. Rather than discuss relaxation for the general difference equation, the discussion here will be limited to the form of difference equation which repreSents the wave equation. The form of the difference equations for the wave equation is easily obtained. The finite difference expression for the two dimensional Laplacian was given in equation (2-8). Substituting this into the two dimensional wave equation (2-1), one Obtains 25 Eff/:2 IE3 +54 ’45: 2 "' Lil/(alga (3-3) Equation (3-3) is the basic difference equation form of the wave equation and will be used later. The subject of the wave equation will now be left for a moment. The following is strictly about relaxation. For the sake of simplicity in understanding relaxation assume equation (3-3) is instead E1+E2+Ea-+E9 ’95) ‘7 “'/0 (3-4) In equation (3-4) so is the value of the field at the central point of a grid system and El, E2, E3, and E4 are the values of the fields at nodes about the central point. A simple example using difference equations like equation (3-4) will suffice to show the principle of relaxation. Assume there are only two nodes in the grid system, as shown in Figure 3.1, and hence only two simultaneous equations. Call the value Of the E's at the two points, 0 O E: E1 0 Figure 3.1 26 El and E2. If the boundary points are all zero as shown, the two simultaneous equations would be, from equation (3-4) £2_+-C>+- O 1152 Figure 4.1 three digit sexadecimal number. The order in which the four groups are punched on the tape is irrelevant. Node 0 corresponds to a zero boundary condition. If there were a normal derivative boundary condition rather than a zero boundary condition, the tape would appear as 020 04K 001 132 indicating node 1 as its own neighbor. This information is read in and assembled into a MISTIC word with nine binary bits per node number and stored consecutively starting at location 17+n. This is done by the packing routine and must be done before the main routine is entered. When the packing routine is through, each node has stored in some memory location 55 the addresses of its surrounding four points. This information is used by the program time and time again in computing both new node values and new values OfJ . The packing routine is given in Appendix II. When the number of nodes is great, writing the interconnections can be quite tedious, especially since they must be written in sexadecimal. Therefore, a program was written to write the interconnection tape auto- matically. This program requires four bits of information about the grid system. It requires the number of points in the grid system and the number of horizontal rows of points. It also must be given the number of points in each horizontal row. Finally, it also requires the number one must add to a point in a row to obtain the point below it in the next lower horizontal row. This number is the same for each individual row, but one must be given for each horizontal row. The interconnections are then written assuming zero boundary conditions. With a slight modification the program can be used with other types of boundary conditions. This program is given in Appendix III. When the iteration program is entered, the first thing it does is to compute a value of J} using equation (3-11), the initial E values, and the interconnection information. Now that the program has a value of.S to 54 use it starts one complete iteration, changing the E values, using equation (3-9) and the interconnection information. While it computes new E values it also computes the numerator and denominator of equation (3-11). Then when it has completed an iteration, it has only to divide the sum corresponding to the numerator, by the sum corresponding to the denominator of equation (3-11) in order to obtain the trial eigenvalue for the next iteration. The iteration process can now begin all over again. After each iteration the absolute value of the new JTis compared with the absolute value of the old I. When the absolute difference Of these two values is less than a predetermined amount called the accuracy constant, the program quite and transfers control back to a master program. Before it transfers control it puts the final eigenvalue in location 16. Location 17 and up contain the correct E values. Thereafter, the print out of this information is easily available. B. Application of the Program The values at the nodes may represent any rectangular component of the electric or magnetic field one wishes. However, the boundary conditions are entirely dependent upon which component is used. These boundary conditions then show up in the interconnections. If a 55 TM mode is to be investigated, the most probable choice of component would be Ez since everywhere on the boundary Ez must be zero no matter what the shape of the cross section. The most probable choice for 3 TE mode would be Hz' In this case, the normal derivative Of Hz must be zero at the boundary. This implies that a node point next to the boundary must have the same value as the point on the boundary. If the cross section of wave guide has any physical symmetry it is usually possible to solve only a partial region, the other points being filled in by symmetry. The boundary conditions on the new boundary between the regions of symmetry have boundary conditions of the same nature as the actual boundaries. If the field has even symmetry with respect to this boundary, the normal derivative of the field must be zero. If the field has odd symmetry with respect to the boundary, the field itself must be zero along the boundary. By using symmetry it is possible to get considerably more accuracy because more points are concentrated in a smaller area of the total cross section. Besides the interconnection tape, the program also requires initial values. Of course the closer the initial values are to a normal mode, the faster the convergence. Convergence will be to the mode which comes closest to the initial distribution Of E values. 56 If the initial values are assumed to be identical, convergence is not necessarily to the dominant mode. Consequently, one must use as initial values, values which are thought to approximate the desired mode. Since the lower modes are usually the ones of interest, one usually does have a general idea as to what the fields will look like. If necessary the desired cross section could be slowly changed from some known shape, each time using for the initial E values on a newly perturbed cross section, the final E values for the unperturbed shape. C. Program Improvements Most likely the greatest improvement with the least amount of change to the program could be accom- plished by using second order or higher difference equations. This would, nevertheless, mean a considerable cut in the maximum number of nodes which could be used. For the second order difference equation eight inter- connection addresses would be needed for each node rather than four. It would probably reduce the maximum number of nodes available from about four hundred at present to less than three hundred. If the cross section to be analyzed was more similar to a circle than to a rectangle, it follows that a similar program written using a cylindrical coordinate 57 system difference equation would give more accurate results than the present program. This would be due to the fact that the cylindrical grid system would more nearly fit the boundary. In order to get an idea of this type of error, the present program was used to compute the characteristic value of the TM01 mode for a circular wave guide. These results are given at the end Of this chapter. More accurate results, using this program could be obtained by using extrapolation. A number of characteristic values could easily be obtained, each from a smaller mesh size. From these values a more accurate characteristic value could be extrapolated. A more classical method of approximating Ki after each iteration was found after the program had been written. It is called the Rayleigh quotient (7, p. 173) and is supposed to give a good approximation to Ki. It is possible that use of the Rayleigh quotient would give faster convergence. It is a weighted average rather than just a simple average. The Rayleigh quotient is shown in equation (4-1). 1 "" EVzE—JX/ /’{c = // ‘7 (4-1) flElJXJj. As a summation equation (4-1) becomes equation (4-2). Obviously, whenever the grid system does not 38 exactly coincide with the boundary some error will result. This error can be made smaller by using Fox's correction 2.. /(L “EEO-V to c —- 2 E2- (4-2) 0 formula (7, p. 65) for points next to a boundary. Fox's formula would be used rather than equation (3-9) and is given in equation (4-3). A typical boundary on which it might be used is shown in Figure 4.2. I I 2E. 29, 2E, 22:. ——-——- + ...—— + .+--—-—-— E 7047) c—(He) 1+"; Ht- (4 3) ° _ 2. 3—+2— —— m E- ’7 ./ 2 . __.=—" l :‘T T h .g7ffir Th 3 I c3 3‘[ .5:l'7h- 0 .7 LI 4 ‘ 5/ ... h 44. Figure 4.2 The more complicated procedures like this one were not employed because it was hOped to keep the program as simple as possible to allow for the maximum number of nodes. 39 D. Results Two results will be discussed. The first shows the accuracy of this method versus other methods and the second shows what type of accuracy can be-expected when the boundary is such that it does not fit the grid system very well. In the first example a simple rectangular wave guide was solved. The dominant mode was investigated by using the Ey component. Ey must be zero on the vertical sides, and the normal derivative of Ey must be zero on the upper and lower boundaries. A grid system of ten points was used as shown in Figure 4.3. E: EL 5": £1 E” 0 60 Et- E.) E? E! 0' a o E}, E, Ea E1 59 0'? L e. e, E9 5, a. X' ~h' L Figure 4.3 The dielectric material is assumed to be free space. The aspect ratio is two to one. The initial values used agreed with the known solution of the TElO mode only in the first digit. This very same problem was solved using the matrix method of Chapter II. If h is given as one meter, which 40 is convenient for comparison purposes, the eigenvalue obtained by the matrix method for the TElo mode is 2 Kc = .26794919. This value is the exact value of the lowest order eigenvalue for this ten by ten matrix. Any other method should approach this value as the correct one. It is not the correct value for the wave guide, in fact, it is in error by some 2% as was mentioned before, but it is the correct eigenvalue for the ten by ten matrix. Therefore, if the iteration method is as accurate, one should get the same answer. The eigenvalue obtained by the iteration method, again for the TE10 mode, is Ki = .2679473. Even greater accuracy is possible since the accuracy constant of the iteration program was not set at its minimum.- It is just as well, for this accuracy is much more than is needed. Greater accuracy would only require more computing time. The answer above was obtained in nineteen iterations. In the first example only ten points were used. In this next example a total of three hundred forty-one points were used. Since a circular wave guide does not fit a rectangular grid system very well it was felt by trying to solve the circular wave guide one might gain insight into the errors involved. An outline of the grid system placed on the circular wave guide is shown in Figure 4.4. The TMOl mode was the mode investigated. The Ez component 41 was used since it is zero everywhere on the boundary. P—Ffi ,7 I r” *— I Figure 4.4 The initial values used somewhat approximated a sine wave, with maximum values toward the center and small values toward the edge. Seventy iterations were required for convergence. This is to be expected since the percentage of boundary points is much lower than in the first example. In this example the cut-off frequencies will be compared rather than the characteristic values. The cut-off frequency of the TM01 mode in a cylindrical wave guide (1, p. 376) of radius a is given in equation ('4-4). . 383 +; -= 0‘ ‘V/U: (4-4) For a cylindrical wave guide with a radius of 11 cm. the cut-off frequency for the TM01 mode is calculated from equation (4-4) to be 2,090 mcs. The eigenvalue computed by the program was 42 KO = .2228. From it the cut-off frequency can be obtained by using equation (1-5). The cut-off frequency as calculated from the iteration program is fo = 2,130 mcs. This answer should be more exact for the polygon approximating the circle than for the cylindrical wave guide, but this answer still only differs from the value for the cylindrical wave guide by 1.8%. This is sufficient accuracy for numerous engineering applications. Slightly greater accuracy could be obtained by using four hundred points rather than three hundred forty-one. CONCLUSION While the program described here is far from perfect it does describe a very practical method of solving this type of problem. Many improvements for the program have been mentioned. A battery of three or four programs using different variations could fit mesh points into almost any type of configuration. The accuracy one can obtain is limited only by the ingenuity of the program writer and most important by the size of the digital computer memory. The greater the memory the greater the number of nodes which can be used. This thesis has dealt exclusively with the two dimensional wave equation. There is no reason why the same procedure can not be extended to the three dimensional wave equation. Again the major problem would be to obtain sufficient memory space. It would then be possible to consider discontinuities in wave guides or wave guides of non-uniform cross section. Finding such things as the resonant frequency of a resonant cavity could be easily done. By using a different type of boundary condition the field patterns of antennas could be studied. 43 44 It has been tacitly assumed that the dielectric medium contains no electric charges and no conduction currents. It is possible by modifying the wave equation that even these effects could be studied by this same procedure. Problems of this type would require computers with quite extensive memories, however. LOCATION 0 APPENDIX I THE ITERATION PROGRAM ORDER KSF " 4268L 5017s J080L SSF 109? 468L 102F 428L 50173 J081L 85F 469L 1023 429L L582L _ L4F 7 L4F L4F L4F NOTES Unpack interconnection addresses Compute first E value using equation (3-9) 45 46 LOCATION ORDER NO ES 10 40891. I 5089L ll 7J9OL 5082L l2 L537L 4213L 13 2213L LSF .J 14 4084L '1 Add to sum representing L49lL numerator and denominator of 15 409lL equation (3—11) LSlOOL l6 4031L 508L 1? J097L 83F 18 3619L 2620L l9 F531L 4031L 20 508L J099L 21 83F 3222L LOCATION 22 23 24 25 26. 27 28 29 30 31 32 33 ORDER 2223L F531L 4031L 509L JO97L 83F 3626L 2627L FSBIL 403lL 509L JO99L 33F 3229L 2230L F531L 4031L L582L 2231L 2632L L484L 2633L L484L 2634L 47 NOTES 48 LOCATION ORDER NOTES 34 L484L 2635L 35 L484L L492L 36 4092L 2637L 37 L584L 40l7F 38 2640L 00F 39 00F 00F __ 4O F54L -7 Reset to compute second and 404L following E values by equation 41 L51L (3-9) L485L 42 401E F537L 43 4037L 2645L 44 00F 00F 45 F586L 4086L 49 LOCATION ORDER NOTES 46 LO3F ‘ 3247L 47 26lL 2648L _ 48 L569L ‘ Reset for new iteration 4012L 49 L577L 4013L 50 L54L L03F 51 404L 0020F 52 46lL L587L 53 2654L 00F 54 4237L L582L 55 4086L F588L 56 40881. 26571. J 57 50921. 7 Compute new I by equation (3-11) 7J90L 50 LOCATION ORDER NOTES 58 5082L 669lL 59 S5F 4083L 6O L582L 409lL 61 L582L 40921. ..4 62 15931:.-1 Check convergence of' S 4094L 63 L783L 4093L 64 L593L L094L 65 4095L L795L 66 L096L 361L 67 2267L L583L 68 40l6F 22F J 69 66831.‘1 Program constants SSF LOCATION 70 71 72 75 74 75 76 77 78 79 so 81 e2 83 e4 e5 86 87 88 89 90 91 92 93 51 ORDER NOTES OOFOOF OOFOOF OOFOOF OOFOOF OOFOOF OOFOOF OOFOOF 2614LOOF OOFOOF OOFOOF 3L3584FLL2048F 00511F002044F OOFOOF 00F00400000000000J (Location of S ) OOFOOF OOlFOOF OOFOOF (Counts E's computed during one iteration) OOl7FOOl7F OOFOOF (iteration count) OOFOOF OOFOOlOOOOOOOOOOOJ OOFOOF OOFOOF OOFOOF 52 LOCATION ORDER NOTES 94 OOFOOF 95 OOFOOF 96 OOFOOlOOOOOJ (Accuracy constant) 97 004095FOOF 98 OOFOOF 99 OOFOO4095F lOO 2231L2632L The following is required before entering the routine: 1. n must be in location 3 2. n must be added to the left hand order in cell 1L 3. location zero (OF) must be cleared 4. n initial conditions must be stored in location 17F through (17 + n - 1) F 5. n interconnections must be stored in locations (17 + n) F through (17 + 2n - 1) F where n = the number of nodes in the grid system APPENDIX II PACKING ROUTINE. LOCATION ORDER NOTES 0 KSF -‘ Pack one word 4215L l 5124L 8012F 2 OO29F L417L 3 40l7L L52L 4 L018L 402L 5 L52OL L023L 6 402OL 3611. ~ ,7 L517L '7 Reset for following words 40l7F 8 L524L 40l7L 53 54 LOCATION ORDER NOTES 9 F57L 407L lO L52L L422L 11 402E L516L 12 4020L F521L 13 4021L LO3F l4 3215L 2615L 15 26lL 22F A 16 OOFOOBF I Program constants 17 OOFOOF 18 OO9FOOF l9 OOFOOF 2O OOFOO3F 21 OOFOOF 22 OO36FOOF 23 OOFOOlF 24 OOFOOF d The following is required before entering the sub-routine: 55 l. n must be stored in 3F 2. 7L must have n added to it before entering 3. n sets of interconnections as shown on page 32 must be in the reader where n = the number of nodes in the grid system LOCATION O wQOm-§\NNH +4 F’ +4 14 F‘ r4 l4 F1 +4 (n -4 ox \n -> v: n) +4 <3 APPENDIX III INTERCONNECTION PROGRAM ORDER NOTES L598FL486L 4086L262L F588L40100F 4088LL086L 3615LF52L 402L262L OOFOOF OOFOOF OOFOOF OOFOOF OOFOOF OOFOOF OOFOOF OOFOOF OOFOOF F589L4089L LO6OOF3622L L5101F2226L 409OLL589L 56 LOCATION 19 20 21 22 25 24 25 26 27 28 29 5o ‘ 31 52 55 54 55 36 57 38 59 4O 41 42 45 ORDER L087L3624L L591LL49OL 0016F2225L L591L2628L 409OL2218L L599FL49OL OOl6F8224F 263OL5091L 0012F2618L 5091L0012F 2623LOOF L57OOFL089L 3633LL591L 409OL264OL L57OOFL089L LO91L3231L L592LLO7OOF 4237L2237L OOFL5F 0012F409OL 264OLOOF L56OOFL089L LO7OIF3648L L5601FL46OOF L089LLO701F 57 NOTES 58 LOCATION ORDER NOTES 44 3645L2648L 45 L592LL4701F 46 0020F4647L 47 L5F2649L 48 2248LL591L 49 L49OL4090L 50 L59OL0016F 51 8224F92131F 52 2653LO0F 55 F592LL493L 54 4O92LL517L 55 L495L4OI7L 56 L524LL493L 57 4024LF594L 58 4094LL06OOF 59 , 3664L2615L 60 OOFOOF 61 OOFOOF 62 OOFOOF 65 OOFOOF 64 F558L4058L 65 4242L0020F 66 4616L464OL 67 LO95L4654L 68 F558LOOZOF 59 LOCATION ORDER NOTES 69 4642LL53OL 7O L493L403OL 71 4633LF543L 72 4043L4245L 73 OO2OF4641L 74 L591L4089L 75 4094LF535L 76 4035LF562L 77 4062LLO99F 78 3679L2615L 79 OFFOOF 80 OOFOOF 81 OOFOOF 82 OOFOOF 83 OOFOOF 84 OOFOOF 85 OOFOOF 86 OOFOOI6F 87 OOFOOZF 88 OOFOOl6F 89 OOFOOF 9O OOFOOF 91 OOFOOF 92 L5100FL5100F 93 OOlFOOF 60 LOCATION ORDER NOTES 94 OOFOOF The following specification tape must immediately follow the program: 9.0992 OOFOOnF where n = total number of nodes OOFOOaF where a = total number of horizontal rows 006OOK OOFOOblF where ba = number of points in OOFOObZF horizontal row number a OOFOObaF OOZOOK OOFOOF where Xa = the number one must OOFOOX1F add to a point in row a to OOFOOXZF obtain the point below it in row 3 (a + l). OOFOOXa_1F OOFOOF 243N REFERENCES Ramo, S. and Whinnery, J.R. Fields and Waves in Modern Radio. John Wiley and Sons, New York. 1953. Spence, R. D. The Propagation of Electromagnetic Waves in Parabolic Pipes. Masters Thesis, Michigan State University. 1942. Cohn, S. B. Properties of Ridge Wave Guide. Proc. Van Bladel, J. and Von Rohr Jr., 0. Semicircular Ridges in Rectangular Waveguides. IRE Trans., vol. Grinter, L. E. (ed.), Numerical Methods of Analysis in Engineering. The Macmillan Co., New York. 1949. Salvadore, M. G. and Baron, M. L. Numerical Methods in Engineering. Prentice-Hall, Inc., New York. 1952. Allen, D. N. deG., Relaxation Methods. McGraw Hill Book Co., Inc., New York. 1954. 61 sun... '93": Q. n 1‘ by.” {1. fflé F 'ijl" P: 9‘3!!- 'E‘Wfl f MICHIGAN STATE UNIVERSITY LIBRARIES 0 3085 0279 3 1293