v u r L:_ i 3‘53 04%? ,Uifii A: Jhaa‘lfi This is to certify that the thesis entitled BOUNDARY SOLUTION OF THE ‘Z‘VE EILUATIOII presented by Elm-{APE J. PEID has been accepted towards fulfillment of the requirements for PhOD. agree in E. E. («We ”MTwA Major professor Date August 21, 1959 0-169 LIBRARY Michigan State University BOUNDARY SOLUTION OF THE WAVE EQUATION By 2,5 RICHARD JF‘AREID 9 1 if, {~_ U) 026.5 A THESIS Submitted to the School for Advanced Graduate Studies of Michigan State University of Agriculture and Applied Science in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Electrical Engineering 1959 ABSTRACT The solution of the wave equation subject to boundary condi- tions existing on irregularly shaped surfaces is obtained by assuming a solution function that is a truncated, n-dimensional Taylor’s series about some point on the boundary or within the region where a solu- tion is required. The assumed series is of degree mi (1 = 1,2,...n) in each co- ordinate (i). The re are then n 17:1 (m1 + 1) coefficients of the series which are chosen to minimize the deviation of the solution function from the prescribed boundary values and also to minimize any residual remaining upon substitution of the as- sumed function into the wave equation. The method of solution has the properties that: (a) the error of the solution can be estimated; (b) boundaries can be of any shape; (c) calculations for the eigenvalues and eigenfunctions of the solution can be made exclusively at coordinates of the boundaries; ii (d) the solution can be made to correspond to a least squares, Chebychev, or other type fit over the region of interest. iii ACKNOWLEDGMENTS The author wishes to acknowledge with gratitude the help and guidance received from his major professor Dr. L. W. Von Tersch. Grateful acknowledgment is also due Dr. G. P. Weeg for his helpful suggestions and assistance. iv I. II. TABLE OF CONTENTS INTRODUCTION .......................... APPROXIMATING SOLUTION FUNCTION .............................. BOUNDARY SOLUTION FOR EIGENVALUES AND EIGENFUNCTIONS .......... EXTENDING THE NUMBER OF VARIABLES ............................. SIMULTANE OUS WAVE EQUATIONS ............ MODIFYING THE SOLUTION FUNCTION .............................. EXAMPLE OF EIGENVALUE AND EIGENFUNCTION DETERMINATION ............. DISCUSSION . . . . . ...................... . 10 15 16 18 20 24 25 I. INTRODUCTION The solution of the electromagnetic wave equation subject to boundary conditions existing on surfaces deviating more than a slight amount from surfaces on which x, r, 9, et cetera, are constant,1 re- quires the use of numerical methods to obtain approximate results. The usual numerical method divides the space over which the equation and boundary conditions are to be satisfied into a mesh and produces at the nodes of the mesh 3 difference equation which ap- prommates the wave equation. This difference equation is obtained by assuming a polynomial expansion of the solution function in each coordinate direction about each particular node. The degree of the polynomial is commonly the same as the order of the wave equation, and the wave equation is used to determine the relation among the coefficients of the polynomial. This method is described by Booth (3) . A more versatile approach to the solution of the wave equa- tion and other partial differential equations has been contributed by 1Methods for small boundary perturbations are given by Slater (1) and by Schelkunoff (2). Milne e_t_a_l. (4), and assumes an n-dimensional polynomial (where n is the number of independent variables) about each node point, and al- lows the degree of the polynomial to be as large as practical from the considerations of required accuracy of solution and of the feasi- bility of calculating the coefficients for the polynomial from the given partial differential equation. All these approaches require the simultaneous solution for the function value at each node subject to the function value satisfying all polynomials that have coefficients determined from the position of this node. An iterative method is usually employed for the solution of these equations. These methods thus require: (a) special calculation of the polynomial coefficients and the function values for given boundary points which do not lie at the nodes of the mesh, or an approximation of the desired boundary that does pass through the nodes; (b) a new determination of the polynomial coefficients and/ or smaller mesh size and then iterating to a solution in order to realize more accuracy in the solution. These will then yield an approximate solution at the nodes of the mesh in terms of a table of values corresponding to the field’s component values. The method to be investigated here is that of assuming a functional form for approximating the solution, since the functional form will allow: (a) (b) (C) calculations giving the error in the solution;2 boundaries to be of any shape since there is no mesh and node system to which they must correspond; calculations for the eigenvalues and eigenfunctions of the solution to be made only at coordinates of the boundaries. The methods previously discussed also permit error analysis. II. APPROXIMATING SOLUTION FUNCTION The wave equation that will be directly considered here is that resulting from Maxwell’s equations for sinusoidal time variations describing one of the coordinate components of the electric or mag- netic field in a charge-free, linear, isotropic medium where the parameters of the medium are constants; i.e., VZE = -w2ueE (1) Equation (1) must be satisfied for some solution function, E, which normally must also satisfy boundary conditions, as in two di- mensions 59—51354) = o ‘ (2) 11 along some curve C1, where n is normal to C1, and E(x,y) = 0 (3) along some curve C2. Since Equation (1) has constant coefficients, if the boundaries C1 and C2 can be formed from a finite number of rectifiable arcs, then, as Bernstein (5) has shown, there ezdsts a solution, E, to Equations. (1), (2), and (3) that is of a class A” (continuous in all derivatives) in the neighborhood of all points on boundaries C1 and CZ. This then allows a Taylor’s series expansion for E about some point on the boundary, and it will be assumed that a finite number of terms of the series will give the soltuion function approximately for some maximum irregularity of the boundaries . The truncated, two-dimensional Taylor’s series, f(x,y) is given by the matrix product mm -= [Yl'[AJ[X] (4) where [Y]' = [1:Yv-nYn-1] ; y = y' - yo and [X]' = [1,x,...,xm-1] ; x = x' - x and [A] is an nxm coefficient matrix. For the definition fyy = [Y] 'nyAHXJ where [Yh'yy = [0,0,2,6y,...,(n-1)(n-2)yn-3] and similarly for fxx' Equation (1) becomes upon substitution of f(x,y) , the E approximating function, 2 _. fxx + fyy + w uef — R1(x,y) in S, the region of solution, where R1(x,y) is the wave equation residual or error function. Also in Equations (2) and (3) is 1 -<*—‘ m1 z... _ ay 2 ’ 1+ dy) 1+ dx along C1,and f=R3(x,y) along 02' A total residual function, RT’ can be defined as RT = J! W1(x, )R12(x,y)da + 5/; W2(x,y)R22(x,y)dl + 5/; W3(x,y)R32(x,y)d1 (5) where W1(x,y), W2(x,y), and W3(x,y) are nonnegative weighting func- tions introduced to allow variance in the degree of satisfaction of Equations (1), (2), and (3). For example, if W1(x,y) is chosen large with respect to W2(x, ) and W3(x,y), then the assumed function will more closely satisfy the wave equation than it will the boundary conditions. This allowance of weighting functions is particularly useful when examining small changes in a section of the boundary, as that boundary section can be relatively heavily weighted, thus eluci- dating its effect upon the solution. The squares of the partial residual functions have been used in Equation (5) in defining the total residual function. This will lead to a solution function fitting the wave equation and boundary condi- tions on a least squares basis. Similarly, for other types of fit, the magnitudes of the partial residual function could be used in Equation (5) or, for a Chebychev type fit, the single maximum mag- nitude of R1(x,y), R2(x,y), or R3(x,y) could be chosen as the total residual function. Of course, this maximum will shift among R1, R2, and R3, and range over the region of solution as the elements of the matrix A are changed. For any of these choices, the total residual function, RT, will be considered as a function of the coefficients in the matrix A; , ..., a ). 1.e., RT(a 1m nm 11, an, a When the least squares basis is used, as in Equation (5), the total residual function will be quadratic in the elements of the matrix A. In order to obtain the absolute minimum value for the total residual function, the requirements for a relative minimum will first be imposed; i.e., z z “T ..., j=1 i=1 daij and the (m.n) X (m.n) matrix r 62% War 75—1-1” ’9a1193nm : (7) 92% £51. L.)anmAa11’ ’c9anm2 - must be positive definite. Since RT is of quadratic form in the aij’s, the elements of Equation (6), 3-??? , will be linear and Equation (6) will be homoge- neous. In addition, since RT is quadratic and always positive, it is guaranteed of having a single minimum. Let the solution of Equation (6) be a set of aij’s, an. The calculations involved in determining the 511’s will be a check that RTGII’ ..., En. :i: e, Em) > RT(;11’ ..., 31]., ..., Enm) for all i and j and for some small e. The finite nature of e for this check can be allowed since .32; is linear and consequently can have no more than one root in the region 51]. :l: h, h g e. The function RT can possibly be reduced to zero when the assumed solution form is the same as the exact solution form. Also, it is possible to obtain a zero for RT by the numerical method of evaluating Equation (5) . This misleading result can occur if the integrals are evaluated in a particular dimension at fewer points than the degree of the polynomial plus one, since the polynomial would then be nonunique. Careful attention to the evaluation of the integrals in Equation (5) will avoid this difficulty. If, in a particular problem, it is required to determine the eigenvalues of Equation (1) , then wzue must be included as one of the independent variables of the function RT. The total residual function will no longer be of quadratic form, however, and likewise, 3%- and Wag-{Ea will not result in linear equations. This must be so, of cause, since for the exact solution there will be an in- the finity of values for wzue. The resulting difficulty of determining the 2113s with wzue as a variable can be overcome by a knowledge of a good approximation to wzue. This approximation is originally introduced as a constant and the best values for the aij’s determined for the resulting linear case. Then the wzue is treated as a vari- able to obtain wzue. III. BOUNDARY SOLUTION FOR EIGENVALUES AND EIGENFUNCTIONS Consider an exact solution for Equations (1), (2), and (3), G(x, ) so that 2 2 V G(x,y) + w ue G(x, ) = 0 (8) for all x and y in the region of solution S and on the boundary C. Let there be another function H(x,y) so that G(x,y)lC = H(x,y)IC , (9) 60:9: c = .9132; IC ’ (1°) 6G§:,y) C = angzy)|C , (11) 623x?! IC = .2313 Ic’ (12) and C (13) 2‘3!) 2‘1!) dgyx IC = &§yx Higher order derivatives may also be equal when evaluated along C, but in general, could be quite different. If all orders of derivatives were equal along C, or along any other curve, or at 10 11 any point then the functions G(x,y) and H(x,y) would be identical in the neighborhood of this equivalence. This requirement is not to be imposed upon H(x, ) at this point. It can be seen from Equations (10) and (11) that H(x,y) will satisfy the normal derivative boundary condition. Of course, this first order derivative equivalence need only be satisfied upon that part of the curve y = f(x) where the normal derivative boundary condition is to apply. If the function H(x,y) is substituted into the wave equation, the resulting nonhomogeneous partial differential is V2 H(x,y) +w2ue H(x,y) = J(x,y). (14) In order to guarantee that the eigenvalues determined from Equations (8) and (14) are the same, it will be sufficient to force the function J (x,y) to negligible values throughout the region of solution. It is here assumed that the function H(x, ) is also forced to fit the proper boundary conditions so that J (x,y) is already negli- gible on the boundary. The fact that J(x,y) must be reduced within the 'region of solution in order to determine the eigenvalues accurately is illus- trated in the following example, although, as will be shown, it is possible to obtain useful approximations to the eigenvalues when the values for J (x,y) within the region of solution are not negligible. Consider the one -dimensional wave equation 12 2 Lgéfl + wzue G(x) = 0 dx subject to the boundary conditions G(x) = 0 at x = £1. The first imagined solution might be a parabola, say H(x) = 1 - x2 which satisfies the boundary conditions but does not satisfy the wave equation at the boundary. Therefore, first choose the second deriva- tive of H(x) to be zero at the boundary, then construct H(x) itself to be zero at the boundary; i.e., 2 Mg) = 12::2 - 12 clx M=4x3-12x+c dx 1 x4-6x2+cx+c 1 2' H(X) Now by choosing c equal to zero and c2 equal to plus five, the 1 function H(x) also satisfies the boundary conditions now as H(x) = x4 - 6x2 + 5. Substituting these results into the wave equation results in 12x2 - 12 + wzue(x4 - 6x2 + 5) = J(x). Now, while J (x) has roots at x equal to plus. or minus one, it can not be made zero over a finite interval for a constant value of 13 wzue. For the determination of possible roots of J (x) 2 12 - 12x2 12 W 116 = _ x -6x2+5 5a-x2 which has a maximum value of plus three at the boundary of the interval of solution and a minimum of 2.4 at x equal to zero. The average value over the interval is determined as The lowest ordered, exact solution of the one-dimensional wave equation for these boundaries is G(x) = cos £215, and the principal eigenvalue determined for this is wzue =(%)2 = 2.467 The deviations of the approximate eigenvalues from this cor- rect value result because the function J (x) has not been sufficiently reduced within the interval of solution. While the previous results may in some cases be adequate or useful approximations, additional requirements which also can be imposed at the boundaries alone will yield the eigenvalues and the eigenfunctions to any desired accuracy. It is noted that the nonhomogeneous part of Equation (14) , J (x,y) corresponds 'to the partial residual function R1(x,y) in 14 Equation (5). However, the preceding development suggests that the integral of R12(x,y) over the region of solution be replaced in Equa- tion (5) so that the total residual function is. now determined as m n . . 2 RT = S: Xwn' [&1+J.R1(3T:Y).l + f W2(X’Y)R22(X’Y)d1 j=0 i=0 l dxl dyJ J(xo,y0) C1 + f W3(x,y)R32(x,y)dl C2 where the summation terms are the Taylor’s series coefficients for the expansion of the function R1(x,y) about the point (x0,y0). Thus, if these coefficients are reduced sufficiently by the reduction of RT to a minimum, the assumed solution function is automatically fitted to the differential equation in the region about (xo,y0). IV. EXTENDING THE NUMBER OF VARIABLES The case of Taylor’s series in two dimension, Equation (4), can be extended to any number of dimensions by repeatedly applying the extension rule that g(x,y,2) = [Zl'lBl where [Z]' = [1,z,...,zp-1] ; z = z' - z0 and [B] is a column matrix of p elements such that each element bi is given as b , [Y1'[A,1[x1 i = 0,1,2,...,p-1 and there will thus be p separate A1 matrices for the three- dimensional case. For the case of n-dimensions such that an mj degree polyno- mial is assumed in the j’th coordinate, there will then be a total of n 1T j=1(mj + 1) elements belonging to A matrices. Equation (6) would be applied to all the elements of the A’s. 15 V. SIMULTANEOUS WAVE EQUATIONS Although the wave function will often satisfy the condition that it is expandable in a Taylor’s series about some point in the region of solution, when a coefficient of the wave equation (11 or e) has discretely different values in the region of solution, this will no longer be possible. Consider the case of Equation (1) where e must have two distinct values within a region and the boundary between the e values is given by some curve C 3. The solution in either of the regions (as determined by C and the original bound- 3 aries) is continuous in all derivatives, and it is desired that the solution in each of these regions give the same values along the curve C Since from Equation (1) if E(x,y) is to 3. be continuous while e has a change in value, then the second de- rivatives of E(x,y) must be discontinuous which the assumed solu- tion type will not provide. For problems such as these it will be necessary to assume multiple solutions, one for each region where derivatives of all orders may exist. Thus for two such regions, a and b, separated 16 17 by a curve C two series will be used, fa(x,y) and fb(x,y), and the 3 residuals within each of these regions calculated. The total re- sidual function for the entire problem is then calculated as RT =f£ Wla(x,y)R1§(x,y)da +124 W1b(X,Y)R1%(X:Y)d3 2 2 + fc W23(x,y)R2a(x,y)dl + J; W2b(x,y)R2b(x,y)d1 1a 1b + Jézaw3a(x,y)R3§(x,y)dl + 42bW3b(x,y)R3§(x,y)dl + 43W4(xyy)[fa(x,y) - £b(x,y)12dl (15) where Cla is the part of the boundary of region a to which the boundary condition of Equation (2) applies, similarly for Clb’ C2a is the part of the boundary of region a to which the boundary condition of Equation (3) applies, similarly for CZb’ and W4(x,y) is a weighting function applied to the criterion of the degree satisfac- tion that both functions fa(x,y) and fb(x,y) have the same values along the curve C3. The coefficients of the functions fa(x,y) and fb(x,y) can now be calculated by applying Equation (6) to the residual function found from Equation (15). VI. MODIFYING THE SOLUTION FUNCTION Since it is often necessary to obtain solutions for the wave equation where it is known that the magnitude of the solution function varies over a considerable range in the region of interest, it is pos- sible to modify the assumed solution function to correspond in form to the suspected general shape of the actual solution. The following is a list of these modifications and the type of problems to which they apply. Partial radial symmetry When there is some apparent degree of radial symmetry in cases of two or more dimensions, a more general assumed solution is possible that corresponds to the Method of Frobenius in one- dimensional problems with nonconstant coefficients. The analogous functions in two dimensions is. given as _ ‘ ' l [190, 991, 92rp2,..., on lrpn'lj T 11, 2 f(r,9) = A r ;.m-1 where the pl’s will be additional variables in the total residual ftmction and will be calculated in the same way as the aij’s; 18 . ll I.I.Il!. .1}: lIIIDl .I. 19 i.e., by evaluating the partial derivatives as in Equations (6) and (7). Of course, the V2 operator must be in the same coordinate system as the assumed, independent variables. Zero at infinity The assumed function can be modified to yield a zero at in- finity by multiplying the assumed Taylor’s series by a weight func- tion, so that the assumed function is. g(x,y) = 6"(XZ + Y2) f(x,y) where f(x,y) is the normal Taylor’s. series as given in Equation (4). The polynomials corresponding to g(x,y) in one dimension are the Hermite polynomials. First approximation to the solution function When an apparently good approximation to the solution function is known in terms of explicit functions, say h(x,y), then the assumed solution g(x,y) = h(x,y) + f(x,y), where f(x,y) is the Taylor’s series, may result in a low degree for. f(x,y) for good results in g(x,y). Thus g(x,y) may be a relatively simple explicit function easily used in following analytical investiga- tions . VII. EXAMPLE OF EIGENVALUE AND EIGENFUNCTION DETERMINATION A complete, digital -computer program to obtain the approxi- mate solution for the wave equation in two dimensions has been prepared for the Computer Laboratory (6). The program requires that the matrix A be square and have an order equal to or less than ten, so that polynomials through the ninth degree are possible. Mixed boundary conditions of 3% = 0 and E = o are possible. The program user must supply the following information: 1. The order of the assumed matrix A up to and including ten. 2. The rectangular coordinates of (a) points in the region at which it is desired to have the assumed function satisfy the wave equation; (11) points on the boundary where the condition E(x,y) = O is to apply; (c) pairs of points which define a normal to the boundary where the. condition éfiaififl is to apply. 3. An initial approximation for the entries in the matrix A and for wzue. The program will then determine the total residual function as defined in Equation (5) with the weight functions all taken as unity. 20 21 This fimction is minimized as a function of the elements of the ma- trix A and of wzue. The term a is not treated as a variable, how- 11 ever, and maintains the initial value it is given at the outset. This is done to prevent the matrix A from degenerating to the zero matrix during the minimization procedure. The center of the coordinate sys- tem therefore must be chosen to coincide with a point within or on the boundary where it is known that the solution function will have a nonzero value. The solution for the following problem has been calculated using this program. Given that the single boundary condition is E(x,y) = 0 along the circular boundary R = 0.5 meters. Find the smallest eigenvalue cor- responding to this boundary condition and the corresponding eigen- function. For this problem the origin for the solution function was chosen to be at the center of the circular boundary, since the func- tion will not be zero there. The computer program was arranged to hold the initial ap- proximation of the eigenvalue constant until the eigenfunction had been fitted relatively closely to the proper solution as indicated by a sufficiently small total residual function. The initial conditions chosen were: and F+0100 +0000 +0495 +0000 +0410 +0000 -0132 +0000 _+0023 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 20 -0495 +0000 +2400 +0000 -2000 +0000 +0640 +0000 -0110 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0410 +0000 -2000 +0000 +1680 +0000 -0530 +0000 +0090 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 -0132 +0000 +0640 +0000 -0530 +0000 +0169 +0000 -0030 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0023— +0000 -0110 +0000 +0090 +0000 -0030 +0000 +0005‘ where the initial AC, was chosen to be the truncated series solution for the square boundary within which the circular boundary in this problem may be considered inscribed. With the above as initial conditions, the eigenvalue and eigen- function were calculated to be and wzue = 23.183 "+0100 +0000 -0595 +0000 +0709 +0000 +0368 +0000 _-0027 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 -0595 +0000 +2700 +0000 -1900 +0000 -1910 +0000 +2640 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0709 +0000 -1900 +0000 +1880 +0000 -0181 +0000 +0742 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0368 +0000 -1910 +0000 -0181 +0000 +0170 +0000 -0038 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 +0000 -0027‘ +0000 +2640 +0000 +0742 +0000 ~0033 ' +0000 +0000_ 23 The total residual function was reduced so that the a1 j’s above are correct to about the third decimal place. At this time the residual function numerically determined to approximate Equa- tion (5) by replacing the integral by unweighted sums was RT = 0.10004. The correct eigenvalue obtained from the analytic solution of the wave equation in cylindrical coordinates is wzue = 23.150 which is in close agreement with the value determined for this prob- lem. VIII. DISCUS SION The methods developed here have been applied exclusively for solutions of the wave equation. These methods can also be used on any other partial differential equation where a solution is known to exist and adequate approximating functions are known. The same ap- plies to ordinary differential equations where the development in Part III obviates the need for integrating throughout the intervals of solu- tion. This should be most valuable in ordinary differential equations with two-point boundary conditions. For the present, application of these methods is somewhat limited by computational speeds available. For the example worked herein, a total calculation time of about 20 hours was necessary in a machine with an add time of 100 microseconds. Thus, it seems for the wave equation, at least, only rough approximations could be obtained if three independ- ent variables were included. The major factor in time-consumption is the inadequacy of the numerical methods used in finding minima of functions of several variables. Two numerical minimization methods were used on the example herein with no significant differences in the calcula- tion times or the results. The methods were that of “steepest de- scen ” and minimization with respect to one variable at a time. 24 IX. REFERENCES CITED Slater, J. C. Microwave Electronics, pp. 80—83. D. Van Nostrand Co., 1950. Schelkunoff, S. A. Applied Mathematics for Engineers and Scientists, p. 298. D. Van Nostrand Co., 1948. Booth, A. D. Numerical Methods, pp. 110—37. Butterworth Scientific Publications, 1955. Milne, W. E., et al. Wright Air Development Center Technical Report 57-556. Mathematics for Digital Computers, Volume I, Multivariate Interpolation, February, 1958. Bernstein, D. L. Existence Theorems in Partial Differential Equations. Annals of Mathematics Studies No. 23. Princeton University Press, 1950. Computer Laboratory. Library Program Gl-M, Wave Equation Solution. Michigan State University Computer Labora- tory, 1959. 25 .nnr. r“ ”g! F. 4;! " “1 t‘ BSS‘I ‘91:;th -' \u' ...H‘f‘nn «- 0:. n .- ‘2‘ JV 4!: ‘ WW 2 3a. n-n‘fifit‘f’ . l .1... 2"-» --‘- I 1 mm:- . "I -' h. F , s l _ ‘ ' *‘u ' s “4