⬇ Download Markdown · ⬇ Download LaTeX · ⬇ original · view original · CSV: table 1, table 2, table 3, table 4, table 5, table 6
vorgelegt von
Dipl.-Inf. Ibrahim Abu El-Khair von der Fakultät IV-Elektrotechnik und Informatik
der Technischen Universität Berlin Zur Erlangung des akademischen Grades Doktor der Ingenieurwissenschaft
-Dr.-Ing.- genehmigte Dissertation Promotionsausschuss:
Vorsitzender:
Prof. Dr. B. Lutterbeck Prof. Dr. G. Mönich Prof. Dr. G. Bärwolff Tag der wissenschaftlichen Aussprache: 24. November 2006 Berlin 2006
D83 i
The objectiveness of this thesis xvii Dedication xvii Acknowledgment 1 1 Principles of EM Radiation 1 1.1 A single Wire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Principle of field detachment shown by small loop and finite length antennas . . . . 4 1.2.1 7 2 Derivation of Time-varying Fields for a Hertzian Dipole . . . . . . . . . . . .
Maxwell’s Equations 2.1 15 The classes of Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.1.1 Maxwell’s Equations in differential form . . . . . . . . . . . . . . . . . . . . . 18 2.1.2 Maxwell’s equations in integral form . . . . . . . . . . . . . . . . . . . . . . . 19 2.1.3 Maxwell’s equations in transform domain . . . . . . . . . . . . . . . . . . . . 20 2.1.4 Introductory to the phasor-notation . . . . . . . . . . . . . . . . . . . . . . . . 22 2.1.5 Duality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Contents ii 3 Tracing of Fieldlines 29 3.1 The longitudinal sectors’ approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2 A synoptical glance on classical methods of field lines representation: . . . . . . . . . 32 3.2.1 A line as the intersection of two surfaces . . . . . . . . . . . . . . . . . . . . . 32 Line density versus flux intensity . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4 E-wall , H-wall and Image theory 35 4.1 Magnetic and electric walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2 Image Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5 From EM Fields to fieldlines 43 6 Start points for line tracing 49 6.1 7 Topology of surface and displacement currents . . . . . . . . . . . . . . . . . . . . . . 49 6.1.1 Discretization of the conducting plane . . . . . . . . . . . . . . . . . . . . . . . 50 6.1.2 Line families . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.1.3 Chronological visualization of the field evolution . . . . . . . . . . . . . . . . 53 6.1.4 The new surface definitions, applied to a straight wire . . . . . . . . . . . . . 54 On the arbitrarily oriented radiators 7.1 57 A provoked solution: Replacement Rules . . . . . . . . . . . . . . . . . . . . . . . . . 58 How does it really work? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.1.1 7.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Rotation Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Contents iii 7.2.1 The Adopted Rotation Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 7.2.2 What are transpose Matrices really good for? . . . . . . . . . . . . . . . . . . . 63 A correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 7.3 To Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 An example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Yet another example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 7.3.1 For an arbitrarily oriented Hertzian Dipole . . . . . . . . . . . . . . . . . . . . 69 magnetic fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Electric Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 7.3.2 For an arbitrarily oriented finite length linear dipole . . . . . . . . . . . . . . . 70 Magnetic Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Electric Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 8 Numerical Approach 73 A minimal requirement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 A well-posed Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Sources of error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 8.0.3 8.1 8.2 Order of Accuracy, Convergence and Stability of numerical Methods . . . . . 76 Why to undertake a numerical approach? . . . . . . . . . . . . . . . . . . . . . . . . . 78 8.1.1 The main reason: Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . 78 8.1.2 Which strategies are used to plot the field lines? . . . . . . . . . . . . . . . . . 79 A system of differential equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 8.2.1 Application: Kepler’s Problem; a test problem . . . . . . . . . . . . . . . . . . 86 Contents iv Kepler’s Laws first . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 9 Results 107 9.1 An omnidirectional radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 9.2 A directional radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 9.2.1 Appendix A passage to 4D visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 115 A On small square loop Antenna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 B On small circular loop Antenna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 B.1 The Magnetic Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 B.2 The Electric Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 C A finite length dipole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 D Mathematical Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 D.1 Taylor Polynomials and Taylor Series . . . . . . . . . . . . . . . . . . . . . . . 149 D.2 Numerical techniques to approximate derivatives. . . . . . . . . . . . . . . . . 150 What does O(h p ) , or Big O mean? . . . . . . . . . . . . . . . . . . . . . . . . . 152 D.3 E Summary, revision and more . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 The Algorithms: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 E.1 One-step methods: fixed step-size . . . . . . . . . . . . . . . . . . . . . . . . . 161 Euler’s Method in Mathematica code . . . . . . . . . . . . . . . . . . . . . . . 161
Contents v Heun’s Method in Mathematica code . . . . . . . . . . . . . . . . . . . . . . . 161 The Runge-Kutta Method RK4 in Mathematica code . . . . . . . . . . . . . . 162 E.2 One-step methods: variable step-size . . . . . . . . . . . . . . . . . . . . . . . 163 The Runge-Kutta-Fehlberg Method RK45 in Mathematica code . . . . . . . . 163 E.3 Multi-step methods: fixed step-size . . . . . . . . . . . . . . . . . . . . . . . . 167 The Adams-Bashforth-Moulton Method ABM in Mathematica code . . . . . 167 The Milne-Simpson Method MS in Mathematica code . . . . . . . . . . . . . 169 The Hamming’s Method in Mathematica code . . . . . . . . . . . . . . . . . . 171 E.4 Multi-step methods: variable step-size . . . . . . . . . . . . . . . . . . . . . . . 173 The Adams-Bashforth-Moulton Variable step-size Method ABM V in Math-
ematica code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 Bibliography 177 Index 177 vii
1 Hertz’s Drawings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1.1 synoptical view on both field types, illustrating Maxwell’s equations . . . . . . . . . 3 1.2 Single wire Radiation, configurations for different possible cases. . . . . . . . . . . . 3 1.3 Single wire Radiation, configurations for different possible cases. . . . . . . . . . . . 4 1.4 A small circular loop antenna. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 A short linear dipole antenna. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.6 Using Biot-savart law to Derive the formula of a Hertzian Dipole . . . . . . . . . . . 8 2.1 The classes of Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Differential flux at a point on closed surface . . . . . . . . . . . . . . . . . . . . . . . . 19 3.1 Explaining the longitudinal sectors’ approach . . . . . . . . . . . . . . . . . . . . . . . 30 3.2 Integration of flux for a given area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.1 The electric fieldlines of a circular loop antenna 4.1(a) and a linear dipole 4.1(b) . . . 36 4.2 The electric fieldlines of a circular loop antenna 4.2(a) and a linear dipole 4.2(b) . . . 36 4.3 The E-wall and the respective surficial electric currents. . . . . . . . . . . . . . . . . . 37 4.4 The H-wall and the respective surficial electric currents. . . . . . . . . . . . . . . . . . 37 List of Figures viii 4.5 The E-Wall for a small circular loop antenna and a linear one. . . . . . . . . . . . . . 39 4.6 The H-Wall for a small circular loop antenna and a linear one. . . . . . . . . . . . . . 39 4.7 E-Wall, source and image for a small circular loop antenna and a linear one. . . . . . 40 4.8 H-Wall, source and image for a small circular loop antenna and a linear one. . . . . . 40 4.9 E-Wall, H-Wall, source and image for a linear radiator. . . . . . . . . . . . . . . . . . . 41 5.1 The field lines and their symmetry surfaces . . . . . . . . . . . . . . . . . . . . . . . . 44 6.1 Border contours with the region of interest and current lines . . . . . . . . . . . . . . 50 6.2 Regions of interest, border contours, sectors and field lines . . . . . . . . . . . . . . . 51 7.1 Rotating Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.2 A Z-oriented Dipole centered at origin, along with its electric symmetry plane! . . . 61 7.3 The first Rotation about Z-axis! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 7.4 The second Rotation about X’-axis! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 7.5 The third Rotation about Z’-axis! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 7.6 Two coordinates systems! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 7.7 Dimensions for a finite length linear dipole! . . . . . . . . . . . . . . . . . . . . . . . . 68 8.1 A block diagram showing a rough configuration of ‘Divide and Conquer Scheme’ . 79 8.2 The processing of an order of tasks form . . . . . . . . . . . . . . . . . . . . . . . . . . 80 8.3 Center of Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 8.4 The Expert System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 8.5 The Interaction among ‘Divide and Conquer’, Expert System, and Display Unit . . . 83 8.6 Poor Performance of Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 List of Figures ix 8.7 π Euler performance for h = 800 , it takes 3100 points, and 0.093 seconds for the result . 93
8.8 π Euler performance at h = 1600 , it takes 4100 points, and 3.171 seconds for the result . 93
8.9 π Euler performance for h = 3200 , it takes 8100 points, and 15.813 seconds for the result 94
result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 8.11 Second Order accuracy is obtained by using the initial derivative at each step to
find a point halfway across the interval, then using the midpoint derivative across the full width of the interval. In the figure solid dots represent final function values, while open dots represent function values that are discarded once their derivatives have been calculated and used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
π 8.13 Heun performance for h =
1600 , it takes 4100 points, and 4.265 seconds for the result 96
π 8.14 Heun performance for h =
3200 , it takes 8100 points, and 15.343 seconds for the result 97
8.15 Fourth-order Runge-Kutta method. In each derivative is evaluated four times: once
at the initial point, twice at trial midpoints, and once at a trial endpoint. From these derivatives the final function value(shown as a filled dot) is calculated. . . . . . . . . 98
the result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
π , and Tol = 10 −13 8.17 Adaptive step size Runge-Kutta-Fehlberg performance for h =
800 it takes 3100 points, and 2.203 seconds for the result . . . . . . . . . . . . . . . . . . . 100 8.18 Integration over [x n , x n+1 ] in Adams-Bashforth method. . . . . . . . . . . . . . . . . 101
seconds for the result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
the result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 π 8.21 Hamming performance for h =
800 , it takes 1800 points, and 0.953 seconds for the result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
List of Figures x 9.1 Field lines and their corresponding planes of symmetry . . . . . . . . . . . . . . . . . 107 9.2 Field lines and their corresponding planes of symmetry . . . . . . . . . . . . . . . . . 108 9.3 Different circulation of inner and outer formations. . . . . . . . . . . . . . . . . . . . 108 9.4 Same circulation of inner and outer formations. . . . . . . . . . . . . . . . . . . . . . . 109 9.5 synoptical view on both field types, illustrating Maxwell’s equations . . . . . . . . . 109 9.6 3-D Radiation patterns of arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 9.7 The directivity of arrays: two elements . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 9.8 The directivity of arrays: four elements . . . . . . . . . . . . . . . . . . . . . . . . . . 112 9.9 The directivity of arrays: six-element broadside array 3D fashion . . . . . . . . . . . 112 9.10 The Visualization of the field lines in 4D fashion . . . . . . . . . . . . . . . . . . . . . 113 11 A small square loop Antenna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 12 A small circular loop Antenna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 13 A finite length dipole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 14 A finite length dipole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 xi
2.1 Overview of Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Some time-domain functions and their phasor-domain equivalents. . . . . . . . . . . 25 3.1 Synoptic view of the three examples above, showing fieldlines as the intersection of
two surfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 7.1 Replacement rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 7.2 The rotated fields components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 7.3 two more dimensions for a finite length . . . . . . . . . . . . . . . . . . . . . . . . . . 68 8.1 The relation between evaluations number and local truncation error . . . . . . . . . . 98 1 Some available Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 xiii
The title of thesis is “On The Visualization of Time-varying complex Electromagnetic Field lines in 3D and 4D Fashions”, which means that either the visualization depicts one instant, taken at a desired time, or depicts the field lines distribution and propagation over an interval of time.
Michael Faraday’s Concept of lines of force: Faraday’s research into electricity and electrolysis was guided by the belief that electricity is only one of many manifestation of the unified forces of nature, which included heat, light, magnetism and chemical affinity. Although this idea was erroneous, it led him into the field of electromagnetism, which was still in its infancy. Inspired by the discovery of Oersted and Ampere, that an electric current produces a magnetic field, and his own idea about conservation of Energy led him to believe that since an electric current could cause a magnetic field, a magnetic field should be able to produce an electric current. He demonstrated this principle of induction in 1831. Faraday expressed the electric current induced in the wire in terms of the number of lines of force that are cut by a wire.
Faraday’s introduction of the concept of lines of force was rejected by most of the mathematical physicists of Europe, however, this descriptive theory of lines of forces moving between bodies with electrical an magnetic properties enabled Clerk Maxwell formulate an exact mathematical theory of electromagnetic waves.
The lines of forces or field lines as we call them today, are the mean ever to demonstrate the intensity and distribution of a field, a field is a vector quantity defined by its magnitude and direction at each point in space. In case of current harmonic distribution flowing in a linear Dipole Antenna, for instance these field lines tend to form closed loops, magnetic circular loops surround the dipole axis, lie in parallel planes right-angled to dipole itself, and electrical kidney-shaped loops in meridian planes that contain the dipole itself. In 1888 Hertz drew an influential series of diagrams to accompany his 1889 paper on dipole radiation, which was translated as “The forces of electric oscillations, treated according to Maxwell’s theory.” Hertz took a great deal of care with these drawings. They have been reprinted innumerable times since their first appearance, often directly from Hertz’s List of Tables xiv originals (see Fig. 1). Merely by using a simple formula that expressed the electric field in term of the magnetic field multiplied by the radius ρ of a cylinder surrounds a z-oriented dipole, he came out with a formula which depends on ρ and sin 2 (θ).
Unfortunately, eversince the technique has been all but abandoned, to be replaced by:
' Tables containing numerical data about the field intensity etc!
' Radiation Patterns, to give an overall pictorial illustration of the behavior of the radiant
object in the far-field.
' Different visualizations depicting contours represent different intensities.
' Or by using many 2D or 3D arrows.
All these methods, though they are certainly useful somewhere else, are not as near as good in representing the field lines as a solid line, viewed as tangential trajectory of the field itself at every point which the field line passes.
Hertz’s approach, is unfortunately restricted to a single Hertzian dipole, and demands rotational of symmetry as a prerequisite, to produce the field lines described above.
Using multi-start-points chosen carefully according to certain scheme, and relying on numerical methods to trace these trajectories, enabled us to produce the desired visualization to our best expectation. The following terms are corner stones in our approach:
At least one Plane of Symmetry A prerequisite to ensure the closure of a field line, to form a closed loop. Here we distinguish between two types of these planes:
' Electric plane of symmetry which is, for convenience, a plane that make a right-angle to axis
of the dipole.
' Magnetic plane of symmetry which is, a plane that contains axis of the dipole or parallel to
it.
The Image Theory The image theory boosts the applicability of this technique to cover general cases comprise arbitrary oriented and located dipoles, the symmetry plane works like sort of a magic mirror, that reflects the images of individuals within a group of radiators, according to prescribed manner. The sum of the fields of sources and images due to the symmetry enable us to treat them as a resultant or net field to be visualized.
List of Tables xv Searching Starting and launch Points A scheme has been developed to guarantee that each field line depicted represents the same amount of Flux ψ or φ. The border lines on the infinite planes of symmetry, moreover, have additional function beside being delimiters between launch points and falling points, they are the place where the Birth/Death Process of those field lines takes place. This is literally the passage to the time-varying visualization of the fields lines over a given time interval, or merely parts of it, reduced to discrete points, within an interval.
The entire work,is intended to be a humble gesture of deference to Michael Faraday, the greatest experimentalist with real clairvoyant powers and perception.
Figure 1: Hertz’s Drawings xvii
This work is dedicated to the shining stars of mine,
Their Highnesses My beloved Parents!
xix
I would like to thank professor G. Mönich, my chief supervisor for his unlimited support and encouragement. Thanks go to professor G. Bärwolff and professor B. Lutterbeck, members of my supervisory board together with professor G. Mönich for their assistance and suggestions over the course of this work and in the preparation of this thesis.
Pretty special thanks go to the dear ones of mine, Milady and sweetheart Jin Wei and Milady Dr. Susanne Bräuning for their absolute and true love, support, encouragement and patience.
Also thanks go to lady Kader Polat from the Residency Office of Berlin for her friendly attitude.
1
In physics, we learned that current, either alternating or direct, had been always associated with its magnetic and electric fields. It is usual to refer to the combined fields as electromagnetic fields.
It was also taught that high frequency electromagnetic currents in a wire antenna, also in turn result in high frequency electromagnetic fields around the antenna, which results in the electromagnetic radiation, in waves that move away from the antenna into the free space at the velocity of light (≈ 3 · 10 8 meter per second.)
Conservation of energy is possibly the most important, and certainly the most practically useful of several conservation laws in physics.
The law states that the total inflow of energy into a system must equal the total outflow of energy from the system, plus the change in the energy contained within the system. In other words, energy can be converted from one form to another, but it cannot be created or destroyed.
This law helps us to understand the mechanism of EM radiation, as both electric and magnetic field, possess a form of energy each. In case of time-varying fields, this energy can either be absorbed or leveled, in accordance, the intensities of those fields experience respectively growth or decay phases.
Considering a lossless and nonconducting medium as free space only, one comes to the conclusion that only two forms of energy are possibly available there, namely the electric and magnetic energy. Therefore, it is obvious that a magnetic field which experiences a growth-phase, can only 1 2 Principles of EM Radiation extort the required energy from a decaying electric field that sets its energy in the free space. The very same is valid for an electric field that experience a growing-phase.
The energy thus, undergoes conversion process and be passed from the electric to magnetic fields and vice versa, as long as, the radiator is fed. The response to the feeding signal is accompanied with a certain delay time however, that prevents the abrupt collapse or growth of these fields.
The electric and magnetic fields, which induced by alternative sources are always coupled, and always occur together as given in Fig. 4.1.
Considering a linear thin wire for the moment we conclude that; for the existence of electromagnetic radiation, either time-varying current or an acceleration or deceleration of charge has to be assumed. We usually refer to currents in context with time-harmonic applications, while charge is most often mentioned in transients. Periodic charge acceleration (or deceleration) or time-varying current is also created when charge is oscillating in time-harmonic motion.
Therefore, whether EM Radiation takes place, one comes to the following cases:
¬ If a charge is not moving, current is not created and there is no radiation.
If charge is moving with a uniform velocity, equals that of light:
• There is no radiation if the wire is straight, and infinite in extent. • There is radiation if the wire is bent, discontinuous, truncated or composed of two segments
connected through a mismatched impedance , as shown in Fig. 1.2
® If charge is oscillating in a time-harmonic motion, then:
• There is no radiation if the oscillation occurs at the velocity of light and wire has straight
infinite extent. • There is radiation if the wire is bent, discontinuous, truncated or composed of two segments
connected through a mismatched impedance , as shown in Fig. 1.2 In all the cases mentioned above, where the EM Radiation occurs, sort of harmony has been disturbed for a while, the response, in form of EM Radiation resembles a very loud protest, or like 1.1 A single Wire 3 pain signals flow through the nerve system, conveying to us that sort of disturbance had taken control.
(a) The electromagnetic field lines of a small circular loop antenna (b) The electromagnetic field lines of zoriented λ 2 antenna Figure 1.1: synoptical view on both field types, illustrating Maxwell’s equations
(a) Two Wires connected through an impendence (b) Bent Wire Figure 1.2: Single wire Radiation, configurations for different possible cases.
1 4
(a) Discontinuous Wire Principles of EM Radiation (b) Truncated Wire Figure 1.3: Single wire Radiation, configurations for different possible cases.
A qualitative understanding of the radiation mechanism may be obtained by considering the creation of the field lines. We will demonstrate that by two examples, the first represents the magnetic field lines of a small circular loop, and later by demonstrating the electric field lines of a finite length linear (short) dipole, as shown in Fig. 1.4 and Fig. 1.5 Fig. 1.4(a) depicts the magnetic field lines, during the first quarter of a time-harmonic cycle of an alternative current up to its peak, the outmost loops were detached before, and they travel outwardly of the center of the loop antenna. The magnetic fields induced by the present current enjoy a growth period, and increase in number and expansion in size.
Fig. 1.4(b) depicts the immediate instant after entering the second quarter and shortly before leaving it, the current starts to decay, and so do the field lines in accordance. Due to the fact that the outmost ones among them, have expanded considerably in space, the decay of the current will influence but the portions which are close to the loop antenna, so that saddle points start to form. The closer the lines are, the stronger the force acts on them, the very closed field lines shrink gradually and proportionally to the decay of the current, till most of them get absorbed in the antenna inductance, while the further ones experience the detaching process, depicted in Fig. 1.4(c), so new kidney-shaped loops are thus created, with a different polarization.
1.2 Principle of field detachment shown by small loop and finite length antennas Fig. 1.4(c) emphasizes the detaching process, for simplicity, just the last loop is considered here, the saddle points coincide at one point and then the formed loops detach, in a manner that reminds us with an “Ameba Splitting In Half During Fission.” Two daughter-loops replace the original one, which undergoes detachment, the outer one moves outwardly and expands its dimension, while the other one keeps shrinking and get closer to the loop antenna, till it gets absorbed within it.
Fig. 1.4(d) depicts the second half of the current-wave, where the polarization takes the opposite sign than before. It shows the instant in which the peak is reached.
While the loop antennas, due to their configuration form closed loops, only currents are expected within them, the linear dipole, however works as depositories for charges, which place themselves on either rod of the dipole in a different polarity. Their distribution is mostly found near the tips of the antenna, once the first quarter of the charging cycle is over. Immediately after that, the antenna undergoes the conducting phase, which reaches its maximum a quarter cycle afterward, then the charging phase starts once over again and reaches its maximum also after another quarter of the cycle. This process continues as long as the antenna is fed.
The external source accelerates the charges and set them in motion and thus occurs the radiation,while at the end of the wire the buildup of charge concentration induces counter field, that decelerates the charges. Therefore, charge acceleration due to an exciting electric field and deceleration due to impedance discontinuities or smooth curves of the wire are mechanisms responsible for EM Radiation.
While both current density J ~ and charge density ρ e appear as sources in terms of Maxwell’s Equations, charge is viewed as a more fundamental quantity, especially for transient fields. Even though this interpretation of radiation is primarily used for transients, it can be used to explain steady-state radiation as well.
Fig. 1.5 depicts the creation of the electric field lines of such type of antennas. Fig. 1.5(a) depicts a totally energized antenna during the first quarter of the cycle, positive charges are assumed at the upper half, and their negative counterparts at the lower, the electric field lines originate from the positive charges and land at their negative counterparts forming a pattern similar to that of the electrostatic dipole, with the exception that these curves move outwardly with time instead. The outmost kidney-shaped loops were as well created during a previous cycle, so they are not of significant interest to us right now.
5 1 6 Principles of EM Radiation Fig. 1.5(b) depicts the situation during the second quarter of the cycle, the conduction phase starts, and charges diminishes, the inmost curves are absorbed by the signal generator, which the larger one get smaller. Detaching points starts to form, and the effects appears clearly on the next lines to the antenna, which detached in two parts, one forms a closed loop that moves outwardly, which the other part keep diminishing, till eventually get swallowed by the signal generator.
Fig. 1.5(c) depicts the start of the third quarter, where the charging phase starts, and the polarity reverses here, so the new electric fields originate from the lower half of the antenna and end at the upper, in the meanwhile the detachment process of the lines induced by the former quarter of the cycle is accomplished, once the third quarter is reached, which depicts in Fig. 1.5(d). The entire process repeat itself from this time on, in a manner similar to that in Fig. 1.5(b) and Fig. 1.5(c) but in a different circulation.
Figure 1.4: A small circular loop antenna.
1.2 Principle of field detachment shown by small loop and finite length antennas
Figure 1.5: A short linear dipole antenna.
determined by Biot-Savart Law
the point P in space, and r is the distance from the current element to the point of interest. The magnetic field is perpendicular to both the direction of the line and the direction from the current element to the point P .
7 1 8 Principles of EM Radiation
Figure 1.6: Using Biot-savart law to Derive the formula of a Hertzian Dipole A couple of simplification can be first done:
Which, since d ` ~ is constant, becomes
Now i(t) was assumed at the start to be a time-varying current given by i(t, r = 0) = I 0 sin(ωt) at the dipole. At distances away from the dipole, Eq.(1.4) cab be written as:
1.2 Principle of field detachment shown by small loop and finite length antennas 9 ~ = ∇ × A, ~ we then extract the value of the z-oriented Magnetic Vector Potential to be: With B
However, there is a retardation effect due to the time required for the EM field to propagate to the
Therefore we obtain:
Similarly we get (1.9)
(1.10)
1 10
(1.11)
and
From Ampere-Maxwell Equation (Maxwell’s equations are described in details in Chapter 2 on page 15)
~ = ε 0 E, ~ then we have where in free space J ~ = 0 and D
(1.14)
(1.15)
(1.16)
(1.17)
1.2 Principle of field detachment shown by small loop and finite length antennas 11
and finally E x is obtained by integrating Eq. (1.19) with respect to the time
(1.20)
1 12 Principles of EM Radiation
and finally E y is obtained by integrating Eq. (1.22) with respect to the time
1.2 Principle of field detachment shown by small loop and finite length antennas 13
x The second term − ∂H
∂y is similarly obtained
(1.25)
1 14 Principles of EM Radiation
and finally E z is obtained by integrating Eq. (1.27) with respect to time
15
Figure 2.1: The classes of Maxwell’s equations A passage of Feynman (1) summarizes how physical problems and phenomena should be seen and understood 2 16 Maxwell’s Equations Mathematicians, or people who have very mathematical minds, are often led astray when “studying” physics because they lose sight of the physics. They say: ‘Look, these differential -the Maxwell equations-are all there is to electrodynamics; it is admitted by the physicists that there is nothing which is not contained in the equations. The equations are complicated, but after all they are only mathematical equations and if I understand them mathematically inside out, I will understand the physics inside out.’ Only it doesn’t work that way. Mathematicians who study physics with that point of view-and there have been many of them-usually make little contribution to physics and, in fact, little to mathematics. They fail because the actual physical situations in the real world are so complicated that it is necessary to have a much broader understanding of the equations.
Thanks a lot, Mr. Feynman, that is very true indeed, and we may add, that a genuine engineer should be a genuine combination of both; a mathematician and a physicist, yet more biased to be a physicist, that is. A physical understanding is a completely unmathematical, imprecise, and inexact thing, but absolutely necessary to understand a certain phenomenon as a whole, instead of fractions or parts.
The unification of electric and magnetic phenomena in a complete mathematical theory was the achievement of the Scottish physicist Maxwell (1850’s). See Table (2.1). In a set of four elegant equations, Maxwell formalized the relationship between electric and magnetic fields. In addition, he showed that a linear magnetic and electric field can be self-reinforcing and must move at a particular velocity, the speed of light. Thus, he concluded that light is energy carried in the form of opposite but supporting electric and magnetic fields in the shape of waves, i.e. self-propagating electromagnetic waves.
In doing this, Maxwell moved physics to a new realm of understanding. Maxwell actually was inspired by Faraday’s experimental work and by the mental picture provided through the “lines of force” that Faraday introduced in developing his theory of electricity and magnetism. By using field theory as the core to electromagnetism, we have moved beyond a Newtonian world-view where objects change by direct contact and into a theory that uses invisible fields. This introduces a type of understanding which can only be described with a type of mathematics that cannot be directly translated into language. In other words, scientists where restricted in talking about electromagnetic phenomenon strictly through the use of a new type of language, one of pure math. Fig. 2.1 offers a possible classification of problems which, Maxwell’s equations deals.
¬ Maxwell’s new theory provides a new description of light, as electromagnetic waves.
2.1 The classes of Maxwell’s equations 17 Electromagnetism represents a sharp change in the way nature is described, i.e. the use of
invisible fields and understanding that can only be communicated with mathematics.
® Light is electromagnetic radiation propagating through space.
¯ The wavelength of the light determines its energy and characteristics.
Electromagnetic radiation is energy that is propagated through free space or through a material medium in the form of electromagnetic waves, such as radio waves, visible light, and gamma rays. The term also refers to the emission and transmission of such radiant energy. The wavelength of the light determines its characteristics. For example, short wavelengths are high energy gamma-rays and x-rays, long wavelengths are radio waves. The whole range of wavelengths is called the electromagnetic spectrum.
In 1887 Heinrich Hertz, a German physicist, provided experimental confirmation of Maxwell’s ideas by producing the first man-made electromagnetic waves and investigating their properties. Subsequent studies resulted in a broader understanding of the nature and origin of radiant energy.
Differential Form
Table 2.1: Overview of Maxwell’s equations 2 18
Maxwell’s Equations
Maxwell’s equations in differential or point, form are:
(2.1)
(2.2)
(2.3) (2.4)
where
electric field (in V /m), magnetic field (in A/m), electric flux density (in C/m 2 ), magnetic flux density (in W b/m 2 ), electric charge density (in C/m 3 ), magnetic charge density (in W b/m 3 ), electric current density (in A/m 2 ), magnetic current density (in V /m 2 ), and all eight quantities are, in general, functions of position ~r and time t.
Presently, there is no experimental evidence to confirm the existence of isolated magnetic charges, therefore, many authors set both ρ m and J ~ m equal to zero at the onset, in equations (2.2) and (2.4) (See Table (2.1)). There are, however, at least two reasons for retaining ρ m and J ~ m : first, symmetry between electric and magnetic is retained, second, equivalent magnetic sources appear in a variety of applications, such as radiation and scattering from aperture antennas or penetrable bodies.
By taking ∇· of (2.1) using (2.3) and the identity ∇ · ∇× ≡ 0, and by assuming the interchangeability of space and time differentiation, we obtain the continuity equation
which expresses conservation of electric charge. Similarly, from (2.2) and (2.4) we obtain the continuity equation 2.1 The classes of Maxwell’s equations
19
Equations (2.1) through (2.4) are necessary but not sufficient for determination of the eight field quantities (six vectors and two scalars) which appear in them. We still must specify the primary sources of the electromagnetic fields, as well as the physical properties of the medium in which the field exists; these physical properties take the form of functional relations among the various field quantities, called constitutive relations, we observe that in vacuo (or free space) the constitu~ = µ 0 ~ ~ = 0 E ~ and B H, where 0 and µ 0 are two constants called the electric tive relations are: D permittivity and the magnetic permeability of free space. In addition to 0 and µ 0 , a material is characterized by its conductivity σ, the term σ E ~ is called the conduction current density, which occurs in response to the impressed (feeding) current. The current density σ E ~ is a current flowing on a nearby conductor due to the electric field E ~ created by source J ~ (3).
Consider a fixed volume v bounded by the closed surface S, whose outward unit normal is ~n as shown in Fig. 2.2.
Figure 2.2: Differential flux at a point on closed surface Integration of ((2.3), (2.4)) and ((2.5), (2.6)) over v, followed by the use of the divergence theorem,
2 20 Maxwell’s Equations yields:
where Q e and Q m are the total electric charge (in C) and magnetic charge (in W b) inside the
~ through the Eq. (2.7) means that the total electric charge contained in v equals the outgoing of D surface S of v. A similar interpretation applies to (2.8), since Q m is zero, the total outgoing flux of ~ through any fixed closed surface is zero. B Eq. (2.9) means that the rate of decrease of the total electric charge Q e inside v equals the amount of electric charge which leaves v in unit time by traveling outward through S; thus (2.3) is obviously a statement of conservation of electric charge. A similar interpretation applies to (2.10); since Q m is zero, the outgoing flux of J m through S is also zero.
Maxwell’s equations ((2.1) to (2.4)) are a system of eight first-order partial differential equations in four independent variables: three space coordinates and time, whose solution is often quite complicated.
It may be advantageous to eliminate the dependence of the field upon one or more of the independent variables by applying Fourier (or Laplace) transform to ((2.1) to (2.4)), solving the resulting equations in the transform domain, and then obtaining the desired field quantities by an inverse transformation.
Obviously, the main advantage of a transform technique with respect to an independent variable is to change the dependence of the equations on that variable from a differential one to an algebraic one; thus, a four-fold Fourier transform can change the differential system ((2.1) to (2.4)) to an algebraic system in the transform domain.
The Fourier transform pair:
2.1 The classes of Maxwell’s equations 21
allows us to transform the electric field from time domain, where the appropriate field vector is ~ E(~r, t), to the frequency domain, where the appropriate field vector is ~ E(~r, ω), and viceversa.
Identical transformation can be applied to all field variables in ((2.1) to (2.4)); the appropriate symbols are listed in Table 2.2.
time-domain symbol frequency-domain symbol
Table 2.2: Symbols If eq. (2.11) and similar formulas are used in ((2.1) to (2.4)) we obtain, on equating the integrands:
(2.13)
(2.15)
(2.14)
(2.16)
similarly, the continuity equations (2.5) and (2.6) become:
(2.17) (2.18)
Observe that (2.13) to (2.18) are formally obtained from (2.1) √ to (2.6) by replacing the differential operator ∂/∂t with the multiplicative factor ω, where = −1 and ω is measured in rad/s.
Example 1. Consider, the sinusoidal (or time-harmonic) electric field 2 22
Maxwell’s Equations (2.19)
with angular velocity ω 0 and initial phase ϕ(~r ). The corresponding field in the frequency domain is obtained by substituting (2.19) into (2.11) and using the integral representation
for the one-dimensional δ-function; we find that
(2.21)
If all fields quantities vary sinusoidally with time, with angular frequency ω 0 , the electric field (2.19) may be written as:
where the asterisk means the complex conjugate, and
~ is measured in should be compared with the field E ~ of (2.21) in the frequency domain. While E ~ ~ ~ ~ times the V /m, E is measured in V s/m. The electric field E is obtained from E by multiplying E time-dependence factor e ω 0 t and taking the real part of the product, and from E ~ via the inverse transform (2.12).
When all fields quantities are written as phasors, Maxwell’s equations assume the form ((2.13)- (2.18)) and are thus indistinguishable from the equations in the frequency domain. Basically, phasors quantities and frequency-domain quantities lead to the same analytical derivations. Unless indicated otherwise, in the following we will use phasors and indicate the angular frequency with ω instead of ω 0 .
Phasor analysis is a useful tool for solving problems involving linear systems in which the excitation is a periodic time function. Many engineering problems are cast in the form of linear integro-differential equations. If the excitation or forcing function varies sinusoidally with time, the 2.1 The classes of Maxwell’s equations use of phasor notation to represent time-dependent variables allows us to convert the integrodifferential equation into a linear algebraic one with no sinusoidal functions, thereby simplifying the solution considerably. After solving for the desired variable, such as voltage or current, conversion from the phasor domain back to the time domain provides the desired result. A sinusoidal waveform hast two attributes, magnitude and phase, and thus sinusoids are natural candidates for representation by phasors. Why might such a representation be useful? One reason is mentioned already, namely it simplifies the description and solution since a complete spatial or temporal waveform is reduced to just a single point represented by the tip of a phasor’s arrow. Thus changes in the waveform are easily documented by the trajectory of the point in the complex plane. In electromagnetism, however, it is a convenient way to decouple the spatial-dependent parts and time-varying ones, and cast the problem in pure spatial-dependent functions instead.
The second reason is that it helps us to visualize how an arbitrary sinusoid maybe decomposed in the sum of a pure sine and pure cosine waveform.
The phasor technique can also be used for analyzing linear systems when the forcing function is any arbitrary (nonsinusoidal) periodic time function, such as a square wave or a sequence of pulses. By expanding the forcing function into a Fourier series of sinusoidal components, we can solve for the desired variable using the phasor analysis for each Fourier component of the forcing function separately.
According to the principle of superposition, the sum of the solutions due to all Fourier components gives the same results as one would obtain had the problem been solved entirely in the time domain without the aid of Fourier representation. The obvious advantage of the Phasor-Fourier approach is simplicity and speed, once routines are used to solve the similar parts using a PC. Moreover, in the case of nonperiodic source functions, such as a single pulse, the function can be expressed Fourier integrals, and a similar application of the principle of superposition can be used as well.
Various forms of the Fourier series description for periodic signals are based on alternate ways of writing a cosine signal. Consider
with amplitude a > 0, frequency ω > 0, and radian phase angle ϕ. (The case of negative amplitude is treated by adding π to ϕ).
23 2 24 Maxwell’s Equations
Equivalence of these expressions can be verified by using Euler’s formula,
Phasors as a e (ω t + ϕ) can be viewed as a vector at the origin of the complex plane having length a and, at any time t, angle (ω t + ϕ). The vector rotates counterclockwise with time, since ω > 0, and the projection on the real axis is described by
For graphical representation, projection on the real (horizontal) axis is inconvenient, and therefore we rotate phasors by π 2 radians and project on the vertical axis. This makes use of the mathematical relationship
The Phasor technique can be fully described as follows:
• Adopt a cosine reference, which means that the forcing function should be as a cosine, e.g.,
• Express time-dependent variables as phasors: which means that any time-varying function
x(t) can be expressed in form x(t) = <{X e (ω t }, where X is a time-independent complexvalued function called the phasor of the instantaneous function x(t).
• Recast the differential/integral equation in a phasor form. (see Table 3.1).
• Solve the phasor-domain equations.
• Find the instantaneous values by multiplying the phasor value Z for instance, by e (ω t) and
then take the real part.
2.1 The classes of Maxwell’s equations 25 The following table provides a summary of some time-domain functions and their phasor equivalents.
Table 2.3: Some time-domain functions and their phasor-domain equivalents.
Example 2. In a medium characterized by σ = 0, µ = µ 0 , 0 , and
~ calculate β and H.
Solution:
~ by making E ~ and H ~ This problem can be solved directly in time domain or using phasors. We find β and H satisfy Maxwell’s four equations.
Method 1 (time domain): Let us solve this problem the harder way- in time domain. It is evident that Gauss’s law for electric fields is satisfied; that is,
From Faraday’s law,
But 2 26
Hence,
It is readily verified that
showing that Gauss’s law for magnetic fields is satisfied. Lastly, from Ampere’s law
But
~ in (2.24) has been substituted: Thus eq. (2.25) becomes where H
~ we have Comparing this with the given E,
or
(2.25)
2.1 The classes of Maxwell’s equations 27 From eq. (2.24),
Method 2 (using phasors:)
Again
or
~ = 0 is satisfied. Notice that ∇ · H
~ in eq. (2.27) into eq. (2.28) gives Substituting H
~ in eq. (2.26), we have Comparing this with the given E or
as obtained before. From eq. (2.27),
(2.28)
2 28 Maxwell’s Equations
as obtained before. It should be obvious that working with phasors provides a considerable simplification compared with working directly in time domain. Also, notice that we have used
because the given E ~ is in sine form and not cosine. We could have used
in which case sine is expressed in terms of cosine and eq. (2.26) would be
or
and we follow the same procedure.
The Maxwell’s equation (2.1) to (2.4) are invariant under the substitutions:
(2.29)
The duality relation (2.29) may be restated by saying that Maxwell’s equations are unchanged when each electric quantity is replaced by the corresponding magnetic quantity, and viceversa.
29
First introduced by H. Hertz to depict the electric fields of an infinitesimal dipole, which has been named after him later. In his paper (2) which was translated as “The forces of electric oscillations, treated according to Maxwell’s theory”, he mentioned very briefly the concept of a field line, as the intersection of a surface of revolution and a meridian plane. He wrote:
“ For the purpose of solving the equation (EM Energy), we will limit ourselves to the special case, but the important case where the distribution of electric field is symmetrical about the z-axis, and hence at every point (P); lies in a meridian plane intersecting the z-axis, and only depends on the z-coordinate of that point, and its distance ρ from the z-axis.” After introducing his vector Π as a function of ρ and z, he defined another function Q = ρ dΠ dρ , then he added:
“The function Q is of importance. The lines in which the surface of revolution Q = const, cut the meridian planes, are lines of electric forces; the construction of the same meridian planes furnishes at every instant an immediate presentation of force distribution.”
3 30 Tracing of Fieldlines Hertz went in rather long ambiguous calculations and definitions of some quantities as m = m λ ,
Unfortunately, I never found a reference that describes his approach, or refer even to it, except of a few lines of Mönich. Though almost all the classics on Electromagnetism have a copy of Hertz’s original drawings. This is really odd.
Let’s instead try to get to the same result, using more familiar terms:
Figure 3.1: Explaining the longitudinal sectors’ approach ~ along the way of electric field E ~ for a z-oriented Let’s consider a tangential infinitesimal segment d S ~ and E ~ are Hertzian as shown in Fig. 3.1. Their cross-product of both yields a zero, as both d S parallel to each other.
3.1 The longitudinal sectors’ approach 31 Expressed the quantities in the spherical coordinates
eq. (3.2) is then can be expressed as:
This means that every term of Eq. (3.3) must be zero. The first two terms vanish due to E ϕ = 0 and dϕ = 0, however, the third term
is the interesting one which yields
or
(3.5)
Using Maxwell-Amperes curl formula, which in phasor-notion is ~ = H ϕ ~e ϕ is valid. Solving Eq. (3.7), yields the following: For a z-oriented Dipole H
(3.8)
Substituting these values in Eq. (3.4) and multiplying each side by r sin ϑ gives the following:
3 32 Tracing of Fieldlines
Let’s now express ρ H ϕ in the Time-domain
Which is the function of the surface of the revolution, we seek.
In the past, the majority of the most impressive plots of electromagnetic fields were achieved by making field lines visible. As the field vector in a couple of cases can be represented as the curl of another vector ~v with only one component, the iso-lines of the latter are the sought for lines. The classic example, even in use today, is the Hertzian Dipole. The procedure used by Hertz is applied in the making of animated films. The considered quantity is the magnetic field strength, directed in azimuth, multiplied by the radius ρ in cylindrical coordinates eq. (3.11). As the method can be used in all rotational symmetric cases, it can be also applied to longer dipoles, even with lumped loads and multiple feeds.
~ Some authors sometimes seem to forget to mention, that in reality displacement current ∂
∂t D lines are found, when ρ ∗ H ϕ iso-lines are depicted.
Changing to cartesian coordinates, arrays of several dipoles can be treated, if they are all parallel, thus producing a vector potential with only one component, say A z . The magnetic lines are given
A generalization in the sense of a three dimensional interpretation of the three examples above leads to the realization, that these lines are always given by the intersection of two surfaces. Their equations can be found in Table 3.1 3.2 A synoptical glance on classical methods of field lines representation: 33 Problem considered Reference
| Problem considered | Reference frame | depicted field ~v | surf.1 xi = k | surf.2U = k | Normal vector 1 | Normal vector 2 | Involved metric factor f | |||
|---|---|---|---|---|---|---|---|---|---|---|
| Some parallel dipoles | cartesian | H ~ | z = k | Az = k | ~ez | ∇Az /\ | ∇Az \ | 1 | ||
| Long wire arbitrary length | cylindrical | d D/dt ~ | ϕ = k | ρHϕ = k | ~eϕ | ∇(ρHϕ ) \ | ∇ρHϕ \ | ρ | ||
| spherical wave modes | spherical | E, ~ H ~ | r = k | Πr = k | ~er | ∇Πr /\ | ∇Πr \ | 1 | ||
Long wire arbitrary length
spherical
wave modes
Table 3.1: Synoptic view of the three examples above, showing fieldlines as the intersection of two surfaces.
The proof of the validity of relations in the tables can be given for all three examples simultaneously:
If the intersection of surface 1 with surface 2 is a field line, then the vectors ~v must be parallel to the cross product of their normal components ~n 1 and ~n 2 respectively.
(3.12)
Line density versus flux intensity Calculating the flux Φ, intended to be represented by the number of lines through a given area F, limited by given values of the coordinates in Fig.4.3(a) gives by using Stokes theorem:
where |d~s | = f · dx. The closed integration path contains two sections parallel to the coordinate x i and two perpendicular to it. The latter do not contribute to the result, the first run along a path of constant U ,which can be taken out of the integral:
3 34
Tracing of Fieldlines (3.14)
If the differences a − b and U 1 − U 2 are kept equal for all areas F the same amount of flux Φ is ensured within the boundaries of F . See fig. 3.2.1.
Figure 3.2: Integration of flux for a given area
35
Magnetic walls and electric walls (aka H-Walls and E-Walls) are merely explanatory mathematical tools, that help us to understand the behavior of the fields. They are so ideal, that it were a real challenge to manufacture H-Wall in reality. They are usually assumed planes or surfaces, which possess infinite dimension, and have total reflection and infinite electric conductivity as well.
Our convention states the following:
• E-Wall
– The only possible component of an Electric Field of a radiator, is the one which is normal
to the surface, while the tangential components vanishes. (See Fig. 4.1) – The only possible components of an Magnetic Field of a radiator, are the ones who are
tangential to the surface, while the normal components vanish.(See Fig. 4.2(b)) – the electric field lines experience an odd symmetry with respect to their surface of sym-
metry. – the magnetic field lines experience an even symmetry with respect to their surface of
symmetry. – Magnetic Fields induce surficial electric currents on the surface of the E-wall. (See Fig.
4.3) – The surface is normal to the axis of the radiator, if the radiator is linear, otherwise it is
parallel to it.
4 36
(a) The electric field lines of a small circular loop antenna, the arrows show an odd symmetry E-wall , H-wall and Image theory (b) The electric field lines of z-oriented λ 2 antenna, the arrows show an odd symmetry.
Figure 4.1: The electric fieldlines of a circular loop antenna 4.1(a) and a linear dipole 4.1(b)
(a) The magnetic field lines of a small circular loop antenna, the arrows show an odd symmetry to the H-Wall (b) The magnetic field lines of z-oriented λ 2 antenna, the arrows show an odd symmetry to the E-Wall.
Figure 4.2: The electric fieldlines of a circular loop antenna 4.2(a) and a linear dipole 4.2(b)
4.1 Magnetic and electric walls
(a) The E-Wall for a small circular loop antenna and the respective surficial electric currents 37 (b) The E-Wall for a z-oriented λ 2 antenna and the respective surficial electric currents Figure 4.3: The E-wall and the respective surficial electric currents.
(a) The H-Wall for a small circular loop antenna and the respective surficial electric currents (b) The H-Wall for a zoriented λ 2 antenna and the respective surficial electric currents Figure 4.4: The H-wall and the respective surficial electric currents.
4 38 E-wall , H-wall and Image theory • H-Wall
– The only possible component of an Magnetic Field of a radiator, is the one which is
normal to the surface, while the tangential components vanish. (See Fig. 4.2) – The only possible components of an Electric Field of a radiator, are the ones which are
tangential to the surface, while the normal components vanish. (See Fig. 4.2(a)) – the electric field lines experience an even symmetry with respect to their surface of sym-
metry. – the magnetic field lines experience an odd symmetry with respect to their surface of
symmetry. – Electric Fields induce surficial magnetic currents on the surface of the H-Wall. (See Fig.
4.3) – The surface is either parallel to the axis of the radiator, or contains it, if the radiator is
linear, otherwise it is normal to it.
E-Walls and H-Walls are surfaces of symmetry, where the fieldlines either start from certain points on them, or land on other points. The start and landing points are separated by lines, which we call Border Contours ‘BC’. The existence of E-Walls and H-Walls as surfaces of symmetry for the radiators, is essential if we were concerned in have closed loops as fieldlines. Due to the symmetry, every half of these loops, is provides by one half of a symmetrically built radiator, which is located at the origin. A dislocated radiator, or more generally an arbitrary oriented and located radiator, or even one with complex construction, that has no symmetry, will certainly fail to produce closed loops as fieldlines. The Image theory, offers a great aid to cope with this very situation, and help us to analyze the performance of a radiator. Every source has its counterpart image, that accounts for the reflected reflection on the surface. These images are not real sources, but imaginary ones, which when combined with the real sources, form an equivalent surface. This equivalent system gives the same radiated field above the surface of symmetry as the actual system itself, for that below this surface or within the field vanishes.
The amount of reflection is generally determined by the respective constitutive parameter below and above surface of symmetry ‘SOS’. For a perfect electric surface (E-Wall) for instance, the entire incident wave will completely reflected, as the fields there vanish. According to the boundary conditions, the tangential components of the electric field must vanish at all points along the interface. The polarization of both the incident field, and that of the reflected wave must satisfy the boundary condition, in order to assume the right direction, in which currents flow in both of 4.2 Image Theory 39 the source and its image.
(a) The E-Wall for a small circular loop antenna (b) The E-Wall for a z-oriented λ 2 antenna Figure 4.5: The E-Wall for a small circular loop antenna and a linear one.
(a) The H-Wall for a small circular loop antenna (b) The H-Wall for a z-oriented λ antenna 2 Figure 4.6: The H-Wall for a small circular loop antenna and a linear one.
This also valid for the (H-Wall). Important is, to observe the Boundary conditions as stated in section 4.1. For a single radiator located at the origin, the source and image coincide. Fig. 9.1 depicts the situation, when an E-Wall is invoked, and Fig. 9.2 depicts the situation, when an H-Wall is invoked.
Now, if the center of the radiator is displaced by a distance z 0 above the surface, an image at the same opposite distance will be produced, the orientation of the currents of both source and image
4 40 E-wall , H-wall and Image theory is shown in Fig. 9.3 showing the E-Wall.
Figure 4.7: E-Wall, source and image for a small circular loop antenna and a linear one.
Figure 4.8: H-Wall, source and image for a small circular loop antenna and a linear one.
4.2 Image Theory 41 Fig. 9.4 shows the H-Wall.
Eventually Fig. 9.5 depicts the situation for an arbitrary oriented dipole.
Radiators that have more complex construction like these, we just handled, could be like-wisely turn into equivalent system of sources and images, just by following the same procedure as we mentioned before.
43
They are vector functions, (1) who vary with respect to both spatial position ~r and time t. They thus can in free space be expressed as
where,
When t varies, each of these vectors above traces out a curve, in general a space curve which varies in magnitude and direction, such a curve is called a hodograph.
~ ˙ For an instance of time t i , every point p i along the trajectories of the numerical integration of D ~ ˙ with respect to the time, can be express as and B
5 44 From EM Fields to fieldlines Here we use a mathematical trick, by discretizing the time interval t, into chosen subintervals, each of which, is fixed-valued, so Eq. (5.2) can be expressed in terms of the coordinates axes only, and one yet can trace the “time-varying factor” out.
(5.3)
Figure 5.1: The field lines and their symmetry surfaces Plotting a series of many field lines, representing different ‘start points’ taken for consecutive subintervals each, will do the trick, to visualize quantities which vary with respect to ~r and t. The difference between Eq. (5.2) and Eq. (5.3) is obvious, the former requires a 4-orthogonal coordinate system, which is not invented yet, while the latter is very content with the familiar 3- orthogonal coordinate system. In addition the discretization is the core of all numerics, simply because the computer has only finite storage capacity, that will have hard times handling continuous intervals, either in space or time, and hence there is no way to represent continuous data, except approximately as a sequence of discrete values. Using the rectangular coordinate system as a standard, has despite the inconvenience, the advantage to sum the vector fields of many radiators ~ ˙ and up, besides its mostly the standard coordinate system in the simulation programs. Writing D
45 ~ ˙ in terms of their components, we come to a situation as shown in Fig. 9.2. B
Note The quantities in Eq. (5.4) lead their counterparts in Eq. (5.5) by Phase displacement of π 2 due to time-harmonic current distribution assumption.
Note Fig. 4.1 depicts the field lines and ‘SOS’ for a linear z-oriented Dipole located at origin, case we have a small xy-loop antenna centered likewise, ‘IES’ and ‘IMS’ interchange their positions, and so also the ‘Electric’ and ‘Magnetic’ field lines, and their corresponding ‘Boundary contours’, ~ ˙ → B ~ ˙ and vice versa. Consider Fig. 4.1(a) and its counterpart Fig. 5.1(b), there an infinitesi- i.e., D ~ ˙ and d B ~ ˙ between two adjacent points p 1 and p 2 along their corresponding loops mal segment d D ~ ˙ are respectively labeled. In terms of their components, they are ~ ˙ and B D
can be rewritten as dot product (·) of their gradients respectively, and a displacement vector d~r =
and
(5.9)
5 46 From EM Fields to fieldlines
The Equations (5.9) through (5.14) are functions of functions, each of which is a 4-dimensional func~ ˙ and B ~ ˙ respectively. These could tion, however they all have their the dimensions of their fields, D be used whenever the field intensity is concerned. So that every field line can be associated with a certain strength, or intensity. To yield another sort of functions, however, i.e., functions with arguments in terms of x, y, z and t, we need to turn the equations (5.9) through (5.14) into dimensionless ones. The trick is done by dividing them each by their corresponding vector-field magnitude. Using another script, to show the new unit vectors:
(5.15)
(5.16)
47 where their corresponding components are
(5.17)
Equivalently we have
(5.18)
These six quantities represent four dimensional differential equation each, these have to be solved simultaneously, for the values of x, y, z and t, where as we mentioned above only fixed values of t are admissible, the algorithm form the values for the points on the integration path in form of p i = p(x, y, z) | t=const , accumulate then in lists, coded them, then use them in plotting.
49
Although start points must be found in magnetic and electric walls as well, the following explanation shall be given for the perfectly conducting plane. Here the surface currents do indeed exist, in the other planes they have to be derived from the field. So the more familiar scenario is invoked for purposes of explanation.
Currents in perfectly conducting surfaces are two dimensional vector fields, here in the x,y-plane. In a time harmonic field each of the two components J x and J y vanishes along a curve for a given time. The intersections of both curves define points of zero current density, called PZD in the following.
Under normal conditions the magnetic field above metallic surfaces shows an elliptical polarization. Due to J ~ = ~n × ~ H, this holds true for the current density J ~ too. Along certain lines - the lines of linear polarization - the polarization ellipses degenerate into straight lines. So the current can vanish only at these lines, PZD can be found only here.
Now the border contour ‘BC’ between starting and landing regions can be determined. Keeping in mind the complex formulation jωD z , the time dependent formulation is given by:
6 50 Start points for line tracing So vanishing values for an arbitrary time t 1 are given by:
(6.2)
In the past, the lines of constant amplitude and constant phase were widely used to characterize the time harmonic electric field above metallic surfaces. They could be easily measured with a monopole probe and a vector voltmeter.
(6.3)
As the arctan function is normally defined over [−π · · · + π], however, the phase angle π has to be added or subtracted depending on the quadrant.
Comparing the equations (6.3) and (6.2) shows that the BC for a given time t 1 can be referred to the phase front:
(6.4)
The tracing of the BC’s can be started from any arbitrary point ‘Start’ in the ‘SOS’ Fig. 6.1, where Ḋ z = 0 holds true. Integrating the magnetic field along the BC with the coordinate s yields the current I, crossing the BC between the start and any other point of the path.
Figure 6.1: Border contours with the region of interest and current lines 6.1 Topology of surface and displacement currents 51 (a) Sector, subdivided into areas launching the same displacement current (b) Fieldlines, embedded in surfaces of constant I or constant ζ, forming horizontal or vertical families (c) Field line, appearing or disappearing at a BC Figure 6.2: Regions of interest, border contours, sectors and field lines
Now the area of interest can be indicated by I max .
Segmentation of the BC between ’Start’ and ’Stop’ into N parts with the same amount of crossing currents ∆I = I max /N divides the total current into equal portions ∆I.
Now the flow lines of Fig.4.3(b) can be started, all parallel to the surface current density vector ~ J.
6 52 Start points for line tracing They usually end in a PZD. The current i enclosed between two lines is given by:
The conduction current I gradually becomes the displacement current, stemming from the plane. The continuity equation requires that the amount of displacement current, launching from the dashed area w sin α∆η in Fig.6.2(a) is equal to the difference between the currents i entering at η 1 and leaving at η 2 :
Eq.(6.8) shows, that the maximum of i occurs in the BC, where Ḋ z = 0.
The next task is to subdivide the area 1, 2, PZD in Fig.6.2(a) into pieces, from which the same amount of displacement current arises. A helpful tool to this end can be a second vector field ~j in the plane, which is parallel to J ~ everywhere, but divergence-less,
Now a new fictive current ĩ, crossing the lines η = const. can be referred to
which does not depend on η.
Extending the term in brackets in eq.(6.8) by the function g defined in eq. (6.9) yields after differentiation:
~ sin αw is independent of η. The first differentiation on the right side leads to zero, as ĩ = g| J| Together with eq. (6.8):
6.1 Topology of surface and displacement currents 53 Rearranging makes a direct integration on either side possible:
(6.10)
Now a re-scaling of the η coordinate can be done in such a way, that the remaining current i(η) = ∆I/g in the sector is normalized ∆I.
Starting from the BC and keeping in mind that g(η = 0) = 1 results in
with the new scaling ζ from 0 at the BC to 1 in the PZD.
Up to now the area of interest has been defined and divided into sectors, starting from a BC. Another division of these sectors into subareas, each launching the same displacement current, led to the start points for field lines. These could be placed in the middle of the subareas or, like done here, in the intersection points of the coordinate lines of I and ζ of the new skew coordinate frame for the area of interest.
Once the displacement currents lines have been calculated, the values of I and ζ can be transferred for all points in the volume of interest above the plane. So I = const. and ζ = const. in Fig.6.2(b) describe bunches of surfaces, in which lines are embedded. So the lines can be assumed again as an intersection of two surfaces, as it has been done in the 2 nd section for the simple cases. So each line can be associated to the values of I and ζ of its start point. On the other hand, two different groups or families can be distinguished: lines, starting from points with constant I lead to a group of lines, we call horizontal family. Starting from points with constant ζ, leads to a vertical family.
The basis of an animated film are frames, i.e. field plots following one another. The question therefore arises as to a line can be referred to another in a plot for an earlier as well as later time. To answer this question, one should regard the ‘birth‘ or ‘death‘ of a line, the limits of its ‘life‘, 6 54 Start points for line tracing which can occur on BCs only. In Fig.6.2(c) a field-line in the vicinity of a BC is shown at three
Consider a BC at two different times, the start time t 0 and an arbitrary time whereby t > t 0 . The results for I max calculated with eq. (6.6) for both times will be different most probably. Proceeding as described above would lead to different ∆I. The pieces however, in which the sector 1,2,PZD has to be divided, must launch the same amount of displacement current at both times. So the remaining current i(η) must not be normalized w.r.t. different ∆I, but only w.r.t. ∆I(t 0 ) for both times t 0 and t, whereas the number N has to be held constant.
As the lines appear or disappear on the BC, a new scale ζ̃ instead of ζ appears expedient, beginning in the PZD with ζ̃ = 0 and ending at the BC.
Keeping the number N for both times is easy to enforce for a closed BC. The maximum of ζ̃,
It is worth to formulate the new surface definition to the rotational symmetric case of a straight filament current with a symmetry to the xy-plane . With
one gets an integration path with constant ρ the path element ds = ρdφ and constant values of H φ . So ρH φ can be taken out of the integral.
where the index 0 denotes a radius where Ḋ z = 0. So, planes of constant φ are the first group.
6.1 Topology of surface and displacement currents 55 The second is given by constant values of ζ or, closer to the original calculation, of g, given by eq. (6.10). With Ḋ z , | ~ J| and the path element dη in the cylindrical coordinates one obtains
due to symmetry. So eq. (6.12) becomes:
57
Dr. Zuhrt H. (5) in his book “Electromagnetic radiation fields”, Springer 1953”, seems to be the only one, sofar, who made an attempt to deal the general case of orientation, he described a rather simple method to derive a formula, unfortunately, that method covers but the Far-Field of a horizontally oriented Hertzian.
Though we are grateful to his revolutionary advance, indeed, it was not at all satisfactory, first for that he considered the far-field case only, and finally, that he restricted it on just a Hertzian. The Hertzian Dipole, though of some scientific interest to understand the behavior of more complicated antennas, is by no means a real antenna that could exist for its own. It is nothing more than a piece of uniform current of infinitesimal length, that is. This chapter presents a technique used to derive the entire field formulas for both arbitrary oriented Hertzian and finite length linear Dipoles, and visualize the electric and magnetic field lines, afterwards. The arbitrary oriented half-wave linear dipole though, will be taken as a standard example for convenience.
One encounters such a dipole in all imaginable orientation, in almost every aspect, that has anything to do with broadcasting, transmitting or receiving signals, either in the free space or under the sea. Sometimes, one finds such a dipole in a stationary variant, like being a part of a complicated radiating object, or also in a transient variant, that change its orientation continuously, just consider a jet or a submarine, tracing certain courses. Now, as it is obviously clear, that such a dipole is significant, how comes that everybody left it over? The reason, obviously, is the difficulty to derive its formulas in the conventional manner, besides, there is always alternatives for everything in life, an arbitrary oriented dipole antenna, can also be replaced be weighted dipoles, that one is accustomed with their formulas and radiation, with the suitable combination of these, one can get the equivalent of an arbitrary oriented. It works, but whence it gets to closed field lines visualization, the entire process will certainly collapse, simply because having an equivalent out of three other dipoles oriented in the direction of the three axes, has no symmetry plane, and thus the enclosure of these lines won’t occur. For this very reason, we were very concerned in deriving 7 58 On the arbitrarily oriented radiators formulas belong to a real arbitrary oriented, not its substitute.
Four years ago, we came out with a simple and very efficient solution that deals with dipoles, which are oriented in the other axes rather the z-axis. Astonishingly this process worked fine with all sort of the dipoles considered, i.e., they do not have to be of linear form. The solution was to develop some replacement rules, applied to the formulas of a Z-oriented antenna, to get in no time new formulas for the desired orientation to either of the other rectangular coordinates, that is, X or Y.
How does it really work?
Let’s explain the process briefly right here. We will furthermore adopt mathematica notation to express, vectors, sets and replacement rules, to keep as close as possible, to the items used in the programs used to visualize the results. Consider first, the following three sequences of the
they are the rotations of each other, in a positive sense of the right hand rule. Fig.7.1 on page 58
Figure 7.1: Rotating Coordinates Now, notice the last symbol of each sequence, and compare this to its own axes in figure 1, amazing, all we need is, in the axes in Fig. 7.1(b), z replaced by X, X by Y, and Y by Z, and in Fig.
7.1 A provoked solution: Replacement Rules 7.1(c) to the right Z replaced by Y, Y by X, and X by Z, right? And exactly these our replacement rules, that we apply to the formulas of a Z-oriented dipole to obtain either an X orientation or a Y orientation, by applying the set of rules belongs to the desired orientation. In other words: To derive formulas for an X oriented antenna from, old formulas derived for a
¬ Replace all X’s in an expression by Y’s or in mathematica notation x : → y
Replace all Y’s in an expression by Z’s or in mathematica notation y : → z
® Replace all Z’s in an expression by X’s or in mathematica notation z : → x Or just follow the following steps instead to derive formulas for an y-oriented antenna from, old formulas derived for a z-oriented one, do the following:
¬ Replace all X’s in an expression by Z’s or in mathematica notation x : → z
Replace all Y’s in an expression by X’s or in mathematica notation y : → x
® Replace all Z’s in an expression by Y’s or in mathematica notation z : → y Quite elementary, isn’t? We though, decided to be completely sure of the correctness of these rules, so we derived, in conventional method, the formulas for the X and Y oriented antennas, and the comparison showed that the results were identical, a matter that, we never doubted actually. The technique, in reality, based on simple and elegant concept, that is, always create a new coordinates system, where the dipole antenna keeps its Z orientation, but in terms of the new belonging to coordinates.
In other words, we simulate a new up-down world for the dipole in which, only the dipole itself, believes that it were Z oriented, to the old coordinates system. That’s true, case it take the role of an X oriented, the X axis of {Y, Z, X} system coincides with the Z axis of {X, Y, Z} system, and the dipole will radiate further, without noticing the drastic changes, of the orientation.
However, though these rules were very effective in time and efforts’ sparing, they were not very satisfactory to us; we wanted to extend them further to cover real arbitrary cases, many methods were considered, and finally we decided for Euler’s rotation theorem, it is the most straightforward approach to our method, and it is besides its simplicity to be understood and programmed a very effective tool too.
59
7 60
On the arbitrarily oriented radiators
According to Euler’s theorem of rotation, any rotation can be described using three angles, namely a triple (ϕ, ϑ, ψ). Each of these rotations can be written and described in terms a rotation matrix, and then a general matrix A, which is the multiplication of all the three matrices, can be constructed. Unfortunately, there are several conventions for Euler angles, depending on the axes about which the rotations are carried out. One can obtain 12 different rotation sequences; only three of them are widely used, in graphics and computer games or simulators.
The so-called “x-convention”, is the most common definition. In this convention, the rotation given by Euler angles (ϕ, ϑ, ψ), where the first rotation is by an angle ϕ about the z-axis, the second is by an angle ϑ about the x-axis, and the third is by an angle ψ about the z-axis (again).
Note, however, that several notational conventions for the angles are in common use. Goldstein (7)(1960, pp. 145-148) and Landau and Lifschitz (1976) (8)use (ϕ, ϑ, ψ), Tuma (1974) (11) says (ψ, ϑ, ϕ) is used in aeronautical engineering in the analysis of space vehicles (but claims that (ϕ, ϑ, ψ) is used in the analysis of gyroscopic motion), while Bate et al. (1971) use (Ω, ι, ω) . Goldstein remarks that continental authors usually use, and warns that left-handed coordinate systems are also in occasional use (ψ, ϑ, ϑ)(Osgood 1937 (10), Margenau (9)). Varshalovich (12)(1988, pp. 21-23) uses the notation (α, β, γ) or (α 0 , β 0 , γ 0 ) to denote the Euler angles, and gives three different angle conventions, none of which corresponds to the x-convention. No matter, we will describe and use, the x- convention, a.k.a ZXZ convention. Let’s call the first rotation about the Z axis R z (ϕ), and the second rotation about the rotated Xaxis(X’) R x (ϑ), and finally the third and last rotation about the rotated Z-axis(Z’) R z (ψ). The definition of R z (ϕ) is given by:
(7.1)
The primed symbols indicate that a rotation has taken place. The next rotation R x (ϑ) that follows the first one immediately, is about the rotated X-axis (X’)by an angle ϑ and given by:
Fig.7.4 depicts its effect:
62 7 On the arbitrarily oriented radiators
The last rotation that follow is about the rotated Z-axis (Z’) by an angle ψ and given by R x (ψ) Rotation Matrix:
Its final effect results in the desired orientation of the dipole, as the depicted in Fig.7.5. Note the double prime of the symbols now. Now we have got a dipole, which orientation coincides with Z”, which match exactly the desired
The three successive rotations can be expressed in terms of one general rotation matrix A, which 7.2 Rotation Matrices 63 is their multiplication, as given by Equation 4.
This can be written in terms of a matrix
(7.5)
where
Rotation matrices are orthogonal matrices, i.e. unit matrices (which are matrices for which M · M T = I) with real coefficients. The module of the determinant of unit matrices is one, among the orthogonal 3 × 3 matrices, only the ones having a positive determinant (+1) are rotation matrices. Notice we used the transpose of a matrix instead of its inverse, simply because, their determinant were the unity. One notes that our replacement-rules were, intuitively and implicitly based on Euler’s angles
principal. For a Z-oriented
Dipole, we set ϕ = ϑ = ψ = 0. The A matrix attributed to this
The transpose of a matrix A is a matrix formed from A by interchanging the rows and columns such that row i of matrix Abecomes column i of the transposed matrix. The transpose of A is 7 64 On the arbitrarily oriented radiators denoted by A T . Each element a ij in A becomes element a ji in A T . The following represents A T the transpose of A:
(7.9)
Rotating a dipole implies rotating its field along with its orientation. Multiplying the rotation
! x ~ = y , will invoke the matrix A by the vector, containing the original coordinates, say X
together different order: H x 0 0 , H y 0 0 , H z 0 0 refer to their corresponding new coordinates,{X 0 , Y 0 , Z 0 } they are by no means, components of the old world of {X, Y, Z}.
Is that a real problem? Certainly not to a human, but for a machine that is programmed, to deal matters sequentially, it is a serious problem. Quite often, one has a sophisticated software written to simulate a real 3D world, like ours, and its by convention that x-component been placed first, followed by y-component and finally z-component, to visualize the fields lines belong to the new orientation, we must rearrange the component within the new vector, that they match the old order. Otherwise, the visualization will make no difference from those the z-oriented one.
A correction The multiplication by A T the transpose of A corrects this behavior and produces a new vector
right here reserved for the corrected fields.
The next installment To Formulas deals the matter of more details, and finally gives expressions for the arbitrary oriented dipole for both a Hertzian and finite length dipole.
7.3
To Formulas 65
Now, in the XYZ system, a space vector ~r from the origin to the point of observation P in the space can be expressed as
and its magnitude is
in the X’Y’Z’ system (see Fig. 7.6), however, the very same vector is defined by
of orientation, we are making a convention by reserving the unprimed symbols for the case of a Z-oriented cases, primed otherwise. ~ = We just agreed that original non-rotated coordinates can be expressed in a vector notation as X
~ 0 = Y 0 , so that every point in space P ~ (x, y, z) can be expressed in Y , the rotated one is X
both coordinates systems components equivalently as
Table # 7.1 summarizes the entire replacement rules used in this work. We will show a sign of 7 66 On the arbitrarily oriented radiators courtesy to mathematica, that helped us doing the visualizations though, by adopting its notation, so the symbol slash dot (/·) will be used as a replacement operator.
In X Y Z System
In X’ Y’ Z’ System
Table 7.1: Replacement rules Applying a set of rules, let’s say {x :→ x 0 y :→ y 0 z :→ z 0 } on either an electric field vector
new field component is indeed a combination of all the three old components, weighted by the corresponding component of the rotation matrix A.
Magnetic field components
Electric field components
Table 7.2: The rotated fields components The next example will show that these components, though completely correct symbolically, their order won’t satisfy a machine, a robot for instance, or another program, that used to work on sets of input data, which come in a certain order.
7.3 To Formulas 67 An example To keep the discussion as simple and obvious as possible, we will consider a magnetic field vector of a Z-oriented linear dipole regardless of its sort. It has no component along its
extract for an X-oriented dipole, and
expressions
alright, for that the H x vanishes, but the order of the components
is not
satisfactory.
Now its been done!
Yet another example Now, we want to extract a Y-oriented Dipole out
of the wellknown Z-oriented formulas, and
The coordinates of any point P (x, y, z) in space can expressed in terms of newly generated coor-
7 68 On the arbitrarily oriented radiators is, once more, indeed expressions for an Y-Oriented dipole alright, for that the H y vanishes, but the order of the components is not satisfactory.
To cover the case of dealing an arbitrary oriented linear finite length radiating object as well, the distance from each tip (please see Fig. # 7.7) has got to be defined as in table 7.3 on page # 68.
Figure 7.7: Dimensions for a finite length linear dipole!
Z-oriented
Arbitrary Oriented
Table 7.3: two more dimensions for a finite length 7.3
To Formulas 69
magnetic fields
The amazing matter, that a new component H z has appeared right from a blue sky! Remember no dipole has a magnetic component along it axis.
In analogous way we go ahead to derive the corresponding expressions of the electric fields:
Electric Fields
(7.16)
(7.17)
7 70
On the arbitrarily oriented radiators
Magnetic Fields
Electric Fields
(7.22)
7.3 To Formulas
71
(7.24)
73
NUMERICAL solution of ordinary differential equations is the most important technique in continuous time dynamics. Since most ordinary differential equations are not soluble analytically, numerical integration is the only way to obtain information about the trajectory. Many different methods have been proposed and used in an attempt to solve accurately various types of ordinary differential equations. However there are a handful of methods known and used universally (i.e., Runge-Kutta, Adams–Bashforth–Moulton and Backward Differentiation Formulae methods). All these discretize the differential system to produce a difference equation or map. The methods obtain different maps from the same differential equation, but they have the same aim; that the dynamics of the map should correspond closely to the dynamics of the differential equation. From the Runge–Kutta family of algorithms come arguably the most well-known and used methods for numerical integration.
Differential equation is an equation involving an unknown function and one or more of its derivatives. The equation is an ordinary differential equation (ODE) if the unknown function depends on only one independent variable.
As for numerical methods, bear the following items in mind:
( Solution of a DE is function in infinite-dimensional space.
( Numerical solutions of a DE is based on finite-dimensional approximation.
( Differential equation is replaced with algebraic equation, called Difference equation whose
approximates that of given DE.
( There is no method, which can deliver the best results, all times.
With a differential equation, we can associate initial conditions or boundary conditions, auxiliary 8 74 Numerical Approach conditions on the unknown function and its derivatives. If these conditions are specified at a single value of the independent variable, usually at start points, they are referred to as initial conditions and the combination of the differential equation and an appropriate number of initial conditions is called an initial value problem (IVP). If these conditions are specified at more than one value of the independent variable,usually at start and end points, they are referred to as boundary conditions and the combination of the differential equation and the boundary conditions is called a boundary value problem (BVP). Partial differential equations are beyond the scope of this text, instead we shall take a brief look at some methods for solving the single first-order ordinary differential equation, written in the standard form:
we expand the scope of eq (8.1) to cover a system of differential equations in section 8.2.
A minimal requirement The functions we shall be considering here when discussing the implemented numerical methods one by one, are assumed to be continuous, since this is a minimal requirement for predictable behavior. The set of formulas, in hand may have derivatives of all orders as well C ∞ (X), and non-stiff equations as well, so it is just right to adopt methods used to solve the IVPs, to solve them numerically.
Functions that are not continuous can skip over point of interest, which is not satisfactory trait when attempting to approximate a solution. More sophisticated understanding of the nature of a problem, leads to a better formulation, in our Test Problem in section 8.2.1 on page 86, we tried to give a comprehensive description of the Kepler’s Problem, which has been chosen, due the fact, that it is very sensitive to the initial conditions, therefore, once the used algorithms succeed there, they will also succeed in solving much well-posed set of formulas.
A well-posed Problem The term well-posed problem stems from a definition given by French mathematician J. S. Hadamard. He believed that mathematical models of physical phenomena should have the properties that:
A good method possesses the following three characteristics:
Ê A solution exists.
Ë The solution is unique ( that is, there is one and only one solution to the differential equa-
tion).
75 Ì The solution depends continuously on the data, in some reasonable topology. Topology is
the study or science of places.
If a problem is well-posed, then it stands a good chance of solution on a computer using a stable algorithm. If it is not well-posed, it needs to be re-formulated for numerical treatment. Typically this involves including additional assumptions, such as smoothness of solution. This process is known as regularization, which means making a function regular or smooth. A smooth function is one that is infinitely differentiable, i.e., has derivatives of all finite orders.
In the course of carrying out a mathematical computation, one has to deal with problem of errors. There are three ways in which errors can enter a calculation.
Sources of error The main sources of error in obtaining numerical solutions to mathematical problems are:
1 The model: its construction usually involves simplifications and omissions;
1 The data: there may be errors in measuring or estimating values;
1 The numerical method : generally based on some approximation (truncation);
1 The representation of numbers: (roundoff) for example, π cannot be represented exactly by
a finite number of digits;
1 The arithmetic: frequently errors are introduced in carrying out operations such as addition
(+) and multiplication (∗) .
Each of these errors is unavoidable in a calculation. Hence, the “problem of errors” is not one of preventing their occurrence. Instead, it is one of controlling their size in order to obtain a final result that is as accurate as possible. This process is called error control. It is concerned with the propagation of errors throughout a computation. For example, we want to be sure that the error that results from performing an arithmetic operation on two or more numbers, which are themselves in error, is within tolerable limits. In addition, the propagation or cumulative effects, of that error in subsequent calculations should also be kept under control.
In this regard Picard theorem is very useful (See page 154).
8 76
Numerical Approach
To quantify the order of accuracy of a numerical approximation the Taylor expansion of y(x) ∈ C r
where ξ ∈ [x n , x n+1 ] and substitute this into the approximate solution x n−1 .
Let us assume that the f ∈ C r , i.e. of the same order of the solution. Definition 1. An approximate method is r-th order accurate if the term in h r in the Taylor expansion of the unknown is correctly reproduced: it is indicated by O(h r ).
The knowledge of the order r is of importance, because it allows us to tell by how much the results are improved when the step is reduced. In general, one prefers methods with a large , since then a reduction of h promises a large gain in accuracy.
Note: The order r depends only on the method and not on the differential equation, provided f satisfies the assumption. Definition 2. A numerical method is convergent if
for all x over the finite interval [x 0 , x N ], i.e. if the sequence of improved values converges to the true value of x. A method which is not convergent is said to be divergent.
We now study what conditions a numerical method must satisfy, if it is to be convergent. One would expect that one such condition would be that it has to be sufficiently accurate representation of the differential system.
The value of the residual E n+k is considered as measure of accuracy:
The leading order deviation is called local truncation error at the n-th step
and the local unit truncation error is given by
77 where it is considered the Cauchy problem
(8.5)
starting at the exact value at time x n , i.e. it is calculated with single-step of integration at the
Definition 3. The total (global) truncation error is the difference between the solution y(x n+1 ) and y n+1 (the solution calculated after n + 1 steps):
The definition of convergence can be formulated using the total error: a numerical method is convergent if
It is natural to expect that the error will accumulate steadily as an integration proceeds, but that is not so. The stability of the problem can cause errors to decrease as well as to increase.
first thought on the appropriate level of accuracy, that might be needed for convergence, is that E n+1 → 0 as h → 0.
In-depth remarks show that this is not going to be enough. The appropriate level is to demand that θ n+k → 0 as h → 0.
First of all some important properties for a numerical method are given:
Definition 4. A numerical method is called consistent ( with the differential system ) if the local truncation error
Although convergence implies consistency, the converse is not true. It can happen that the difference system, produced by applying a numerical method, suffers an in-built instability which persists even in the limit as h → 0 and prevents convergence.
Definition 5. A method is said to be stable if a small deviation from the true solution does not tend to grow as the solution is iterated.
Let us consider a form of stability, which is concerned with the stability of the difference system in the limit
discretized problem (i.e. difference equation generated by the method) and let {y n , n = 0, 1, · · · , N } and 8 78 Numerical Approach {y n ∗ , n = 0, 1, · · · , N } be the resulting perturbed solutions. Then if there exist positive constant S and h 0
then the method is said to be zero-stable.
It is a characteristic of non zero-stable methods, that decreasing the step-size actually makes matters worse, i.e. the error grows at an increasing pace. Any error due to discretization and round-off could be interpreted as being equivalent to perturbing the problem. The zero-stability is thus a requirement that the difference system be likewise insensitive to perturbations. If the inequality (8.9) is not satisfied, then no acceptable solution will be produced.
One may ask, what bother with numerical methods, IVPs and stuff ? Especially, whence we already, have the solutions in form of neat mathematical formulas?
Capturing a fieldline means tracing it as close as possible, as if it were a route of course which we strictly follow. On the contrary of fields, no solid evidence was ever delivered either to affirm or even to exclude the existence of ‘fieldlines’ as were proposed by Faraday. Taking refuge in mathematics and observing ‘mother nature’ while performing varieties of phenomena of all kinds, may imply that existence of such ‘a thing’ is likely to be true. We know for sure that fields have certain values and directions all over the space, who surrounds their origins. These field are described by Differential Equations, which once solvable possess at least one solution each, which can be graphically represented. Otherwise a family of solutions or no solution at all. These solutions can be drafted by curves, no matter whether they are plane or space curves, they reveal us a great deal information about the behavior of these solutions by merely just a glance. Therefore, taking the existence of the ‘fieldlines’ for possible is really advantageous. By Means of mathematics, Hertz was able to trace the fieldlines of an ideal theoretical dipole, which has been an event in his times, these fieldlines of his, turn out to form ‘closed loops’, pertaining certain direction every half a period. Accordingly, bearing that in mind left the way open to conclude that such unseen lines which connect certain points in the domain of an oscillating radiator could be much more than ‘just likely to occur’, they could as good as undiscovered yet fact.
The reason is quite simple; namely, the analytical solutions to DE’s (once available) deal with continuous functions; i.e. functions that depend continuously on the independent variables. A computer, however, has only finite storage capacity, and hence there is no way to represent continuous 8.1 Why to undertake a numerical approach?
data, except approximately as a sequence of discrete values, this what the process of discretization all about.
The best passage to is to pave ground by a pictorial mean, this is, to show what our visualization is all about, and thereafter whether an analytical solution could or could not possibly cope?
Obviously we are solely interested in visualizing the propagated fieldlines radiated by either a dipole antenna, or an array of antennas regardless of the location, orientation and configuration, our algorithm must deliver good results.
The methodology of implementation can be summarized by means of a simple flow chart, as shown in Fig.8.1.
This is known as Divide and Conquer Strategy. It consists in breaking a problem into simpler subproblems, usually of similar type, next to solve these subproblems, possibly in parallel once desired , and finally to amalgamate the obtained results into a solution to the problem. Thus the strategy displays two aspects; the first one breaks the problem into subproblems, the second one merges the partial results into the total results. Though, this sounds rather elementary and evidently, the process of ‘Dividing’ and ‘Merging’ could turn into real cumbersome and tedious matter, if the choice was wrong.
Figure 8.1: A block diagram showing a rough configuration of ‘Divide and Conquer Scheme’ Summary: The task is divided in three subtasks, (see Fig.8.1, each of which processes and manipulates certain expressions, performs numerical integration operation, labels the outcomes, and stores them in form of sets of sets of numerical data, each of which is a 3-dimensional representation of a point in space. These data, thereafter deposited in local storages, where only one provider 79
8 80 Numerical Approach has access to a certain storage, where the labeled and packed data are ready to be delivered to the next stage, to start working. Each stage is activated, once it obtains the required data, and the order to start.
Order of tasks Center of Information
Analysis & Classification Electric Fields Magnetic Fields Expert System Divide & conquer
Scheme Numerical Processing
Unit Display Unit Figure 8.2: The processing of an order of tasks form In our case, we are interested in a ‘genetic model’ that stands for a family of programs, which are capable to solve families of various problems, Fig.8.2. This model itself is based on three subprograms, which adopt Centralized Control/Decentralized Execution strategy in addition, in hope, it is a valid tenet of excellent performance, to ensure a frictionless flow of information, and to reduce the possibility of conflicts between the active processes. Fig.8.2 through Fig.8.5.
8.1 Why to undertake a numerical approach?
Specify the Grouping of the Field Lines and their inducing
radiators.
81 Specify the Sort of
the Field Lines and their inducing
radiators.
Identify, Label,
& Check up
Unit Magnetic Time Instance? Interval?
Currents Phasing Intensity Electric Radiators Number Grouping Location Orientation Field Lines Number Distribution Grouping SOS
& Arrows Number Grouping Location Orientation Expert System Figure 8.3: Center of Information As we mentioned, these three subprograms work on data delivered to them by special providers, once called, they manipulate these data locally, label and encoded them, and send the encoded labeled data to next stage for further manipulation, and their own to the ‘display unit’.
8 82 Numerical Approach Sets of rules Data Bank Basic Formulas Constructors
& Modifiers Modified formulas Data Bank Graphical Symbols Modified graphics Formulas processing
| Unit | Unit | ||
|---|---|---|---|
| Vectors | Starting | Vectors | |
| Points | |||
| Boundary | Finder | Surficial | Field |
| Contours | Currents | Lines | |
& Conquer Scheme Numerical Processing Unit Display Unit Figure 8.4: The Expert System The core of such a model is an expert system, which possesses a “Data Bank” and “Sets of Rules”, the data bank contains basic symbolic formulas for a variety of radiators, usually located at origin, and also auxiliary graphical symbols to represent each radiator and the surfaces of symmetry ‘SOS’. Once an ‘order of tasks’ is made, this activates the expert-system, such an ‘order of tasks’ contains information about the radiators, the fields to be visualized and the interval of time; as for the radiators, the number of the radiators, their types, current distribution, location, and orientation.
8.1 Why to undertake a numerical approach?
Expert System 83 Divide & Conquer
Scheme Numerical
Unit Splitter
& Selector Algorithms Providers of BC Starting Points
BC Determination Sets of Points Labeling & Packing Unit Storage Surficial Currents Determination Providers
Sets of Points Labeling & Packing Unit Storage Sets of Points Display Unit Providers Field Lines Determination Display Unit Display Unit Figure 8.5: The Interaction among ‘Divide and Conquer’, Expert System, and Display Unit The fields are characterized by their sorts, ‘electric’, ‘magnetic’ or ‘both’, their number and distribution within a contour, and the number of these contours. Eventually one has the choice between ‘Snapshot visualization’ and ‘clip visualization’, the former considers the field lines distribution in a certain instance of time, while the latter over selected points of an interval. These are all used 8 84 Numerical Approach as inputs in the header of the main program. The ‘Expert System’ checks these inputs, abstracts the required formulas, constructs their images, both symbolically and graphically, and stores the final expressions in a local warehouse. In this warehouse, the process of computing three distinguished ‘unit vectors’ for the contours, surficial currents in the ‘SOS’ and eventually the field lines, is accomplished. Determining the start points for the boundary contours on the ‘SOS’ follows in a manner matching the requirements of their desired number (spatially) and occurrence (time varying). These points get labeled and packed into a special package, and given free to a certain provider, which deliver them to one of the three subprograms, namely the first one. Having the unit vectors as differential equations and the start points as initial values, the program solves them numerically, and constructs the results, into packages of sets of sets of 3-dimensional x, y, z points, which can be first sent to the display facility, and afterwards manipulated to choose start points among them, for the next stage namely the surficial currents on the ‘SOS’. The second subprogram starts where the first has ceased, and after drafting the required currents, it delivers start points for the final stage, namely the plotting of field lines.
The terminology used in (8.1) is misleadingly simple, because y could be a vector-valued function. Thus, if we are working in R N , and x is permitted to take on any real value, then the domain and range of the function f which defines a differential equation and the solution to this equation are given by
Since we might be interested in time values that lie only in some interval [a, b], we sometimes
specific problems, it is often convenient to focus, not on the vector-valued function f and y, but on individual components. Thus, instead of writing a differential equation system in the form (8.1), we can write coupled equations for the individual components:
A differential equation for which f is a function of y only, is said to be ‘autonomous.’ Some equations arising in physical modeling are more naturally expressed in one form or the other, but we
8.2 A system of differential equations 85 ‘emphasize’ that is always possible to write a ‘non-autonomous’ equation in an equivalent ‘autonomous’ form. All we need to do to change the formulation is to introduce an additional component y N +1 into the y vector, and ensure that this can always maintain the same value as x, by 0 associating it with the differential equation y
N +1 = 1. Thus, the modified system is
A system of differential equations alone does not generally define a unique solution, and it is necessary to add to the formulation of the problem a number of addition conditions. These are either “boundary conditions”, if further information is given at two or more values of x, or “initial conditions”, if all components of y are specified at a single value of x. If the value of y(x 0 ) = x 0 is given, then the pair of equations
is known as ‘initial value problem IVP’. Our main interest is with exactly this problem, where the aim is to obtain approximate values of y(x) for specific values of x, usually with x > x 0 , corresponding to the prediction of future states of a differential equation system.
Note that for an N -dimensional system, the individual components of an initial-value vector need to be given specific values. Thus, we might write
For many naturally occurring phenomena, the most appropriate form in which to express a differential equation is as a high order system. For example, an equation might be of the form
of the motion of physical systems subject to forces are the equation system of the form
.. .. . .
(8.14)
8 86 Numerical Approach where the equations, though second order, do have the advantages of being autonomous, (that is
To write (8.13) in what will become our standard first order system form, we can introduce additional components y N +1 , y N +2 , . . . , y 2N . The differential equation (8.14) can now be written as first order system
.. . .. .
.. . .. .
Historically, the first description of the orbits of the Solar System’s planets was empirical, given by the three Kepler’s laws. These empirical laws were demonstrated by Newton, as a simple 2-body application of the theory of universal gravitation. But, in reality, the Solar System is a n-body system and not a 2-body one, since the orbit of each planet is not only determined by the Sun’s attraction. Every planet is also influenced by the other planets. So, considering the Solar System a 2-body system is a very rough estimation of the reality.
This problem can be overridden in two ways:
@ first it is possible to use computers to simulate to a certain extent a n-body problem; • @ a second solution is to consider that the planets’ orbits can be described by Kepler’s laws, • applying to the orbits some corrections due to the fact that planets interact gravitationally.
8.2 A system of differential equations 87 Kepler’s Laws first The planets’ orbits around the Sun are described by three laws first formulated by Johannes Kepler (1571-1630), a great astronomer of the sixteenth century. These empirical laws have been initially formulated on the basis of astronomical observations of Mars’ orbit, made by Tycho Brahe (1546- 1601). Later, Isaac Newton (1642-1727) demonstrated that these three laws can be obtained as an application of the law of universal gravitation in an approximated case (a 2-body problem where the bodies are point-like, the Sun is still and the planets don’t interact ). To formulate this demonstration, Newton invented the infinitesimal calculus. Here are Kepler’s three laws and their main consequences:
@ Kepler’s first law The orbits of the planets are ellipses with the Sun at one focus. • @ Kepler’s second lawThe line joining the planet to the Sun sweeps equal areas in equal • times. • @ Kepler’s third lawThe ratio of the squares of the revolutionary period and the cube of
p 2 p 2
1 2 the semi-major axis is constant for all planets: A 1 3 = A 2 3 .
Now the first law implies the following:
W the orbit of a planet is planar, this means that it all lies on a single plane; • W the orbit is also closed and periodic; • W the distance between Sun and planet is not constant during the orbit.
• A very important consequence of the second law is that the planet’s speed will not be constant during the orbit. In fact, due to the elliptical shape of the orbit (first law) the line joining the Sun and the planet will not be a constant distance. Since, as the second law says, the area swept by this line must be constant, this means that the speed of the motion will vary during the orbit. So, as a consequence of the first and second Kepler’s laws, the planet will go faster near the perihelion(the point where the object is nearest to the sun) and slow down when it comes near the aphelion (the point of this orbit where the object is farthest).
The third law implies that the period of a planet increases rapidly with the radius of its orbit: the more the planet is distant from the sun, the bigger is its major semi-axis, and the longer is its period P. This means that the farther planets of the solar system are much slower than the inner planets, and therefore have much longer years (the time it takes them to make a complete orbit around the Sun).
8 88 Numerical Approach The simplest formulation of the 2-body problem: Newton’s solutions:
A first approximated example of the 2-body problem is the case when the smallest of the two bodies has a negligible mass. In this case, the heaviest body can be considered as having an infinite mass, and so being still, while the lighter body evolves around it. Newton showed that in the hypothesis formulated above, using the second principle of dynamics, and knowing that the two bodies are attracted by each other by the force of universal attraction F ~ = −γ m r 2 M e ~ r , where,
It is possible to determine the orbit of the lighter object which comes near the Sun from an infinite distance. The orbits resulting from this calculations are four curves called conics, represented on the drawing on the right.
The orbits found by Newton as a solution to the 2-body problem are very different from each other: some of them are closed and periodic, meaning that the body will periodically continue following that trajectory (circle and ellipse) while the others are opened orbits (parabola and hyperbola).
But, in which of the possible orbits determined by Newton the body will lie? This depends on the initial binding energy of the two bodies, or in other words the energy required to bring the two bodies at an infinite distance.
In particular, the binding energy E is a measure of how much a body is tied to another : it corresponds to the work necessary to bring the two bodies at an infinite distance. This energy can be negative if the two bodies are gravitationally bound (as in the case of closed orbits, such as the circle or the ellipse) or positive if the two bodies are not bound, and the orbit is opened. The value of the binding energy E is inversely proportional and of opposite sign of the major axes. This means that a very little semi-axis corresponds to a strongly negative energy and therefore to a very bound orbit.
This principle, applied to the solar system, marked the end of the ptolemaic system. In fact, being the force exerted by a body proportional to its mass, it seemed impossible that the Earth could make the Sun orbit around it.
Even if the two assumptions made in this case (the fact that there are only 2 bodies, isolated from the universe, and the assumption that the heaviest of the two bodies is still) seem fairly restrictive, they allow us to give a first explanation of why Kepler’s three laws describe, in a first approximation, the dynamics of the solar system.
In the reality, it is difficult to imagine a body having an infinite mass. If the mass of the lighter body is not negligible (or the two bodies have comparable masses) the problem is a little more complicated. In fact, in this case, it is necessary to define the center of mass of the system that will be the point around which the two bodies will evolve.
8.2 A system of differential equations 89 For the restricted 2-body problem, some important simplifications can be introduced in the equations expressing Newton’s laws. • The larger body may be assumed to remain at rest, and its position may then be chosen
as the origin of the coordinate system. • The motion of the smaller body must lie in the plane determined by its initial velocity
direction and the position of the larger body. So the system is only 2-dimensional. For obvious reasons we exclude the case in which the initial velocity vector points exactly at the larger body. • Finally, we can rescale the space and time coordinates to make the constant γ m equal to
Formulation Solving the famous Kepler’s Problem is an application to the previous approach, this problem describes the motion of a single planet about a heavy sun. By this we mean that, although the sun exerts a gravitational attraction on the planet, we regard the corresponding attraction of the planet on the sun as negligible, and the that the sun will be treated as being stationary. This approximation to the physical system can be interpreted in another way: even though both bodies are in motion about their center of mass, the motion of the planet relative to the sun can be modeled using the simplification we just described. We also make a further assumption, that the motion of the planet is confined to a planet.
Let y 1 (x) and y 2 (x) denote rectangular coordinates centered at the sun, specifying at time x the position of the planet. Also let y 3 (x) and y ( x) denote the components of velocity in y 1 and y 2 directions, respectively. If M denotes the mass of the sun, γ the gravitational constant and m the mass of the planet, then attractive force on the planet will have the magnitude
Resolving this force in the coordinate directions, we find that the components of acceleration of the planet, due to this attraction, are
8 90 Numerical Approach where the negative sign denotes the inward direction of the acceleration. We can write the equations of motion:
(8.16)
By adjusting the scales of the variables, the factor γ M can be removed from the formulation, and we arrived at the the equations
(8.17)
The solutions of this system are known to be conic sections, that is ellipses, parabolas or hyperbolas, if we ignore the possibility that the trajectory is a straight line directed either towards or away from the sun.
Considering the ellipse case only, we find that all points on the trajectory lie on the ellipse
= 1 where e ≥ 0 is the eccentricity of an ellipse on which orbit lies.
This very problem well be our test to evaluate the performance of the implemented algorithms. The eight implemented algorithms are divided evenly between two major methods: W One-Step Method
R The fixed step size algorithms:
8.2 A system of differential equations 91
@ Euler @ Heun @ Runge-Kutta of order tour or RK4 P The variable step size algorithm: : Runge-Kutta-Fehlberg or RKF45 W Multi-step Methods R
The fixed step size algorithms: @ Adams-Bashforth-Moulton @ Milne-Simpson @ Hamming P The variable step size algorithm: :
Adams variable step size The one-step methods are referred to as memoryless. That is, the values of Y i (x) for x before x j do not directly affect the values of Y i (x) for x after x j .
In the multi-step methods information at more than one previous point is used to estimate the solution at the next point. That is, the points at x k , x k−1 , x k−2 , . . . could make the approximation of y(x k+1 ), which could make the process more efficient.
Performance
8 One Step Methods = First: Euler You notice the poor performance of this method to solve an IVP (8.1), taking the scalar form
for explanatory purpose, we recall that basic idea of Euler is to use the tangent line to the actual solution curve as an estimate of the curve itself(see Fig. 8.6), with the though that provided we don’t project too far along the tangent line on a particular step, these two won’t drift too far apart. In reality, this turn out to be expecting pretty much. Even when extremely small step sizes are used, over a large number of steps the error starts to accumulate and the two solutions are part company! The general rule of this type is the Euler update formula:
8 92 Numerical Approach which tells us how to compute y 1 once we know y 0 , and then how to compute y 2 once we know
Where the actual solution curve, as in Fig. 8.6 is concave-up, its tangent line will underestimate the vertical coordinate of the next point.
For each step the eq. (8.18) should include the correction O(h 2 ) which refers to its local Truncation
Besides, there are actually two kinds of Truncation Errors:
Local - occurs at each step!
propagational - due to the error at previous steps!
Their sum makes the global truncation error!
There are many reasons that Euler’s method is not recommended for practical use, among them, J the method is not accurate compared with other methods run at the same step-size!
J neither is it very stable!
8.2 A system of differential equations
8 94
= Second:Heun Heun’s method has two versions, a.k.a modified Euler, or Runge-Kutta second order
which use the midpoint method, i.e. take a “trial” step to the midpoint of the interval. Then use
8.2 A system of differential equations 95 the value of both x and y at that midpoint to compute the “real” step across the whole interval. (See Fig. 8.11)
The other form is:
Both versions have the Truncation Error of O(h 2 ). They have two evaluations per step as you notice, while Euler has but one, accordingly they deliver rather satisfactory results!
8 96
8.2 A system of differential equations
= Third:Runge-Kutta Runge-Kutta methods requires f our evaluation per a step h.
See Fig. 8.15 for the graphical interpretation.
An excellent discussion of the pitfalls in constructing a good Runge-Kutta algorithm is given in (17) and (? ).
8 98 Numerical Approach
It has a truncation error of O(h 4 ) then. Butcher (16) has established the relationship between the number of evaluations per step and the local truncation error shown in Table.
Table 8.1: The relation between evaluations number and local truncation error While for many scientific users, fourth order Runge-Kutta is the α and ω, especially if you combine it with an adaptive step size, but keep in mind Runge-Kutta is just a field-horse and no race-horse, so it’s not the best option though, once speed and high accuracy are required.
The Predictor-corrector methods or Bulirsch-Stoer method proved to be very much more efficient for problems where very high accuracy is a requirement. These methods are Arabians, real high-strung racehorses. So use Runge-Kutta whenever:
d you don’t know any better, or d you have an intransigent problem where Predictor-Corrector is falling, or d you have a trivial problem where computational efficiency is of no concern.
8.2 A system of differential equations
This will be superior to the midpoint methods Heun and Modified Euler, if at least twice as large a step is possible with RK4 for the same accuracy, and more accurate than Euler with one-quarter of the step size h. Is that so? The answer is: often, perhaps even usually, but surely not always. This means that high order does not always imply high accuracy, however for sure it is superior to both Euler and Heun in a sense of practice of science rather than as a statement of strict mathematics.
Adaptive Step Size A good ODE integrator should exert some adaptive control over its own Progress, making frequent changes in its step size. Usually the purpose of this step size control is to achieve some predetermined accuracy Tolerance through estimated truncation error in the solution with minimum computational efforts. Many small steps should tiptoe through treacherous terrain, while a few great strides should speed through smooth uninteresting countryside. The resulting gain in efficiency can sometimes be factors of ten, a hundred, or more, depends on the problem in hand.
= Fourth:Runge-Kutta-Fehlberg RKF45 This technique consists of using a Runge-Kutta method with local truncation error of order five
to estimate the local error in a Runge-Kutta method of order four,
8 100 Numerical Approach where
See (17) chapter 5 for details. An advantage of this method is that only six evaluations of f are required per step. Arbitrary Runge-Kutta methods of order four and five used together require (see Table 8.1) at least four evaluations of f for the fourth-order method an an additional six for the f ifth-order method for a total of at least ten functional evaluations. The value of q is determined at the n th step for two purposes:
¬ To reject, if necessary, the initial choice of h at the n th step and repeat the calculation using 8.2 A system of differential equations 101
To predict an appropriate choice of h for the (n + 1) st step.
Finally it has the local truncation error O(h 6 ). 8 Multi Steps Methods All the previous methods use the information from one previous point to compute the successive point; that is, only the initial point (x 0 , y 0 ) is used to compute (x 1 , y 1 ) and in general, y n is used to compute y n+1 . After several points have been found, it is feasible to use several prior points in the calculation. These methods are not self-starting; four ini-
A desirable feature of a multi-step method is that the local truncation error can be determined and a correction term can be included, which improve the accuracy of the answer at each step. Also, it is possible to determine if the step size small enough to obtain an accurate value for y n+1 , yet large enough so that unnecessary and time-consuming calculations are eliminated, or even introduce an adaptive step size, to increase the speed of the calculations considerably.
Using the combinations of predictor, an explicit method and corrector, an implicit method requires only two function evaluations of f (x, y) per step. Besides boosts the power of the algorithm, to cope with the stiff-equations, because the implicit methods can handle this issue fairly enough. = Fifth: Adams-Bashforth-Moulton The Adams-Bashforth-Moulton method is a multi-step method derived from the fundamental theorem of calculus:
Figure 8.18: Integration over [x n , x n+1 ] in Adams-Bashforth method.
8 102 Numerical Approach The predictor uses the Lagrange polynomial approximation for f (x, y(x)) based on the points
(8.22). This process produces the Adams-Bashforth Explicit Predictor:
The corrector is developed similarly. The value p n+1 just computed can now be used. A second Lagrange polynomial for f (x, y(x)) is constructed, which is based on the points (x n−2 , y n−2 ),
then integrated over [x n , x n+1 ] producing the Adams-Moulton Implicit Corrector
Fig. 8.18 shows the nodes for the Lagrange polynomials that are used in developing formulas (8.24) and (8.23), respectively.
According to Mathews (18) and Burden (17) the local truncation error error for both the predictor and corrector is of the order O(h 5 )
8.2 A system of differential equations
= Sixth: Milne-Simpson Another popular Predictor-Corrector scheme. Its predictor is based on
The predictor uses the Lagrange polynomial approximation for f (x, y(x)) based on the points
the Milne Explicit Predictor
The value p n+1 can be used to derive the corrector, another Lagrange polynomial which based on
structed, and the integration is taken over [x n−1 , x n+1 ], the results is the familiar Simpson Implicit Corrector
8 104
According to Mathews (18) and Burden (17) also the local truncation error error for both the predictor and corrector is of the order O(h 5 )
This method, however in practice has less stability than the former one! = Seventh: Hamming Almost has same performance as Adams-Bashforth-Moulton in this example it has been even 1.34418 times faster, these two horses are running neck-and-neck in most application.
The Hamming Explicit Predictor over the interval [x n−3 , x n+1 ] :
and The Hamming Implicit Corrector over the interval [x n−2 , x n+1 ] :
8.2 A system of differential equations 105 Mathews (18) suggests that the following the following stringent inequalities should be used for the step size h.
= Eighth: Adaptive Step Size Adams-Bashforth-Moulton Due to its good performance and stability, we decided to make it our standard program in all applications, and to integrate an adaptive step size process within this algorithm. The results, either in Kepeler’s Problem or field lines visualization were the best and fastest ever.
8 106
107
Fields of a half-wave dipole can be computed almost as fast as those of the Hertzian dipole (15), whereby the former is of greater relevance. So in the following presentation is made used of this type of dipole only.
Figure 9.1: Field lines and their corresponding planes of symmetry
The first example considers a single radiator, located at the origin, both linear λ/2 dipole and circular small loop antenna are present, the figures 9.1 through 9.5. Fig. 9.1 establishes the parallel 9 108 Results attitude of the magnetic field lines of a loop antenna Fig. 9.1(a) and the electric field lines of a linear antenna Fig. 9.1(b), as well as their magnetic and electric surfaces of symmetry, they both form the familiar kidney-shaped closed loops.
Figure 9.2: Field lines and their corresponding planes of symmetry
Figure 9.3: Different circulation of inner and outer formations.
9.1 An omnidirectional radiation
109
Figure 9.4: Same circulation of inner and outer formations.
Figure 9.5: synoptical view on both field types, illustrating Maxwell’s equations Fig.9.5 is devoted to show the Maxwellian equations pictorially by a synoptic view of both fields.
9 110
Results
Once more than one radiator is involved, the effects of directivity appears clearly in a way that old omnidirectional fashion of radiation is no longer valid. Omnidirectional antennas have a radiation pattern that is donut shaped with the antenna at the center of the donut. This means that with the antenna oriented vertically, the signal coverage is equal in all directions in the horizontal plane. Directional antennas have a radiation pattern that is more focused than omnidirectional antennas. The coverage area is limited to a conical area of various widths depending on the type of directional antenna. Fig.9.6 is of great relevance, once we decide for a scheme to plot the field lines, where they are so dense, and discard them otherwise.
(a) Two-element end-fire (c) five-element end-fire (b) four-element end-fire (d) sixelement broadside Figure 9.6: 3-D Radiation patterns of arrays 9.2 A directional radiation 111 The first example considers two elements, a quarter wavelength apart from each other, and fed as an endfire array. They develop low directivity only, as shown in Fig.9.6(a), the three dimensional radiation pattern. The formation of its electric and magnetic field lines is given in Fig.9.7(a), Fig.9.7(b) respectively, to improve the visibility in Fig.9.7(a), only one half of the surface of symmetry is shown.
(a) The electric field lines of twoelement end-fire array (b) The magnetic field lines of twoelement end-fire array Figure 9.7: The directivity of arrays: two elements In the second example a four-element endfire array is shown, which may be compared with a short Yagi-Uda-Antenna or the active region of a medium sized logarithmic periodic dipole antenna (LPD). The pattern shows higher directivity (see Fig.9.6(b), Fig.9.8(a) and Fig.9.8(b) ). The BC, the current lines and the displacement (field) lines in Fig.9.8 are restricted to the main beam. In the moment of the snapshot the inner members of the vertical families have already detached from the first dipole, whereas detachment for the outer members has not yet taken place.
9 112
(a) The electric field lines of fourelement endfire array Results (b) The magnetic field lines of four-element endfire array Figure 9.8: The directivity of arrays: four elements The Figures 9.6(d) and 9.9 consider the electric field launched by a six- element broadside array. In 9.9 two lines families are heading in each main lobe, and one into each side-lobe.
Figure 9.9: The directivity of arrays: six-element broadside array 3D fashion
The pattern of the four element array in Fig.9.6(b) reveals a null backwards. Therefore a subdivision of the BC, which leads to displacement current lines starting from the negative x-axis, is 9.2 A directional radiation 113 to be avoided. This would represent an almost vanishing portion of the displacement current, Fig.9.8(a).
The next example shall demonstrate that lines can be traced over a given time interval. Therefore a coverage of the full range from the positive to the negative x-axis is necessary. A five-element endfire array with a medium back-lobe in Fig.9.6(c) allows this coverage and lines can be started from the negative x-axis. In Fig.9.10(b) the inner BC is calculated for a certain moment, the outer for an eighth of a period later. The displacement currents in Fig.9.10(a) of both moments are shown (magenta color for the later moment) synoptically, this is an application of the 4D visualization, i.e., the spatial and time dependency.
(a) Displacement lines of a five-element endfire array at two moments with a time difference of 45 degrees (b) Two boundary contours for two different times of a five-element endfire array Figure 9.10: The Visualization of the field lines in 4D fashion
115
The field of an electrically small loop antenna are dependent of the area of loop, but are independent of its shape. In Fig. 11, a model of xy square antenna centered at the origin. The area of the loop is assumed to be (∆`) 2 and it carries a uniform current I 0 .
Figure 11: A small square loop Antenna The loop may be viewed as four segments, which each represents a Hertzian dipole, each parallel pair carries opposite directed current. Due to the tiny dimensions of that loop, relative to a point Appendix 116 P in space, the distance vectors from the centers of the four segments become almost parallel, and can be expressed as
(1)
The individual magnetic vector potential contributions due to the four segments of the current loop are given by
(2)
We need to define another four functions, which can be expressed in a Maclaurin series about ∆`, these functions are:
A On small square loop Antenna 117
(3)
Using the expression (4), each term can be expanded as in Eq. (5), which are very good approximations by considering just the first two terms of Eq. (4):
(4)
(5)
Appendix 118
(6)
The terms of Eq. (6) can be expressed in the rectangular coordinates
(7)
The magnetic fields are determined according to Eq. (8)
where
(8)
(9)
A On small square loop Antenna 119 Therefore,
and (10)
(11)
(12)
Appendix 120
Accordingly adding Eq. (13) to Eq. (14) yields H z
The electric fields are obtained by the potential theory
A On small square loop Antenna 121 ~ = ! 0? Where Let’s check up whether ∇ · A
(17)
The first term of Eq. (17) is
and the second term of Eq. (17) is
~ = 0. Indeed adding Eq. (18) to Eq. (19) yields that ∇ · A Eq. (16) is thus reduced to
(20)
Appendix 122
(21)
and
(22)
and finally
Any antenna used for transmission has the same goal. By generating an alternating current in a conductor, a rapidly changing electrical field is induced in free space. An electrical field cannot exist without a magnetic field component. These fields propagate until reaching another conductor - in which an alternating current similar to the original is induced. The only difference being the much lower energy level.
To make the transfer of energy from transmitter to antenna, and from antenna to space most efficient, the antenna is made resonant - in other words the alternating current flowing within the antenna swings to and fro at a rate that matches the arrival of each peak of current from the transmitter. It’s like pushing a child on a swing - if you push at the same rate that the swing is moving, you get bigger swings for your effort...
A magnetic loop antenna is so called because unlike a conventional wire dipole the small transmitting loop aims to generate the magnetic field rather than the electrical one. The result is the same once the radio wave is produced. Both the magnetic and electrical components propagate as before. The significance of generating a magnetic field lies in the tendency of magnetic fields to be captured by nearby conductors, but to be less affected by nearby ground. This may be because the magnetic field falls away more rapidly with distance than the electric field. This doesn’t mean that small loops prefer being near the ground - but it does seem to suggest that they may be happier than the equivalent ’electrical’ antenna at the same hight. Very low dipoles tend to radiate straight up!
Symmetric loop antennas have a plane of symmetry running along the feed and through the loop. Planar loop antennas lie in a single plane which also contains the conductors of the feed.
Three-dimensional loop antennas have wire which runs in all of the x, y, and z directions (in a rectangular Cartesian system). By definition they are not planar. They may, however, be symmetric about planes which contain the feed.
In the case of vanishingly small loops, the traditional calculation assumes that the current ele~ 0 is the same everywhere around the loop perimeter, as can be seen diagrammatically in ment I 0 d` Fig. 12. In this case, the radiation along any loop diameter arrives from an oppositely directed but parallel element of current after a short time delay, which puts in a phase shift so that the radiation contributions do not entirely cancel.
Appendix 124
Even a single, straight wire conductor has an electrical and magnetic field around it so long as the current flowing in it varies over time. However by coiling this conductor, the magnetic lines of force are concentrated (this is the basis of an electromagnet). Small transmitting loop antennas, are therefor small (single turn) coil antennas! Multiple turn small transmitting loop (STL) antennas have been successfully built and used.
The resonance of a dipole antenna is determined by its size - the smallest resonant length is a half wavelength long (slightly smaller in fact because radio frequency travels more slowly in a conductor than in free space ). It is this ’half-wavelength limitation’ that gives the magnetic loop its sense of appeal, or even dare I say, mystery?
A magnetic loop antenna is comprised of a (usually but not necessarily) round conductor of a total length perhaps one quarter or more commonly for short wave one fifth or one eight of a wavelength. Loops made from conductors longer than a quarter wave-length, probably behave more like delta loop antennas.
Since frequency multiplied by wavelength always equals the speed of light, the higher the oper-
B On small circular loop Antenna 125 ating frequency, the smaller the half wave antenna can be. This isn’t so small. The magnetic loop, or small transmitting loop (STL) antenna appears to offer some performance without needing that mystical half-wavelength measurement. Indeed the whole operation of magnetic loops generates considerable debate. The pivot of this argument is the possibility that a small magnetic loop might be as effective as a half wave wire antenna. It is surprisingly difficult to establish one way or the other as to whether this is in fact the case.
Unlike the approach used before, in deriving Time-varying expression for a linear Dipole antenna, we will use the phasor notation this time, because it is much simpler. The shape of the antenna is no matter, it could be instead, a rectangular, square, triangle, ellipse or any symmetrical figure, to simplify the calculations in analysis and construction.
~ of a circular loop of radius a Fig. 12, at the observation Firstly the Magnetic Vector Potential MPV A point in space P has to be found, where R a. Due to the cylindrical symmetry, one deduces that ~ is independent of the azimuthal angle ϕ 0 . The position vectors of P and the integration point Q, A ~ 0 are defined as the following: which lies in the middle of a current-element I 0 d`
(24) (25)
Therefore we have,
Applying the binomial theorem
(28)
Appendix 126 ~ is The infinitesimal segment of the loop d 0 is defined as d 0 = adϕ 0 , therefore in a thin wire A given by:
The current I ~ e (x 0 , y 0 , z 0 ) can be expressed in its three cartesian components
Assuming a very thin loop, enables us to assume that I x = I x = I 0 , where I 0 is constant, which yields:
using:
expansion:
(31)
B On small circular loop Antenna 127 ~ Using the previous relations in A
Accordingly we have
which can be expressed in the cartesian coordinates x, y, z
Similarly we obtain
(36)
Appendix 128 The magnetic field components can now be obtained using the relation:
where
(37)
(39)
(40)
B On small circular loop Antenna 129
(41)
We evaluate the first term of Eq. (41)
The second term of Eq. (41) yields
Appendix 130 Adding (42) to (43) results in H z
~ = −ω A ~ − ∇ϕ, to derive the electric fields expressions, This time we will use the potential theory E
Let’s examine whether the second term of Eq. (45) yields a zero?
(45)
B On small circular loop Antenna 131 and
The electric fields can be now easily obtained
All we need is to substitute for Eq. (35) and Eq. (36) in their counterparts in Eq. (47) to get the results:
Appendix 132
(49)
(50)
Hertzian dipole and small loops are explanatory tools to explain the behavior of real antennas used in practice. In Fig. ?? we have a case of a real antenna, a center-fed finite length dipole of length `. For convenience a sinusoidal current distribution is assumed, that is:
Appendix 134
(51)
Assuming a very thing wire, implies further simplifications, that is x 0 = y 0 = 0, and from any point on the dipole to a space point P , the distance R can be expressed as:
(52)
From Maxwell’s equations we know that
therefore, in free space we have
Accordingly equation (53) can be expressed as
Using Euler’s Identity
and (51) in (53), the x-component of the magnetic field yields:
(55)
C A finite length dipole 135
be rearranged as:
The term in (57) can considerably simplified if we tailored a suitable differentiation that exactly matches it, therefore consider the the following differentiation:
(58)
Appendix 136 The first term in (58) yields:
(59)
The second term in (58) yields:
and eventually the third term yields:
Therefore their sum yields:
(62)
C A finite length dipole
For further simplifications, of (62) we introduce two more dimensions R 1 and R 2 , which represent each the distance between the space point P and the upper and lower tips of the dipole, as shown in Fig. 14, these are defined as:
Appendix 138
Again finding a suitable differentiation for this term simplifies its integration considerably, regard this term:
The first term in (66) yields:
(67)
The second term in (66) yields:
and eventually the third term yields:
(68)
C A finite length dipole 139
together they make the sum
(70)
We, thus can use (70) in (65) and get:
Appendix 140
(73)
Thus
For integrand in (75), the differentiation ∂z ∂ 0
C A finite length dipole
Thus
We set (64), (71), (74) and (77) in (56) to evaluate H x
Appendix 142
C A finite length dipole 143 Which simplifies to:
Therefore H x can be given as:
(78)
Appendix 144
(79)
and
Deriving expressions for the electric field in free space follows by using Ampere-Maxwell’s curl equation:
and
C A finite length dipole
145
(81)
and
(82)
Finally, we Evaluate the z-component E z
(83)
Appendix 146
(84)
Using (84) and (85) in (83), we obtain E z C A finite length dipole 147
D Mathematical Notes
149 Taylor Series are of fundamental importance in numerical analysis. They are the most basic tool for talking about the approximation of functions. Consider a function f (x) that is smooth – when we say “smooth”, what we mean is that its derivatives exist and are bounded (for the following discussion, we need f to have (n + 1) derivatives). We would like to approximate f (x) near the point x = x 0 , and we can do it as follows:
where
is the nth order Taylor polynomial of f about x 0 , and
is the remainder term or truncation error. The point ξ(x) in the error term lies somewhere between the points x 0 and x. If we look at the infinite sum ( i.e. let n → ∞), then the resulting infinite sum is called the Taylor series of f (x) about x = x 0 . This result is also know as Taylor’s Theorem.
Remember that we assumed that f (x) is smooth (in particular, that its derivatives up to order (n + 1) exist and are finite). That means that all of the derivatives appearing in P n and R n are bounded. Therefore, there are two ways in which we can think of the Taylor polynomial P n (x) as
since as (x − x 0 ) gets small, the error term R n (x) goes to zero (n is considered fixed and all terms depending on n are thus constant). Therefore, the approximation improves when x gets closer and closer to x 0 .
more and more terms in the series. When n is increased, the factorial in the denominator of the error term will eventually dominate the (x − x 0 ) n+1 term (regardless of how big (x − x 0 ) is), and thus drive the error to zero.
In summary, we have two ways of improving the Taylor polynomial approximation to a function: by evaluating it at points closer to the point x 0 ; and by taking more terms in the series.
Appendix 150
Note
Difference formula:
That is, forward-difference is used for f 0 (x 0 ) and backward-difference is used for f 0 (x n ).
Q: How accurate is the difference formula?
By Taylor expansion,
Q: How to get more accurate approximation? Using interpolation polynomial. Given data {(x i , y i )} ni=0 , we construct the interpolation polynomial p n (x). Without loosing gener-
Error analysis:
Define D Mathematical Notes 151
Then it can be proved that
If the truncation error is O(h n ) , we say the method has order of accuracy n.
Consider n = 1 : Lagrange form:
Now based on the above error analysis we may use this formula for any x ∈ [x 0 , x 1 ] and trunca-
Truncation error for p 1 (x) :
formula:
Truncation error for central-difference: O(h 2 )
Appendix 152
Note that for our purpose p will always be a positive integer. The notation is occasionally known as the Landau Notation.
Definition 7. (the rigorous Landau O notation) Let g : R → R . If there exists a constant C and > 0 such that
then we say that
Notice that the equality in Equation (88) is a slight abuse of the notation, as O(x p ) is not a function, but rather a class of functions. The pedantic should write g(x) ∈ O(x p ) but since nobody wants to admit pedantry the equality sign is used universally.
Consider the initial value problem for a system of n ordinary differential equations of first order:
on the interval [a, b] where the function f (t, y) :
is continuous in t and Lipschitz continuous in y, that is,
for some positive constant M . Under these conditions, problem (89) admits one and only one
there are many numerical methods in service to solve (89). Methods may be explicit or implicit, one-step or multi-step, actually just few out of many possible methods, they we won’t bother to mention.
We would first recall few methods to warm the reader up and to introduce some notations.(See Table (1))
D Mathematical Notes 153 Explicit Euler: Implicit Euler: Trapezoidal Method:
a) Predictor: b) Corrector:
Simple, Low Order Methods Error
Modified Euler Method Error
Classical Runge-Kutta Method Error
Table 1: Some available Methods Notes: Things to bear in mind ¬ w i denotes the numerical solution, while either y(t i ) or y i stands for the exact solution.
An ODE solver needs to produce a numerical solution w i of system (89) which converges to
stable at working step size.
® In each method, the weights(the coefficients that multiply function f ) add to 1.
¯ The accuracy of(the exponent on h) equals the number of f evaluations.
° Runge-Kutta methods achieve higher accuracy by retaining the one-step form but sacrificing
linearity. One-step form makes it easy to change the step size but it makes it difficult to estimate the local errors.
± Linear methods combine values of w i+1 , w i , w i−1 , . . . in a linear way to achieve higher ac-
curacy, but sacrifice the one-step format (which is easier to understand and to implement).
Appendix 154 Linearity allows for a simple local error estimate, but it makes it difficult to change the step size.
² If the stability of an explicit method requires an unduly small step size, we say that the
problem is stiff and use an implicit method. We will address the stiffness issue early enough here.
³ inasmuch as the best numerical methods we do have today, the entire numerics, regardless
of all sincere efforts, is still in the sweet days of innocence and juvenility, having precious amateurs who may cope with many classes of solutions, at times, does not imply the claim of making professionals out of them, even if they were the best we had had sofar.
´ There is no method, which can deliver the best results, all times.
• Solution of a DE is function in infinite-dimensional space. • Numerical solutions of a DE is based on finite-dimensional approximation. • Differential equation is replaced with algebraic equation, called Difference equation
whose approximates that of given DE.
Theorem 1. Picard.
L|u − v| for a ≤ x ≤ b, all u, v. Then there exists a unique solution y(x) to the initial value problem
Picard’s Theorem also holds for systems of equations like
.. . with initial conditions
where each f i satisfies a Lipschitz condition in each y j . Note that any higher order initial value problem
converted (in many different ways) to a first order initial value problem.
Stiffness and Numerical Stability • Stiffness of ODE: A stiff equation is a differential equation for which certain numerical
methods for solving the equation are numerically unstable, unless the step size is taken to be extremely small. Or one can say, the equations are still once certain implicit methods, in particular backward differentiation methods, perform much better than explicit ones. It has
D Mathematical Notes 155 proved difficult to formulate a precise definition of stiffness. However, Given system (89), consider the n × n Jacobian matrix
We assume that the n eigenvalue λ 1 , . . . , λ n of the J matrix have negative real parts, <λ j < 0, and ordered as follows:
W Stiffness ratio The stiffness ratio of the system y = f (t, y) is the positive number 0
where the eigenvalues of the Jacobian matrix of the system satisfy the relation
The phenomenon of stiffness appears under various aspects
– A linear constant coefficient system is stiff if all of its eigenvalues have negative real
parts and the stiffness ratio is large. – Stiffness occurs when stability requirements, rather than those of accuracy, constrain
the step length. – Stiffness occurs when some components decay more rapidly than others. We will
comeback to this point soon. – A system is said to be stiff in a given interval I containing t if in I the neighboring
solution curves approach the solution curve at rate which is very large in comparison with the rate at which the solution varies in that interval.
• Numerical stability: A numerical stability is a property of numerical algorithms. It de-
scribes how errors in the input data propagate through the algorithm. In a stable method, the errors due to the approximations get damped out as the computation proceeds. In an unstable method, any errors in processing get magnified as the calculation proceeds. Unstable methods quickly generate garbage and are useless for numerical processing.
– The numerical stability of a method together with the condition number defines how
good a result one can get when using approximate methods to calculate a certain mathematical problem. The condition number associated with a numerical problem is a measure of that quantity’s amenability to digital computation, that is, how well-posed the problem is. A problem with a low condition number is said to be well-conditioned, while a problem with a high condition number is said to be ill-conditioned.
Appendix 156 – Sometimes a single calculation can be achieved in several ways, all of which are alge-
braically identically in terms of ideal real or complex numbers, but in practice yield different results as they have different levels of numerical stability. One of the common tasks of numerical analysis is to try to select algorithms which are robust that is to say, have good numerical stability in a wide range of situations. – Definition: Given an algorithm f (x), with x the input data and ε the error in the input
data, ∗ one says that the algorithm is “numerically stable” for the absolute error if
∗ and numerically stable for the relative error if
∗ One says that an algorithm is “numerically unstable” for the absolute error if
∗ and numerically unstable for the relative error if
– Absolute Error Given some value a and an approximation b of a, the “absolute error”
is
– Absolute Error
where the vertical bars denote the absolute value.
More on the Phenomenon of Stiffness Many time dependent phenomena which can be modeled via systems of differential equations are of such a nature that significantly different time scales occur, i.e., that rapidly and slowly varying processes appear simultaneously. While the long-term behavior of such a system is governed by the solution components connected with the slow time scale, transient phenomena are described by rapidly varying components. Stability of such a system means that these strongly varying components do not grow but die out quickly. Despite the stable behavior of such a dynamical system, the numerical simulation by standard techniques often leads to serious difficulties and numerical instabilities. In the direction field of the corresponding system of differential equations, widely differing directions occur even for points which are close together. Conventional discretization methods based on locally frozen directions get into trouble in such a situation. Even small discretization errors may lead to the evaluation of significantly wrong directions; if such a direction is frozen and used for integrating forward up to the next mesh point, the numerical solution obtained at this point will be dramatically inaccurate. Such an instability can only be suppressed if extremely small step-sizes are used, which is out of question due to efficiency reasons. Algorithms which are robust against such a difficulty have been well known for a long time; originally, the design and understanding of such methods was strongly based on simple models, leading to concepts like A-stability.
D Mathematical Notes 157 Jacobian matrix The Jacobian matrix is the matrix of all first-order partial derivatives of a vectorvalued function. Its importance lies in the fact that it represents the best linear approximation to a differentiable function near a given point. In this sense, the Jacobian is akin to a derivative of a multivariate function.
Suppose F : R n → R m is a function from Euclidean n-space to Euclidean m-space. Such a func-
derivatives of all these functions (if they exist) can be organized in an m × n matrix, the Jacobian matrix of F, as follows:
This matrix is denoted by :
The i th row of this matrix is given by the gradient of the function y i for i = 1, . . . , m.
If p is a point in R n and F is differentiable at p, then its derivative is given by J F (p) (and this is the easiest way to compute said derivative). In this case, the linear map described by J F (p) is the best linear approximation of F near the point p, in the sense that
for x close to p.
Example 3. The Jacobian matrix of the function F : R 3 → R 4 with components:
is:
Appendix 158 Example 4. What does Stiffness, look like? Significant difficulties can occurs the standard numerical techniques are applied to approximate a DE when the exact solution contains terms of the form e λ t , where λ ∈ C, with negative real parts. This very term decays to zero with increasing t. The problem ist particularly acute when the exact solution consists of a steady-state-term that does not grow significantly with t, together with a transient term that decays rapidly to zero, that is; it has a transient region whose behavior is on a different scale from that outside this transient region. A physical example of a stiff system involves chemical reaction rates, where typically the convergence to a final solution can be quite rapid. The transient portion will dominate the calculations and produce meaningless results.
An important characteristic of a stiff system is that the equations are always stable, meaning that they converge to a solution. The following example clarifies this characteristic:
The analytical solution is
For any negative |J(t)|, the solutions y and p will eventually converge, demonstrating that the equation is stable. If |J(t)|, is large, the convergence will be quite rapid. During this transient period, a small step size might be necessary for numerical stability. Using the following stiff system,
the analytical solution to this system is
The solutions u and v have a second term e −1000t that has a negligible effect on the solution for t > 0 but can restrict the step size of a numerical solver, it is who causes stiffness. You must usually reduce this step size to limit the numerical instability, h 1000.
We have now enough knowledge, that make us claim that the Formulas we are working on, are non-stiff ones, they don’t contain any negative real part at all, therefore, our decision to adopt D Mathematical Notes techniques used to solve non-stiff DE, based on IVPs was a good one! The two sets of formulas below, consider the simplest cases we treated at all, the linear Dipole Antennas either a Hertzian or a finite Length, we studied pretty much complicated cases actually, that comprise arrays of dipoles, of every possible orientation, location and configuration.
159
E The Algorithms:
E.1 One-step methods: fixed step-size 161 Euler’s Method in Mathematica code • Having IVP Equations System
• INPUT the initial values t 0 , x 0 , y 0 , y 0 , z 0 , a fixed step-size h and the number of the calculated
• Build an iterative loop as follows
Heun’s Method in Mathematica code • Having IVP Equations System
Appendix 162 • INPUT the initial values t 0 , x 0 , y 0 , y 0 , z 0 , a fixed step-size h and the number of the calculated
• Build an iterative loop as follows
The Runge-Kutta Method RK4 in Mathematica code • Having IVP Equations System
• INPUT the initial values t 0 , x 0 , y 0 , y 0 , z 0 , a fixed step-size h and the number of the calculated
• Build an iterative loop as follows
E The Algorithms: 163
E.2 One-step methods: variable step-size The Runge-Kutta-Fehlberg Method RK45 in Mathematica code • Having IVP Equations System
Appendix 164 • INPUT the initial values t 0 , x 0 , y 0 , y 0 , z 0 , a fixed step-size h and the number of the calculated
• Build an iterative loop as follows
E The Algorithms: 165
Appendix 166
E
The Algorithms: 167
The Adams-Bashforth-Moulton Method ABM in Mathematica code • Having IVP Equations System
• INPUT the initial values t 0 , x 0 , y 0 , y 0 , z 0 , a fixed step-size h and the number of the calculated
• Run a starter sub-algorithm RK4
• Now start the Explicit-Implicit Predictor-Corrector main algorithm:
Appendix 168
Explicit Adams-Bashforth Predictor
Implicit Adams-Moulton Corrector
E The Algorithms: 169 The Milne-Simpson Method MS in Mathematica code • Having IVP Equations System
• INPUT the initial values t 0 , x 0 , y 0 , y 0 , z 0 , a fixed step-size h and the number of the calculated
• Now start the Explicit-Implicit Predictor-Corrector main algorithm:
Explicit Milne’s Predictor
Appendix 170
Implicit Simpson’s Corrector
E The Algorithms: 171 The Hamming’s Method in Mathematica code • Having IVP Equations System
• INPUT the initial values t 0 , x 0 , y 0 , y 0 , z 0 , a fixed step-size h and the number of the calculated
• Now start the Explicit-Implicit Predictor-Corrector main algorithm:
Explicit Hamming’s Predictor Appendix 172
Implicit Hamming’s Corrector
E
The Algorithms: 173
The Adams-Bashforth-Moulton Variable step-size Method ABM V in Mathematica code • Having IVP Equations System
• INPUT the initial values t 0 , x 0 , y 0 , y 0 , z 0 , a fixed step-size h and the number of the calculated
• Now start the Explicit-Implicit Predictor-Corrector main algorithm:
Appendix 174 Explicit Adams-Bashforth Predictor
Implicit Adams-Moulton Corrector
175
[1] Richard Feynmann, The Feynmann Lectures on Physics, Volume 2, Benjamin Cummings 1963,
ISBN: 0-201-02117-X
[2] H. Hertz, Die Kräfte electrischer Schwingungen, behandelt nach der Maxwell’schen Theorie, An-
nalen der Physik und Chemie 1889
[3] Warren L. Stutzman, Gary A. Thiele Antenna theory and design John Wiley & Sons,Inc, New
York, 1981
[4] P. Lorraine, D.P. Corson, F. Lorraine, Electromagnetic Fields and Waves W.H. Freeman and
company, New York
[5] Zuhrt H. “Electromagnetic radiation fields”, Springer 1953”
[6] Arfken, G. Mathematical Methods for Physicists, 3rd ed. Orlando, FL: Academic Press, pp. 198-
200, 1985.
[7] Goldstein, H. “The Euler Angles” and “Euler Angles in Alternate Conventions. ” §4-4 and Appendix
B in Classical Mechanics, 2nd ed. Reading, MA: Addison-Wesley, pp. 143-148 and 606-610, 1980.
[8] Landau, L. D. and Lifschitz, E. M. Mechanics, 3rd ed. Oxford, England: Pergamon Press, 1976.
[9] Margenau, H. and Murphy, G. M. The Mathematics of Physics and Chemistry, 2 vols. Princeton,
NJ: Van Nostrand, 1956-64.
[10] Osgood, W. F. Mechanics. New York: Macmillan, 1937.
[11] Tuma, J. J. Dynamics. New York: Quantum Publishers, 1974.
[12] Varshalovich, D. A.; Moskalev, A. N.; and Khersonskii, V. K. “Description of Rotation in
Terms of the Euler Angles. Ӥ1.4.1 in Quantum Theory of Angular Momentum. Singapore: World Scientific, pp. 21-23, 1988.
[13] MATHWORLD Website http://mathworld.wolfram.com/EulerAngles.html
[14] G. Mönich Correlation between Gain and Wave Propagation in the Aperture Region of long trans-
mitting Yagi-Uda -Antennas Int. URSI Symposium 1980, Munich Aug. 26. to 29.
176 Bibliography [15] C.A. Balanis, Antenna Theory, Analysis and Design John Wiley and Sons, New York
[16] Butcher, J. C., “Numerical Methods for Differential Equations” John Wiley and Sons 2003.
[17] R. L. Burden. J. D. Faires, Numerical Analysis. 5th Ed. PWS Publishing Company Boston 1993.
[18] J. H. Mathews, Numerical Methods Using MATLAB, 3rd Ed, 1999, Prentice Hall. Inc.
[19] L. F. Shampine, Numerical Solution of Ordinary Differential Equations, Chapman and Hall, 1994.
177
Big O notation, 152 Biot-Savart Law, 7
continuity equation, 18
density
| Accuracy, 150 | detaching |
|---|---|
| algorithms, 161 | -points, 6 |
| -Multi-step methods, 167 | -process, 5 |
| Adams-Bashforth-Moulton Method, 167 | dipole |
| Hamming’s Method, 171 | -finite length, 4 |
| Milne-Simpson Method, 169 | -Hertzian , xiv |
| variable Adams-Bashforth-Moulton Method, | -radiation, xiii |
| 173 | duality, 28 |
| -One-step methods, 161 | |
| Euler’s Method , 161 | EM |
| Heun’s Method , 161 | -Radiation, 2 |
EM Euler’s Method , 161 - Radiation, 2 Heun’s Method , 161 - experimental confirmation, 17 Runge-Kutta Method , 162 - spectrum, 17 Runge-Kutta-Fehlberg Method , 163 excitation, 22 antenna
linear dipole , xiii field - lines of force, xiii
detachment, 4 field lines, xiii - retardation effect, 9 - Radiation Patterns, xiv frequency domain, 22 - the image theory, xiv antennas H. Hertz, 29
the image theory, 35 hodograph, 43 arbitrarily oriented
finite length, 70 Jacobian matrix, 157
electric, 70 magnetic, 70 loop - Hertzian, 69 - kidney-shaped , 4
Electric, 69 - small, 4 Magnetic, 69 - Radiators, 57 Maxwell’s
Index 178 - Ampere Equation, 10 - classes and forms, 15 - differential Equations, 18 - Equations, 15 - integral Equations, 19 - integro-differential equations, 17, 22
nonperiodic, 23 Numerical
phasor
notation, 22 Picard, 154 plane
meridian , xiii polarization, 4
Radiation
single wire
charge, 2 - time-harmonic applications, 2 - time-varying current, 2 superposition, 23 symmetry
surfaces, 38
Taylor, 149 The longitudinal sectors’ approach, 29 time domain function, 23
walls