On the visualization of time-varying complex
electromagnetic field lines in 3Dand 4D
fashions
vorgelegt von
Dipl.-Inf.
Ibrahim Abu El-Khair
von der Fakult¨
at IV-Elektrotechnik und Informatik
der Technischen Universit¨
at Berlin
Zur Erlangung des akademischen Grades
Doktor der Ingenieurwissenschaft
-Dr.-Ing.-
genehmigte Dissertation
Promotionsausschuss:
Vorsitzender: Prof. Dr. B. Lutterbeck
1. Berichter: Prof. Dr. G. M¨
onich
2. Berichter: Prof. Dr. G. B¨
arwolff
Tag der wissenschaftlichen Aussprache: 24. November 2006
Berlin 2006
D83
i
Contents
The objectiveness of this thesis xvii
Dedication xvii
Acknowledgment 1
1 Principles of EM Radiation 1
1.1 AsingleWire ......................................... 2
1.2 Principle of field detachment shown by small loop and finite length antennas . . . . 4
1.2.1 Derivation of Time-varying Fields for a Hertzian Dipole . . . . . . . . . . . . 7
2 Maxwell’s Equations 15
2.1 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
ii Contents
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 Magneticandelectricwalls ................................. 35
4.2 ImageTheory ......................................... 38
5 From EM Fields to fieldlines 43
6 Start points for line tracing 49
6.1 Topology of surface and displacement currents . . . . . . . . . . . . . . . . . . . . . . 49
6.1.1 Discretization of the conducting plane . . . . . . . . . . . . . . . . . . . . . . . 50
6.1.2 Linefamilies...................................... 53
6.1.3 Chronological visualization of the field evolution . . . . . . . . . . . . . . . . 53
6.1.4 The new surface definitions, applied to a straight wire . . . . . . . . . . . . . 54
7 On the arbitrarily oriented radiators 57
7.1 A provoked solution: Replacement Rules . . . . . . . . . . . . . . . . . . . . . . . . . 58
Howdoesitreallywork? .............................. 58
7.1.1 Limitations ...................................... 59
7.2 RotationMatrices....................................... 60
Contents iii
7.2.1 The Adopted Rotation Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
7.2.2 What are transpose Matrices really good for? . . . . . . . . . . . . . . . . . . . 63
Acorrection...................................... 64
7.3 ToFormulas .......................................... 65
Anexample...................................... 67
Yetanotherexample ................................. 67
7.3.1 For an arbitrarily oriented Hertzian Dipole . . . . . . . . . . . . . . . . . . . . 69
magneticfields .................................... 69
ElectricFields..................................... 69
7.3.2 For an arbitrarily oriented finite length linear dipole . . . . . . . . . . . . . . . 70
MagneticFields.................................... 70
ElectricFields..................................... 70
8 Numerical Approach 73
Aminimalrequirement ............................... 74
Awell-posedProblem................................ 74
Sourcesoferror.................................... 75
8.0.3 Order of Accuracy, Convergence and Stability of numerical Methods . . . . . 76
8.1 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
8.2 A system of differential equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
8.2.1 Application: Kepler’s Problem; a test problem . . . . . . . . . . . . . . . . . . 86
iv Contents
Kepler’sLawsfirst.................................. 87
Formulation...................................... 89
Performance...................................... 91
9 Results 107
9.1 An omnidirectional radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
9.2 Adirectionalradiation....................................110
9.2.1 A passage to 4Dvisualization............................112
Appendix 115
A OnsmallsquareloopAntenna ...............................115
B On small circular loop Antenna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
B.1 TheMagneticFields .................................128
B.2 TheElectricFields ..................................131
C Afinitelengthdipole.....................................133
D MathematicalNotes......................................149
D.1 Taylor Polynomials and Taylor Series . . . . . . . . . . . . . . . . . . . . . . . 149
D.2 Numerical techniques to approximate derivatives. . . . . . . . . . . . . . . . . 150
What does O(hp), or Big O mean?.........................152
D.3 Summary, revision and more . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
E TheAlgorithms:........................................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 RK4in 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 ABMV in Math-
ematicacode................................173
Bibliography 177
Index 177
vii
List of Figures
1 Hertz’sDrawings....................................... 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
viii List of Figures
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 RotatingCoordinates..................................... 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 Twocoordinatessystems! .................................. 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 CenterofInformation..................................... 81
8.4 TheExpertSystem ...................................... 82
8.5 The Interaction among ‘Divide and Conquer’, Expert System, and Display Unit . . . 83
8.6 PoorPerformanceofEuler.................................. 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
8.10 Euler performance for h=π
6400 , it takes 14,000 points, and 49.032 seconds for the
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
havebeencalculatedandused. ............................... 95
8.12 Heun performance for h=π
800 , it takes 3100 points, and 2.219 seconds for the result 96
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
8.16 Runge-Kutta performance for h=π
800 , it takes 3100 points, and 3.172 seconds for
theresult............................................ 99
8.17 Adaptive step size Runge-Kutta-Fehlberg performance for h=π
800 , and Tol = 10−13
it takes 3100 points, and 2.203 seconds for the result . . . . . . . . . . . . . . . . . . . 100
8.18 Integration over [xn, xn+1]in Adams-Bashforth method. . . . . . . . . . . . . . . . . 101
8.19 Adams-Bashforth-Moulton performance for h=π
800 , it takes 1800 points, and 1.281
secondsfortheresult.....................................103
8.20 Milne-Simpson performance for h=π
800 , it takes 1800 points, and 1.265 seconds for
theresult............................................104
8.21 Hamming performance for h=π
800 , it takes 1800 points, and 0.953 seconds for the
result ..............................................105
8.22 Adaptive Step Size Adams-Bashforth-Moulton performance for h=π
800 , and Tol =
10−13 it takes 1000 points, and 0.609 seconds for the result . . . . . . . . . . . . . . . 106
xList of Figures
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 3Dfashion . . . . . . . . . . . 112
9.10 The Visualization of the field lines in 4Dfashion .....................113
11 AsmallsquareloopAntenna ................................115
12 A small circular loop Antenna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
13 Afinitelengthdipole.....................................133
14 Afinitelengthdipole.....................................137
xi
List of Tables
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
twosurfaces........................................... 33
7.1 Replacementrules....................................... 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 SomeavailableMethods...................................153
xiii
The objectiveness of this thesis
The title of thesis is “On The Visualization of Time-varying complex Electromagnetic Field lines in 3D
and 4DFashions”, 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 accom-
pany 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
xiv List of Tables
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 sin2(θ).
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 2Dor 3Darrows.
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 ex-
perimentalist with real clairvoyant powers and perception.
Figure 1: Hertz’s Drawings
xvii
Dedication
This work is dedicated to the shining stars of mine,
Their Highnesses My beloved Parents!
xix
Acknowledgement
I would like to thank professor G. M¨
onich, my chief supervisor for his unlimited support and
encouragement. Thanks go to professor G. B¨
arwolff and professor B. Lutterbeck, members of my
supervisory board together with professor G. M¨
onich 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¨auning 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
Chapter 1
Principles of EM Radiation
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 electro-
magnetic radiation, in waves that move away from the antenna into the free space at the velocity of
light (≈3·108meter 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
21 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 accompa-
nied 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.
1.1 A single Wire
Considering a linear thin wire for the moment we conclude that; for the existence of electro-
magnetic 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 dis-
turbed 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 z-
oriented λ
2antenna
Figure 1.1: synoptical view on both field types, illustrating Maxwell’s equations
(a) Two Wires connected through an im-
pendence
(b) Bent Wire
Figure 1.2: Single wire Radiation, configurations for different possible cases.
41 Principles of EM Radiation
(a) Discontinuous Wire (b) Truncated Wire
Figure 1.3: Single wire Radiation, configurations for different possible cases.
1.2 Principle of field detachment shown by small loop and finite length
antennas
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 5
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 them-
selves 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 radi-
ation,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 ~
Jand charge density ρeappear 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.
61 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.
(a) The growth of the field lines during one quarter
of a cycle
(b) The decaying of the field lines
(c) The detachment process (d) The growth of the field lines during the last quar-
ter of a cycle
Figure 1.4: A small circular loop antenna.
1.2 Principle of field detachment shown by small loop and finite length antennas 7
(a) The growth of the field lines during one
quarter of a cycle
(b) The decaying of the field lines and detachment
process
(c) detached loops (d) The growth of the field lines during the last
quarter of a cycle
Figure 1.5: A short linear dipole antenna.
1.2.1 Derivation of Time-varying Fields for a Hertzian Dipole
In Fig. 1.6, a time-varying current i(t) = I0sin(ωt)flows in a z-oriented infinitesimal piece of wire
d~
`, in a spatial point P(x, y, z)distanced by r=px2+y2+z2, the magnetic flux Density ~
Bcan be
determined by Biot-Savart Law
~
B=µ0i(t)
4π d~
`×~er
r2!(1.1)
where ~er=x
|~r|~ex+y
|~r|~ey+z
|~r|~ezis the unit vector along the line joining the current element and
the point Pin space, and ris 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.
81 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:
~er
r2=−∇(1
r)
~
B=µ0i(t)
4πd~
`×−∇1
r(1.2)
From the vector identity ∇×(ϕ~
C) = ϕ(∇× ~
C) + ~
C×∇(ϕ)Equation (1.2) can be written as
~
B=µ0i(t)
4πr ∇×d~
`+∇× µ0i(t)d~
`
4πr (1.3)
Which, since d~
`is constant, becomes
~
B=∇× µ0i(t)d~
`
4πr (1.4)
Now i(t)was assumed at the start to be a time-varying current given by i(t, r = 0) = I0sin(ωt)at
the dipole. At distances away from the dipole,
Eq.(1.4) cab be written as:
~
B=∇× µ0I0d`
4πr sin(ωt)~ez(1.5)
1.2 Principle of field detachment shown by small loop and finite length antennas 9
With ~
B=∇× ~
A, we then extract the value of the z-oriented Magnetic Vector Potential to be:
~
Az=I0d`
4πr sin(ωt)(1.6)
However, there is a retardation effect due to the time required for the EM field to propagate to the
point P, where r6= 0, i.e., the time varying part becomes sin(ωt −ωr
c) = sin(ωt −βr), therefore
i(t, r) = I0
sin(ωt −βr), and Eq.(1.5) becomes
~
B=∇× µ0I0d`
4πr sin(ωt −βr)~ez(1.7)
where ~
B=µ0~
H, then we get
~
H=∇× I0d`
4πr sin(ωt −βr)
| {z }
Az
~ez(1.8)
~
H=Hx~ex+Hy~ey+Hz~ez=
~ex~ey~ez
∂
∂x
∂
∂y
∂
∂z
0 0 Az
Therefore we obtain:
Hx=∂
∂y (Az)
Hy=−∂
∂x(Az)
Hz= 0 (1.9)
Hx=I0d`
4π
∂
∂y sin(ωt −βr)
r
=I0d`
4π
∂
∂y −β r cos(ωt −βr)∂r
∂y −sin(ωt −βr)∂r
∂y
r2!
=I0d`
4π
∂
∂y −β r cos(ωt −βr)y
r−sin(ωt −βr)y
r
r2
=−I0d`
4πy
r3sin(ωt −βr) + βy
r2cos(ωt −βr)(1.10)
Similarly we get
10 1 Principles of EM Radiation
Hy=−I0d`
4π
∂
∂x sin(ωt −βr)
r
=−I0d`
4π
∂
∂x −β r cos(ωt −βr)∂r
∂x −sin(ωt −βr)∂r
∂x
r2!
=−I0d`
4π
∂
∂x −β r cos(ωt −βr)x
r−sin(ωt −βr)x
r
r2
=I0d`
4πx
r3sin(ωt −βr) + βx
r2cos(ωt −βr)(1.11)
and
Hz= 0 (1.12)
From Ampere-Maxwell Equation (Maxwell’s equations are described in details in Chapter 2 on page
15)
∇× ~
H=~
J+∂~
D
∂t (1.13)
where in free space ~
J= 0 and ~
D=ε0~
E,then we have
∇× ~
H=∂~
D
∂t =ε0
∂~
E
∂t (1.14)
∂~
E
∂t =1
ε0∇× ~
H(1.15)
~
E=Z∂~
E
∂t dt =~
Ex~ex+~
Ey~ey+~
Ez~ez(1.16)
∂~
D
∂t =
~ex~ey~ez
∂
∂x
∂
∂y
∂
∂z
~
Hx~
Hy0
=~ex(−∂
∂z Hy)
| {z }
∂Dx
∂t
+~ey(∂
∂z Hx)
| {z }
∂Dy
∂t
+~ez(∂
∂xHy−∂
∂y Hx)
| {z }
∂Dz
∂t
(1.17)
1.2 Principle of field detachment shown by small loop and finite length antennas 11
∂Dx
∂t =−∂Hy
∂z
=−I0d`
4π
∂
∂z x
r3sin(ωt −βr) + βx
r2cos(ωt −βr)
=−I0d`
4π x
r3
∂
∂z (sin(ωt −βr)) + sin(ωt −βr)∂
∂z x
r3
+βx
r2
∂
∂z (cos(ωt −βr)) + cos(ωt −βr)∂
∂z x
r2!
=−I0d`
4π
∂
∂z x
r3sin(ωt −βr) + βx
r2cos(ωt −βr)
=−I0d`
4π −βx
r3
z
r(cos(ωt −βr)) −3x z
r5cos(ωt −βr)
+ββx
r2
z
r(sin(ωt −βr)) −2x z
r4cos(ωt −βr)!
=−I0d`
4π 3x z
r5cos(ωt −βr)−β3x z
r4cos(ωt −βr) + β2x z
r3sin(ωt −βr)!
=I0d`
4π 3x z
r5+β3x z
r4cos(ωt −βr)−β2x z
r3sin(ωt −βr)!(1.18)
From Eq. (1.18) we can derive ∂Ex
∂t
∂Ex
∂t =I0d`
4π ε 3x z
r5+β3x z
r4cos(ωt −βr)−β2x z
r3sin(ωt −βr)!(1.19)
and finally Exis obtained by integrating Eq. (1.19) with respect to the time
Ex=Z∂Ex
∂t dt
=I0d`
4π ε ω 3x z
r5+β3x z
r4sin(ωt −βr) + β2x z
r3cos(ωt −βr)!(1.20)
Similarly we can obtain the formula for ∂Dy
∂t =∂
∂z Hx
12 1 Principles of EM Radiation
∂Dy
∂t =∂Hx
∂z
=−I0d`
4π
∂
∂z y
r3sin(ωt −βr) + βy
r2cos(ωt −βr)
=−I0d`
4π y
r3
∂
∂z (sin(ωt −βr)) + sin(ωt −βr)∂
∂z y
r3
+βy
r2
∂
∂z (cos(ωt −βr)) + cos(ωt −βr)∂
∂z y
r2!
=−I0d`
4π
∂
∂z y
r3sin(ωt −βr) + βy
r2cos(ωt −βr)
=−I0d`
4π −βy
r3
z
r(cos(ωt −βr)) −3y z
r5cos(ωt −βr)
+ββy
r2
z
r(sin(ωt −βr)) −2y z
r4cos(ωt −βr)!
=−I0d`
4π 3y z
r5cos(ωt −βr)−β3y z
r4cos(ωt −βr) + β2y z
r3sin(ωt −βr)!
=I0d`
4π 3y z
r5+β3y z
r4cos(ωt −βr)−β2y z
r3sin(ωt −βr)!(1.21)
From Eq. (1.21) we can derive ∂Ey
∂t
∂Ey
∂t =I0d`
4π ε 3y z
r5+β3y z
r4cos(ωt −βr)−β2y z
r3sin(ωt −βr)!(1.22)
and finally Eyis obtained by integrating Eq. (1.22) with respect to the time
Ey=Z∂Ey
∂t dt
=I0d`
4π ε ω 3y z
r5+β3y z
r4sin(ωt −βr) + β2y z
r3cos(ωt −βr)!(1.23)
After obtaining the x-component and y-component of the electric field ~
E, we go ahead to derive
the z-component, where ∂Dz
∂t =∂Hy
∂x −∂Hx
∂y
1.2 Principle of field detachment shown by small loop and finite length antennas 13
We start here with the first term ∂Hy
∂x
∂Hy
∂x =I0d`
4π
∂
∂x x
r3sin(ωt −βr) + βx
r2cos(ωt −βr)
=I0d`
4π −βx
r3
x
rcos(ωt −βr) + r2−3x2
r5sin(ωt −βr)
+β−βx
r2
x
r(sin(ωt −βr)) −r2−2x2
r4cos(ωt −βr)!
=I0d`
4π r2−3x2
r5−β2x2
r4sin(ωt −βr)
+βr2−2x2
r4−x2
r4cos(ωt −βr)!
=I0d`
4π r2−3x2
r5−β2x2
r4sin(ωt −βr)
+βr2−3x2
r4cos(ωt −βr)!(1.24)
The second term −∂Hx
∂y is similarly obtained
−∂Hx
∂y =I0d`
4π
∂
∂y y
r3sin(ωt −βr) + βy
r2cos(ωt −βr)
=I0d`
4π −βy
r3
y
rcos(ωt −βr) + r2−3y2
r5sin(ωt −βr)
+β−βy
r2
y
r(sin(ωt −βr)) −r2−2y2
r4cos(ωt −βr)!
=I0d`
4π r2−3y2
r5−β2y2
r4sin(ωt −βr)
+βr2−2y2
r4−y2
r4cos(ωt −βr)!
=I0d`
4π r2−3y2
r5−β2y2
r4sin(ωt −βr)
+βr2−3y2
r4cos(ωt −βr)!(1.25)
14 1 Principles of EM Radiation
Adding the first term (1.24) to the second term (1.25) yields ∂Dz
∂t
∂Dz
∂t =I0d`
4π r2−3 (x2+y2)
r5−β2(x2+y2)
r4sin(ωt −βr)
+βr2−3 (x2+y2)
r4cos(ωt −βr)!
=I0d`
4π 3z2−r2
r5−β2z2−r2
r3sin(ωt −βr) + β3z2−r2
r4cos(ωt −βr)!(1.26)
From Eq. (1.26) we can get ∂Ez
∂t
∂Ez
∂t =I0d`
4πε 3z2−r2
r5−β2z2−r2
r3sin(ωt −βr) + β3z2−r2
r4cos(ωt −βr)!(1.27)
and finally Ezis obtained by integrating Eq. (1.27) with respect to time
Ez=Z∂Ez
∂dt
=I0d`
4πε ω −3z2−r2
r5−β2z2−r2
r3cos(ωt −βr) + β3z2−r2
r4sin(ωt −βr)!(1.28)
16 2 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 ad-
mitted 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 situa-
tions 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.
2.1 The classes of Maxwell’s equations
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. Sub-
sequent studies resulted in a broader understanding of the nature and origin of radiant energy.
Differential Form Integral Form Remarks
∇· ~
DHs~
D·d~s =Rvρedv Gauss’s law
∇· ~
BHs~
B·d~s = 0 Gauss’s law for magnetism
∇× ~
E=−∂~
B
∂t H`~
E ·d~
`=−∂
∂t Rs~
B·d~s Faraday’s law
∇× ~
H=~
J+∂~
D
∂t H`~
H·d~
`=Rs(~
J+∂~
D
∂t )·d~s Ampere’s circuit law
Table 2.1: Overview of Maxwell’s equations
18 2 Maxwell’s Equations
2.1.1 Maxwell’s Equations in differential form
Maxwell’s equations in differential or point, form are:
∇× ~
H=~
Je+∂~
D
∂t ,(2.1)
∇× ~
E=−~
Jm−∂~
B
∂t ,(2.2)
∇· ~
D=ρe,(2.3)
∇· ~
B=ρm,(2.4)
where
~
E=electric field (in V/m),
~
H=magnetic field (in A/m),
~
D=electric flux density (in C/m2),
~
B=magnetic flux density (in Wb/m2),
ρe=electric charge density (in C/m3),
ρm=magnetic charge density (in Wb/m3),
~
Je=electric current density (in A/m2),
~
Jm=magnetic current density (in V/m2),
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 ρmand ~
Jmequal 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 ρmand ~
Jm: first, symmetry
between electric and magnetic is retained, second, equivalent magnetic sources appear in a va-
riety 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 interchange-
ability of space and time differentiation, we obtain the continuity equation
∇· ~
Je+∂ρe
∂t = 0 (2.5)
which expresses conservation of electric charge. Similarly, from (2.2) and (2.4) we obtain the con-
tinuity equation
2.1 The classes of Maxwell’s equations 19
∇· ~
Jm+∂ρm
∂t = 0 (2.6)
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-
tive relations are: ~
D=0~
Eand ~
B=µ0~
H,where 0and µ0are two constants called the electric
permittivity and the magnetic permeability of free space. In addition to 0and µ0, a material is char-
acterized by its conductivity σ, the term σ~
Eis called the conduction current density, which occurs
in response to the impressed (feeding) current. The current density σ~
Eis a current flowing on a
nearby conductor due to the electric field ~
Ecreated by source ~
J(3).
2.1.2 Maxwell’s equations in integral form
Consider a fixed volume vbounded 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,
20 2 Maxwell’s Equations
yields:
Is
~
D·~n dS =Zv
ρedv =Qe,(2.7)
Is
~
B·~n dS =Zv
ρmdv =Qm,(2.8)
Is
~
Je·~n dS =−dQe
dt ,(2.9)
Is
~
Jm·~n dS =−dQm
dt ,(2.10)
where Qeand Qmare the total electric charge (in C) and magnetic charge (in Wb) inside the
volume v, respectively.
Eq. (2.7) means that the total electric charge contained in vequals the outgoing of ~
Dthrough the
surface S of v. A similar interpretation applies to (2.8), since Qmis zero, the total outgoing flux of
~
Bthrough any fixed closed surface is zero.
Eq. (2.9) means that the rate of decrease of the total electric charge Qeinside vequals the amount of
electric charge which leaves vin 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 Qm
is zero, the outgoing flux of Jmthrough S is also zero.
2.1.3 Maxwell’s equations in transform domain
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
~
E(~r, ω) =
∞
Z
−∞
~
E(~r, t)e−ωtdt, (2.11)
~
E(~r, t) = 1
2π
∞
Z
−∞
~
E(~r, ω)eωtdω, (2.12)
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 ~
E~
H~
D~
Bρeρm~
Je~
Jm
frequency-domain symbol ~
E~
H~
D~
B ρeρm~
Je~
Jm
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:
∇× ~
H=~
Je+ω ~
D, (2.13)
∇× ~
E=−~
Jm−ω ~
B, (2.14)
∇· ~
D=ρe,(2.15)
∇· ~
B=ρm,(2.16)
similarly, the continuity equations (2.5) and (2.6) become:
∇· ~
Je+ω ρe= 0 (2.17)
∇· ~
Jm+ω ρm= 0 (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 =√−1and ωis measured in rad/s.
Example 1. Consider, the sinusoidal (or time-harmonic) electric field
22 2 Maxwell’s Equations
~
E(~r, t) = E0cos(ω0t+ϕ(~r)) (2.19)
with angular velocity ω0and initial phase ϕ(~r). The corresponding field in the frequency domain is obtained
by substituting (2.19) into (2.11) and using the integral representation
δ(ω) = 1
2π
∞
Z
−∞
e−ωtdt (2.20)
for the one-dimensional δ-function; we find that
~
E(~r, ω) = πE0(~r)[eϕ(~r)δ(ω−ω0) + e−ϕ(~r)δ(ω+ω0)] (2.21)
2.1.4 Introductory to the phasor-notation
If all fields quantities vary sinusoidally with time, with angular frequency ω0, the electric field
(2.19) may be written as:
~
E(~r, t) = 1
2[~
E(~r)eω0t+~
E∗(~r)eω0t] = <{~
E(~r)eω0t},(2.22)
where the asterisk means the complex conjugate, and
~
E(~r) = E0eϕ(~r)(2.23)
should be compared with the field ~
Eof (2.21) in the frequency domain. While ~
Eis measured in
V/m,~
Eis measured in V s/m. The electric field ~
Eis obtained from ~
Eby multiplying ~
Etimes the
time-dependence factor eω0tand taking the real part of the product, and from ~
Evia 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 23
use of phasor notation to represent time-dependent variables allows us to convert the integro-
differential 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
x(t) = acos(ω t +ϕ)
with amplitude a > 0, frequency ω > 0, and radian phase angle ϕ. (The case of negative amplitude
is treated by adding πto ϕ).
24 2 Maxwell’s Equations
•Trigonometric:x(t) = acos(ϕ) cos(ω t)−asin(ϕ) sin(ω t)
•Complex exponential:x(t) = a
2e ϕ e ω t +a
2e− ϕ e− ω t
•Phasor real part:x(t) = <{a e(ω t +ϕ)}
Equivalence of these expressions can be verified by using Euler’s formula,
e ϑ = cos(ϑ) + sin(ϑ)
Phasors as a e(ω t +ϕ)can be viewed as a vector at the origin of the complex plane having length
aand, at any time t, angle (ω t +ϕ). The vector rotates counterclockwise with time, since ω > 0,
and the projection on the real axis is described by
<{a e(ω t +ϕ)}=<{acos(ω t +ϕ) + sin(ω t +ϕ)}
=acos(ω t +ϕ)
For graphical representation, projection on the real (horizontal) axis is inconvenient, and therefore
we rotate phasors by π
2radians and project on the vertical axis. This makes use of the mathematical
relationship
<{a e(ω t +ϕ)}=={a e(ω t +ϕ+π
2)}
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.,
replace x(t) = asin(ω t +ϕ)by x(t) = acos(ω t +ϕ−π
2)
•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 Xis a time-independent complex-
valued 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 Zfor 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 equiv-
alents.
z(t)Z
a cos(ω t)⇐⇒ a
a cos(ω t +ϕ)⇐⇒ a eϕ
a cos(ω t +β r +ϕ)⇐⇒ a e(β r +ϕ)
a e−α x cos(ω t +β r +ϕ)⇐⇒ a e−α x e(β r +ϕ)
a sin(ω t)⇐⇒ a e−π
2
a sin(ω t +ϕ)⇐⇒ a e(ϕ−π
2)
d
dt(z1(t)) ⇐⇒ ω Z1
d
dt(a cos(ω t +ϕ)) ⇐⇒ ω a eϕ
Rz1(t)dt ⇐⇒ 1
ω Z1
Rasin(ω t +ϕ)dt ⇐⇒ 1
ω a e(ϕ−π
2)
Table 2.3: Some time-domain functions and their phasor-domain equivalents.
Example 2. In a medium characterized by σ= 0,µ=µ0,0, and
~
E=E0sin(108t−βz)~eyV/m
calculate βand ~
H.
Solution:
This problem can be solved directly in time domain or using phasors. We find βand ~
Hby making ~
Eand ~
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,
∇· ~
E=∂Ey
∂y = 0
From Faraday’s law,
∇× ~
E=−µ∂~
H
∂t =⇒~
H=−1
µZ(∇× ~
E)dt
But
26 2 Maxwell’s Equations
∇× ~
E=
~ex~ey~ez
∂
∂x
∂
∂y
∂
∂z
0Ey0
=−∂Ey
∂z ~ex+∂Ey
∂x ~ez
=E0βcos(108t−βz)~ex+ 0
Hence,
~
H=−E0β
µZcos(108t−βz)dt ~ex
=−E0β
µ108sin(108t−βz)~ex(2.24)
It is readily verified that
∇· ~
H=∂Hx
∂x = 0
showing that Gauss’s law for magnetic fields is satisfied. Lastly, from Ampere’s law
∇× ~
H=σ~
E+∂~
E
∂t =⇒~
E=1
Z(∇× ~
E)dt (2.25)
because σ= 0.
But
∇× ~
H=
~ex~ey~ez
∂
∂x
∂
∂y
∂
∂z
Hx0 0
=−∂Hx
∂z ~ey−∂Hx
∂y ~ez
=E0β2
µ108cos(108t−βz)~ey+ 0
where ~
Hin (2.24) has been substituted: Thus eq. (2.25) becomes
~
E=−E0β2
µ108Zcos(108t−βz)dt ~ey
=−E0β2
µ1016 sin(108t−βz)~ey
Comparing this with the given ~
E, we have
E0β2
µ1016 =E0
or
β=±108√µ =±108pµ0·40=±108(2)
c=±108(2)
3×108
=±2
3
2.1 The classes of Maxwell’s equations 27
From eq. (2.24),
~
H=±E0
60πsin(108t±2z
3)~exA/m
Method 2 (using phasors:)
~
E=={~
Eeωt}=⇒~
E=E0e−βz~ey(2.26)
where ω= 108.
Again
∇· ~
E=∂Ey
∂y = 0
∇× ~
E=−ωµ ~
H=⇒~
H=∇× ~
E
−ωµ
or
~
H=1
−ωµ ∂Ey
∂z ~ex=−E0β
ωµ e−βz~ex(2.27)
Notice that ∇· ~
H= 0 is satisfied.
∇× ~
H=ω ~
E=⇒~
E=∇× ~
H
ω (2.28)
Substituting ~
Hin eq. (2.27) into eq. (2.28) gives
~
E=1
ω
∂Hx
∂z ~ey=E0β2e−βz
ω2µ ~ey
Comparing this with the given ~
Ein eq. (2.26), we have
E0=E0β2
ω2µ
or
β=±ω√µ =±2
3
as obtained before. From eq. (2.27),
~
H=±E0(2/3)e±βz
108(4π×10−7)~ex=±E0
60πe±βz~ex
28 2 Maxwell’s Equations
~
H=={~
Heωt}
=±E0
60πsin(108t±βz)~exA/m
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
~
A=={~
Aeωt}
because the given ~
Eis in sine form and not cosine. We could have used
~
A=<{~
Aeωt}
in which case sine is expressed in terms of cosine and eq. (2.26) would be
~
E=E0cos(108t−βz −π/2)~ey=<{~
Eeωt}
or ~
E=E0e−βz−π
2~ey=− E0e−βz~ey
and we follow the same procedure.
2.1.5 Duality
The Maxwell’s equation (2.1) to (2.4) are invariant under the substitutions:
~
E → ~
H,~
E → −~
H,
~
D → ~
B,~
B → −~
D,
(2.29)
~
Je→~
Jm,~
Jm→ − ~
J,
ρe→ρm, ρm→ −ρe,
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
Chapter 3
Tracing of Fieldlines
3.1 The longitudinal sectors’ approach
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 symmetri-
cal 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.”
30 3 Tracing of Fieldlines
Hertz went in rather long ambiguous calculations and definitions of some quantities as m=m
λ,
n=π
T,Π = E ` sin(m r−n t)
rand in a flash, quite suddenly he eventually presented this Qas
Q=E ` m cos(m r −n t)−sin(m r −n t)
m r sin2ϑ(3.1)
Unfortunately, I never found a reference that describes his approach, or refer even to it, except of
a few lines of M¨
onich. 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
Let’s consider a tangential infinitesimal segment d~
Salong the way of electric field ~
Efor a z-oriented
Hertzian as shown in Fig. 3.1. Their cross-product of both yields a zero, as both d~
Sand ~
Eare
parallel to each other.
d~
S×~
E= 0 (3.2)
3.1 The longitudinal sectors’ approach 31
Expressed the quantities in the spherical coordinates
~
E=Er~er+Eϑ~eϑ+Eϕ~eϕ
d~
S=dr ~er+r dϑ ~eϑ+rsin ϑdϕ ~eϕ
eq. (3.2) is then can be expressed as:
d~
S×~
E= (rdϑ Eϕ−rsin ϑ dϕ Eϑ)~er+ (rsin ϑ dϕ Er−dr Eϕ)~eϑ
+ (dr Eϑ−rdϑ Er)~eϕ= 0 (3.3)
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
dr Eϑ−rdϑ Er= 0 (3.4)
is the interesting one which yields
dr Eϑ=rdϑ Er
or Er
Eϑ
=rdϑ
dr (3.5)
Using Maxwell-Amperes curl formula,
∇× ~
H=∂~
D
∂t (3.6)
which in phasor-notion is
∇× ~
H=ωε0~
E(3.7)
For a z-oriented Dipole ~
H=Hϕ~eϕis valid. Solving Eq. (3.7), yields the following:
Er=1
ωε0
1
rsin ϑ
∂(Hϕsin ϑ)
∂ϑ
Eϑ=−1
ωε0
1
r
∂(r Hϕ)
∂r (3.8)
Eϕ= 0
Substituting these values in Eq. (3.4) and multiplying each side by rsin ϑgives the following:
∂(
const
z }| {
rsin ϑHϕ)
∂ϑ dϑ +∂(
const
z }| {
rsin ϑHϕ)
∂r dr = 0 (3.9)
32 3 Tracing of Fieldlines
which means that rsin ϑHϕ=ρ Hϕis constant, where ρ=rsin ϑwhere
Hϕ=β I d`
4π r sin ϑ1 + 1
βre−βr (3.10)
Let’s now express ρ Hϕin the Time-domain
ρHϕ=<(Hϕeωt)
=βI d`
4πsin2ϑ1
βr cos(ωt −β r)−sin(ωt −βr)(3.11)
Which is the function of the surface of the revolution, we seek.
3.2 A synoptical glance on classical methods of field lines representa-
tion:
3.2.1 A line as the intersection of two surfaces
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 ∂~
D
∂t 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 Az. The magnetic lines are given
by Az= const. due to ~
H=∇× ~
A.
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 Reference depicted surf.1 surf.2 Normal Normal Involved
considered frame field ~v xi=k U =kvector 1 vector 2 metric factor f
Some
parallel cartesian ~
H z =k Az=k ~ez∇Az/|∇Az|1
dipoles
Long wire
arbitrary cylindrical d~
D/dt ϕ =k ρHϕ=k ~eϕ∇(ρHϕ)
|∇ρHϕ|ρ
length
spherical
wave spherical ~
E, ~
H r =kΠr=k ~er∇Πr/|∇Πr|1
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 simultane-
ously:
If the intersection of surface 1with surface 2is a field line, then the vectors ~v must be parallel to
the cross product of their normal components ~n1and ~n2respectively.
~v =∇×(~ei
U
f) = U∇×~ei
1
f
| {z }
0
+∇U×(~ei
1
f)(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:
Φ = Zb
aZ2
1
~v ·d~
F=Zb
aZ2
1∇×(~ei
U
f)·d~
F=I(~ei
U
f)·d~s (3.13)
where |d~s|=f·dx. The closed integration path contains two sections parallel to the coordinate xi
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:
35
Chapter 4
E-wall , H-wall and Image theory
4.1 Magnetic and electric walls
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.
36 4 E-wall , H-wall and Image theory
(a) The electric field lines of a small circular
loop antenna, the arrows show an odd sym-
metry
(b) The electric field lines of z-oriented λ
2
antenna, the arrows show an odd symme-
try.
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 sym-
metry to the H-Wall
(b) The magnetic field lines of z-oriented λ
2an-
tenna, 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 37
(a) The E-Wall for a small
circular loop antenna and
the respective surficial
electric currents
(b) The E-Wall for a z-oriented λ
2antenna
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 an-
tenna and the respective surficial electric
currents
(b) The H-Wall for a z-
oriented λ
2antenna and
the respective surficial
electric currents
Figure 4.4: The H-wall and the respective surficial electric currents.
38 4 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.
4.2 Image Theory
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 cir-
cular loop antenna
(b) The E-Wall for a z-oriented λ
2antenna
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
λ
2antenna
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 situ-
ation, 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 z0above the surface, an image at the
same opposite distance will be produced, the orientation of the currents of both source and image
40 4 E-wall , H-wall and Image theory
is shown in Fig. 9.3 showing the E-Wall.
(a) The E-Wall for a small cir-
cular loop antenna, above is the
source and below is the image
(b) The E-Wall for a z-oriented λ
2
antenna, the source and image
Figure 4.7: E-Wall, source and image for a small circular loop antenna and a linear one.
(a) The H-Wall for a small circular
loop antenna, above is the source
and below is the image
(b) The H-Wall for a z-oriented λ
2an-
tenna, the source and image
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.
(a) The E-Wall for a z-oriented λ
2
antenna, the source and image
(b) The H-Wall for a z-oriented λ
2
antenna, the source and image
Figure 4.9: E-Wall, H-Wall, source and image for a linear radiator.
43
Chapter 5
From EM Fields to fieldlines
What are the EM Fields?
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
~
E=~
E(t,~r),~
D=~
D(t,~r)
~
H=~
H(t,~r),~
B=~
B(t,~r)
~r =~r(x, y, z)(5.1)
where,
~
D=ε0~
E,
~
B=µ0~
H
When tvaries, 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.
The right side of both “Maxwell’s curl Equations”, ∇ × ~
E=−∂~
B
∂t and ∇ × ~
H=∂~
D
∂t , include time
derivative terms, that is, the rate of change of the “displacement current” ˙
~
D=∂~
D
∂t with respect to
time and ˙
~
B=∂~
B
∂t the rate of change of the “magnetic flux density” with respect to time.
For an instance of time ti, every point pialong the trajectories of the numerical integration of ˙
~
D
and ˙
~
Bwith respect to the time, can be express as
pi=p(ti, xi, yi, zi)|ti=const (5.2)
44 5 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.
pi=p(ti, xi, yi, zi)|ti=const =p(xi, yi, zi)|ti=const (5.3)
(a) The electric field lines loops of a z-oriented single
linear Dipole lie in a meridian yz-plane, notice that
the Boundary contour ˙
Dz= 0 lies in xy-plane. Even
symmetry to ‘IMS’, and odd symmetry to ‘IES’, ~n =
{0,0,1}.
(b) The magnetic field lines loops of a z-oriented single
linear Dipole lie in an equatorial xy-plane, notice that
the Boundary contour ˙
By= 0 lies in xz-plane. Even
symmetry to ‘IES’, and odd symmetry to ‘IMS’, ~n =
{0,1,0}.
Figure 5.1: The field lines and their symmetry surfaces
Plotting a series of many field lines, representing different ‘start points’ taken for consecutive subin-
tervals 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 coor-
dinate 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 con-
tinuous 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
up, besides its mostly the standard coordinate system in the simulation programs. Writing ˙
~
Dand
45
˙
~
Bin terms of their components, we come to a situation as shown in Fig. 9.2.
˙
~
D=˙
Dx~ex+˙
Dy~ey+˙
Dz~ez
˙
~
B=˙
Bx~ex+˙
By~ey+˙
Bz~ez(5.4)
~
E=Ex~ex+Ey~ey+Ez~ez
~
H=Hx~ex+Hy~ey+Hz~ez(5.5)
Note The quantities in Eq. (5.4) lead their counterparts in Eq. (5.5) by Phase displacement of π
2due
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’,
i.e., ˙
~
D → ˙
~
Band vice versa. Consider Fig. 4.1(a) and its counterpart Fig. 5.1(b), there an infinitesi-
mal segment d˙
~
Dand d˙
~
Bbetween two adjacent points p1and p2along their corresponding loops
˙
~
Dand ˙
~
Bare respectively labeled. In terms of their components, they are
d˙
~
D=d˙
Dx~ex+d˙
Dyd~ey+d˙
Dz~ez
d˙
~
B=d˙
Bx~ex+d˙
By~ey+d˙
Bzd~ez(5.6)
The components d˙
Dx, d ˙
Dy, d ˙
Dyand d˙
Bx, d ˙
By, d ˙
Byeach represent spatial displacement, they each
can be rewritten as dot product (·)of their gradients respectively, and a displacement vector d~r =
dx ~ex+dy ~ey+dz ~ez
d˙
Dx=∇(˙
Dx)·d~r, d ˙
Dy=∇(˙
Dy)·d~r, d ˙
Dz=∇(˙
Dz)·d~r (5.7)
d˙
Bx=∇(˙
Bx)·d~r, d ˙
By=∇(˙
By)·d~r, d ˙
Bz=∇(˙
Bz)·d~r (5.8)
and
d˙
Dx=∇(˙
Dx)·d~r =∂˙
Dx
∂x dx +∂˙
Dx
∂y dy +∂˙
Dx
∂z dz
= Φ1(˙
Dx(t, x, y, z)) (5.9)
46 5 From EM Fields to fieldlines
d˙
Dy=∇(˙
Dy)·d~r =∂˙
Dy
∂x dx +∂˙
Dy
∂y dy +∂˙
Dy
∂z dz
= Φ2(˙
Dy(t, x, y, z)) (5.10)
d˙
Dz=∇(˙
Dz)·d~r =∂˙
Dz
∂x dx +∂˙
Dz
∂y dy +∂˙
Dz
∂z dz
= Φ3(˙
Dz(t, x, y, z)) (5.11)
Similarly
d˙
Bx=∇(˙
Bx)·d~r =∂˙
Bx
∂x dx +∂˙
Bx
∂y dy +∂˙
Bx
∂z dz
= Ψ1(˙
Bx(t, x, y, z)) (5.12)
d˙
By=∇(˙
By)·d~r =∂˙
By
∂x dx +∂˙
By
∂y dy +∂˙
By
∂z dz
= Ψ2(˙
By(t, x, y, z)) (5.13)
d˙
Bz=∇(˙
Bz)·d~r =∂˙
Bz
∂x dx +∂˙
Bz
∂y dy +∂˙
Bz
∂z dz
= Ψ3(˙
Bz(t, x, y, z)) (5.14)
The Equations (5.9) through (5.14) are functions of functions, each of which is a 4-dimensional func-
tion, however they all have their the dimensions of their fields, ˙
~
Dand ˙
~
Brespectively. These could
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,zand t, we need to turn the equations (5.9) through (5.14) into di-
mensionless ones. The trick is done by dividing them each by their corresponding vector-field
magnitude. Using another script, to show the new unit vectors:
d˙
~
D=d˙
~
D
|d˙
~
D|
=d˙
Dx~ex+d˙
Dy~ey+d˙
Dz~ez(5.15)
d˙
~
B=d˙
~
B
|d˙
~
B|
=d˙
Bx~ex+d˙
By~ey+d˙
Bz~ez(5.16)
47
where their corresponding components are
d˙
Dx=d˙
Dx
|d˙
~
D|
=f1(t, x, y, z),d˙
Dy=d˙
Dy
|d˙
~
D|
=f2(t, x, y, z),
d˙
Dz=d˙
Dz
|d˙
~
D|
=f3(t, x, y, z)(5.17)
Equivalently we have
d˙
Bx=d˙
Bx
|d˙
~
B|
=g1(t, x, y, z),d˙
By=d˙
By
|d˙
~
B|
=g2(t, x, y, z),
d˙
Bz=d˙
Bz
|d˙
~
B|
=g3(t, x, y, z)(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
tare admissible, the algorithm form the values for the points on the integration path in form of
pi=p(x, y, z)|t=const, accumulate then in lists, coded them, then use them in plotting.
49
Chapter 6
Start points for line tracing
Although start points must be found in magnetic and electric walls as well, the following explana-
tion 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.
6.1 Topology of surface and displacement currents
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 Jxand Jyvanishes along a curve for a given
time. The intersections of both curves define points of zero current density, called PZD in the
following.
Jx= 0,Jy= 0
Under normal conditions the magnetic field above metallic surfaces shows an elliptical polariza-
tion. Due to ~
J=~n×~
H, this holds true for the current density ~
Jtoo. 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ωDz, the time dependent formulation is given by:
˙
Dz=<(jωDz) cos(ωt)−=(jωDz) sin(ωt)(6.1)
50 6 Start points for line tracing
So vanishing values for an arbitrary time t1are given by:
tan(ωt1) = <(jωDz)
=(jωDz)−→ −ωt1= arctan=(Ez)
<(Ez)±nπ (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.
|Ez|2==2(Ez) + <2(Ez) arg(Ez) = arctan=(Ez)
<(Ez)(±π)(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 t1can be referred to
the phase front:
arg(Ez) = −ωt1(±π)(6.4)
6.1.1 Discretization of the conducting plane
The tracing of the BC’s can be started from any arbitrary point ‘Start’ in the ‘SOS’ Fig. 6.1, where
˙
Dz= 0 holds true. Integrating the magnetic field along the BC with the coordinate syields 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 ar-
eas launching the same displace-
ment current
(b) Fieldlines, embedded in
surfaces of constant I or con-
stant ζ, 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
I=
s1
Z
Start
~
Hd~s =
s1
Z
Start
~n ×~
Jd~s (6.5)
Now the area of interest can be indicated by Imax.
Imax =
Stop
Z
Start
~
Hd~s =
Stop
Z
Start
~n ×~
Jd~s (6.6)
Segmentation of the BC between ’Start’ and ’Stop’ into N parts with the same amount of crossing
currents ∆I=Imax/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.
52 6 Start points for line tracing
They usually end in a PZD. The current ienclosed between two lines is given by:
i=|~
J|·sin α·w(6.7)
The conduction current Igradually becomes the displacement current, stemming from the plane.
The continuity equation requires that the amount of displacement current, launching from the
dashed area wsin α∆ηin Fig.6.2(a) is equal to the difference between the currents ientering at η1
and leaving at η2:
∂i
∂η ·∆η=−∂
∂η (|~
J|sin αw)∆η=wsin α∆η·˙
Dz(6.8)
Eq.(6.8) shows, that the maximum of ioccurs in the BC, where ˙
Dz= 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 ~
jin
the plane, which is parallel to ~
Jeverywhere, but divergence-less,
~
j=g(x, y)·~
J(6.9)
Now a new fictive current ˜
i, crossing the lines η=const. can be referred to
˜
i=|~
j(η)|·sin α·w= ∆I
which does not depend on η.
Extending the term in brackets in eq.(6.8) by the function gdefined in eq. (6.9) yields after differ-
entiation:
−∂
∂η g|~
J|sin αw1
g·∆η=−1
g
∂
∂η g|~
J|sin αw·∆η+1
g2
∂g
∂η (g|~
J|sin αw)·∆η
The first differentiation on the right side leads to zero, as ˜
i=g|~
J|sin αw is independent of η.
Together with eq. (6.8):
˙
Dzwsin α∆η=1
g2
∂g
∂η (g|~
J|sin αw)·∆η
6.1 Topology of surface and displacement currents 53
Rearranging makes a direct integration on either side possible:
˙
Dz
|~
J|=1
g
∂g
∂η ln g=
η
Z
η=0
˙
Dz
|~
J|dη (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
ζ= (1 −1/g)(6.11)
with the new scaling ζfrom 0 at the BC to 1 in the PZD.
6.1.2 Line families
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 Iand ζof the new skew coordinate
frame for the area of interest.
Once the displacement currents lines have been calculated, the values of Iand ζ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 2nd section for the simple cases. So each
line can be associated to the values of Iand ζof its start point. On the other hand, two different
groups or families can be distinguished: lines, starting from points with constant Ilead to a group
of lines, we call horizontal family. Starting from points with constant ζ, leads to a vertical family.
6.1.3 Chronological visualization of the field evolution
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‘,
54 6 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
different times t1, t2, t3. For t1< t2< t3a line appears, for t1> t2> t3a line disappears.
Consider a BC at two different times, the start time t0and an arbitrary time whereby t>t0. The
results for Imax 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(t0)for both
times t0and t, whereas the number Nhas 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.
˜
ζ=∆I(t)
∆I(t0)·1
g=∆Imax(t)
∆Imax(t0)·1
g0<˜
ζ < ∆I(t)
∆I(t0
Keeping the number Nfor both times is easy to enforce for a closed BC. The maximum of ˜
ζ,
˜
ζmax = ∆I(t)/∆I(t0)determines whether lines appear (˜
ζmax >1) or disappear (˜
ζmax <1)
6.1.4 The new surface definitions, applied to a straight wire
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
I=Z~
Hd~s =
φ0
Z
φ=0
Hφ(˙
Dz= 0) ·ρdφ
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.
I=ρHφ˙
Dz=0
φ0
Z
φ=0
dφ =ρ0Hφ(ρ0)φ0(6.12)
where the index 0 denotes a radius where ˙
Dz= 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 ˙
Dz,|~
J|and the path element dη in the cylindrical coordinates one obtains
˙
Dz=1
ρ
∂
∂ρρHφ|J|=Hφdη =dρ
due to symmetry. So eq. (6.12) becomes:
ln g=
ρ
Z
ρ=ρ0
1
ρHφ
∂
∂ρρHφdρ = lnρHφ
ρ
ρ0
g=ρHφ(ρ)
ρ0Hφ(ρ0)=ρHφ(ρ)·const.
57
Chapter 7
On the arbitrarily oriented radiators
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 hori-
zontally 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 com-
plicated 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 lin-
ear 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 any-
thing 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 com-
plicated 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 every-
thing 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
58 7 On the arbitrarily oriented radiators
formulas belong to a real arbitrary oriented, not its substitute.
7.1 A provoked solution: Replacement Rules
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
coordinates system, {X, Y, and Z},{Y, Z, and X},{Z, X, and Y }. At once, one will notice that
they are the rotations of each other, in a positive sense of the right hand rule. Fig.7.1 on page 58
depicts the situation.
(a) X Y Z (b) Y Z X (c) Z X Y
Figure 7.1: Rotating Coordinates
Now, notice the last symbol of each sequence, and compare this to its own axes in figure 1, amaz-
ing, 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 59
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 Xoriented antenna from, old formulas derived for a
Zoriented one, do the following:
¬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.
7.1.1 Limitations
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 straightfor-
ward approach to our method, and it is besides its simplicity to be understood and programmed
a very effective tool too.
60 7 On the arbitrarily oriented radiators
7.2 Rotation Matrices
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 con-
structed. 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.
7.2.1 The Adopted Rotation Matrix
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. Gold-
stein (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 (Ω, ι, ω). Gold-
stein 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 Rz(ϕ), and the second rotation about the rotated X-
axis(X’) Rx(ϑ), and finally the third and last rotation about the rotated Z-axis(Z’) Rz(ψ). The
definition of Rz(ϕ)is given by:
Rz(ϕ) =
cos(ϕ) sin(ϕ) 0
−sin(ϕ) cos(ϕ) 0
0 0 1
(7.1)
7.2 Rotation Matrices 61
Figure 7.2: A Z-oriented Dipole centered at origin, along with its electric symmetry plane!
Applying this Rotation Matrix Rz(ϕ)to a Z-Oriented linear dipole 7.2, it results in that, the X-axis
and Y to rotate clockwisely about the Z-axis by an angle ϕ, Fig.7.3 depicts its effect:
Figure 7.3: The first Rotation about Z-axis!
The primed symbols indicate that a rotation has taken place.
The next rotation Rx(ϑ)that follows the first one immediately, is about the rotated X-axis (X’)by
an angle ϑand given by:
Rx(ϑ) =
1 0 0
0 cos(ϑ) sin(ϑ)
0−sin(ϑ) cos(ϑ)
(7.2)
Fig.7.4 depicts its effect:
62 7 On the arbitrarily oriented radiators
Figure 7.4: The second Rotation about X’-axis!
The last rotation that follow is about the rotated Z-axis (Z’) by an angle ψand given by Rx(ψ)
Rotation Matrix:
Rz(ψ) =
cos(ψ) sin(ψ) 0
−sin(ψ) cos(ψ) 0
0 0 1
(7.3)
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
orientation, we had in mind, i.e., expressed in terms of Euler triple (ϕ, ϑ, ψ).
Figure 7.5: The third Rotation about Z’-axis!
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.
A=Rz(ψ)·Rz(ϑ)·Rz(ϕ)(7.4)
This can be written in terms of a matrix
A=
a11 a12 a13
a21 a22 a23
a31 a32 a33
(7.5)
where
a11 = cos(ψ) cos(ϕ)−cos(ϑ) sin(ϕ) sin(ψ)
a12 = cos(ψ) sin(ϕ) + cos(ϑ) cos(ϕ) sin(ψ)
a13 = sin(ψ) cos(ϑ)(7.6)
a21 =−sin(ψ) cos(ϕ)−cos(ϑ) sin(ϕ) cos(ψ)
a22 =−sin(ψ) sin(ϕ) + cos(ϑ) cos(ϕ) cos(ψ)
a23 = cos(ψ) sin(ϑ)(7.7)
a31 = sin(ϑ) sin(ϕ)
a32 =−sin(ϑ) cos(ϕ)
a33 = cos(ϑ)(7.8)
Rotation matrices are orthogonal matrices, i.e. unit matrices (which are matrices for which M·
MT =I) with real coefficients. The module of the determinant of unit matrices is one, among the
orthogonal 3×3matrices, 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 Amatrix attributed to this
orientation is A=1 0 0
0 1 0
0 0 1 , to obtain an orientation in the x-Axis, we just set ϕ= 0, ϑ =3π
2and ψ =
3π
2instead, to obtain A=0 1 0
0 0 1
1 0 0 . The y-orientation is easily obtained, likewise, by setting ϕ=
π
2, ϑ =π
2and ψ = 0, and thus A=0 0 1
1 0 0
0 1 0 .
7.2.2 What are transpose Matrices really good for?
The transpose of a matrix Ais a matrix formed from Aby interchanging the rows and columns
such that row iof matrix Abecomes column iof the transposed matrix. The transpose of Ais
64 7 On the arbitrarily oriented radiators
denoted by AT. Each element aij in Abecomes element aji in AT. The following represents AT
the transpose of A:
AT=
a11 a21 a31
a12 a22 a32
a13 a23 a33
(7.9)
Rotating a dipole implies rotating its field along with its orientation. Multiplying the rotation
matrix Aby the vector, containing the original coordinates, say ~
X= x
y
z!, will invoke the
replacement rules, ~
X0=A·~
X= x0
y0
z0!these when applied to a field vector, e.g., the magnetic
vector ~
H=Hx~ex+Hy~ey+Hz~ez, or in vector notation ~
H= Hx
Hy
Hz!will result in a vector field
~
H0=
H0
x0
H0
y0
H0
z0
, these expressions though completely correct in content and symbols, they have al-
together different order: H0
x0, H0
y0, H0
z0refer to their corresponding new coordinates,{X0, Y 0, Z0}
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 ATthe transpose of Acorrects this behavior and produces a new vector
~
H=Hx~ex+Hy~ey+Hz~ez, or as vector ~
H=AT·
H0
x0
H0
y0
H0
z0
=
Hx
Hy
Hz
, so the calligraphic text is
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
7.3 To Formulas
Now, in the XYZ system, a space vector ~r from the origin to the point of observation Pin the space
can be expressed as
~r =x~ex+y~ey+z~ez(7.10)
and its magnitude is
r=|~r|=px2+y2+z2(7.11)
in the X’Y’Z’ system (see Fig. 7.6), however, the very same vector is defined by
~r0=x0~e0
x+y0~e0
y+z0~e0
z(7.12)
Figure 7.6: Two coordinates systems!
Thought they both obviously have the very same magnitude r=r0=px2+y2+z2, we have
e.g.,x6=x0, y 6=y0, z 6=z0,~ex6=~e0
x,~ey6=~e0
y,&~ez6=~e0
z. To distinguished between the different cases
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=
X
Y
Z!, the rotated one is ~
X0= X0
Y0
Z0!, so that every point in space ~
P(x, y, z)can be expressed in
both coordinates systems components equivalently as x0
x0
z0!=A· x
y
z!or x
x
z!=AT· x0
y0
z0!.
Table #7.1 summarizes the entire replacement rules used in this work. We will show a sign of
66 7 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
x=a11x0+a21y0+a31x0x0=a11x+a12y+a13z
x=a12x0+a22y0+a32x0x0=a21x+a22y+a23z
x=a13x0+a23y0+a33x0x0=a31x+a32y+a33z
Table 7.1: Replacement rules
Applying a set of rules, let’s say {x:→x0y:→y0z:→z0}on either an electric field vector
~
E=
Ex
Ey
Ez
or a magnetic field vector ~
H=
Hx
Hy
Hz
, rotates the fields and induces
~
E / · {x:→x0y:→y0z:→z0}=~
E0= E0
x
E0
y
E0
z!
and ~
H / · {x:→x0y:→y0z:→z0}=~
H0= H0
x
H0
y
H0
z!respectively. Table #7.2 shows that every
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
H0
x=a11Hx+a12Hy+a13HzE0
x=a11Ex+a12Ey+a13Ez
H0
x=a21Hx+a22Hy+a23HzE0
x=a21Ex+a22Ey+a23Ez
H0
x=a31Hx+a32Hy+a33HzE0
x=a31Ex+a32Ey+a33Ez
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 vec-
tor of a Z-oriented linear dipole regardless of its sort. It has no magnetic component along its
axis or Hz= 0, so the magnetic field vector reduces to ~
H=
Hx
Hy
0
, we furthermore want to
extract expressions for an X-oriented dipole, and have, therefore, constructed a rotation matrix
A=
0 1 0
0 0 1
1 0 0
and its transpose AT=
0 0 1
1 0 0
0 1 0
. The coordinates of any point P(x, y, z)in space
can expressed in terms of newly generated coordinates by x0
y0
z0!= 010
001
100!· x
y
z!= y
z
x!in
space can expressed in terms of newly generated coordinates by y
z
x!, and so the magnetic field
that yields ~
H0= 010
001
100!· Hx
Hy
0!= Hy
Hz
0!, this is indeed expressions for an X-Oriented dipole
alright, for that the Hxvanishes, but the order of the components is not satisfactory.
Now multiply ~
H0by ATcorrects the order, as we desire it: ~
H=
Hx
Hy
Hz
=
0 0 1
1 0 0
0 1 0
·
Hy
Hz
0
=
0
Hy
Hz
Now its been done!
Yet another example
Now, we want to extract a Y-oriented Dipole out of the wellknown Z-oriented formulas, and
have, therefore, constructed a rotation matrix A=
0 0 1
1 0 0
0 1 0
and its transpose AT=
0 1 0
0 0 1
1 0 0
.
The coordinates of any point P(x, y, z)in space can expressed in terms of newly generated coor-
dinates by x0
y0
z0!= 001
100
010!· x
y
z!= z
x
y!in space can expressed in terms of newly generated
coordinates by z
x
y!, and so the magnetic field that yields ~
H0= 001
100
010!· Hx
Hy
0!= Hz
Hx
0!, this
68 7 On the arbitrarily oriented radiators
is, once more, indeed expressions for an Y-Oriented dipole alright, for that the Hyvanishes, but
the order of the components is not satisfactory.
Now multiply ~
H0by ATcorrects the order, as we desire it: ~
H=
Hx
Hy
Hz
=
0 1 0
0 0 1
1 0 0
·
Hz
Hx
0
=
Hx
0
Hz
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
R1=qx2+y2+ (z−`
2)2R0
1=qx02+y02+ (z0−`
2)2
R2=qx2+y2+ (z+`
2)2R0
2=qx02+y02+ (z0−`
2)2
Table 7.3: two more dimensions for a finite length
7.3 To Formulas 69
7.3.1 For an arbitrarily oriented Hertzian Dipole
magnetic fields
Hx=I0d`
4π −a11y0+a21x0
r03+ β −a11y0+a21x0
r02!e− β r0(7.13)
Hy=I0d`
4π −a12y0+a22x0
r03+ β −a12y0+a22x0
r02!e− β r0(7.14)
Hz=I0d`
4π −a13y0+a23x0
r03+ β −a13y0+a23x0
r02!e− β r0(7.15)
The amazing matter, that a new component Hzhas 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
Ex=I0d`
4π ω ε 3z0(a11x0+a21y0) + a31(3z02−r02)
r05
+ β 3z0(a11x0+a21y0) + a31(3z02−r02)
r04
−β2z0(a11x0+a21y0) + a31(z02−r02)
r03!e− β r0(7.16)
Ey=I0d`
4π ω ε 3z0(a12x0+a22y0) + a32(3z02−r02)
r05
+ β 3z0(a12x0+a22y0) + a32(3z02−r02)
r04
−β2z0(a12x0+a22y0) + a32(z02−r02)
r03!e− β r0(7.17)
70 7 On the arbitrarily oriented radiators
Ez=I0d`
4π ω ε 3z0(a13x0+a23y0) + a33(3z02−r02)
r05
+ β 3z0(a11x0+a21y0) + a33(3z02−r02)
r04
−β2z0(a13x0+a23y0) + a33(z02−r02)
r03!e− β r0(7.18)
7.3.2 For an arbitrarily oriented finite length linear dipole
Magnetic Fields
Hx=I0(−a11y0+a21x0)
4π(x02+y02) e− β R0
1+e− β R0
2−2 cos β `
2e− β r0!(7.19)
Hy=I0(−a12y0+a22x0)
4π(x02+y02) e− β R0
1+e− β R0
2−2 cos β `
2e− β r0!(7.20)
Hz=I0(−a13y0+a23x0)
4π(x02+y02) e− β R0
1+e− β R0
2−2 cos β `
2e− β r0!(7.21)
Electric Fields
Ex=η I0
4π(x02+y02) e− β R0
1+e− β R0
2−2 cos β `
2e− β r0∗
a11x02+a21y02−a31x02+y02!(7.22)
7.3 To Formulas 71
Ey=η I0
4π(x02+y02) e− β R0
1+e− β R0
2−2 cos β `
2e− β r0∗
a12x02+a22y02−a32x02+y02!(7.23)
Ez=η I0
4π(x02+y02) e− β R0
1+e− β R0
2−2 cos β `
2e− β r0∗
a13x02+a23y02−a31x02+y02!(7.24)
73
Chapter 8
Numerical Approach
NUMERICAL solution of ordinary differential equations is the most important technique in con-
tinuous 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 dynam-
ics 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 deriva-
tives. 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
74 8 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:
y0=f(x, y)(8.1)
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 prop-
erties 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:
1The model: its construction usually involves simplifications and omissions;
1The data: there may be errors in measuring or estimating values;
1The numerical method : generally based on some approximation (truncation);
1The representation of numbers: (roundoff) for example, πcannot be represented exactly by
a finite number of digits;
1The 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).
76 8 Numerical Approach
8.0.3 Order of Accuracy, Convergence and Stability of numerical Methods
To quantify the order of accuracy of a numerical approximation the Taylor expansion of y(x)∈Cr
around xnis considered, i.e.
yT(xn+h) = y(xn) + hy0(xn) + h2
2y00(xn) + ···+hr
r!y(r)(ξ)(8.2)
where ξ∈[xn, xn+1]and substitute this into the approximate solution xn−1.
Let us assume that the f∈Cr, i.e. of the same order of the solution.
Definition 1. An approximate method is r-th order accurate if the term in hrin the Taylor expansion of
the unknown is correctly reproduced: it is indicated by O(hr).
The knowledge of the order ris 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
hpromises a large gain in accuracy.
Note: The order rdepends only on the method and not on the differential equation, provided f
satisfies the assumption.
Definition 2. A numerical method is convergent if
lim
h→0yn=y(x)
for all xover the finite interval [x0, xN], 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 En+kis considered as measure of accuracy:
The leading order deviation is called local truncation error at the n-th step
En+k=
k
X
j=0
αjy(xn+j)−hφf(y(xn+k, y(xn+k−1,··· , y(xn), xn;h)(8.3)
and the local unit truncation error is given by
θn+k=En+k
h(8.4)
77
where it is considered the Cauchy problem
(u0=f(y)
u0=y(xn)(8.5)
starting at the exact value at time xn, i.e. it is calculated with single-step of integration at the
position (tn+kh)
Definition 3. The total (global) truncation error is the difference between the solution y(xn+1)and yn+1
(the solution calculated after n+ 1 steps):
en+1 =|y(xn+1)−yn+1|(8.6)
The definition of convergence can be formulated using the total error: a numerical method is convergent if
lim
h→0max
n=0,1,··· ,N |en|= 0 (8.7)
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 En+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
lim
h→0θn +k= 0 (8.8)
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→0and 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
as htends to zero.
Definition 6. Let {δn, n = 0,1,··· , N}and {δ∗
n, n = 0,1,··· , N}be any two perturbations of the
discretized problem (i.e. difference equation generated by the method) and let {yn, n = 0,1,··· , N}and
78 8 Numerical Approach
{y∗
n, n = 0,1,··· , N}be the resulting perturbed solutions. Then if there exist positive constant Sand h0
such that, for all h∈(0, h0]:
|yn−y∗
n| ≤ S whenever |δn−δ∗
n| ≤ , 0≤n≤N(8.9)
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.
8.1 Why to undertake a numerical approach?
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 mathe-
matics 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 graph-
ically 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.
8.1.1 The main reason: Discretization
The reason is quite simple; namely, the analytical solutions to DE’s (once available) deal with con-
tinuous functions; i.e. functions that depend continuously on the independent variables. A com-
puter, however, has only finite storage capacity, and hence there is no way to represent continuous
8.1 Why to undertake a numerical approach? 79
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.
8.1.2 Which strategies are used to plot the field lines?
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 sub-
problems, usually of similar type, next to solve these subproblems, possibly in parallel once de-
sired , 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 ev-
idently, 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 ma-
nipulates 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 representa-
tion of a point in space. These data, thereafter deposited in local storages, where only one provider
80 8 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? 81
Specify the Sort of
the Field Lines
and their inducing
radiators.
Specify the
Grouping of
the Field Lines
and their inducing
radiators.
Identify,
Label,
&
Check up
Unit
Magnetic Electric
Time
Interval?
Instance?
Currents
Phasing
Intensity
Radiators
N
umber
Grouping
Location
Orientation
SOS
&
Arrows
N
umber
Grouping
Location
Orientation
Field Lines
Number
Distribution
Grouping
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’.
82 8 Numerical Approach
Sets of rules
Constructors
&
Modifiers
Data Bank
Graphical
Symbols
Data Bank
Basic
Formulas
Modified formulas Modified graphics
Formulas processing
center
Unit
Vectors
Boundary
Contours
Unit
Vectors
Surficial
Currents
Unit
Vectors
Field
Lines
Starting
Points
Finder
Divide
&
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 ori-
gin, 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 orienta-
tion.
8.1 Why to undertake a numerical approach? 83
Providers
of BC
Starting
Points
Expert System
Divide & Conquer
Scheme
Splitter
&
Selector
BC Determination
Surficial Currents
Determination
Field Lines
Determination
Numerical
Unit
Algorithms
Sets of Points
Labeling &
Packing Unit
Sets of Points
Storage Providers
Sets of Points
Labeling &
Packing Unit
Storage Providers
Display Unit 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 distri-
bution 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
84 8 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 distin-
guished ‘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 con-
structs 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.
8.2 A system of differential equations
The terminology used in (8.1) is misleadingly simple, because ycould be a vector-valued function.
Thus, if we are working in RN, and xis permitted to take on any real value, then the domain and
range of the function fwhich defines a differential equation and the solution to this equation are
given by
f:R×RN→RN,
y:R→RN
Since we might be interested in time values that lie only in some interval [a, b],we sometimes
consider problems in which y: [a, b]→RN, and f: [a, b]×RN→RN.When dealing with
specific problems, it is often convenient to focus, not on the vector-valued function fand 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:
y0
1(x) = f1(x, y1, y2, . . . , yN)
y0
1(x) = f1(x, y1, y2, . . . , yN)
.
.
..
.
. (8.10)
y0
N(x) = fN(x, y1, y2, . . . , yN)
A differential equation for which fis a function of yonly, is said to be ‘autonomous.’ Some equa-
tions 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 ‘au-
tonomous’ form. All we need to do to change the formulation is to introduce an additional com-
ponent yN+1 into the yvector, and ensure that this can always maintain the same value as x, by
associating it with the differential equation y0
N+1 = 1.Thus, the modified system is
y0
1(x) = f1(yN+1, y1, y2, . . . , yN)
y0
1(x) = f1(yN+1, y1, y2, . . . , yN)
.
.
..
.
. (8.11)
y0
N(x) = fN(yN+1, y1, y2, . . . , yN)
y0
N+1(x)=1.
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 yare specified at a single value of x. If the value of y(x0) = x0is
given, then the pair of equations
y0(x) = f(x, y(x)), y(x0) = y0(8.12)
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 > x0,corre-
sponding 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
y0= [η1η2. . . ηN]T.
For many naturally occurring phenomena, the most appropriate form in which to express a dif-
ferential equation is as a high order system. For example, an equation might be of the form
y(n)=φx, y, y0, y00, . . . , y(n−1),(8.13)
with initial values given for y(x0), y0(x0), y00(x0), . . . , y(n−1)(x0).Especially important in modeling
of the motion of physical systems subject to forces are the equation system of the form
y00
1(x) = f1(x, y1, y2, . . . , yN)
y00
1(x) = f1(x, y1, y2, . . . , yN)
.
.
..
.
. (8.14)
y00
N(x) = fN(x, y1, y2, . . . , yN),
86 8 Numerical Approach
where the equations, though second order, do have the advantages of being autonomous, (that is
fis a function not of x, but of yonly) and without y0
1, y0
2, . . . , y0
Noccurring amongst the arguments
of f1, f2, . . . , fN.
To write (8.13) in what will become our standard first order system form, we can introduce addi-
tional components yN+1, yN+2, . . . , y2N.The differential equation (8.14) can now be written as first
order system
y0
1(x) = yN+1
y0
2(x) = yN+2
.
.
..
.
.
y0
N(x) = y2N
y0
N+1(x) = f1(x, y1, y2, . . . , yN),(8.15)
y0
N+2(x) = f1(x, y1, y2, . . . , yN),
.
.
..
.
.
y0
2N(x) = fN(x, y1, y2, . . . , yN),
8.2.1 Application: Kepler’s Problem; a test problem
Description
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
the semi-major axis is constant for all planets: p2
1
A3
1=p2
2
A3
2.
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).
88 8 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 M
r2~er, where,
~er=~r
|r|=x~ex+y ~ey+z ~ez, and r =px2+y2+z2
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 Eis a measure of how much a body is tied to another : it cor-
responds 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 equa-
tions 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
1. For example, to interpret these calculations in terms of the motion, in our solar system,
of a planet at approximately one earth-orbit radius, the unit of length would correspond
to 7.2∗1010 m(approximately 45 million miles), and the unit of time to 1.7∗106seconds,
approximately 19 days.
Formulation
Solving the famous Kepler’s Problem is an application to the previous approach, this problem de-
scribes 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 approxima-
tion 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 y1(x)and y2(x)denote rectangular coordinates centered at the sun, specifying at time xthe
position of the planet. Also let y3(x)and y(x)denote the components of velocity in y1and y2
directions, respectively. If Mdenotes the mass of the sun, γthe gravitational constant and mthe
mass of the planet, then attractive force on the planet will have the magnitude
γ M m
y2
1+y2
2
Resolving this force in the coordinate directions, we find that the components of acceleration of
the planet, due to this attraction, are
−γ M y1
(y2
1+y2
2)3/2
−γ M y2
(y2
1+y2
2)3/2
90 8 Numerical Approach
where the negative sign denotes the inward direction of the acceleration. We can write the equa-
tions of motion:
dy1
dx =f1(y1, y2, y3, y4) = y3,
dy2
dx =f2(y1, y2, y3, y4) = y4,
dy3
dx =f3(y1, y2, y3, y4) = −γ M y1
(y2
1+y2
2)3/2,(8.16)
dy4
dx =f4(y1, y2, y3, y4) = −γ M y2
(y2
1+y2
2)3/2,
By adjusting the scales of the variables, the factor γ M can be removed from the formulation, and
we arrived at the the equations
dy1
dx =f1(y1, y2, y3, y4) = y3,
dy2
dx =f2(y1, y2, y3, y4) = y4,
dy3
dx =f3(y1, y2, y3, y4) = −y1
(y2
1+y2
2)3/2,(8.17)
dy4
dx =f4(y1, y2, y3, y4) = −y2
(y2
1+y2
2)3/2,
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
(y1+e)2+y2
2
1−e2= 1 where e≥0is 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 Yi(x)for xbefore xjdo
not directly affect the values of Yi(x)for xafter xj.
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 xk, xk−1, xk−2, . . . could make the approximation
of y(xk+1), which could make the process more efficient.
Performance
The initial value y(0) = [1,0,0,1]T, defines the solution y(x) = [cos(x),sin(x),−sin(x),cos(x)]T.
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:
yn+1 =yn+h·f(xn, yn), xn+1 ≡xn+h, y(x0) = y0(8.18)
92 8 Numerical Approach
which tells us how to compute y1once we know y0, and then how to compute y2once we know
y1, and so on.
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.
Figure 8.6: Poor Performance of Euler
For each step the eq. (8.18) should include the correction O(h2)which refers to its local Truncation
error τi=|yexact −ycomputed|at each step. The global Error τis one degree smaller, i.e. τ=O(h).
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,
Jthe method is not accurate compared with other methods run at the same step-size!
Jneither is it very stable!
8.2 A system of differential equations 93
(a) An ellipse represents the
orbit of the planet
(b) A 3Dcurve represents the
velocity of the planet
(c) A 3Dcurve represents the ac-
celeration of the planet
Figure 8.7: Euler performance for h=π
800 , it takes 3100 points, and 0.093 seconds for the result
(a) An ellipse represents the
orbit of the planet
(b) A 3Dcurve represents the
velocity of the planet
(c) A 3Dcurve represents the
acceleration of the planet
Figure 8.8: Euler performance at h=π
1600 , it takes 4100 points, and 3.171 seconds for the result
94 8 Numerical Approach
(a) An ellipse represents the
orbit of the planet
(b) A 3Dcurve represents
the velocity of the planet
(c) A 3Dcurve represents
the acceleration of the
planet
Figure 8.9: Euler performance for h=π
3200 , it takes 8100 points, and 15.813 seconds for the result
(a) An ellipse represents the
orbit of the planet
(b) A 3Dcurve represents
the velocity of the planet
(c) A 3Dcurve represents
the acceleration of the
planet
Figure 8.10: Euler performance for h=π
6400 , it takes 14,000 points, and 49.032 seconds for the
result
=
Second:Heun Heun’s method has two versions, a.k.a modified Euler, or Runge-Kutta second
order
k1=hf(xn, yn)
k2=hf(xn+1
2h, yn+1
2k1)
yn+1 =yn+k2+O(h3)
(8.19)
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 xand yat that midpoint to compute the “real” step across the whole interval.
(See Fig. 8.11)
Figure 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.
The other form is:
k1=hf(xn, yn)
k2=hf(xn+2
3h, yn+2
3k1)
yn+1 =yn+1
4k1+3
4k2
(8.20)
Both versions have the Truncation Error of O(h2).They have two evaluations per step as you
notice, while Euler has but one, accordingly they deliver rather satisfactory results!
96 8 Numerical Approach
(a) An ellipse represents the or-
bit of the planet
(b) A 3Dcurve represents the
velocity of the planet
(c) A 3Dcurve represents the
acceleration of the planet
Figure 8.12: Heun performance for h=π
800 , it takes 3100 points, and 2.219 seconds for the result
(a) An ellipse represents the or-
bit of the planet
(b) A 3Dcurve represents the
velocity of the planet
(c) A 3Dcurve represents the ac-
celeration of the planet
Figure 8.13: Heun performance for h=π
1600 , it takes 4100 points, and 4.265 seconds for the result
8.2 A system of differential equations 97
(a) An ellipse represents the or-
bit of the planet
(b) A 3Dcurve represents the
velocity of the planet
(c) A 3Dcurve represents the
acceleration of the planet
Figure 8.14: Heun performance for h=π
3200 , it takes 8100 points, and 15.343 seconds for the result
=
Third:Runge-Kutta Runge-Kutta methods requires four evaluation per a step h.
See Fig. 8.15 for the graphical interpretation.
k1=hf(xn, yn)
k2=hf(xn+1
2h, yn+1
2k1)
k3=hf(xn+1
2h, yn+1
2k2)
k4=hf(xn+h, yn+k1)
yn+1 =yn+1
6(k1+ 2k2+ 2k3+k4) + O(h5)
(8.21)
An excellent discussion of the pitfalls in constructing a good Runge-Kutta algorithm is given in
(17) and (?).
98 8 Numerical Approach
Figure 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.
It has a truncation error of O(h4)then. Butcher (16) has established the relationship between the
number of evaluations per step and the local truncation error shown in Table.
Evaluations per Step 2 3 4 5 ≤n≤7 8 ≤n≤9 10 ≤n
Best Possible
Local Truncation Error O(h2)O(h3)O(h4)O(hn−1)O(hn−2)O(hn−3)
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:
dyou don’t know any better, or
dyou have an intransigent problem where Predictor-Corrector is falling, or
dyou have a trivial problem where computational efficiency is of no concern.
8.2 A system of differential equations 99
(a) An ellipse represents the or-
bit of the planet
(b) A 3Dcurve represents the
velocity of the planet
(c) A 3Dcurve represents the
acceleration of the planet
Figure 8.16: Runge-Kutta performance for h=π
800 , it takes 3100 points, and 3.172 seconds for the
result
This will be superior to the midpoint methods Heun and Modified Euler, if at least twice as large a
step is possible with RK4for 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 solu-
tion 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
ˆyn+1 =yn+16
135k1+6656
12825k3+28561
56430k4−9
5k5+2
55k6,
to estimate the local error in a Runge-Kutta method of order four,
yn+1 =yn+25
216k1+1408
2565k3+2197
4104k4−1
5k5,
100 8 Numerical Approach
where
k1=h f(xn, yn),
k2=h f(xn+1
4h, yn+1
4k1),
k3=h f(xn+3
8h, yn+3
32k1+9
32k2),
k4=h f(xn+12
13h, yn+1932
2197k1−7200
2179k2+7296
2179k3),
k5=h f(xn+h, yn+439
216k1−8k2+3680
513 k3−845
4104k4),
k6=h f(xn+1
2h, yn−8
27k1+ 2k2−3544
2565k3+1859
4104k4−11
240k5),
Error Estimate |ˆyn+1 −yn+1|
Step size qualifier q=Tol h
2|ˆyn+1 −yn+1|1/4
0≤hmin ≤h≤hmax
(a) An ellipse represents the or-
bit of the planet
(b) A 3Dcurve represents the
velocity of the planet
(c) A 3Dcurve represents the ac-
celeration of the planet
Figure 8.17: Adaptive step size Runge-Kutta-Fehlberg performance for h=π
800 , and Tol = 10−13
it takes 3100 points, and 2.203 seconds for the result
See (17) chapter 5for details. An advantage of this method is that only six evaluations of fare
required per step. Arbitrary Runge-Kutta methods of order four and five used together require (see
Table 8.1) at least four evaluations of ffor the fourth-order method an an additional six for the
fifth-order method for a total of at least ten functional evaluations. The value of qis determined
at the nth step for two purposes:
¬To reject, if necessary, the initial choice of hat the nth step and repeat the calculation using
8.2 A system of differential equations 101
h=q h, and
To predict an appropriate choice of hfor the (n+ 1)st step.
Finally it has the local truncation error O(h6).
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 (x0, y0)is used to compute (x1, y1)
and in general, ynis used to compute yn+1.After several points have been found, it is feasi-
ble to use several prior points in the calculation. These methods are not self-starting; four ini-
tial points (x0, y0),(x1, y1),(x2, y2)and (x3, y3)must be given in advance in order to generate the
points {(xn, yn) : n≥4}.
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 yn+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:
yn+1 =yn+Zxn+1
xn
f(x, y(x))dx (8.22)
(a) The four nodes for the Adams-Bashforth
predictor (extrapolation is used)
(b) The four nodes for the Adams-Moulton
corrector (interpolation is used)
Figure 8.18: Integration over [xn, xn+1]in Adams-Bashforth method.
102 8 Numerical Approach
The predictor uses the Lagrange polynomial approximation for f(x, y(x)) based on the points
(xn−3, yn−3),(xn−2, yn−2),(xn−1, yn−1)and (xn, yn). It is integrated over the interval [xn, xn+1]in
(8.22). This process produces the Adams-Bashforth Explicit Predictor:
pn+1 =yn+h
24 (−9fn−3+ 37fn−2−59fn−1+ 55fn)(8.23)
The corrector is developed similarly. The value pn+1 just computed can now be used. A sec-
ond Lagrange polynomial for f(x, y(x)) is constructed, which is based on the points (xn−2, yn−2),
(xn−1, yn−1),(xn, yn),and the new point (xn+1, fn+1)=(xn+1, f(xn+1, pn+1)).This polynomial is
then integrated over [xn, xn+1]producing the Adams-Moulton Implicit Corrector
yn+1 =yn+h
24 (fn−2−5fn−1+ 19fn+ 9fn+1)(8.24)
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(h5)
τn+1(h) = 251
720y(5)(µn)h5µn∈(xn−3, xn+1)(L.T.E for the predictor)
τn+1(h) = −19
720y(5)(µn)h5µn∈(xn−2, xn+1)(L.T.E for the corrector)
8.2 A system of differential equations 103
(a) An ellipse represents the orbit
of the planet
(b) A 3Dcurve represents the ve-
locity of the planet
(c) A 3Dcurve represents the ac-
celeration of the planet
Figure 8.19: Adams-Bashforth-Moulton performance for h=π
800 , it takes 1800 points, and 1.281
seconds for the result
=
Sixth: Milne-Simpson Another popular Predictor-Corrector scheme. Its predictor is based on
integration of f(x, y(x)) over the interval [xn−3, xn+1].
yn+1 =yn−3+Zxn+1
xn−3
f(x, y(x))dx (8.25)
The predictor uses the Lagrange polynomial approximation for f(x, y(x)) based on the points
(xn−3, fn−3),(xn−2, fn−2),(xn−1, fn−1)and (xn, fn)over the interval [xn−3, xn+1]. This produces
the Milne Explicit Predictor
pn+1 =yn−3+4h
3(2fn−2−fn−1+ 2fn)(8.26)
The value pn+1 can be used to derive the corrector, another Lagrange polynomial which based on
the points (xn−1, fn−1),(xn, fn),and the new point (xn+1, fn+1) = (xn+1, f(xn+1, pn+1)) is con-
structed, and the integration is taken over [xn−1, xn+1], the results is the familiar Simpson Implicit
Corrector
yn+1 =yn−1+h
3(2fn−1−4fn+fn+1)(8.27)
104 8 Numerical Approach
(a) An ellipse represents the orbit
of the planet
(b) A 3Dcurve represents the ve-
locity of the planet
(c) A 3Dcurve represents the ac-
celeration of the planet
Figure 8.20: Milne-Simpson performance for h=π
800 , it takes 1800 points, and 1.265 seconds for
the result
According to Mathews (18) and Burden (17) also the local truncation error error for both the pre-
dictor and corrector is of the order O(h5)
τn+1(h) = 28
90y(5)(µn)h5µn∈(xn−3, xn+1)(L.T.E for the predictor)
τn+1(h) = −1
90y(5)(µn)h5µn∈(xn−2, xn+1)(L.T.E for the corrector)
This method, however in practice has less stability than the former one!
=
Seventh: Hamming Almost has same performance as Adams-Bashforth-Moulton in this exam-
ple 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 [xn−3, xn+1] :
pn+1 =yn−3+4h
3(2fn−2−fn−1+ 2fn)(8.28)
and
The Hamming Implicit Corrector over the interval [xn−2, xn+1] :
yn+1 =yn−2+ 9yn
8+3h
8(−fn−1+ 2fn+ 2fn+1)(8.29)
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.
h < 0.75
|fy(x, y)|(Adams-Bashforth-Moulton)
h < 0.45
|fy(x, y)|(Milne-Simpson)
h < 0.69
|fy(x, y)|(Hamming)
(a) An ellipse represents the or-
bit of the planet
(b) A 3Dcurve represents the
velocity of the planet
(c) A 3Dcurve represents the ac-
celeration of the planet
Figure 8.21: Hamming performance for h=π
800 , it takes 1800 points, and 0.953 seconds for the
result
=
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.
106 8 Numerical Approach
(a) An ellipse represents the or-
bit of the planet
(b) A 3Dcurve represents the veloc-
ity of the planet
(c) A 3Dcurve represents the accel-
eration of the planet
Figure 8.22: Adaptive Step Size Adams-Bashforth-Moulton performance for h=π
800 , and Tol =
10−13 it takes 1000 points, and 0.609 seconds for the result
107
Chapter 9
Results
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.
(a) The magnetic field lines of a small
circular loop antenna, and its mag-
netic plane of symmetry
(b) The electric field lines of z-
oriented λ
2antenna, and its electric
plane of symmetry
Figure 9.1: Field lines and their corresponding planes of symmetry
9.1 An omnidirectional radiation
The first example considers a single radiator, located at the origin, both linear λ/2dipole and
circular small loop antenna are present, the figures 9.1 through 9.5. Fig. 9.1 establishes the parallel
108 9 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.
(a) The electric field lines of a small
circular loop antenna, and its electric
plane of symmetry
(b) The magnetic field lines of z-
oriented λ
2antenna, and its magnetic
plane of symmetry
Figure 9.2: Field lines and their corresponding planes of symmetry
(a) The magnetic field lines of a
small circular loop antenna
(b) The electric field lines of z-
oriented λ
2antenna, and its mag-
netic plane of symmetry
Figure 9.3: Different circulation of inner and outer formations.
9.1 An omnidirectional radiation 109
(a) The electric field lines of a
small circular loop antenna
(b) The magnetic field lines of z-
oriented λ
2antenna, and its mag-
netic plane of symmetry
Figure 9.4: Same circulation of inner and outer formations.
(a) The electromagnetic field lines of a
small circular loop antenna
(b) The electromagnetic field lines of
z-oriented λ
2antenna
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.
110 9 Results
9.2 A directional radiation
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 (b) four-element end-fire
(c) five-element end-fire (d) six-
element
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 two-
element end-fire array
(b) The magnetic field lines of two-
element 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.
112 9 Results
(a) The electric field lines of four-
element endfire array
(b) The magnetic field lines
of four-element endfire ar-
ray
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 3Dfashion
9.2.1 A passage to 4Dvisualization
The pattern of the four element array in Fig.9.6(b) reveals a null backwards. Therefore a subdi-
vision 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 4Dvisualization,
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 4Dfashion
115
Appendix
A On small square loop Antenna
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 (∆`)2and it carries a uniform current I0.
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
116 Appendix
Pin space, the distance vectors from the centers of the four segments become almost parallel, and
can be expressed as
R1≃r+∆`
2sin ϑsin ϕ
R2≃r−∆`
2sin ϑcos ϕ
(1)
R3≃r−∆`
2sin ϑsin ϕ
R4≃r+∆`
2sin ϑcos ϕ
The individual magnetic vector potential contributions due to the four segments of the current
loop are given by
~
A1=µ0I0∆`
4π
e−β(r+∆`
2sin ϑsin ϕ)
r+∆`
2sin ϑsin ϕ~ex
~
A2=µ0I0∆`
4π
e−β(r−∆`
2sin ϑcos ϕ)
r−∆`
2sin ϑcos ϕ~ey
(2)
~
A3=µ0I0∆`
4π
e−β(r−∆`
2sin ϑsin ϕ)
r−∆`
2sin ϑsin ϕ(−~ex)
~
A4=µ0I0∆`
4π
e−β(r+∆`
2sin ϑcos ϕ)
r+∆`
2sin ϑcos ϕ(−~ey)
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
f1=e−β(r+∆`
2sin ϑsin ϕ)
r+∆`
2sin ϑsin ϕ
f2=e−β(r−∆`
2sin ϑcos ϕ)
r−∆`
2sin ϑcos ϕ
f3=e−β(r−∆`
2sin ϑsin ϕ)
r−∆`
2sin ϑsin ϕ(3)
f4=e−β(r+∆`
2sin ϑcos ϕ)
r+∆`
2sin ϑcos ϕ
Using the expression (4), each term can be expanded as in Eq. (5), which are very good approxi-
mations by considering just the first two terms of Eq. (4):
f=f(0) + f0(0)∆`+1
2f00(0)(∆`)2+. . .
(4)
f0(0) = ∂f
∂∆`∆`=0
, f00(0) ∂2f
∂∆`2∆`=0
, . . .
f1=1
r−∆`
2β
r+1
r2sin ϑsin ϕe−βr
f2=1
r+∆`
2β
r+1
r2sin ϑcos ϕe−βr
f3=1
r+∆`
2β
r+1
r2sin ϑsin ϕe−βr (5)
f4=1
r−∆`
2β
r+1
r2sin ϑcos ϕe−βr
118 Appendix
We can now simply combine the x-directed and y-directed terms, which yields:
Ax=|~
A1+~
A3|
=µ0I0∆`
4π
1
r−∆`
2β
r+1
r2sin ϑsin ϕ−
1
r−∆`
2β
r+1
r2sin ϑsin ϕ!e−βr
=−µ0I0(∆`)2
4π β
r+1
r2sin ϑsin ϕ!e−βr
(6)
Ay=|~
A2+~
A4|
=µ0I0∆`
4π
1
r+∆`
2β
r+1
r2sin ϑcos ϕ−
1
r+∆`
2β
r+1
r2sin ϑcos ϕ!e−βr
=µ0I0(∆`)2
4π β
r+1
r2sin ϑcos ϕ!e−βr
The terms of Eq. (6) can be expressed in the rectangular coordinates
Ax=−µ0I0(∆`)2
4πe−βr β
r+1
r2
px2+y2
r
y
px2+y2
=−µ0I0(∆`)2
4πy e−βr β
r2+1
r3(7)
Ay=µ0I0(∆`)2
4πe−βr β
r+1
r2
px2+y2
r
x
px2+y2
=µ0I0(∆`)2
4πx e−βr β
r2+1
r3
Az= 0.
The magnetic fields are determined according to Eq. (8)
~
H=1
µ0∇× ~
A=Hx~ex+Hy~ey+Hz~ez(8)
where
∇× ~
A=
70
∂Az
∂y −∂Ay
∂z
~ex+
∂Ax
∂z −0
∂Az
∂x
~ey+∂Ay
∂z −∂Ax
∂y ~ez(9)
A On small square loop Antenna 119
Therefore,
Hx=−1
µ0
∂Ay
∂z
=−I0(∆`)2
4πx∂
∂z β
r2+1
r3e−βr!
=−I0(∆`)2
4πx β
r2+1
r3∂e−βr
∂z +e−βr ∂
∂z β
r2+1
r3!
=−I0(∆`)2
4πx β
r2+1
r3−β z
re−βr+e−βr β 2z
r4+3z
r5!
=−I0(∆`)2
4πx β2z
r3−β z
r4−β 2z
r4−3z
r5!e−βr
=−I0(∆`)2
4πx β2z
r3−β 3z
r4−3z
r5!e−βr
=I0(∆`)2
4π 3z x
r5+β 3z x
r4−β2z x
r3!e−βr (10)
Hy=−1
µ0
∂Ax
∂z
=I0(∆`)2
4πy∂
∂z β
r2+1
r3e−βr!
=I0(∆`)2
4πy β
r2+1
r3∂e−βr
∂z +e−βr ∂
∂z β
r2+1
r3!
=I0(∆`)2
4πy β
r2+1
r3−β z
re−βr+e−βr β 2z
r4+3z
r5!
=I0(∆`)2
4πy β2z
r3−β z
r4−β 2z
r4−3z
r5!e−βr
=I0(∆`)2
4πy β2z
r3−β 3z
r4−3z
r5!e−βr
=−I0(∆`)2
4π 3z y
r5+β 3z y
r4−β2z y
r3!e−βr (11)
and
Hz=1
µ0∂Ay
∂x −∂Ax
∂y (12)
120 Appendix
the first term 1
µ0
∂Ay
∂x yields
1
µ0
∂Ay
∂x =I0(∆`)2
4π
∂
∂x βx
r2+x
r3e−βr!
=I0(∆`)2
4π βx
r2+x
r3∂e−βr
∂x +e−βr ∂
∂x βx
r2+x
r3!
=I0(∆`)2
4π βx
r2+x
r3−β x
re−βr+e−βr β r2−2x2
r4+r2−3x2
r5!
=I0(∆`)2
4π r2−3x2
r5+β r2−3x2
r4+β2x2
r3!e−βr (13)
the second term −1
µ0
∂Ax
∂y yields
−1
µ0
∂Ax
∂y =I0(∆`)2
4π
∂
∂y βy
r2+y
r3e−βr!
=I0(∆`)2
4π βy
r2+y
r3∂e−βr
∂y +e−βr ∂
∂y βy
r2+y
r3!
=I0(∆`)2
4π βy
r2+y
r3−β y
re−βr+e−βr β r2−2y2
r4+r2−3y2
r5!
=I0(∆`)2
4π r2−3y2
r5+β r2−3y2
r4+β2y2
r3!e−βr (14)
Accordingly adding Eq. (13) to Eq. (14) yields Hz
Hz=1
µ0∂Ay
∂z −∂Ax
∂y
=I0(∆`)2
4π r2−3 (x2+y2)
r5+β r2−3 (x2+y2)
r4+β2(x2+y2)
r3!e−βr
=I0(∆`)2
4π 3z2−r2
r5+β 3z2−r2
r4+β2z2−r2
r3!e−βr (15)
The electric fields are obtained by the potential theory
~
E=−ω ~
A+∇(∇· ~
A)
ωµ0ε(16)
A On small square loop Antenna 121
Let’s check up whether ∇· ~
A!
= 0? Where
∇· ~
A=∂Ax
∂x ~ex+∂Ay
∂y ~ey+>0
∂Az
∂z ~ez(17)
The first term of Eq. (17) is
∂Ax
∂x =−µ0I0(∆`)2
4πy∂
∂x β
r2+1
r3e−βr!
=−µ0I0(∆`)2
4πy β
r2+1
r3∂e−βr
∂x +e−βr ∂
∂x β
r2+1
r3!
=−µ0I0(∆`)2
4πy β
r2+1
r3−β x
re−βr+e−βr β −2x
r4+−3x
r5!
=µ0I0(∆`)2
4π −3x y
r5−β 3x y
r4+β2x y
r3!e−βr (18)
and the second term of Eq. (17) is
∂Ay
∂y =µ0I0(∆`)2
4πx∂
∂y β
r2+1
r3e−βr!
=µ0I0(∆`)2
4πx β
r2+1
r3∂e−βr
∂y +e−βr ∂
∂y β
r2+1
r3!
=µ0I0(∆`)2
4πx β
r2+1
r3−β y
re−βr+e−βr β −2y
r4+−3y
r5!
=µ0I0(∆`)2
4π 3x y
r5+β 3x y
r4−β2x y
r3!e−βr (19)
Indeed adding Eq. (18) to Eq. (19) yields that ∇· ~
A= 0.
Eq. (16) is thus reduced to
~
E=−ω ~
A=−ωAx~ex−ωAy~ey
=Ex~ex+Ey~ey(20)
=⇒
122 Appendix
Ex=− ωAx
= ωµ0
I0(∆`)2
4πy β
r2+1
r3!e−βr
=β
√µ0ε0
µ0
I0(∆`)2
4πy β
r2+1
r3!e−βr
= βrµ0
ε0
I0(∆`)2
4πy β
r2+1
r3!e−βr
= β η I0(∆`)2
4πy β
r2+1
r3!e−βr (21)
and
Ey=− ωAy
=− ωµ0
I0(∆`)2
4πx β
r2+1
r3!e−βr
=−β
√µ0ε0
µ0
I0(∆`)2
4πx β
r2+1
r3!e−βr
=− βrµ0
ε0
I0(∆`)2
4πx β
r2+1
r3!e−βr
=− β η I0(∆`)2
4πx β
r2+1
r3!e−βr (22)
and finally
Ez= 0.(23)
B On small circular loop Antenna
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 conduc-
tor - 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 trans-
mitting 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 zdirections (in a rect-
angular 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-
ment I0~
d`0is the same everywhere around the loop perimeter, as can be seen diagrammatically in
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.
124 Appendix
Figure 12: A small circular loop Antenna
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.
Firstly the Magnetic Vector Potential MPV ~
Aof a circular loop of radius aFig. 12, at the observation
point in space Phas to be found, where Ra. Due to the cylindrical symmetry, one deduces that
~
Ais independent of the azimuthal angle ϕ0. The position vectors of Pand the integration point Q,
which lies in the middle of a current-element I0~
d`0are defined as the following:
~r =~
r0+~
R(24)
~
r0=acos ϕ0~ex+asin ϕ0~ey(25)
~r =x~ex+z ~eztaken at ϕ0= 0 for simplicity (26)
Therefore we have,
~
R=~
r0−~
r0
=px2+z2−2a x cos ϕ0+a2
=pr2−2a r sin ϑcos ϕ0+a2
=rr1−2a
rsin ϑcos ϕ0+ (a
r)2
The term a
r2can be neglected due to the assumption ra, thus we have
R≃rr1−2a
rsin ϑcos ϕ0(27)
Applying the binomial theorem
(a+b)n=n an−1b+n(n−1)
2! an−2b2+n(n−1)(n−2)
3! an−3b3+. . .
on Eq. (27), further simplification can be obtained:
R≃rr1−a
rsin ϑcos ϕ0(28)
126 Appendix
The infinitesimal segment of the loop d`0is defined as d`0=adϕ0, therefore in a thin wire ~
Ais
given by:
~
A(x, y, z) = µ0
4πZe
~
Ic(x0, y0, z0)e−βR
Rd`0(29)
The current ~
Ie(x0, y0, z0)can be expressed in its three cartesian components
~
Ie(x0, y0, z0) = Ix(x0, y0, z0)~ex+Iy(x0, y0, z0)~ey+Iz(x0, y0, z0)~ez
The current ~
Ie(x0, y0, z0)in Fig. 12 flows in ~eϕ=−sin ϕ0~ex+ cos ϕ0~eydirection, therefore it is
~
Ie=−Ixsin ϕ0~ex+Iycos ϕ0~ey
Assuming a very thin loop, enables us to assume that Ix=Ix=I0, where I0is constant, which
yields:
~
A=µ0I0a
4πZ2π
0
e−β(r−asin ϑcos ϕ0)
(r−asin ϑcos ϕ0)dϕ0~eϕ
=−µ0I0a
4πZ2π
0
sin ϕ0e−β(r−asin ϑcos ϕ0)
(r−asin ϑcos ϕ0)dϕ0~ex
+µ0I0a
4πZ2π
0
cos ϕ0e−β(r−asin ϑcos ϕ0)
(r−asin ϑcos ϕ0)dϕ0~ey(30)
For a small loop, the function f=e−β(r−asin ϑcos ϕ0)
(r−asin ϑcos ϕ0)can be expanded in Maclaurin Series about a= 0
using:
f=f(0) + f0(0)a+1
2!f00(0)a2+···+1
(n−1)!fn−1an−1+. . .
where f0=∂f
∂a |a=0, f00 =∂2f
∂a2|a=0, . . . , fm=∂mf
∂am|a=0 .We consider only the first two terms of the
expansion:
f(0) = e−βr
r
f0(0) = aβ
r+1
r2sin ϑcos ϕ0e−βr
r
f≃1
r+aβ
r+1
r2sin ϑcos ϕ0e−βr
r(31)
B On small circular loop Antenna 127
Using the previous relations in ~
A
~
A=µ0I0a
4πZ2π
0 −sin ϕ0e−βr
r1
r+aβ
r+1
r2sin ϑcos ϕ0dϕ0
+ cos ϕ0e−βr
r1
r+aβ
r+1
r2sin ϑcos ϕ0dϕ0!(32)
and observing that the integration R2π
0sin ϕ0cos ϕ0= 0 the Eq. (32) can be further simplified to
~
A=Aϕ~eϕ=µ0I0a
4πZ2π
0
cos ϕ0e−βr
r 1
r+aβ
r+1
r2sin ϑcos ϕ0!dϕ0
=µ0I0a
4π
e−βr
rsin ϑZ2π
0 1
r+aβ
r+1
r2cos2ϕ0!dϕ0~eϕ
=µ0I0a
4
e−βr
r 1
r+aβ
r+1
r2!sin ϑ ~eϕ
≃µ0I0a2
4
e−βr
rβ
r+1
r2sin ϑ ~eϕ(33)
Accordingly we have
Ax=−sin ϕAϕ
Ay= cos ϕAϕ(34)
Az= 0
which can be expressed in the cartesian coordinates x, y, z
Ax=−a
22µ0I0e−βr β
r+1
r2
px2+y2
r
y
px2+y2
=−a
22µ0I0e−βr y
rβ
r+1
r2(35)
Similarly we obtain
Ay=a
22µ0I0e−βr β
r+1
r2
px2+y2
r
x
px2+y2
=a
22µ0I0e−βr x
rβ
r+1
r2(36)
128 Appendix
The magnetic field components can now be obtained using the relation:
~
H=Hx~ex+Hy~ey+Hz~ez=1
µ0∇× ~
A(37)
where
∇× ~
A=
70
∂Az
∂y −∂Ay
∂z
~ex+
∂Ax
∂z −0
∂Az
∂x
~ey+ ∂Ay
∂x −∂Ax
∂y !~ez(38)
B.1 The Magnetic Fields
Hx=−1
µ0
∂Ay
∂z
=−a
22I0x∂
∂z β
r2+1
r3e−βr!
=−a
22I0x β
r2+1
r3∂e−βr
∂z +e−βr ∂
∂z β
r2+1
r3!
=−a
22I0x β
r2+1
r3−β z
re−βr−e−βr 2β z
r4+3z
r5!
=a
22I0 3x z
r5+β 3x z
r4−β2x z
r3!e−βr (39)
Hy=1
µ0
∂Ax
∂z
=−a
22I0y∂
∂z β
r2+1
r3e−βr!
=−a
22I0y β
r2+1
r3∂e−βr
∂z +e−βr ∂
∂z β
r2+1
r3!
=−a
22I0y β
r2+1
r3−β z
re−βr−e−βr 2β z
r4+3z
r5!
=a
22I0 3y z
r5+β 3y z
r4−β2y z
r3!e−βr (40)
B On small circular loop Antenna 129
Hz=1
µ0∂Ay
∂x −∂Ax
∂y (41)
We evaluate the first term of Eq. (41)
1
µ0
∂Ay
∂x =a
22I0 ∂
∂x xβ
r2+1
r3e−βr!
=a
22I0 ∂
∂x β x
r2+x
r3e−βr!
=a
22I0 β x
r2+x
r3∂e−βr
∂x +e−βr ∂
∂x β x
r2+x
r3!
=a
22I0 −β x
re−βr β x
r2+x
r3+e−βr β r2−2x2
r4+r2−3x2
r5!
=a
22I0 r2−3x2
r5+β r2−3x2
r4+β2x2
r3e−βr!(42)
The second term of Eq. (41) yields
−1
µ0
∂Ax
∂y =a
22I0 ∂
∂y yβ
r2+1
r3e−βr!
=a
22I0 ∂
∂y β y
r2+y
r3e−βr!
=a
22I0 β y
r2+y
r3∂e−βr
∂y +e−βr ∂
∂y β y
r2+y
r3!
=a
22I0 −β x
re−βr β y
r2+y
r3+e−βr β r2−2y2
r4+r2−3y2
r5!
=a
22I0 r2−3y2
r5+β r2−3y2
r4+β2y2
r3e−βr!(43)
130 Appendix
Adding (42) to (43) results in Hz
Hz=1
µ0∂Ay
∂x −∂Ax
∂y
=a
22I0r2−3x2+r2−3y2
r5+β r2−3x2+r2−3y2
r4+β2x2+y2
r3e−βr
=a
22I02r2−3 (x2+y2)
r5+β 2r2−3 (x2+y2)
r4+β2x2+y2
r3e−βr
=a
22I03z2−r2
r5+β 3z2−r2
r4+β2z2−r2
r3e−βr (44)
This time we will use the potential theory ~
E=−ω ~
A− ∇ϕ, to derive the electric fields expressions,
where according to Lorentz Gauge ϕ=−1
ωµ0ε0∇(∇· ~
A), thus we have
~
E=−ω ~
A+1
ωµ0ε0∇(∇· ~
A)(45)
Let’s examine whether the second term of Eq. (45) yields a zero?
∇· ~
A=∂Ax
∂x +∂Ay
∂y +0
∂Az
∂z
=⇒
∂Ax
∂x =−a
22µ0I0y βr
r2+1
r3∂e−βr
∂x +e−βr ∂
∂x β
r2+1
r3!
=−a
22µ0I0y β
r2+1
r3−β x
re−βr+e−βr −β 2x
r4−3x
r5!
=−a
22µ0I0 β2x y
r3−β 3x y
r4−3x y
r5!e−βr
=a
22µ0I0 −β2x y
r3+β 3x y
r4+3x y
r5!e−βr
B On small circular loop Antenna 131
and
∂Ay
∂y =a
22µ0I0x βr
r2+1
r3∂e−βr
∂y +e−βr ∂
∂y β
r2+1
r3!
=a
22µ0I0x β
r2+1
r3−β y
re−βr+e−βr −β 2y
r4−3y
r5!
=a
22µ0I0 β2x y
r3−β 3x y
r4−3x y
r5!e−βr
=⇒
∂Ax
∂x +∂Ay
∂y =∇· ~
A= 0.(46)
The electric fields can be now easily obtained
B.2 The Electric Fields
~
E=−ω ~
A=−ω Ax~ex−ω Ay~ey=Ex~ex+Ey~ey(47)
All we need is to substitute for Eq. (35) and Eq. (36) in their counterparts in Eq. (47) to get the
results:
Ex=−ω Ax
=ω a
22µ0I0yβ y
r2+y
r3e−βr
=µ0ωa
22I0yβ y
r2+y
r3e−βr
=µ0
β
√µ0ε0a
22I0yβ y
r2+y
r3e−βr
=βrµ0
ε0a
22I0yβ y
r2+y
r3e−βr (48)
132 Appendix
Ey=−ω Ay
=−ω a
22µ0I0xβ x
r2+x
r3e−βr
=−µ0ωa
22I0xβ x
r2+x
r3e−βr
=−µ0
β
√µ0ε0a
22I0xβ x
r2+x
r3e−βr
=−βrµ0
ε0a
22I0xβ x
r2+x
r3e−βr (49)
and
Ez=−ω
>0
Az= 0.(50)
C A finite length dipole
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:
iz(x0, y0, z0) =
~ezI0sin β`
2−z00≤z0≤`
2
~ezI0sin β`
2+z0−`
2≤z0≤0
Figure 13: A finite length dipole
For a z-oriented dipole, like this, the magnetic vector potential takes the orientation of the feeding
current, and can be written as:
~
A=A ~ez
which in any point Pin space can be expressed as:
134 Appendix
~
A=µ I0
4πZ0
−`
2
sin β`
2−z0e−βR
Rdz0
+ sin β`
2+z0e−βR
Rdz0~ez(51)
Assuming a very thing wire, implies further simplifications, that is x0=y0= 0, and from any
point on the dipole to a space point P, the distance Rcan be expressed as:
R=p(x−x0)2+ (y−y0)2+ (z−z0)2
=p(x)2+ (y)2+ (z−z0)2(52)
From Maxwell’s equations we know that
∇· ~
B= 0 =⇒~
B=∇× ~
A
therefore, in free space we have
~
B=µ0~
H
~
H=1
µ0∇× ~
A(53)
=⇒
∇× ~
A=
∂Az
∂y −0
∂Ay
∂z
~ex+
0
∂Ax
∂z −∂Az
∂x
~ey+
*0
∂Ay
∂x −∂Ax
∂y
~ez
=∂Az
∂y ~ex+∂Az
∂x ~ey(54)
Accordingly equation (53) can be expressed as
~
H=1
µ0
∂Az
∂y ~ex+1
µ0
∂Az
∂x ~ey=Hx~ex+Hy~ey(55)
Using Euler’s Identity
sin β`
2±z0=eβ(`
2±z0)+e−β(`
2±z0)
2
and (51) in (53), the x-component of the magnetic field yields:
C A finite length dipole 135
Hx=I0
4π(∂
∂y Z0
−`/2
(eβ(`
2+z0)−e−β(`
2+z0))
2
e−βR
Rdz0
+Z`/2
0
(eβ(`
2−z0)−e−β(`
2−z0))
2
e−βR
Rdz0)
=I0
8π(eβ`/2Z0
−`/2
∂
∂y e−jβ(R−z0)
Rdz0
−e−β`/2Z0
−`/2
∂
∂y e−jβ(R+z0)
Rdz0
+eβ`/2Z`/2
0
∂
∂y e−jβ(R+z0)
Rdz0
−e−β`/2Z`/2
0
∂
∂y e−jβ(R−z0)
Rdz0)(56)
Now we examine the terms in (56) closely, the term eβ`/2R`/2
0
∂
∂y e−jβ(R+z0)
Rdz0for instance can
be rearranged as:
eβ`/2Z`/2
0
∂
∂y e−jβ(R+z0)
Rdz0=eβ`/2Z`/2
01
R
∂
∂y e−jβ(R+z0)+e−jβ(R+z0)∂
∂y 1
R
=y eβ`/2Z`/2
0−β e−jβ(R+z0)
R2−e−jβ(R+z0)
R3dz0(57)
The term in (57) can considerably simplified if we tailored a suitable differentiation that exactly
matches it, therefore consider the the following differentiation:
∂
∂z0e−jβ(R+z0)
R(R+z0−z)=∂
∂z0R−1(R+z0−z)−1e−jβ(R+z0)
=e−jβ(R+z0) 1
(R+z0−z)
∂R−1
∂z0+1
R
∂
∂z0(R+z0−z)−1
−β ∂
∂z0
(R+z0)
R(R+z0−z)!(58)
136 Appendix
The first term in (58) yields:
1
(R+z0−z)
∂R−1
∂z0=z−z0
R3(R+z0−z)(59)
The second term in (58) yields:
1
R
∂
∂z0(R+z0−z)−1=−1
R2
1
(R+z0−z)(60)
and eventually the third term yields:
−β
R(R+z0−z)
∂
∂z0(R+z0) = −β1
R2(61)
Therefore their sum yields:
∂
∂z0e−jβ(R+z0)
R(R+z0−z)=e−jβ(R+z0) z−z0
R3(R+z0−z)
−1
R2(R+z0−z)−β 1
R2!
=e−jβ(R+z0) −1
R3−β 1
R2!(62)
We, thus can use (62) in (57) and get:
C A finite length dipole 137
eβ`/2Z`/2
0
∂
∂y e−jβ(R+z0)
Rdz0=y eβ`/2Z`/2
0−β e−jβ(R+z0)
R2−e−jβ(R+z0)
R3dz0
=y eβ`/2Z`/2
0−β e−jβ(R+z0)
R2−e−jβ(R+z0)
R3dz0
=y eβ`/2Z`/2
0
∂
∂z0e−jβ(R+z0)
R(R+z0−z)dz0
=y eβ`/2"e−jβ(R+z0)
R(R+z0−z)
`/2
0
(63)
For further simplifications, of (62) we introduce two more dimensions R1and R2, which represent
each the distance between the space point Pand the upper and lower tips of the dipole, as shown
in Fig. 14, these are defined as:
R1=rx2+y2+ (z−`
2)2
R2=rx2+y2+ (z+`
2)2
Figure 14: A finite length dipole
Using (52), R1and R2in (63) we get:
138 Appendix
eβ`/2Z`/2
0
∂
∂y e−jβ(R+z0)
Rdz0=y eβ`/2 e−jβ(R+z0)
R1(R1+`
2−z)−e−jβr
r(r−z)!(64)
In a similar manner, eβ`/2R0
−`/2
∂
∂y e−jβ(R−z0)
Rdz0can also be computed:
eβ`/2Z0
−`/2
∂
∂y e−jβ(R−z0)
Rdz0=eβ`/2Z0
−`/21
R
∂
∂y e−jβ(R−z0)+e−jβ(R−z0)∂
∂y 1
R
=y eβ`/2Z0
−`/2−β e−jβ(R−z0)
R2−e−jβ(R−z0)
R3dz0(65)
Again finding a suitable differentiation for this term simplifies its integration considerably, regard
this term:
∂
∂z0−e−jβ(R−z0)
R(R−z0+z)=−∂
∂z0R−1(R−z0+z)−1e−jβ(R−z0)
=−e−jβ(R−z0) 1
(R−z0+z)
∂R−1
∂z0+1
R
∂
∂z0(R−z0+z)−1
−β ∂
∂z0
(R−z0)
R(R−z0+z)!(66)
The first term in (66) yields:
1
(R−z0+z)
∂R−1
∂z0=z−z0
R3(R−z0+z)(67)
The second term in (66) yields:
1
R
∂
∂z0(R−z0+z)−1=−1
R2
1
(R−z0+z)(68)
and eventually the third term yields:
C A finite length dipole 139
−β
R(R−z0+z)
∂
∂z0(R−z0) = β1
R2(69)
together they make the sum
∂
∂z0−e−jβ(R−z0)
R(R−z0+z)=−e−jβ(R−z0) z−z0
R3(R−z0+z)
−1
R2(R−z0+z)−β 1
R2!
=e−jβ(R−z0) −1
R3−β 1
R2!(70)
We, thus can use (70) in (65) and get:
eβ`/2Z0
−`/2
∂
∂y e−jβ(R−z0)
Rdz0=y eβ`/2Z0
−`/2−β e−jβ(R−z0)
R2−e−jβ(R−z0)
R3dz0
=y eβ`/2Z0
−`/2−β e−jβ(R−z0)
R2−e−jβ(R−z0)
R3dz0
=y eβ`/2Z0
−`/2
∂
∂z0e−jβ(R−z0)
R(R−z0+z)dz0
=y eβ`/2"e−jβ(R−z0)
R(R−z0+z)
0
−`/2
=y eβ`/2 e−βr
r(r+z)+e−β(R2+`
2)
R2(R2+`
2+z)!(71)
The term −e−β`/2R`/2
0
∂
∂y e−jβ(R−z0)
Rdz0can also be computed:
−e−β`/2Z`/2
0
∂
∂y e−jβ(R−z0)
Rdz0=−e−β`/2Z`/2
01
R
∂
∂y e−jβ(R−z0)+e−jβ(R−z0)∂
∂y 1
R
=y e−β`/2Z`/2
0β e−jβ(R−z0)
R2+e−jβ(R−z0)
R3dz0(72)
140 Appendix
The differentiation ∂
∂z0 e−β(R−z0)
R(R−z0+z)!fits perfectly for the integrand in (72)
∂
∂z0e−jβ(R−z0)
R(R−z0+z)=∂
∂z0R−1(R−z0+z)−1e−jβ(R−z0)
=e−jβ(R−z0) 1
(R−z0+z)
∂R−1
∂z0+1
R
∂
∂z0(R−z0+z)−1
−β ∂
∂z0
(R−z0)
R(R−z0+z)!
=e−jβ(R−z0) 1
R3+β 1
R2!(73)
Thus
e−β`/2Z`/2
0
∂
∂y e−jβ(R−z0)
Rdz0=y e−β`/2Z`/2
0β e−jβ(R−z0)
R2+e−jβ(R−z0)
R3dz0
=y e−β`/2Z`/2
0
∂
∂z0e−jβ(R−z0)
R(R−z0+z)dz0
=y e−β`/2"e−jβ(R−z0)
R(R−z0+z)
`/2
0
=y e−β`/2 −e−βr
r(r+z)+e−β(R1−`
2)
R1(R1−`
2+z)!(74)
The last term −e−β`/2R0
−`/2
∂
∂y e−jβ(R+z0)
Rdz0can also be computed:
−e−β`/2Z0
−`/2
∂
∂y e−jβ(R+z0)
Rdz0=−e−β`/2Z0
−`/21
R
∂
∂y e−jβ(R+z0)+e−jβ(R+z0)∂
∂y 1
R
=y e−β`/2Z0
−`/2β e−jβ(R+z0)
R2+e−jβ(R+z0)
R3dz0(75)
For integrand in (75), the differentiation ∂
∂z0 −e−β(R+z0)
R(R+z0−z)!fits perfectly:
C A finite length dipole 141
∂
∂z0−e−jβ(R+z0)
R(R+z0−z)=−∂
∂z0R−1(R+z0−z)−1e−jβ(R+z0)
=e−jβ(R−z0) 1
R3+β 1
R2!(76)
Thus
−e−β`/2Z0
−`/2
∂
∂y e−jβ(R+z0)
Rdz0=y e−β`/2Z0
−`/2
∂
∂z0e−jβ(R−z0)
R(R+z0−z)dz0
=y e−β`/2"e−jβ(R−z0)
R(R+z0−z)
0
−`/2
=y e−β`/2 −e−βr
r(r−z)+e−β(R2−`
2)
R2(R2−`
2−z)!(77)
We set (64), (71), (74) and (77) in (56) to evaluate Hx
142 Appendix
Hx=I0
8π"y eβ`/2e−β(R2+`
2)
R2(R2+`
2+z)−e−βr
r(r+z)
+y eβ`/2e−β(R1+`
2)
R1(R1+`
2−z)−e−βr
r(r−z)
+y e−β`/2e−β(R1−`
2)
R1(R1−`
2+z)−e−βr
r(r+z)
+y e−β`/2e−β(R2−`
2)
R2(R2−`
2−z)−e−βr
r(r−z)#
=I0
8π"y eβ`/2e−βR2·e−β `
2
R2(R2+`
2+z)+e−βR1·e−β `
2
R1(R1+`
2−z)
+y e−β`/2e−βR1·e−β `
2
R1(R1−`
2+z)+e−βR2·e−β `
2
R2(R2−`
2−z)
−y eβ`/2e−βr
r(r+z)+e−βr
r(r−z)
−y e−β`/2e−βr
r(r+z)+e−βr
r(r−z)#
=I0
8π"y e−βR1R1+
`
2−
z+R1−
`
2+
z
R1(R1+`
2−z)(R1−`
2+z)
+y e−βR2R2+
`
2−
z+R2−
`
2+
z
R2(R2+`
2+z)(R2−`
2−z)
+y e−βr 2
r2−z22 cos β `
2#
C A finite length dipole 143
Which simplifies to:
Hx=I0
8πy"e−βR1 2
R2
1−`
2−z2!
+e−βR2 2
R2
2−`
2+z2!
−2
r2−z22 cos β `
2e−βr#
=I0
8π"y e−βR2
R2R2+`
2+z+e−βR1
R1R1+`
2−z
+e−βR2
R2R2−`
2−z+e−βR1
R1R1−`
2+z!
−y e−βr e−β`1
r(r+z)+1
r(r−z)
+eβ`1
r(r+z)+1
r(r−z)!#
=I0
8πy"e−βR1 1
R1R1+`
2−z+1
R1R1−`
2+z!
+e−βR2 1
R2R2+`
2+z+1
R2R2−`
2−z!
− r+
z+r−
z
r(r+z)(r−z)e−βreβ` +e−β`!#
Therefore Hxcan be given as:
Hx=I0
4π
y
x2+y2 e−βR1+e−βR2−2 cos β`
2e−βr!(78)
Similarly Hy=−1
µ0
∂
∂x can be also derived:
144 Appendix
Hy=−I0
4π
x
x2+y2 e−βR1+e−βR2−2 cos β`
2e−βr!(79)
and
Hz= 0 (80)
Deriving expressions for the electric field in free space follows by using Ampere-Maxwell’s curl
equation:
∇× ~
H=
0
~
J+ω0~
E
=⇒
~
E=1
ω0∇× ~
H
=Ex~ex+Ey~ey+Ez~ez
and
∇× ~
H=
0
∂Hz
∂y −∂Hy
∂z
~ex+
∂Hx
∂z −
0
∂Hz
∂x
~ey+ ∂Hy
∂x −∂Hx
∂y !~ez
Accordingly,
C A finite length dipole 145
Ex=−1
ω0
∂Hy
∂z
=−I0
4πω0
x
x2+y2
∂
∂z e−βR1+e−βR2−2 cos β`
2e−βr!
=−I0
4πω0
x
x2+y2 −β e−βR1∂R1
∂z −β e−βR2∂R2
∂z + 2β cos β`
2e−βr ∂r
∂z !
=I0
4πω0
x
x2+y2β e−βr z−`
2
R1
e−β(R2−r)+z+`
2
R2
e−β(R2−r)−2z
rcos β`
2!
=I0
4πω0
x
x2+y2β e−βr z−`
2
R1
e−β(R2−r)+z+`
2
R2
e−β(R2−r)
−2z
rcos β`
2!(81)
and
Ey=1
ω0
∂Hx
∂z
=−I0
4πω0
y
x2+y2
∂
∂z e−βR1+e−βR2−2 cos β`
2e−βr!
=−I0
4πω0
y
x2+y2 −β e−βR1∂R1
∂z −β e−βR2∂R2
∂z + 2β cos β`
2e−βr ∂r
∂z !
=I0
4πω0
y
x2+y2β e−βr z−`
2
R1
e−β(R2−r)+z+`
2
R2
e−β(R2−r)−2z
rcos β`
2!
=I0
4πω0
y
x2+y2β e−βr z−`
2
R1
e−β(R2−r)+z+`
2
R2
e−β(R2−r)
−2z
rcos β`
2!(82)
Finally, we Evaluate the z-component Ez
Ez=1
ω0 ∂Hy
∂x −∂Hx
∂y !(83)
146 Appendix
Evaluation of the first term 1
ω0
∂Hy
∂x of (83) yields:
1
ω0
∂Hy
∂x =−1
ω0
I0
4π
∂
∂x x
x2+y2e−βR1+e−βR2−2 cos β`
2e−βr!
=I0
4πω0 x
x2+y2−βe−βR1∂R1
∂x −βe−βR2∂R2
∂x
+ 2β cos β`
2e−βr∂r
∂x
+e−βR1+e−βR2−2 cos β`
2e−βr∂
∂xx
x2+y2!
=I0
4πω0 x
x2+y2β−e−βR1
R1−e−βR2
R2
+ 22
rcos β`
2e−βr
+y2−x2
x2+y22e−βR1+e−βR2−2 cos β`
2e−βr!(84)
The second term −1
ω0
∂Hx
∂y of (83) yields:
−1
ω0
∂Hx
∂y =−1
ω0
I0
4π
∂
∂y y
x2+y2e−βR1+e−βR2−2 cos β`
2e−βr!
=I0
4πω0 y
x2+y2−βe−βR1∂R1
∂y −βe−βR2∂R2
∂y
+ 2β cos β`
2e−βr∂r
∂y
+e−βR1+e−βR2−2 cos β`
2e−βr∂
∂y y
x2+y2!
=I0
4πω0 y
x2+y2β−e−βR1
R1−e−βR2
R2
+ 22
rcos β`
2e−βr
+x2−y2
x2+y22e−βR1+e−βR2−2 cos β`
2e−βr!(85)
Using (84) and (85) in (83), we obtain Ez
C A finite length dipole 147
Ez=I0
4πω0
β
>1
x2+y2
x2+y2−e−βR1
R1−e−βR2
R2
+ 22
rcos β`
2e−βr
+
*0
y2−x2+x2−y2
x2+y22e−βR1+e−βR2−2 cos β`
2e−βr!
=−I0β
4πω0e−βR1+e−βR2−2 cos β`
2e−βr(86)
D Mathematical Notes 149
D Mathematical Notes
D.1 Taylor Polynomials and Taylor Series
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 fto have (n+ 1) derivatives). We would like to approximate f(x)near the
point x=x0, and we can do it as follows:
f(x) = Pn(x)
|{z}
Taylor polynomial
+Rn(x)
|{z}
remainder term
,
where
Pn(x) = f(x0) + f0(x0)(x−x0) + f00(x0)
2! (x−x0)2+···+f(n)(x0)
n!(x−x0)n
is the nth order Taylor polynomial of fabout x0, and
Rn(x) = f(n+1)(ξ(x))
(n+ 1)! (x−x0)n+1
is the remainder term or truncation error. The point ξ(x)in the error term lies somewhere between
the points x0and 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=x0. 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 Pnand Rnare
bounded. Therefore, there are two ways in which we can think of the Taylor polynomial Pn(x)as
an approximation of f(x):
1. First of all, let us fix n. Then, we can improve the approximation by letting xapproach x0,
since as (x−x0)gets small, the error term Rn(x)goes to zero (nis considered fixed and all
terms depending on nare thus constant). Therefore, the approximation improves when x
gets closer and closer to x0.
2. Alternatively, we can think of fixing x. Then, we can improve the approximation by taking
more and more terms in the series. When nis increased, the factorial in the denominator of
the error term will eventually dominate the (x−x0)n+1 term (regardless of how big (x−x0)
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 x0; and by taking more terms in the series.
150 Appendix
D.2 Numerical techniques to approximate derivatives.
Note
f0(x0) = lim
x→0
f(x0−h)−f(x0)
h
Difference formula:
f(x0−h)−f(x0)
h≈f0(x0), h 6= 0, h small.
Forward-difference formula: h > 0.
Backward-difference formula: h < 0.
That is, forward-difference is used for f0(x0)and backward-difference is used for f0(xn).
Q: How accurate is the difference formula?
By Taylor expansion,
f(x0+h) = f(x0) + h f0(x0) + 1
2h2f00(ξ), ξ in between x0and x0+h
=⇒f(x0+h)−f(x0)
h−f0(x0) = 1
2h f00(ξ)
Truncation error: 1
2h f00(ξ) = O(h)
Order of accuracy: 1.
Q: How to get more accurate approximation?
Using interpolation polynomial.
Given data {(xi, yi)}n
i=0, we construct the interpolation polynomial pn(x). Without loosing gener-
ality, assume x0< dots < xn. Then ∀x∈[x0, xn],
p0(x)≈f0(x)
Error analysis:
Define
D Mathematical Notes 151
Φn(x) =
n
Y
i=0
(x−xi)
Then it can be proved that
f0(x)−p0(x) = (O(hn),if Φ0
n(x)6= 0
O(hn+1),if Φ0
n(x) = 0
If the truncation error is O(hn), we say the method has order of accuracy n.
Consider n= 1 :Lagrange form:
p1(x) = x−x1
x0−x1
f(x0) + x−x0
x1−x0
f(x1)
=⇒f0(x)≈p0
1(x) = −f(x0)
h+f(x1)
h=f(x1)−f(x0)
h
Now based on the above error analysis we may use this formula for any x∈[x0, x1]and trunca-
tion error is O(h)if Φ0
1(x)6= 0 and O(h2)if Φ0
1(x) = 0.
Q: What point in [x0, x1]makes Φ0
1(x) = 0?
Only mid point x= (x0+x1)/2
Φ0
1(x) = (x−x0)+(x−x1)
=⇒x= (x0+x1)/2
Truncation error for p1(x):
¬O(h)for any x6= (x0+x1)/2.In particular, truncation error is O(h)for x0and x1.
O(h2)for any x= (x0+x1)/2.
For x= (x0+x1)/2we can write x0=x−1
2hand x1=x+1
2h. We thus get the central-difference
formula:
f0(x)≈f(x+1
2h)−f(x−1
2h)
h.
Truncation error for central-difference: O(h2)
152 Appendix
What does O(hp), or Big O mean?
O(x)is read as “ Big O of x” or “ order x”, and similarly O(xp)is “ Big O of xp” or “ order xp”.
Note that for our purpose pwill 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 Cand
> 0such that
|g(x)| ≤ C|xp|,∀|x|< (87)
then we say that
g(x) = O(xp)(88)
Notice that the equality in Equation (88) is a slight abuse of the notation, as O(xp)is not a function, but
rather a class of functions. The pedantic should write g(x)∈O(xp)but since nobody wants to admit
pedantry the equality sign is used universally.
D.3 Summary, revision and more
Consider the initial value problem for a system of nordinary differential equations of first order:
y0=f(t, y), y(t0) = y0,(89)
on the interval [a, b]where the function f(t, y) :
f:R×Rn→Rn,(90)
is continuous in tand Lipschitz continuous in y, that is,
|f(t, y)−f(t, x)|< M|y−x|,(91)
for some positive constant M. Under these conditions, problem (89) admits one and only one
solution, y(t), t ∈[a, b].
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
Simple, Low Order Methods Error
Explicit Euler: wi+1 =wi+h f(ti, wi)O(h)
Implicit Euler: wi+1 =wi+h f(ti+1, wi+1)O(h)
Trapezoidal Method: wi+1 =wti+h
2[f(ti, wi) + f(ti+1, wi+1)] O(h2)
Modified Euler Method Error
a)Predictor: wp(i+1) = wi+h f(ti, wi)
b)Corrector: ˜wi+1 =wi+h
2f(ti, wi) + f(ti+1, wp
i+1)O(h2)
Classical Runge-Kutta Method Error
set k1:k1=h f(ti, y(ti))
set k2:k2=h f(ti+1
2h, y(ti) + 1
2k1)
set k3:k3=h f(ti+1
2h, y(ti) + 1
2k2)
set k4:k4=h f(ti+h, y(ti) + k3)
wi+1 =wi+1
6[k1+ 2k2+ 2k3+k4]O(h4)
Table 1: Some available Methods
Notes: Things to bear in mind
¬widenotes the numerical solution, while either y(ti)or yistands for the exact solution.
An ODE solver needs to produce a numerical solution wiof system (89) which converges to
the exact solution y(ti)as h↓0,with h=ti+1 −ti=b−a
Nfor each ti∈[a, b], and remains
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 fevaluations.
°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 wi+1, wi, wi−1, . . . in a linear way to achieve higher ac-
curacy, but sacrifice the one-step format (which is easier to understand and to implement).
154 Appendix
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.
Let f(x, y)be continuous for a≤x≤band all y, and satisfy the Lipschitz condition |f(x, u)−f(x, v)| ≤
L|u−v|for a≤x≤b, all u, v. Then there exists a unique solution y(x)to the initial value problem
y0(x) = f(x, y(x)) for a≤x≤b, y(a) = A.
Picard’s Theorem also holds for systems of equations like
y01(x) = f1(x, y1, y2, . . . , yn), y1(a) = A1,
.
.
. with initial conditions .
.
.
y0n(x) = fn(x, y1, y2, . . . , yn), yn(a) = An,
where each fisatisfies a Lipschitz condition in each yj. Note that any higher order initial value problem
y(n)=g(x, y, y0, . . .,y(n−1)),a≤x≤b,y(a) = A1,y0(a) = A2, . . . , y(n−1)(a) = Ancan always be
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×nJacobian matrix
J=∂yf(t, y) = ∂fi
∂yjfor each i, j = 1, . . . , n
We assume that the neigenvalue λ1, . . . , λnof the Jmatrix have negative real parts, <λj<0,
and ordered as follows:
<λn≤ ··· ≤ <λ2≤ <λ1<0
W
Stiffness ratio The stiffness ratio of the system y0=f(t, y)is the positive number
r=<λn
<λ1
where the eigenvalues of the Jacobian matrix of the system satisfy the relation
y0=λy, <λ < 0
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 Icontaining tif in Ithe 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 un-
stable 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 math-
ematical 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.
156 Appendix
–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 xthe input data and εthe error in the input
data,
∗one says that the algorithm is “numerically stable” for the absolute error if
·x−(x+)≃f(x)−f(x+)
∗and numerically stable for the relative error if
·x−(x+)
x≃f(x)−f(x+)
f(x).
∗One says that an algorithm is “numerically unstable” for the absolute error if
·x−(x+)<< f(x)−f(x+)
∗and numerically unstable for the relative error if
·x−(x+)
x<< f(x)−f(x+)
f(x).
– Absolute Error Given some value aand an approximation bof a, the “absolute error”
is
:= |a−b|
– Absolute Error
η:= |a−b|
|a|
where the vertical bars denote the absolute value.
More on the Phenomenon of Stiffness Many time dependent phenomena which can be mod-
eled 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 sys-
tem 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 corre-
sponding 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 signif-
icantly 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 vector-
valued 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:Rn→Rmis a function from Euclidean n-space to Euclidean m-space. Such a func-
tion is given by m real-valued component functions, y1(x1, ..., xn), . . . ym(x1, ..., xn). The partial
derivatives of all these functions (if they exist) can be organized in an m×nmatrix, the Jacobian
matrix of F, as follows:
∂y1
∂x1··· ∂y1
∂xn
.
.
.....
.
.
∂ym
∂x1··· ∂ym
∂xn
This matrix is denoted by :
JF(x1, . . . , xn)or by ∂(y1, . . . , ym)
∂(x1, . . . , xn)
The ith row of this matrix is given by the gradient of the function yifor i= 1, . . . , m.
If pis a point in Rnand Fis differentiable at p, then its derivative is given by JF(p)(and this is the
easiest way to compute said derivative). In this case, the linear map described by JF(p)is the best
linear approximation of F near the point p, in the sense that
F(x)≈F(p) + JF(p)·(x−p).
for xclose to p.
Example 3. The Jacobian matrix of the function F:R3→R4with components:
y1=x1
y2= 5x3
y3= 4(x2)2−2x3
y4=x3sin(x1)
is:
JF(x1, x2, x3) =
1 0 0
0 0 5
0 8x2−2
x3cos(x1) 0 sin(x1)
158 Appendix
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:
y0(t) = J(t)y(t) + p0(t)
The analytical solution is
y(t) = (A−p(0)) eJ(t)+p(t)
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,
u0(t)
v0(t)=998 1998
−998 −1998·u(t)
v(t)with initial condition u(0)
v(0)=1
0
the analytical solution to this system is
u(t) = 2e−t
|{z}
steady −e−1000t
| {z }
transient
v(t) = −e−t
|{z}
steady
+e−1000t
| {z }
transient
The solutions uand vhave a second term e−1000tthat has a negligible effect on the solution for t > 0but
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, h1000.
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 159
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
afinite Length, we studied pretty much complicated cases actually, that comprise arrays of dipoles,
of every possible orientation, location and configuration.
E The Algorithms: 161
E The Algorithms:
E.1 One-step methods: fixed step-size
Euler’s Method in Mathematica code
•Having IVP Equations System
Euler[h , t0, x0, y0, z0, n ] = Block [{k1, l1, m1},
fx[t , x , y , z ] := Form1
fy[t , x , y , z ] := Form2
fz[t , x , y , z ] := Form3
•INPUT the initial values t0, x0, y0, y0, z0,a fixed step-size hand the number of the calculated
points n.
t=t0;x=x0;y=y0;z=z0;pts ={{x, y, z}};
•Build an iterative loop as follows
Do[
k1=h fx[t, x, y, z]; l1=h fy[t, x, y, z]; m1=h fz[t, x, y, z];
x=x+k1;y=y+l1;z=z+m1;pts =AppendTo{pts, {x, y, z}};
t=t+h, {i, 1, n}
];
Print[” ”]; Print[”Euler0s methotd is used with ”, n, ”subintervals.”];
]
Heun’s Method in Mathematica code
•Having IVP Equations System
Heun[h , t0, x0, y0, z0, n ] = Block [{k1, k2, l1, l2, m1, m2},
fx[t , x , y , z ] := Form1fy[t , x , y , z ] := Form2fz[t , x , y , z ] := Form3
162 Appendix
•INPUT the initial values t0, x0, y0, y0, z0,a fixed step-size hand the number of the calculated
points n.
t=t0;x=x0;y=y0;z=z0;pts ={{x, y, z}};
•Build an iterative loop as follows
Do[
k1=h fx[t, x, y, z]; l1=h fy[t, x, y, z]; m1=h fz[t, x, y, z];
k2=h fx[t+2
3h, x +2
3k1, y +2
3l1, z +2
3m1];
l2=h fy[t+2
3h, x +2
3k1, y +2
3l1, z +2
3m1];
m2=h fz[t+2
3h, x +2
3k1, y +2
3l1, z +2
3m1];
x=x+1
4(k1+ 3k2); y=y+1
4(l1+ 3l2); z=z+1
4(m1+ 3m2);
pts =AppendTo{pts, {x, y, z}};t=t+h, {i, 1, n}
];
Print[” ”]; Print[”Heun0s methotd is used with ”, n, ”subintervals.”];
]
The Runge-Kutta Method RK4in Mathematica code
•Having IVP Equations System
RK4[h , t0, x0, y0, z0, n ] = Block [{k1, k2, k3, k4, l1, l2, l3, l4, m1, m2, m3, m4},
fx[t , x , y , z ] := Form1fy[t , x , y , z ] := Form2fz[t , x , y , z ] := Form3
•INPUT the initial values t0, x0, y0, y0, z0,a fixed step-size hand the number of the calculated
points n.
t=t0;x=x0;y=y0;z=z0;pts ={{x, y, z}};ptts ={{t, x, y, z, h}};
•Build an iterative loop as follows
E The Algorithms: 163
Do[
k1=h fx[t, x, y, z]; l1=h fy[t, x, y, z]; m1=h fz[t, x, y, z];
k2=h fx[t+1
2h, x +1
2k1, y +1
2l1, z +1
2m1];
l2=h fy[t+1
2h, x +1
2k1, y +2
2l1, z +1
2m1];
m2=h fz[t+1
2h, x +1
2k1, y +2
2l1, z +1
2m1];
k3=h fx[t+1
2h, x +1
2k2, y +1
2l2, z +1
2m2];
l3=h fy[t+1
2h, x +1
2k2, y +1
2l2, z +1
2m2];
m3=h fz[t+1
2h, x +1
2k2, y +1
2l2, z +1
2m2];
k4=h fx[t+h, x +k3, y +l3, z +m3];
l4=h fy[t+h, x +k3, y +l3, z +m3];
m4=h fz[t+h, x +k3, y +l3, z +m3];
x=x+1
6(k1+ 2k2+ 2k3+k4);
y=y+1
6(l1+ 2l2+ 2l3+l4);
z=z+1
6(m1+ 2m2+ 2m3+l4);
pts =AppendTo{pts, {x, y, z}};
ptts =AppendTo{ptts, {t, x, y, z, h}};
t=t+h, {i, 1, n}
Print[” ”];
Print[”The Runge−Kutta methotd is used with ”, n, ”subintervals.”];
]
E.2 One-step methods: variable step-size
The Runge-Kutta-Fehlberg Method RK45 in Mathematica code
•Having IVP Equations System
RK45[h , t0, x0, y0, z0, n ] = Block [{k1, k2, k3, k4, k5, k6, l1, l2, l3, l4, l5, l6, m1, m2, m3, m4, m5, m6},
fx[t , x , y , z ] := Form1fy[t , x , y , z ] := Form2fz[t , x , y , z ] := Form3
164 Appendix
•INPUT the initial values t0, x0, y0, y0, z0,a fixed step-size hand the number of the calculated
points n.
hmax = 64h;hmin =1
64h;Tol = 10−10;t=t0;x=x0;y=y0;
z=z0;h=hmax;pts ={{x, y, z}};
•Build an iterative loop as follows
Do[
k1=h fx[t, x, y, z]; l1=h fy[t, x, y, z]; m1=h fz[t, x, y, z];
k2=h fx[t+1
4h, x +1
4k1, y +1
4l1, z +1
4m1];
l2=h fy[t+1
4h, x +1
4k1, y +1
4l1, z +1
4m1];
m2=h fz[t+1
4h, x +1
4k1, y +1
4l1, z +1
4m1];
k3=h fx[t+3
8h, x +3
32k1+9
32k2, y +3
32l1+9
32l2, z +3
32m1+9
32m2];
l3=h fy[t+3
8h, x +3
32k1+9
32k2, y +3
32l1+9
32l2, z +3
32m1+9
32m2];
m3=h fz[t+3
8h, x +3
32k1+9
32k2, y +3
32l1+9
32l2, z +3
32m1+9
32m2];
k4=h fx[t+12
13h, x +1932
2197k1−7200
2197k2+7200
2197k3,
y+1932
2197l1−7200
2197l2+7200
2197l3, z +1932
2197m1−7200
2197m2+7200
2197m3]
l4=h fy[t+12
13h, x +1932
2197k1−7200
2197k2+7200
2197k3,
y+1932
2197l1−7200
2197l2+7200
2197l3, z +1932
2197m1−7200
2197m2+7200
2197m3]
m4=h fz[t+12
13h, x +1932
2197k1−7200
2197k2+7200
2197k3,
y+1932
2197l1−7200
2197l2+7200
2197l3, z +1932
2197m1−7200
2197m2+7200
2197m3]
E The Algorithms: 165
k5=h fx[t+h, x +439
216k1−8k2+3680
513 k3−845
4104k4,
y+439
216l1−8l2+3680
513 l3−845
4104l4,
z+439
216m1−8m2+3680
513 m3−845
4104m4]
l5=h fy[t+h, x +439
216k1−8k2+3680
513 k3−845
4104k4,
y+439
216l1−8l2+3680
513 l3−845
4104l4,
z+439
216m1−8m2+3680
513 m3−845
4104m4]
m5=h fz[t+h, x +439
216k1−8k2+3680
513 k3−845
4104k4,
y+439
216l1−8l2+3680
513 l3−845
4104l4,
z+439
216m1−8m2+3680
513 m3−845
4104m4]
k6=h fx[t+1
2h, x −8
27k1+ 2k2−3544
2565k3+1859
4104k4−11
40k5,
y−8
27l1+ 2l2−3544
2565l3+1859
4104l4−11
40l5,
z−8
27m1+ 2m2−3544
2565m3+1859
4104m4−11
40m5]
l6=h fy[t+1
2h, x −8
27k1+ 2k2−3544
2565k3+1859
4104k4−11
40k5,
y−8
27l1+ 2l2−3544
2565l3+1859
4104l4−11
40l5,
z−8
27m1+ 2m2−3544
2565m3+1859
4104m4−11
40m5]
m6=h fz[t+1
2h, x −8
27k1+ 2k2−3544
2565k3+1859
4104k4−11
40k5,
y−8
27l1+ 2l2−3544
2565l3+1859
4104l4−11
40l5,
z−8
27m1+ 2m2−3544
2565m3+1859
4104m4−11
40m5]
σ1=N
"Abs 1
360 k1−128
4275 k3−75240
360 k4+1
50 k5+2
55 k6;
h#
σ2=N
"Abs 1
360 l1−128
4275 l3−75240
360 l4+1
50 l5+2
55 l6;
h#
σ3=N
"Abs 1
360 m1−128
4275 m3−75240
360 m4+1
50 m5+2
55 m6;
h#
σ=N
qσ2
1+σ2
2+σ2
3;
166 Appendix
If[σ6= 0,
δ= 0.84 Tol
σ1
4;
If[δ > 4, h = 4h, h =δ h];
If[δ < .1, h =.1h, h =.7δ h];
If[h>hmax, h =.313hmax];
If[h<hmin, Break[ ]];
];
If[δ≤Tol,
x=x+25
216k1+1408
2565k3+2197
4104k4−1
5k5;
y=y+25
216l1+1408
2565l3+2197
4104l4−1
5l5;
z=z+25
216m1+1408
2565m3+2197
4104m4−1
5m5;
pts =AppendTo{pts, {x, y, z}};
];
t=t+h,
{i, 1, n}
Print[” ”];
Print[”The Runge−Kutta−Fehlberg methotd is used with ”, n, ”subintervals.”];
]
E The Algorithms: 167
E.3 Multi-step methods: fixed step-size
The Adams-Bashforth-Moulton Method ABM in Mathematica code
•Having IVP Equations System
ABM[h , t0, x0, y0, z0, n ] = Block [{F1, F2, F3, F4, F5, F6, F7, F8, F9, F10, F11},
fx[t , x , y , z ] := Form1fy[t , x , y , z ] := Form2fz[t , x , y , z ] := Form3
•INPUT the initial values t0, x0, y0, y0, z0,a fixed step-size hand the number of the calculated
points n.
t=t0;x=x0;y=y0;z=z0;pts ={{x, y, z}};ptts ={{t, x, y, z, h}};
•Run a starter sub-algorithm RK4RK4[h, t, x, y, z, 3];
•Now start the Explicit-Implicit Predictor-Corrector main algorithm:
Do[
t=t+h;
F0=fx[ptts [[i−3,1]], ptts[[i−3,2]], ptts[[i−3,3]], ptts[[i−3,4]], ptts[[i−3,5]]] ;
F1=fx[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F2=fx[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F3=fx[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
F4=fy[ptts [[i−3,1]], ptts[[i−3,2]], ptts[[i−3,3]], ptts[[i−3,4]], ptts[[i−3,5]]] ;
F5=fy[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
168 Appendix
F6=fy[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F7=fy[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
F8=fz[ptts [[i−3,1]], ptts[[i−3,2]], ptts[[i−3,3]], ptts[[i−3,4]], ptts[[i−3,5]]] ;
F9=fz[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F10 =fz[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F11 =fz[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
Explicit Adams-Bashforth Predictor
p1=x+h
24 (−9F0+ 37F1−59F2+ 55F3) ;
p2=x+h
24 (−9F4+ 37F5−59F6+ 55F7) ;
p3=x+h
24 (−9F8+ 37F9−59F10 + 55F11) ;
Implicit Adams-Moulton Corrector
x=x+h
24 (F1−5F2+ 19F3+ 9fx[t, p1, p2, p3, p4]) ;
y=x+h
24 (F5−5F6+ 19F7+ 9fy[t, p1, p2, p3, p4]) ;
z=x+h
24 (F9−5F10 + 19F11 + 9fz[t, p1, p2, p3, p4]) ;
pts =AppendTo{pts, {x, y, z}};
ptts =AppendTo{ptts, {t, x, y, z, h}},{i, 4, n}
Print[” ”];
Print[”The Adams−Bashforth−Moulton methotd is used with ”, n, ”subintervals.”];
]
E The Algorithms: 169
The Milne-Simpson Method MS in Mathematica code
•Having IVP Equations System
ABM[h , t0, x0, y0, z0, n ] = Block [{F1, F2, F3, F4, F5, F6, F7, F8, F9, F10, F11},
fx[t , x , y , z ] := Form1fy[t , x , y , z ] := Form2fz[t , x , y , z ] := Form3
•INPUT the initial values t0, x0, y0, y0, z0,a fixed step-size hand the number of the calculated
points n.
t=t0;x=x0;y=y0;z=z0;pts ={{x, y, z}};ptts ={{t, x, y, z, h}};
p1old = 0; p2old = 0; p3old = 0; xold = 0; yold = 0; zold = 0;
•Run a starter sub-algorithm RK4RK4[h, t, x, y, z, 3];
•Now start the Explicit-Implicit Predictor-Corrector main algorithm:
Do[
t=t+h;
F1=fx[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F2=fx[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F3=fx[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
F5=fy[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F6=fy[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F7=fy[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
F9=fz[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F10 =fz[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F11 =fz[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
Explicit Milne’s Predictor
170 Appendix
p1new =ptts[[i−3,2]] + 4h
3(2F1−F2+ 2F3) ;
p2new =ptts[[i−3,3]] + 4h
3(2F5−F6+ 2F7) ;
p3new =ptts[[i−3,4]] + 4h
3(2F9−F10 + 2F11) ;
p1mod =p1new +28 (xold −p1old )
29 ;
p2mod =p2new +28 (yold −p2old )
29 ;
p3mod =p3new +28 (zold −p3old )
29 ;
F0=fx[t, p1mod , p2mod , p3mod ];
F4=fy[t, p1mod , p2mod , p3mod ];
F8=fz[t, p1mod , p2mod , p3mod ];
Implicit Simpson’s Corrector
x=ptts[[i−1,2]] + h
3(F2+ 4F3+F0) ;
y=ptts[[i−1,3]] + h
3(F6+ 4F7+F4) ;
z=ptts[[i−1,4]] + h
3(F10 + 4F11 +F8) ;
pts =AppendTo{pts, {x, y, z}};ptts =AppendT o{ptts, {t, x, y, z, h}};
p1old =p1new ;p2old =p2new ;p3old =p3new ;xold =x;yold =y;zold =z;
F1=F2;F2=F3;F3=fx[t, x, y, z]; F5=F6;F6=F7;
F7=fy[t, x, y, z]; F9=F10;F10 =F11;F11 =fz[t, x, y, z],{i, 4, n}
Print[” ”];
Print[”The Adams−Bashforth−Moulton methotd is used with ”, n, ”subintervals.”];
]
E The Algorithms: 171
The Hamming’s Method in Mathematica code
•Having IVP Equations System
Ham[h , t0, x0, y0, z0, n ] = Block [{F1, F2, F3, F4, F5, F6, F7, F8, F9, F10, F11},
fx[t , x , y , z ] := Form1fy[t , x , y , z ] := Form2fz[t , x , y , z ] := Form3
•INPUT the initial values t0, x0, y0, y0, z0,a fixed step-size hand the number of the calculated
points n.
t=t0;x=x0;y=y0;z=z0;pts ={{x, y, z}};ptts ={{t, x, y, z, h}};
p1old = 0; p2old = 0; p3old = 0; c1old = 0; c2old = 0; c3old = 0;
•Run a starter sub-algorithm RK4RK4[h, t, x, y, z, 3];
•Now start the Explicit-Implicit Predictor-Corrector main algorithm:
Do[
t=t+h;
F1=fx[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F2=fx[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F3=fx[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
F5=fy[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F6=fy[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F7=fy[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
F9=fz[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F10 =fz[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F11 =fz[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
Explicit Hamming’s Predictor
172 Appendix
p1new =ptts[[i−3,2]] + 4h
3(2F1−F2+ 2F3) ;
p2new =ptts[[i−3,3]] + 4h
3(2F5−F6+ 2F7) ;
p3new =ptts[[i−3,4]] + 4h
3(2F9−F10 + 2F11) ;
p1mod =p1new +112 (c1old −p1old )
121 ;
p2mod =p2new +112 (c2old −p2old )
121 ;
p3mod =p3new +112 (c3old −p3old )
121 ;
F0=fx[t, p1mod , p2mod , p3mod ];
F4=fy[t, p1mod , p2mod , p3mod ];
F8=fz[t, p1mod , p2mod , p3mod ];
Implicit Hamming’s Corrector
c1new =1
8(9 ptts[[i, 2]] −ptts[[i−2,2]]) + 3h
8(−F2+ 2F3+F0);
c2new =1
8(9 ptts[[i, 3]] −ptts[[i−2,3]]) + 3h
8(−F6+ 2F7+F4);
c3new =1
8(9 ptts[[i, 4]] −ptts[[i−2,4]]) + 3h
8(−F10 + 2F11 +F8);
x=c1new +9
121 (p1new −c1new ) ;
y=c2new +9
121 (p2new −c2new ) ;
z=c3new +9
121 (p3new −c3new ) ;
pts =AppendTo{pts, {x, y, z}};
ptts =AppendTo{ptts, {t, x, y, z, h}};
p1old =p1new ;p2old =p2new ;p3old =p3new ;c1old =c1new ;c2old =c2new ;c3old =c3new ;
F1=F2;F2=F3;F3=fx[t, x, y, z]; F5=F6;
F6=F7;F7=fy[t, x, y, z]; F9=F10;F10 =F11;F11 =fz[t, x, y, z],{i, 4, n}
Print[” ”];
Print[”The Hamming0s methotd is used with ”, n, ”subintervals.”];
]
E The Algorithms: 173
E.4 Multi-step methods: variable step-size
The Adams-Bashforth-Moulton Variable step-size Method ABMV in Mathematica code
•Having IVP Equations System
ABMV [h , t0, x0, y0, z0, n ] = Block [{F1, F2, F3, F4, F5, F6, F7, F8, F9, F10, F11},
fx[t , x , y , z ] := Form1fy[t , x , y , z ] := Form2fz[t , x , y , z ] := Form3
•INPUT the initial values t0, x0, y0, y0, z0,a fixed step-size hand the number of the calculated
points n.
hmax = 64h;hmin =1
64h;Tol = 10−10;t=t0;x=x0;y=y0;z=z0;
h=hmax;pts ={{x, y, z}};ptts ={{t, x, y, z, h}};
•Run a starter sub-algorithm RK4RK4[h, t, x, y, z, 3];
•Now start the Explicit-Implicit Predictor-Corrector main algorithm:
Do[
t=t+h;
F0=fx[ptts [[i−3,1]], ptts[[i−3,2]], ptts[[i−3,3]], ptts[[i−3,4]], ptts[[i−3,5]]] ;
F1=fx[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F2=fx[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F3=fx[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
F4=fy[ptts [[i−3,1]], ptts[[i−3,2]], ptts[[i−3,3]], ptts[[i−3,4]], ptts[[i−3,5]]] ;
F5=fy[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F6=fy[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F7=fy[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
F8=fz[ptts [[i−3,1]], ptts[[i−3,2]], ptts[[i−3,3]], ptts[[i−3,4]], ptts[[i−3,5]]] ;
F9=fz[ptts [[i−2,1]], ptts[[i−2,2]], ptts[[i−2,3]], ptts[[i−2,4]], ptts[[i−2,5]]] ;
F10 =fz[ptts [[i−1,1]], ptts[[i−1,2]], ptts[[i−1,3]], ptts[[i−1,4]], ptts[[i−1,5]]] ;
F11 =fz[ptts [[i, 1]], ptts[[i, 2]], ptts[[i, 3]], ptts[[i, 4]], ptts[[i, 5]]] ;
174 Appendix
Explicit Adams-Bashforth Predictor
p1=x+h
24 (−9F0+ 37F1−59F2+ 55F3) ;
p2=x+h
24 (−9F4+ 37F5−59F6+ 55F7) ;
p3=x+h
24 (−9F8+ 37F9−59F10 + 55F11) ;
Implicit Adams-Moulton Corrector
c1=x+h
24 (F1−5F2+ 19F3+ 9fx[t, p1, p2, p3, p4]) ;
c2=x+h
24 (F5−5F6+ 19F7+ 9fy[t, p1, p2, p3, p4]) ;
c3=x+h
24 (F9−5F10 + 19F11 + 9fz[t, p1, p2, p3, p4]) ;
σ=19 h
270 p(c1−p1)2+ (c2−p2)2+ (c3−p3)2;
If[σ6= 0,
δ= 0.84 Tol
σ1
4;
If[δ > 4, h = 4h, h =δ h];
If[δ < .1, h =.1h, h =.7δ h];
If[h > hmax, h =.313hmax];
If[h<hmin, Break[ ]];
];
If[δ≤Tol,
x=c1;y=c2;z=c3;
pts =AppendTo{pts, {x, y, z}};
ptts =AppendTo{ptts, {t, x, y, z, h}}
],{i, 4, n}
];
Print[” ”];
Print[”The V araible−step Adams−Bashforth−Moulton method is used with ”, n, ”subintervals.”];
]
175
Bibliography
[1] Richard Feynmann, The Feynmann Lectures on Physics, Volume 2, Benjamin Cummings 1963,
ISBN: 0-201-02117-X
[2] H. Hertz, Die Kr¨afte 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¨
onich 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
Index
Index
Accuracy, 150
algorithms, 161
- Multi-step methods, 167
Adams-Bashforth-Moulton Method, 167
Hamming’s Method, 171
Milne-Simpson Method, 169
variable Adams-Bashforth-Moulton Method,
173
- One-step methods, 161
Euler’s Method , 161
Heun’s Method , 161
Runge-Kutta Method , 162
Runge-Kutta-Fehlberg Method , 163
antenna
- linear dipole , xiii
- lines of force, xiii
field lines, xiii
- Radiation Patterns, xiv
- the image theory, xiv
antennas
- the image theory, 35
arbitrarily oriented
- finite length, 70
electric, 70
magnetic, 70
- Hertzian, 69
Electric, 69
Magnetic, 69
- Radiators, 57
Big O notation, 152
Biot-Savart Law, 7
continuity equation, 18
density
- current ~
J, 5
- electric charge ρe, 5, 18
- electric current Je, 18
- magnetic charge ρm, 18
- magnetic current Jm, 18
- magnetic flux ~
B, 7
detaching
- points, 6
- process, 5
dipole
- finite length, 4
- Hertzian , xiv
- radiation, xiii
duality, 28
EM
- Radiation, 2
- experimental confirmation, 17
- spectrum, 17
excitation, 22
field
- detachment, 4
- retardation effect, 9
frequency domain, 22
H. Hertz, 29
hodograph, 43
Jacobian matrix, 157
loop
- kidney-shaped , 4
- small, 4
Maxwell’s
178 Index
- Ampere Equation, 10
- classes and forms, 15
- differential Equations, 18
- Equations, 15
- integral Equations, 19
- integro-differential equations, 17, 22
nonperiodic, 23
Numerical
- approximate derivatives , 150
- Approach, 73
- Phenomenon of Stiffness, 156
- Stability, 154
- Stiffness, 154
phasor
- notation, 22
Picard, 154
plane
- meridian , xiii
polarization, 4
Radiation
- single wire, 2
Replacement Rules, 58
Rotation Matrix, 60
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
- E-wall (electric), 35
- H-wall (magnetic), 35
