
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00004682/00001
Material Information
 Title:
 Multiscale discretization of electricfield equations
 Creator:
 Huang, ShuJen
 Publication Date:
 2005
 Language:
 English
 Physical Description:
 xii, 83 leaves : ill. ; 29 cm.
Subjects
 Subjects / Keywords:
 Altitude ( jstor )
Approximation ( jstor ) Charge density ( jstor ) Clouds ( jstor ) Conductivity ( jstor ) Electric fields ( jstor ) Lightning ( jstor ) Modeling ( jstor ) Three dimensional modeling ( jstor ) Thunderstorms ( jstor ) Dissertations, Academic  Mathematics  UF ( lcsh ) Mathematics thesis, Ph. D ( lcsh )
 Genre:
 bibliography ( marcgt )
theses ( marcgt ) nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph. D.)University of Florida, 2005.
 Bibliography:
 Includes bibliographical references.
 General Note:
 Printout.
 General Note:
 Vita.
 Statement of Responsibility:
 by ShuJen Huang.
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
 Resource Identifier:
 003320039 ( ALEPH )
AA00004682_00001 ( sobekcm ) 770716514 ( OCLC )

Downloads 
This item has the following downloads:

Full Text 
MULTISCALE DISCRETIZATION OF
ELECTRICFIELD EQUATIONS
By
SHUJEN HUANG
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
2005
Copyright 2005
by
ShuJen Huang
To DengShan,
for now and forever.
ACKNOWLEDGMENTS
First of all, I would like to express my gratitude to my supervisory committee
chair, Professor William W. Hager. Without his assistance, this dissertation could
not have been completed. In addition, I appreciate his encouragement throughout
this period of the research. His confidence in me always strengthens my beliefs in
my ability to complete this research.
Second, I would also like to thank Dr. Jay Gopalakrishnan, Dr. Shari Moskow,
Dr. Sergei S. Pilyugin, and Dr. Vladimir A. Rakov for serving on my supervisory
committee. Their valuable suggestions have been very helpful to my research.
Third, thanks go to my officemates (Dr. Soonchul Park, Beyza, Hongchao, and
Sukanya), and all colleagues and friends in the Department of Mathematics at the
University of Florida. Their company alleviated the stress and frustration of this
time.
Last, but not least, I wish to express my special thanks to my family: to my
husband, DengShan, for his love, patience, and encouragement; to our daughter,
Marie, for being a glorious joy to us; to my mother for her support and assistance
and for helping me to take care of my daughter throughout the dissertation
writing process; to my parentsinlaw for their wholehearted understanding and
encouragement during my studies; and to my sister for her unstopping support
and encouragement. They have always stood by me. Without their support, this
dissertation could not have been completed successfully.
TABLE OF CONTENTS
page
ACKNOWLEDGMENTS ............................ iv
LIST OF TABLES ....... .... ............. ...... vii
LIST OF FIGURES ................................ viii
KEY TO ABBREVIATIONS AND SYMBOLS ................ x
ABSTRACT .. .. .. xiii
CHAPTER
1 INTRODUCTION ... .... ......... ..... .... 1
2 LITERATURE REVIEW .......................... 4
2.1 Dynamic Cloud Models ........ 4
2.2 Approaches for Lowering Charge due to Lightning ......... 6
2.2.1 Lightning Model with Bulk Flashes ............ 6
2.2.2 Lightning Model with Explicit Lightning Channels ..... 7
2.2.3 Hager's M odel ...... ........... ...... 10
3 THE M ODEL ......... ................ ........ 12
3.1 Formulation of the Equations ................... 12
3.2 Discretized Equations ........................ 14
3.3 Important Role and Approximation for Conductivity Function 22
3.4 Discussion: d = (1 e())/).Y ...... ....... 26
4 STABILITY AND ERROR ANALYSIS .................. 28
4.1 Stability ... .. . .... 28
4.2 Error Analysis .................. .......... 31
5 CONVERGENCE ANALYSIS ....................... 38
5.1 MaxNorm Bound for Matrix C1 = (A B) .. 38
5.2 Estimates for A 1 ........................... 40
5.3 MaxNorm Bound for Matrix C1D ................ 46
5.3.1 MaxNorm for Matrix C1D in 1Dimensional Space 47
5.3.2 MaxNorm for Matrix C1D in 2 and 3Dimensional Spaces 52
5.3.3 Numerical Results for Spectral Radius of C1D ...... 57
6 NUMERICAL EXPERIMENTS .......
6.1 Experiment 1 .............
6.2 Experiment 2 .............
6.3 Experiment 3 .............
6.4 Experiment 4 ..............
7 CONCLUSION ................
APPENDIX ....................
A APPROXIMATION FOR CONDUCTIVITY
A.1 Approximation from Figure 9.6 [27] .
A.2 Approximation from Figure 33 ..
B FISHPAK ...................
REFERENCES ...................
BIOGRAPHICAL SKETCH ............
FUNCTION
LIST OF TABLES
Table page
31 Coefficients of the polynomial p(z) .... 25
51 Spectral radius for C1D in one dimensional space ... 57
52 Spectral radius for C1D in two dimensional space ... 58
53 Spectral radius for C1D in three dimensional space ... 58
61 Breakdown time comparison for different values of w ... 62
62 Relative error in the potential at t = 65 seconds using the maximum
norm for CrankNicholson scheme (pA = 1/2). ...... 67
63 Relative error in the potential at t = 65 seconds using the maximum
norm for Backward Euler scheme (pA = 1). ...... 67
64 Estimates for the constant C generated by Equation (6.7) ... 71
65 Estimates for constants p and C generated by Equation (6.10) 73
A1 Four points chosen from Figure 9.6 [27] .... 76
A2 Thirteen points chosen from Figure 33 . 77
LIST OF FIGURES
Figure page
31 Path PXPj connects the node Pi to the node Pj ... 19
32 Directed graph, G(M), for a tridiagonal matrix M with nonzero en
tries on its diagonal, superdiagonal, and subdiagonal. 21
33 Electrical conductivity a and corresponding relaxation time T = Ea1,
where E = 8.8512 F m1, versus altitude under a variety of geo
physical conditions. LL, low latitude, "wavy"; MLPS, midlatitude
typical night (highlatitude, quite); AZTDN, auroral zone, typical
disturbed night; MLD, midlatitude day, quite; MHL, midhigh
latitude, typical of ~ 100 measurements; REP, relativistic electron
(energy from a few MeV to 10 MeV) precipitation event (unusual);
PCA, polar cap absorption event (an unusually large flux of ener
getic, ~ 100 MeV) solar protons within the polar cap). Adapted
from Hale [8]. (Source: reprinted with permission from Rakov and
Uman [21], published in 2003, copyright Vladimir A Rakov and
Martin A Uman, published by Cambridge University Press) 24
sin2
41 Graph of f(x) = sn+ )2, when a = 1.05. ........... 30
(a + cos X)2
sin x
42 Graph of f(x) = cosx)2, when a = 1. .............. 31
(a + cos X)2
61 Cylindrical shell of the volume integral f dV 60
62 Potentials (0(0, 0, z), 0 < z < 100 km) in two domains, Q1 and Q2 64
63 Potentials (0(0, 0, z), 4.6 km < z < 5.4 km) near the negative charge
(5 km) in two domains, Q1 and 2 . 64
64 Potentials (0(0, 0, z), 9.6 km < z < 10.4 km) near the positive charge
(10 km) in two domains, Q2 and 2 . 65
65 Potentials (O(x, 0, 5) and O(x, 0, 10), 25 km < x < 25 km) along the
xaxis at the altitude 5 km and 10 km .. ..... 65
66 Maximum norm of potential versus time, where At = 1/8 s and p = 0 68
A1 Electric conductivity as a function of altitude approximated from curve
B in Figure 9.6 [27] . 76
A2 Electric conductivity as a function of altitude approximated from Fig
ure 33 . . 78
KEY TO ABBREVIATIONS AND SYMBOLS
The list shown below gives a brief description of the major mathematical
symbols defined in this dissertation.
E electric field
J current density
a conductivity
E permittivity
potential
Uz
A = e0(At)
y
1 e/(At)
d =  7
7
= ( e 'Y(A t)) O_._
or
S= {(x,y,z) :0 < z < D, x\ < D\, yl\ < Dy},
where Dx, Dy, and D, are some prescribed values.
O, = {(x, y, z): Ix < 50 km, y < 50 km, 0 < z < 100 km}
Q2 = {(x, y, z) : x < 25 km, y < 25 km, 0 < z < 100 km}
2
A1 (A1
/ 0
B d2
\
1
2
di
0
d3
A2
Am / mxm
A, + I I
A2 = I A1+I I
MT2 X M2
(B,
B2 = Bi
A2 = Ai
A2+I
A = hI
1 B2
(A2
A = A2
m2XM/ 2x
m2 xm2
I
A2+I
i m3xm3
2)
m3Xm3
m3Xm3
~ )mxm
d3
/ mxm
h
C = A B
2
D = AA + (I p)B
2
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
MULTISCALE DISCRETIZATION OF
ELECTRICFIELD EQUATIONS
By
ShuJen Huang
May 2005
Chair: William W. Hager
Major Department: Mathematics
We develop stable finitedifference approximations for equations describing the
electric field between the surface of the Earth and the ionosphere. In this scheme,
we introduce a parameter p in the time discretization. For example, 1 = 0 and
pI = 1 give us forward and backward difference methods respectively, and p = 1/2 is
the CrankNicholson scheme. The approximations are unconditionally stable when
1/2 < p < 1 and unstable when 0 < p < 1/2. Maxnorm error estimates for the
discrete approximations are established. We show that the error is O(At)2 + 0(h2)
if p = 1/2, and O(At) + 0(h2) otherwise. We conduct numerical experiments to
verify our analysis.
CHAPTER 1
INTRODUCTION
Most people are afraid of lightning because it can destroy buildings and
numerous systems critical to daily life. Moreover, it can kill people. Loss of US
property by lightning damage is a huge expense each year.
The physical processes involved in lightning are the focus of intensive research.
As the particles, such as graupel and ice, within a cloud grow in size and increase
in number, some become charged through collisions. Smaller particles tend to
acquire positive charge, while larger particles acquire more negative charge. These
particles tend to separate under the influence of updrafts and gravity, until the
upper portion of the cloud obtains a net positive charge and the lower portion of
the cloud becomes negatively charged. This separation of charge produces huge
electrical potential difference within the cloud as well as between the cloud and the
ground. As a result of these electric fields, a flash occurs, moving charge between
positive and negative regions of a thunderstorm.
Uman [26] gave a detailed history of early lightning research. In the mid
eighteenth century, Benjamin Franklin performed the first scientific study of
lightning. He designed an experiment that conclusively proved the electric nature
of lightning. Little significant progress was made in understanding the properties
of lightning until the late nineteenth century when photographic and spectroscopic
tools became available for lightning research. Lightning current measurements
were first proposed by Pockels [18, 19, 20]. He analyzed the magnetic field induced
by lightning currents to estimate the amount of current. Wilson [31, 32] was the
first researcher to use the electric field measurements to estimate the structure of
thunderstorm charges involved in lightning discharges. He won the Nobel Prize
for inventing the Cloud Chamber and made a major contribution to our present
understanding of lightning. Lightning research continued at a steady pace until the
late 1960s when the research became particularly active. This increased interest
was motivated by
The danger of lightning to aircraft or spacecraft,
Solidstate electronics used in computers and other devices,
Improved measurement and observational capabilities due to advancing
technology.
Most lightning research is done by physicists, chemists, meteorologists, and
electrical engineers. Hager [5, 6, 7] was the first mathematician using Maxwell's
equations to develop a threedimensional mathematical electrical model to sim
ulate a lightning discharge. His discharge model [6] was obtained by discretizing
Maxwell's equations to obtain a relation between the potential field and current
density due to the motion of charged particles moving under the influence of wind.
Spatial derivatives in his equation were approximated by using volume elements in
space, while the temporal derivatives were estimated by a backward Euler scheme
in time. Since conductivity is very large in the region where the electric field
reaches the breakdown threshold, he evaluated the solution limit as the conductiv
ity tends to infinity in the breakdown region. In his model [6], the output was the
electric field as a function of time, and the inputs were currents generated by the
flow of charged particles within the thundercloud under the influence of wind.
This dissertation is based on Hager's mathematical model. Some improve
ments are made compared to Hager's earlier work. For example,
The equations are discretized in a different way allowing us to solve them
by fast fourier transform (FFT) rather than by Cholesky factorization.
Consequently, finer meshes can be employed and larger domain can be
modeled. For example, numerical results reported by Hager et al. [7], the top
3
of the model domain was 48 kilometers, we measure the top of the domain at
100 kilometers (ionosphere).
* The huge increase in conductivity with altitude is taken in account in our
scheme.
CHAPTER 2
LITERATURE REVIEW
Numerous studies in lightning from different aspects have been reported in
the past few decades. This review is focused on numerical models of lightning
discharges.
2.1 Dynamic Cloud Models
Atmospheric scientists, Helsdon and Farley [10] and Ziegler et al. [33, 34, 35],
considered the electromicrophysical processes for the charge exchange between
interacting particles within the cloud domain to obtain dynamic cloud models.
They divided water and ice hydrometeors into five categories: cloud water, cloud
ice, rain, snow, and highdensity precipitating ice (graupel/hail). Each of the five
categories of hydrometeor was allowed to have a charge associated with it. This
charge was then transported with the species of hydrometeor, according to the
transport equation
Q = V VQ + (UtQpa) + V KmVQ 6Q (2.1)
&t pa 9z 9t )inter
where Q was the charge density (in coulombs per cubic meter) carried by hydrom
eteor category (negative or positive), V was the velocity vector, Pa was the air
density, Ut was the massweighted mean terminal fall speed of the hydrometeor
class, Km was the nonlinear eddy coefficient, and (6Q/6t)inter was the charge ex
changed between categories of hydrometeor due to interactions. The cloud model of
Helsdon and Farley [10] was twodimensional; while the models of Ziegler, Straka,
MacGorman, Dye, and Ray [33, 34, 35] were established in threedimensional space.
The first term on the righthand side of Equation (2.1) was the advection term.
The second term was the fallout term, and it was equal to zero for cloudwater
and cloudice charge. The third term was the eddymixing term, and the last term
represented the charge exchanged between any two classes of hydrometeor during
an interaction between those two classes.
Helsdon and Farley's model [10] also considered the presence of small ions.
The growth and decay equation governing the number concentration of small ions
was
= V (nl,2V ni,2/ 1,2E KmVnl,2) + G anin2 + SRC SINK,
where n1,2 was the concentration of small ions with subscript 1 representing
positive ions and subscript 2 representing negative ions, p was the smallion
mobility (in meters squared per volt per second, pressure dependent), E was the
twodimensional electric field vector, G was the iongeneration rate due to cosmic
ray ionization (in ions per cubic meter per second, dependent on height), a was
the ionion recombination coefficient (1.6 x 1012 m3/s), SRC was the source term
for small ions from processes other than cosmicray ionization, and SINK was the
sink term for small ions from processes other than recombination. Ziegler et al. [34]
excluded the contribution from free ions in their model, because they thought this
contribution was very small within the cloud.
With the various types of charges accounted for in the model, at each grid
point a total charge density can be expressed as
p = e(n n2) + Q,
where p was the total charge density (in coulombs per cubic meter) and e was the
electronic charge (1.6 x 1019 C). The first term represented the net charge due to
the difference in smallion concentrations, and it was equal to zero in the model of
Ziegler et al. [34] because they ignored the contribution from free ions. The second
term represented the net charge on the five categories of hydrometeor.
With the total charge determined, they applied Poisson's equation
to obtain the scalar electric potential in volts, where e was the electrical permit
tivity of air (8.8592 x 1012 C2 / ( N m2)). Then the electricfield vector can be
solved from the potential q using Gauss's law
E = Vr.
2.2 Approaches for Lowering Charge due to Lightning
There are different approaches for lowering charge due to lightning. The neu
tralization of charge by lightning in the models with bulk flashes and with explicit
lightning channels are discussed in Section 2.2.1 and Section 2.2.2, respectively.
The approach used in the model of Hager et al. is discussed in Section 2.2.3.
2.2.1 Lightning Model with Bulk Flashes
Rawlins [22] was the first scientist to add a simple lightning parameterization
to a numerical cloud model. When the electric field reached the discharge threshold
of 500 kV/m everywhere within the model domain, he simply reduced the net
charge by a constant fraction everywhere (arbitrarily to 30% of the prestrike
values). This was the simplest representation of the neutralization of charge by
lightning.
On the other hand, Takahashi [25] suggested a more complex parameterization
of a lightning flash for his axisymmetric, twodimensional model. His parameteri
zation neutralized 20 Coulombs of charge by removing that amount in each of the
two oppositely charged regions whenever the maximum electric field had the value
of 340 kV/m. To neutralize the charge, first he found the maximum positive charge
density and the maximum negative charge density. Next, he determined the small
est region about each of the two maxima that contained 20 Coulombs. Finally, he
subtracted positive charge equally from all charged particles in the positive region;
and similarly, subtracted negative charge equally from all charged particles in the
negative region.
Ziegler and MacGorman [33] also applied a simple parameterization, but used
a threedimensional model to treat the bulk effect of lightning on storm charge in
a time step. They neutralized a fraction of the net charge at all grid points where
the magnitude of the net charge density exceeded 0.5 nC/m3 to imitate the effect of
lightning releasing the charge from its channels. The above bulkflash schemes are
relatively easy to implement but lack realism.
2.2.2 Lightning Model with Explicit Lightning Channels
Most lightning discharge models with explicit channel treat only a one
dimensional or unbranched channel. The branching of lightning channels would
be very difficult to treat analytically because of the lack of symmetry and the
incomplete understanding of the processes.
2.2.2.1 The Model of Helsdon et al.
Helsdon et al. [9, 10, 11, 12] estimated both the geometry and charge distri
bution of an intercloud lightning flash in a twodimensional Storm Electrification
Model (SEM). Since then, their parameterization has been extended to a three
dimensional numerical cloud model. They applied the concept of bidirectional
leader [13]: The parameterized lightning propagated bidirectionally (initially
parallel and antiparallel to the electric field) from the point of initial breakdown.
Their parameterization produced a single, unbranched channel that traced
the electricfield line from an initial point. Initiation, propagation direction, and
termination of the discharge were computed using the magnitude and direction
of the electric field vector as the determining criteria. The charge redistribution
associated with lightning was approximated by assuming that the channel remained
electrically neutral over its total length.
An initial criterion: a threshold of electric field was chosen with a value of
400 kV/m. The channel was extended in both directions along the field line until
the ambient electricfield magnitude fell below a certain threshold (150 kV/m)
at the locations of the channeltermination points. They assumed that the linear
charge density at a grid point, P, along the channel was proportional to the
difference between the potential at the point where the discharge originated, and
the potential at P. The linear charge density can be presented by
Qp = k(4p 4o),
where Qp was the charge density at P, and 4p and o0 were the potentials
at P and the initiation point of the discharge, respectively. The value of this
proportionality constant k controlled the amount of charge transferred by the
discharge. Charge distribution at each end of the channel was adjusted in order
to maintain net charge neutrality of the flash, consistent with the theory proposed
by Kasemir [13]. They extended the channel by four grid points at each end,
and adjusted the distribution of linear charge density on the extensions, to force
the integral of lightning charge density for the whole channel to be zero. In this
extended region, they assumed that the charge density decreased like eax2, where
x is the distance from the channel.
2.2.2.2 MacGorman's Model
MacGorman et al. [15] suggested a lightning parameterization that was
considered an extension of the Helsdon et al. [12] parameterization in conjunction
with some of the bulklightning parameterization methods presented by Ziegler and
MacGorman [33]. MacGorman et al. [15] developed a parameterization to enable
cloud models to simulate the location and structure of individual lightning flashes
by using the conceptual model of MacGorman et al. [14] and Williams et al. [30].
Their parameterization [15] of lightning flash structure proceeded in two stages.
First, Helsdon et al. [12] provided the initial development. From initiation at a grid
point, a flash traced the electricfield line outward in both parallel and antiparallel
directions, until the magnitude of the ambient electric field at each end fell below
some certain threshold value. Secondly, if one end of the channel reached ground,
the parameterization terminated at that end, but allowed the other end to continue
developing.
Charge estimation and neutralization were parameterized by applying the
technique proposed by Ziegler and MacGorman [33], except that Ziegler and
MacGorman neutralized charge at all grid points having p(i, j, k) > Pi (where
p(i, j, k) was the net charge density at the grid point (i, j, k) and pi was the
minimum Ip(i, j, k) for all grid points to be involved in lightning beyond initial
propagation) throughout the storm, but their parameterization neutralized charge
only at such grid point within a single localized flash.
2.2.2.3 Mansell's Model
Mansell et al [16] proposed a lightning parameterization derived from
the dielectric breakdown model that was developed by Niemeyer et al. [17] and
Wiesmann and Zeller [29] to simulate electric discharges. They extended the
dielectric breakdown model to a threedimensional domain to represent more
realistic electric field and thunderstorm dynamics.
In their work, the stochastic lightning model (SLM) was an application of
the WiesmannZeller model to simulate bidirectional discharges in the regions
of varying net charge density (e.g., in an electrified thunderstorm). Procedures
for simulating lightning flashes in the thunderstorm model were as follows. A
flash occurred when the magnitude of the electric field exceeded the initiation
threshold E1iit anywhere in the model domain. The lightning initiation point was
chosen randomly from among all the points where the magnitude of the electric
field is greater than 0.9Eint. Both decisions for choosing the initiation threshold
and the initiation point were suggested by MacGorman et al. [15]. Positive and
negative parts of the flash were propagated independently so that up to two new
channel segments (positive and negative) could be added at each step. Both ends
had default initial propagation thresholds of 0.75Enit. For flash neutrality, they
assumed that the channel structure would maintain overall charge neutrality
as long as neither end reached the ground, as suggested by Kasemir [13]. But,
for computational simplicity, their parameterization maintained nearneutrality
(within 5%) by a technique of adjusting the reference potential to the growth of the
lightning structure instead of adjusting the reference potential of the channel.
2.2.3 Hager's Model
Hager et al. [5, 6, 7] proposed a threedimensional lightningdischarge model
that produced bidirectional, axisymmetric IC and CG flashes. The model gen
erated the discharge region, charge transfer, and detailed charge rearrangement
associated with the flash.
Their approach to lightning was quite different from those in Section 2.2.1
and Section 2.2.2. Their breakdown model was based on Maxwell's equations.
They assumed that current due to transport of charge under the influence of wind
was known. They obtained an equation governing the evolution of the electric
potential under the assumption that the time derivative of the magnetic field can
be disregarded. After integrating this equation over boxes and approximating
derivatives by finite differences, they obtained an implicit system of difference
equations describing the evolution of the electric field. Their approach to lightning
was to let the conductivity tend to infinity wherever the electric field reached the
breakdown threshold. This approach appeals to our basic conception of nature:
When the electric field reaches breakdown threshold, conductivity becomes very
large as a plasma forms.
When the electric field reaches the breakdown threshold, the electric potential
changes instantaneously everywhere within the thundercloud. The Inverse Matrix
Modification Formula [4] was applied to evaluate this change:
after = ,before AiU(UTA1U)1UTbefore, (2.2)
where 4 before was the electric potential before discharge, after was the electric
potential after discharge, A was the discrete Laplacian, and U was a matrix with a
+1 and 1 in each column corresponding to the arcs associated with the breakdown.
There were no parameters in Equation (2.2) besides the electric potential before
discharge. This was consistent with experimental observations: The charge is
controlled predominately by a single parameter: the local electrostatic field. This
was observed in experiments reported by Williams et al. [30].
CHAPTER 3
THE MODEL
Hager et al. [6] were the first to use Maxwell's equations to develop a three
dimensional electric model for a thunderstorm. We first establish the model based
on Maxwell's equations.
3.1 Formulation of the Equations
According to Hager et al. [6], the curl of the magnetic field strength H is given
based on Maxwell's equations:
9E
Vx H = e + E + J, (3.1)
where e is the permittivity, a is the conductivity, E is the electric field, and J is the
current density associated with the influence of the wind. In this model, we have
the following assumptions:
Assumption 1: The time derivative of the magnetic field is zero, i.e.,
B B_
=0 VxE 0
> E=V0
Therefore E is the gradient of a potential q.
Assumption 2: Let EB be the breakdown field strength. Then the electric
field magnitude is always less than or equal to the breakdown threshold EB
That is, El < EB.
Assumption 3: When the electric field reaches the breakdown threshold EB at
some point, the conductivity tends to infinity in a small neighborhood at that
point. That is, if IE(x, y, z) = EB, then a(x, y, z) = +oc.
In this dissertation, we do not implement the lightning discharge described in
Section 2.2.3. Instead we focus on integrating the equations up to the discharge
initiation.
We take the divergence of Equation (3.1) to get
0E
eV* + VaE + V J = 0. (3.2)
Boundary conditions: In this model, we consider a bounded domain Q between
the surface of the earth and ionosphere. Let us consider boundary conditions on
9fO:
= 0 on the surface of the earth, since we treat the earth as a perfect
conductor.
= 1 on the ionosphere. 41 = 250,000 V (Adlerman and Williams [1]).
On the sides of the domain Q, assume that the potential 0 is the "fair field"
potential. That is, the current and the time derivative of the electric field are
both zero under fairweather conditions. Hence, the potential satisfies
V aE = 0. (3.3)
Assume that both the potential and the conductivity a depend only on z
on the sides of the domain Q, and since E = V0, we have
a or =0, (3.4)
from Equation (3.3). In order to obtain the potential function on the sides of
the domain Q, the Equation (3.4) is solved with these boundary conditions:
the potential is zero on the surface of Earth and 41 on the ionosphere.
Therefore, we discover that
(z) = dz,
O /z) (Z)
where
01
f dz
and h is the height of the domain Q (h = 100 km for numerical experiments
in Chapter 6). If 01 = 250000 (see [1]), c = 2.852108 x 109.
Finally, using the transformation O(x, y, z) = j(z) + O(x, y, z), we have the
boundary condition: 4 = 0 on 0Q.
Since E = Vq and 0 = ( + 0, Equation (3.2) becomes
e AV)+ UA4'+ua +VJ =0,
at OZ
where A denotes the Laplacian operator defined by A = V V, and az denotes the
derivative of a with respect to z.
Then our model is the partial differential equation with the boundary condi
tion as follows:
PDE: A+ 'AO + _ + =0
at e e 9z e
BC: P = 0 on O
Domain: Q = {(x, y, z) : 0 < z < Dz, Ix < Dx, yl < Dy},
where D,, Dy, and D, are prescribed values. In our computation, we take D, =
100 km. For reference,
Qi = {(x,y, z) : lx< 50 km, y < 50 km, 0 < z < 100 km}. (3.5)
3.2 Discretized Equations
We demonstrate the discretization process for a uniform mesh in a rectangular
coordinate system. Consider the differential equation
A9 aO a + (3.6)
At e E 95z e E
The solution of Equation (3.6) is
An+ = e(At)An e(St+) 7z + ) ds, (3.7)
where superscript n denotes the time level tn, 7 = c/e, and 7 = oz/6 is the
derivative of a/e with respect to z. First, we approximate the integral part of
Equation (3.7), denoted as I:
I = I (st+l) 7Z + ds.
The following steps are conducted.
Step 1: Both ft/Oz and V J are evaluated at time 5(tn + tn+1), which is
denoted as time level n + I. Hence we have the following approximation
1 eY(At) / \n+i (V. J)n+)
e ( ( ) ) + ( /J)
Define
d 1 ey(At)7
y
Then we can rewrite the above approximation as follows :
(/n+ d eY(At) (V J)n+
I Oz 7 e
Step 2: Use (1 p)(VI/Oz)" + t(I09/Oz)"+11 to approximate the term
(ap/a8z)"n+. Hence,
[ (r ) (Y +1] 1 e(At) (V. J)+
I d (1 /)(z 1+/" ) J + 
1O ( 9Z1z 7
where 0 < p < 1.
So Equation (3.7) is then approximated by the following equation
Aon+1 =
AAO d (1 () )n+
1 e^() (V J)(+.
, (3.8)
where A = e"^(), 7 = a/e, and 0 < u < 1.
Now we consider Equation (3.6) but ignore the force term J.
_0 = 7A 7 (3.9)
Similarly, the solution of Equation (3.9) is
an+1 = (At)Aon e_(st+1)y^z ds, (3.10)
and it can be approximated by
Ao++=AAA d (1A) z + /Z (3.11)
Next, the spatial part of Equation (3.11) can be discretized by the finite difference
method. Define
hX = 2D hy D= and h, =
m+1 m+1 m+1
where m is a positive integer. The mesh point (xi, yj, Zk) is defined as
xi = Dx + ihx,
yj = DA, + jhy,
Zk = khz,
where 0 < i < m + 1, 0 < j < m + 1, and 0 < k < m + 1. Let Oi,j,k denote
the discrete approximation to o at the point (xi, yj, zk). For domain Q1: using the
uniform mesh size for x, y, and z (that is, Ax = Ay = Az = h), we have the
approximation
SWi,j,k = i,j,k+l Oij,k1
2h
for the partial derivative O//oz at the point (xi, yj, Zk). The second partial
derivative 0291/0x2 at the point (xi, yj, zk) has an approximation denoted by
OijPi,k:
2,jk = _i +l,j,k 2~ ,j,k + Vil,j,k
8 Wi~j,k 2  ^
Similarly, we use approximations 19ij,k and 930ij,k to the second partial deriva
tives 92?/Oy2 and 02/0z2 :
a2q, Oij+1,k 21)0ijk + Oijlk
2' Oi,j,k h 2
and
i,,k i,j,k+1 20i,j,k + i,j,kl
)3,,k h2
respectively. We define the discrete Laplacian operator Ah as follows:
3
A ,j,k (21 ijk
l=1
After replacing A by Ah and O9i/Oz by the discrete azo4 in Equation (3.11), we
obtain the following finite difference scheme.
Ahvzn k As j., =zh dk n 9,,p (3.12)
A ij, + AkA ,,3k = (1 P) h29z,3,k + P ? ] (3,12)
Equation (3.12) yields a matrixvector form
A n+1 AA n" = h [(1 p)B "n + tB"'++] (3.13)
where 4 denotes the vector with components i,j,k. The matrix A is symmetric
and positive definite and the structure of the matrix A is a block tridiagonal
matrix in the following form:
A2 + Im2 Im2
Im2 A2 + Im2 Im2
1
A = h (3.14)
Im2 A2 + Im2 Im2
Im2 A2 + In,2
where Im2 is the m2 x m2 identity matrix and A2 is an m2 x m2 block tridiagonal
matrix in the form
Ai + Im Im
Im A, + Im Im
where Im is an m x m identity matrix and A1 is a tridiagonal matrix with 2's on its
diagonal and 1's on its superdiagonal and subdiagonal. The matrix A is a diagonal
matrix with m2 blocks on the diagonal, and each block has Ai on its diagonal. The
matrix B is a diagonal block matrix in the following form:
B
B = (3.15)
where the diagonal block B1 is an m x m tridiagonal matrix in the form
0 di
d2 0 d2
d3 0 d3
Varga [28] gives the definition for an irreducible matrix and some theorems to
determine the irreducibility of a matrix. Definition 3.2.1, Definition 3.2.2, Theorem
3.2.3, Definition 3.2.4, and Theorem 3.2.5 are adopted from Varga [28].
Definition 3.2.1 For n > 2, an n x n matrix M is reducible if there exists an n x n
permutation matrix P such that
PMPT = (Mi M1,2
0 M2,2
where Mi,1 is an r x r submatrix and M2,2 is an (n r) x (n r) submatrix, where
1 < r < n. If no such permutation matrix exists, then M is irreducible.
A directed graph is an alternative way to check irreducibility for a matrix. Let
M be any n x n matrix, a finite directed graph G(M) can be constructed in the
following way: Consider any n distinct points (nodes) P1, P2, , Pn in the plane.
For every nonzero entry mij of the matrix, we connect a path directed from the
node Pi to the node Pj, as shown in Figure 31.
P P
Figure 31: Path PiPj connects the node Pi to the node Pj
20
Definition 3.2.2 A directed graph is strongly connected if, for any ordered pair of
nodes Pi and Pj, there exists a directed path
PiPl, PIIP2 12 Pr1 Pr=j
connecting Pi to Pj.
Theorem 3.2.3 An n x n matrix M is irreducible if and only if its directed graph
G(M) is strongly connected.
Definition 3.2.4 An n x n matrix M = (mij) is diagonally dominant if
n
Imi,i E Imij (3.16)
j=iljji
for all 1 < i < n. An n x n matrix M is irreducibly diagonally dominant zf M is
irreducible and diagonally dominant, with strict inequality in Inequality (3.16) for at
least one i.
Theorem 3.2.5 If M is an n x n irreducibly diagonally dominant matrix, then the
matrix M is nonsingular.
Now, we rearrange Equation (3.13) to yield
(A B) "+1= (AA + h(1 M)B n. (3.17)
Let C = A hMB and D = AA + h(1 ,.)B. Then we have
Cxn+l = DW",
where the matrix C is a block tridiagonal matrix given by
/ Cm m I(m
Im 2 Cm2 Im2
C = (3.18)
h2
where Cm2 is an m2 x m2 tridiagonal block matrix with m x m matrices
/ 6 1 dih
1 + d2 6 1 hd2
(3.19)
on the diagonal block, and with negative m x m identity matrices on both its
superdiagonal and subdiagonal blocks.
Lemma 3.2.6 Let M be any nx n tridiagonal matrix. If all entries on the diagonal,
superdiagonal, and subdiagonal of M are nonzero, then M is irreducible.
1 2 3 n1 n
Figure 32: Directed graph, G(M), for a tridiagonal matrix M with nonzero entries
on its diagonal, superdiagonal, and subdiagonal.
Proof: Since the tridiagonal matrix M has nonzero entries on its diagonal,
superdiagonal, and subdiagonal, the directed graph of the matrix M, G(M), is
shown in Figure 32. From Figure 32, it is easy to see that G(M) is strongly
connected, and therefore M is irreducible by Theorem 3.2.3. *
Lemma 3.2.7 The matrix C in (3.18) is irreducible.
Proof: First, we want to show that each diagonal block Cm2 of the matrix C is
irreducible. By Lemma 3.2.6, each diagonal block (3.19) of Cm2 is irreducible;
hence, there exists a path for any two vertices associated with a diagonal block
(3.19). Suppose we are given vertices i < j in two distinct diagonal blocks; the
subdiagonal blocks of Cm2 provides a path between i and i + m and between i + m
and i + 2m etc. Eventually, we reach a vertex k in the block containing j. And
since the vertices k and j are in the same irreducible diagonal block, there exists a
path from k to j. Hence, there exists a path from i to j.
Now the matrix C has irreducible diagonal blocks Cm2. The structure of C
is similar as the structure of Cm2. Similarly, we can show that the matrix C is
irreducible. n
Therefore C = A hpB is irreducibly diagonally dominant by Definition 3.2.4
and Lemma 3.2.7. Based on Theorem 3.2.5, C is nonsingular. Therefore, we have
the iteration
r+1 = C'D ",
where C = A /B, D = AA + (1 p)B, and C1D is the amplification matrix.
3.3 Important Role and Approximation for Conductivity Function
The conductivity function a varies by many orders of magnitude between
the surface of the Earth and the ionosphere. Figure ?? and Figure 33 show
that a grows more than the factor 1010 S/m from the surface of the Earth to the
ionosphere, which leads the differential Equation (3.6) to be a very stiff equation.
For a stiff equation
dy
d ky, (3.20)
where k is a very large positive number, its solution is y(t) = ekty(O) and y(t)
tends to zero quickly. Let us consider the following three schemes:
1. Euler explicit scheme for Equation (3.20):
yn+1 y TL
=^ ky".
At
Therefore, yn = (1 kAt)"y(0). When k is very large, (1 kAt) < 1 and
y" + +oo as n increases. Hence Euler explicit scheme does not work.
2. Euler implicit scheme for Equation (3.20)
yn+l __ n1
y+ kyn+l.
At
Therefore, y" = ) y(0) and 0 < < 1 when k is large. Hence
yn > 0 as n + 00. The scheme works.
3. The exact formula:
n+1 = kAt n
y =e y.
Therefore, yn+1 = (ekIt)n y(O). The factor (ekAt)" tends to zero much
faster than )when k is large. This scheme works better than Euler
implicit scheme.
Therefore, the above explanation demonstrates that the exact formula is the best
method to solve a stiff equation. The conductivity function, a, plays an important
role in the Equation (3.6). That is the reason that the exact solution (3.7) is
chosen in Section 3.2.
At the beginning of the research, we use Figure 9.6 in [27] to find an approxi
mation for the conductivity function o (see Appendix A.1).
For 0 < z < 27.828 km:
o'(z) = 1013.40+0.1600z0.005175z2 +0.00009077z3
For z > 27.828 km:
a(z) = 1012.88+.06718z
Then we have
a, = (0.1600 0.01035z + 0.0002723z2)(ln 10) for 0 < z < 27.828
(7
aI = 0.1546 for z > 27.828
a
and
S() = (0.01035 + 0.0005446z)(lnl0) ,for 0 < z < 27.828
Oz or
O I = 0 for z > 27.828
Figure 9.6 [27] only provides conductivity data up to 27.8 km (about 100,000
ft). In order to produce more accurate results, we use Figure 1.3 in the book
120/ /
110 AZTDN Night Nominal"' Day
100 A'
90 I ,. .
MLTN  
so 
 ', r*" REP Nov '69 PCA. "ight
~ i M'HL
10 
1014 10"1 1012 1011 10so 10. 104 107 108 10" 10,4 103 10"2 10 10
Conductivity, S m1
1 I I i i t l l1 lI
10s 102 101 111 101 102 103 10 10 5 106 10' 10 104 1010 1011
Relaxation time, s
Figure 33: Electrical conductivity or and corresponding relaxation time T = Eoi,
where e = 8.8512 F m1, versus altitude under a variety of geophysical conditions.
LL, low latitude, "wavy"; MLPS, midlatitude typical night (highlatitude, quite);
AZTDN, auroral zone, typical disturbed night; MLD, midlatitude day, quite;
MHL, midhighlatitude, typical of 100 measurements; REP, relativistic electron
(energy from a few MeV to 10 MeV) precipitation event (unusual); PCA, polar cap
absorption event (an unusually large flux of energetic, ~ 100 MeV) solar protons
within the polar cap). Adapted from Hale [8]. (Source: reprinted with permission
from Rakov and Uman [21], published in 2003, copyright Vladimir A Rakov and
Martin A Uman, published by Cambridge University Press)
Lightning: physics and effects [21]. In the book [21], Rakov and Uman give a
figure adapted from Hale [8] that estimates conductivity at altitude up to 120 km.
Figure 1.3 in [21] is reproduced in Figure 33 of this thesis, with permission of the
authors. In Appendix A.2, we give the following approximation based on Figure
33:
a = I0p(z),
where p(z) is a quartic polynomial
4
p(z) = azi, (3.21)
i=0
where z is in kilometers and ai are coefficients of the polynomial p(z) as shown in
Table 31: Coefficients of the polynomial p(z)
i a,
0 1.316 x 101
1 6.351 x 102
2 8.682 x 104
3 8.674 x 106
4 5.781 x 108
Table 31. Then we have
a_ = (0.06351 + 0.001736z 0.00002602z2 + 0.0000002312z3) (In 10)
and
() = (0.001736 0.00005204z + 0.0000006936z2)(ln 10).
For the approximations estimated from Figure 33, we can show that and
_ (g) are bounded in the domain Q1 (defined in (3.5)) as follows:
az < 0.4793
26
and
S(
3.4 Discussion: d = yz(1 e'(t))/y
In this section, we discuss the parameter d, where
d = (1 e(A~)) 7, (3.22)
7 = l/e, and 7, = az/e. The parameter d is a critical factor when we discuss
maxnorm of some important matrices in Chapter 4.
Lemma 3.4.1 For d as defined in (3.22), we have
d = O(At),
while d < 0.4793 in Q, for any At.
Proof: For any At, we have
d (1 e(t)) 7
7
y
< 0.4793.
We also have
1 eY(At) z
d =
= 1 ey^A L s
= e7 At,
where 0 < < 1. Since e7'Y < 1, we have
d < z At = O(At)
E
Lemma 3.4.2 For d as defined in (3.22) and y7 e C2, we have
di d+l = O(At h).
Proof:
di d = (1 e(At)) ez (1 e1 +(At)) ( A
v0 i 0 i+1
+ (1 ei(At) ) (1 e 7i(At)) 
( i ri+1
= (1 e'Y(At)) )  + (eAt) eA(At) )
i [i+l + Oi+1
By the mean value theorem, we have
d+1 = (1 e (t)) Oz z)h + e(t) ( (At)) h ( +
where (E)) is the derivative of ( ) with respect to z at some zi between z, and
zi+,, =i and 7= are the value and the derivative of 7 at some = between zi and zi+1.
Since Q is bounded and a E C2, both z and A(z) are bounded, and we have
Id d,+I < (1 e'(Y ) ) ) h + )(At) (z) h
S (At) (z)ih + 5(At) (z)+,h h
= O(At h).
CHAPTER 4
STABILITY AND ERROR ANALYSIS
4.1 Stability
This section uses Fourier analysis in our study of stability. We use Fourier
analysis on the grid of integers, Z or hZ, where hZ is defined by hZ = {hm : m E
Z}. Strikwerda [23] gives the Fourier transform as follows:
Definition 4.1.1 If v is a grid function in 2 defined for all integers m, its Fourier
transform is given by
) = eimvm
Sm==oo
for ( E [7r,7r], and (7r) = (7r). If the spacing between the grid points is h, by
changing variables we have
()= m0oo
for e E [ir/h, 7r/h].
Next, we consider ) a grid function whose spacing between the grid points is
h. Then we have
1.
1 eimhO(O.+l 20m + Om1)h
V2=s e
000
S m=(oos
= (2 cos 0 2) 1 eimhmh
v m= oo
= (2 cos 0 2) m,
2.
1 oo
1 eimh (m+ Om1
M=00
= 2isin0 1 ei mh
m=oo
(2i sin 90) .
Theorem 4.1.2 The finite difference scheme, Equation (3.12), is unconditionally
stable when 1/2 < p < 1 and is unstable when 0 < p < 1/2 .
Proof: Take Fourier transform of Equation (3.12) to obtain
3
S(2 cos 0 2)n+1
j=1
3
= Ak (2 cos 0j 2)" dkh(1 p)i sin 03 n dkhpi sin 03 +1.
j=1
Sn+i 2Ak E(COS Oj 1) dkh(1 P)isin 03,n
2 E(cos O 1) + dkhpi sin 03
Define the amplification factor as
2Ak E(cos Oj 1) dkh(1 p)isin 03
2 E(cos 0j 1) + dkhpi sin 03
Then the magnitude of p is given by
ip2 4A' (E(cos O 1))2 + dh2 (1 )2 sin2 03
4 (E(cos Oj 1))2 + d2kh22 sin2 03
For 1/2 < p < 1, we have p2 < 1, since Ak = ek(t0) < 1. Hence the scheme
is unconditionally stable. For 0 < p < 1/2, we want to show \p\2 < 1 + K(At),
where K is a constant, does not hold for any Oj E [7r, r]. If pl2 < 1 + K(At), then
we have
4A2 (E(cos 0j 1))2 + dh2(1 ( )2 sin2 03
4 (E(cos 0j 1))2 + dh22 sin2 03 
4 4A ( (cos 0j 1))2 + h2(1 2 sin2 3
< 4 (E(cos0 1) )2 + 4K(At) ( cos j 1)) 2
+dkh2p2 sin2 03 + K(At)dih2p2 sin2 03
4= 4(A2 1 K(At)) ( (cos o0 1))2 < dh2(K(At)P2 + 2z 1) sin2 03
4 4(A2 1 K(At)) < d h2(K(At)L2 + 2p 1) sin 2 3 2
((cos0j 1))
4(1 + K(At)) A ) > d2h2(l
/
/
//
/
/
,/
/
 2p K(At)2j2) Sn 3 
(E(cos 0oj 1))2
sin2 x
Figure 41: Graph of f(x) = cosx)2', when a = 1.05.
(a + cos X)2
Now consider the function
sin2 03 sin2 03
(E(cos 0j 1))2 (a + cos 03)2'
where a = cos 01 + cos 02 3, and 5
sin2 x
(a + cos )2'
31
where 7r < x < 7r, and 5 < a < 1. If a # 1, the graph of f(x) is shown in
Figure 41. We can find the maximum of f(x) is a/(a2 1)3/2. When a tends
to 1, the maximum of f(x) approaches to infinity. If a = 1, the graph of f(x) is
shown in Figure 42. We can see the supremum of f(x) is infinity. Hence,
300s
250
20J
y15:
10
I \
3 2 1 0
1 2 3
x
Figure 42: Graph of
S sin2 x 1
f () = (,, ) when a = 1.
(a + cos )2
sin2 03
sup = +00.
7/h
Therefore, it is impossible to have the inequality
(4.1)
4(1 + K(At)) A2) > d h2(1 2p K(At)j?) sin23 03
(E(cos0j 1))2
In other words, the scheme is unstable when 0 K< p < 1/2. U
4.2 Error Analysis
In this section, the error generated from the discretizing process is analyzed.
Lemma 3.4.1 and Lemma 4.2.1 are applied in the error analysis. Lemma 3.4.1 gives
the estimate for d, and Hager [6] mentions the error estimate for integration as
described in Lemma 4.2.1.
Lemma 4.2.1 Let Wk,([a, b]) denote the space of functions defined on the interval
[a, b] with derivatives through order k essentially bounded. If g E W1'oo and
f E W2',,, then we have the estimate
b f (z)g (z)dz f (c) g(z)dz (  a I f'IIlg'I +  (b 4 a ) 3,
where  denotes the essential supremum and c is the midpoint of the interval
[a, b].
Proof:
jf(z)g(z)dz f (c) I g(z)dz
Ja Ja
= (f (z) f (c))g(z)dx
a
= b ( f'z(s)ds) g(z)dz
= ( [f '(C) + / f"(Tr)dr ds) g(z)dz
Ja (f J
< fb z f'(c)ds) g(z)dz + b ( s f"(r)drds) g(z)dz
SI+ 12,
where
Ib = f z f'(c)ds) (g(c) + j g'(s)ds) dz
< f (z c)f'(c)g(c)dz + (z c)f'(c) ( g'(s)ds) dz
a a c
J b
< f'Illjg'I (z c)2dz
and
12 <_ IIf"I (j Z s dTds) dz
(b a)3li11,H l.
24
This proves the lemma. n
Note that 7 is considered small in Lemma 4.2.2, Lemma 4.2.3, and Theo
rem 4.2.4.
Lemma 4.2.2 If we use the finite difference scheme (3.12) to discretize the integral
equation (3.10), then the truncation error is
O(At)m + O(h2),
where m = 3 if p = 1/2 and m = 2 if p 6 1/2,
Proof: Define R1 as follows:
where d =
I
S= eY(s+1)zds z e(st+')ds
ft. 8z 09z
tn l n7z
=I e7 tn+1) yds d Pl
et ?Stn )l ) rhA"ds
= 7z(1 e7(at))/7. Apply Lemma 4.2.1 to obtain the following estimate:
= **eT+ \ i) ] et+l
R1 = e st y)7P ds 7z f e(stl)ds
< C(At)3.
(4.2)
Define R2 as follows:
R2=d \z
I
d (1p) +P  
9z (9z
where 0 < p < 1. By Taylor's expansion, we have
+ T2 } 9
(At 2 (a
2 Qz) it
0z (tn+I) =
az
z (t.+) +
(At (a~ )
2} Oz ]
1 (At2
(t+ +
2 2 2
(9Z t
Oz (t. ) + (1F ) 2
+2 (1 )
' (f)2++/t
t)
( ,)2 (n+ ) +
This implies that if p = 1/2,
(1 p) (z
9z n
+ ( +1
z a
and if p # 1/2,
1&,n 2
19Z
+pI
(19 +
= O(At).
Therefore
R21 = O(At)m,
where m = 3 if p = 1/2, and m = 2 if po 1/2.
Let R3 be defined as follows:
AA" (Ah+ Ak Ah) k.
i ,,,k k j,k) 
80
9z (t, =
00(t.+ )
19Z 2
Hence
(z tn+ )
[ tz)
= O(At)2,
(4.3)
\z
R3 = A"n+l
/I ( )]
(1 ) z
8zp n
By Taylor's expansion, we have
JR31 = O(h2). (4.4)
Combining (4.2), (4.3), and (4.4), we see that when the solution to (3.10) is
substituted in (3.12), we obtain
RI + R2 + R3 = O(At)m + 0(h2),
where m = 3 if p = 1/2 and m = 2 if p 1/2. m
Now if we add the force term J back to Equation (3.10), then the correspond
ing finite difference scheme is
A ijk = kAhI,'j,k 2(1 2 )9zij,k + IUijk]
1 e(t) (V j)n+ (4.5)
7k E
Lemma 4.2.3 If we use the finite difference scheme (4.5) to discretize Equation
(3.7), then it has the truncation error
O(At)m + O(h2),
where m = 3 if p = 1/2 and m = 2 if p 0 1/2.
Proof: We do not include the force term J in Lemma 4.2.2. Now we only need
to focus on the error caused by estimating the J term when the finite difference
scheme (4.5) is applied to Equation (3.7). In the scheme (4.5). we approximate the
V J term in (3.7) by the V J term in (4.5). By Lemma 4.2.1, we have
e t+ ds e)f t+)ds = 0(At)3.
Combining this estimate with Lemma 4.2.2, we find that when the solution to
(3.7) is substituted in (4.5), we obtain
O(At)' + O(h2),
where m = 3 if p = 1/2, and m = 2 if p # 1/2. The finite difference scheme
(4.5) yields a matrixvector relation of the form
Ax"p+ AA 2" = h [(1 j)BqIn + Bn+l] 1 e I, (4.6)
7 e
where In corresponds to the term V J of (4.5) evaluated at time level n + 5.
The next theorem gives us the maxnorm error estimate when the solution to
(3.7) is substituted into Equation (4.6). Define the maxnorm for a vector x as
ixoo = max xil,
1
where x = (X1, x2,.. ., x)T, and define the maxnorm for a matrix M as
n
in
loo = max E mj ,
j=1
where M = (mi,j) is an n x n matrix.
Theorem 4.2.4 Let On denote the continuous solution of (3.7) evaluated at
t = nAt and at the mesh points of the grid {(xi, yj, Zk)}. Let %n be the vector with
components ?Ii,j,k in Equation (4.6). Then we have
1 On E7" 0 = O(At)m + O(h2),
where m = 2 if p = 1/2 and m = 1 if pM 1/2.
Proof: Insert en into Equation (4.6) and we get
h 1 ep(At) I
A0+i = [(1 p)BEn + MB En+] I" + Tn, (4.7)
2y 6
where Tn is the truncation error and
ITn.1oo < O(At)m + O(h2),
where m = 3 if p = 1/2, and m = 2 if p 0 1/2. Then by subtracting Equation (4.7)
from Equation (4.6), we have
,Pn+l n+1 = C1D ('n _ a) C1Tn, (4.8)
where C = A AB and D = AA + h(1 tp)B. Following from Equation (4.8), it
is found
n1
r" = M" (o0 _0) ZMIC1 n_1, (4.9)
1=0
where M = C1D. Note that we prove IMIoo = 1 + O(At) in Section 5.3 for 0 <
pi < 1 if 7 is small. Therefore, IMoo can be written as IMIoo < 1 + cAt < ecat.
Taking the maxnorm of Equation (4.9), then we have
11 e"no < e"^ct I1 6o11o + IIC1lloo 11.r1 111oo .
1=0
Note that IJC'1 oo is bounded by a constant which is independent of h (see
Lemma 5.1.2). Consequently, if I0 011oo = O(h2), then
1n 9Onl = O(At)m + O(h2),
where m= 2 ifp= 1/2, and m = 1 if p 1/2. *
CHAPTER 5
CONVERGENCE ANALYSIS
5.1 MaxNorm Bound for Matrix C1 = (A ,/B)
In this section we obtain a maxnorm bound for C1 = (A Lp1B)
where A and B are defined in Equation (3.14) and Equation (3.15), respectively.
Lemma 5.1.1 from Varga's Matrix Iterative Analysis [28] is used to show that
the matrix C is invertible and the maxnorm of the matrix C1 is bounded by a
constant which is independent of h.
Lemma 5.1.1 If M = (mj) is a real irreducibly diagonally dominant n x n matrix
with m,j <_ 0 for i 7 j, and mi,j > 0 for all 1 < i < n, then M1 > 0, i.e., all
elements of M1 are positive.
Lemma 5.1.2 If h is sufficiently small, then the matrix C is invertible. Moreover,
IIC11 o is bounded by a constant which is independent of h.
Proof: The proof of this result is the same proof used in [6] for a slightly different
matrix. Note that all the elements in the matrix C satisfy the conditions in
Lemma 5.1.1 except the elements on the subdiagonal of C. The elements on the
subdiagonal have the form 1 + pjd, where d is a bounded positive number from
Section 3.4, and 0 < p < 1. In order to use Lemma 5.1.1, we must satisfy the
condition
1 + s d < 0. (5.1)
2
Accordingly, if h is sufficiently small, C1 > 0, i.e., the elements of C1 are positive
by Lemma 5.1.1. Hence, IC'Ioo = FIC11Hoo, where 1 is the vector with every
component equals to 1. Furthermore, IICI[ _< C FG j F if F > 1.
Note that the matrix C is the coefficient matrix corresponding to the expres
sion
A ,,,k dkyOz2f,j,k, (5.2)
where Ah is the discrete Laplacian operator, z, is the discrete 9O/O9z, and dk =
(1 ek ) (A7/7)). Let
Nh = {(i,j,k) : (xi, yj, Zk) E }.
Next, we show that IC11loo < 2a2, where a = hi, and I is the maximum value
of the index i associated with points in the set Nh. Let us define O,j,k = (hi)2.
Observe that
S[h2(i + 1)2 2h2i2 + h (i 1)2] = 2 > 1.
h2
Moreover, define wij,k = (hi)2 for (i, j, k) outside Nh. If (i, j,k) E N', we choose
wi,j,k to satisfy
Ahwi,j,k dkPQzIw,i,k = 0. (5.3)
Then from Equation (5.3), we obtain
Wi,j,k = Wi+lj,k + il,j,k + Wi,j+l,k + 2i,jl,k + 1 + hdkP Wi,j,k+l
+ 1 dk) Wij,kl] (5.4)
for (i, j, k) E Nh. If h is sufficiently small, the condition (5.1) is satisfied, that is,
1 hdkp > 0. Hence, by Equation (5.4), we have
max lWi,j,k < max Iwij,k,
N4h fONh
where ONh is the set of boundary points of Nh. Hence,
ui,jd,kl (hi)2 (h)2 = a2,
for all (i, j, k) E Nh. If i = i,j,k i,j,k+ Wi,j,k, then we have
i,j,k = 0 if (i, j, k) e ON,
Ah ij,k dk p^z i,j,k > 1 if (i, j, k) E Nh.
Let T, E, and fl denote the vectors formed from iij,k, Oi,j,k, and wi,j,k for (i, j, k) e
Nh. Take F = CF > 1 to obtain
llC_1ll _< llC1Fllo = 11,11.o
= Ii8 + llo
S Ill lloo + ll0 l,
< a2 + a2 = 2a2,
where a is a constant which is independent of h. *
Remark: As noted in the condition (5.1), the matrix C will satisfy the conditions
of Lemma 5.1.1 if h is sufficiently small. In fact, using the estimate for d from
Section 3.4, we obtain the condition of h < 4.173 kilometers when p = 1 for the
domain Q1.
5.2 Estimates for A~1
In this section, we will get some preliminary results for finding the maxnorm
of the amplification matrix C1D.
Lemma 5.2.1 Let A1 be the n x n symmetric tridiagonal matrix with 2's on the
diagonal and with 1 's on the superdiagonal and on the subdiagonal. Then the eth
row of the inverse matrix of A1 has the form
(A) [ni+1,2(ni+l), ,i(ni+l),i(ni),i(ni1), ,2i,i].
n +1
Proof: Note that we compute the (i, j)th entry of the matrix A1 'A by the
formula
(AiA1)i,,j = (Aj1)i,k(Ai)kj
k=1
= (a1),jjajj, + (a1)i,jaj,j + (a1)i,j+laj+,j,
where (a1)i,j denotes the (i, j)th entry of the matrix A1'. For j = i, we have
(Ai Al),,i = [(i 1)(n i + 1)(1) + i(n i + 1)(2) + i(n i)(1)] = 1.
For j 7 i, we consider two cases:
1
(A11A)i,j = 1 [(j1)(ni+1)(1)+j(ni+1)(2)+(j+1)(ni+1)(1)] = 0,
ifj= 1,2,.. ,i1, and
1
(Ail1A)i', =  [i(ni(m1))(1)+i(nim)(2)+i(ni(m+1))(1)] = 0,
if j = i + m, where m = 1, 2, ,n i. Therefore, A'A = I.
Remark: Dr. J. Gopalakrishnan suggests to use the finite element method to
find Ai1 Consider the twopoint value problem:
U"= f on (0, 1),
U() = U(1) =0.
Then the Green's function G(x, s) is given by
G(x, (1 s)x if x < s,
(1 x)s if x > s.
Let 0 = Xo < xi < . < Xn+i = 1 with xi+1 xa = 1/(n + 1) for all i = 0,1,...,n.
Define
a(u, v) = u'(x)v'(x)dx,
42
and S = {v : v(x) is continuous on [0,1], linear on each interval [xi, xi+,] for all
i = 0,1,... n, and v(0) = v(1) = 0}. The discrete Green's function Gs(x) E S is
defined by
a(G', v) = v(xj) for all v E S,
(5.5)
for any grid vertex xj. Then
G'(x) := G(x, xj) =
(1 Xj)x if x
(1 x)xj if x>Xj.
G(xi, x) i(n + 1 j) if <
(n+ 1)2 i
G(xi, xj) = j(n+ 1i) i>j
(n+ 1)2 '
Let {}.on+1 be a basis of S, where
x) = {
1 if x = xa,
0 ifx=xj (j i)
Then the discrete Green function can be written as
G's(x) = (1 x,) k+
k=0O
= (1 j) kxk +
k=1
n+1
k=j+1
n )
k=j+l
since xo = 0 and xn+l = 1. Thus
OkXk ) Xj
n+1
k=j+l
k nE+1
\A=j+1
3
(1 L X)SkXk+ X
( k=
Ok~ k) ) I, )
k=1
j n
= L[a(Ok,0)(1 xj)Xk] + [a(k, i)Xj ( Xk)]
k=l k=j+l
= (a(qi, 01), a(i, 2), .., a( On))
(1 Xj)Xi1
(1 Xj)z2
(1 Xj)Xk
:j(1 .k Ik l)
Xj(1 Xk+2)
\ Xj(1 xn)
(5.6)
By Equation (5.5), we know
(where = 1 if i ,) == j, and(j) = = 0 otherwise. Note thatj
where 5,j = 1 if i = j, and ,,j = 0 otherwise. Note that
2(n+ 1) ifj=i
a(i, Oj) = (n + 1) ifj=ilorj=i+l
0 otherwise
Let G be the matrix with the (i,j)th entry G(xi, xj). Equation (5.6) can be written
as
(n + 1)A(i,). G(.,j) = 6ij,
where Ai(i, ) is the ith row of the matrix A1 and G(,j) is the jth column of the
matrix G. Then, we have
A1= (n + 1)G.
a(Gj, )
= a(
44
Lemma 5.2.2 Let A1 be the matrix described in Lemma 5.2.1. Let T be the n x n
tridiagonal matrix with O's on the diagonal, with 1's on the superdiagonal, and with
1's on the subdiagonal. Then IIAI'TIoo = n 1.
Proof: Since
(AIIT)i,j = E(Aj1)i,(W)kj
k=1
= (a1)i,jtjl,j + (a1)ij+ltj+l,j
= (a)il (a )i,j+1
we have
(AiT)i,j <
Then, the infinity norm for
IIAr1'TIK
2(n+li) ifj=l,2,.i 1
n + 2i 1i)
if 3 = i
2, if j =i+1,i + 2, ,n1
n+1
2i
2i if j = n
n+1'
the matrix Ai'T is
n
= maxE (A1'T)i,3
j=l
[2(n i + ~1) n+2i1 2i ]
= max (i 1)+ + (n i)
i n+1 n+1 n+1 J
S 4i2 + 4i(n + 1) 2(n+ 1) 12i (n + 1) (5.7)
Smax +  (5.7
Now focus on the term inside brackets of (5.7):
4i2 + 4i(n + 1) 2(n + 1) 12i (n + 1)I
+
n+1 n+1
[4i2 4i(n + 1) + (n + )2] + (n + )2 2(n + 1) 2i (n + 1)
n+1 n+1
n=+ (n 1) (5.8)
< nl,
for all 1 < i < n, where i and n are both integers. The maximum value of Equation
(5.8) is n 1 when
[2i (n + 1)]2 + 12i (n + 1) = 0. (5.9)
Now we consider two cases:
Case 1: n is an odd number. The condition (5.9) is satisfied if we choose
i = (n + 1)/2. Putting this i into Equation (5.7), we obtain n 1.
Case 2: n is an even number. The condition (5.9) is satisfied if we choose
i = n/2. Putting this i into Equation (5.7), we obtain n 1. a
Lemma 5.2.3 Let fi = e+(At) e7(), and 7 e C2. Then Il = O(At h) and
I +1i = O(At. h2).
Proof: Both qualities are proved by the mean value theorem:
(1)
fi g7.+l(At) _g7.(At)
where 7, and 7Y are the value and the derivative of 7 at some z between zi and
zij+. Hence,
If, = O(At. h).
(2)
f, f+ = e7i+(At) e7i(At) (ei+2(At) e,+(At
= ei(At)(Y(At)h) e+(At)(i+(At)h)
= (At)h (e +(t)1 e t)
= (At)h (e (At)(.fy(At))y=h + e Yi) I
where 7, and Y7 are the value and the derivative of 7 at some z between z, and
zi+l, 7i+l, and Y+1 are the value and the derivative of 7 at some Z+l between
zi+1 and zi+2, and =7, =7, and 7 are the value, the first derivative, and the second
derivative of 7 at some zi between zi and zi+2. Hence,
Ifi li+1i (At)h ((J,)2(At)h + 'h)
= O(At.h2).
5.3 MaxNorm Bound for Matrix C1D
The maxnorm bound for the matrix C1D, where C = A 2/hB, and
D = AA + (1 p)B are discussed in this section. Only the case p = 1 is discussed
for the convergence analysis. Note that
C1D = (A B) (AA)
= (A B AA+ (A ) (AA AA)
S(A B A B A+ A ( B) BA
+ A B (AA AA)
=A+h I A1B ABA+ I A'B AAAAA)
2 ( 2 /\ 2
Hence,
<1 + Ah II ABI +, IIAI(AA AA)oo
21 IIAIBII 1 IIAIBlo
whenever /,lA1'Boo < 1, and where A = max A,. Our goal in this section is to
show that ICDlloo = 1 + O(At).
5.3.1 MaxNorm for Matrix C1D in 1Dimensional Space
The preliminary results from Section 5.2 are used to show that IIC1DIk =
1 + O(At) in one dimensional space.
Lemma 5.3.1 Let Ax be the matrix described in Lemma 5.2.1. Let B1 be the n x n
tridiagonal matrix with O's on the diagonal, with did's on the superdiagonal, and with
di's on the subdiagonal. Then
h1Al'Blloo = O(At).
Proof: The matrix B1 has the form
0 di
d2 0 d2
d3 0 d3
d4 0 d4
where each d, = O(At). We can write the matrix B1 as
0 di
0 0 d2
di 0 d3
d2 0 d4
We will compute Ai B' and A B" respectively.
(1) compute Ai'B':
(AIB'),,j = dj1 ((a1)i,j1 (a1)i,+l)
where (a1)i, is the (i, j)th entry of A,1. Hence,
0
2(n + i 1)
(A,B')i,j = +1
n + 2i 1
2n+1
2."I^i
I +
Therefore, for any i,
SI(A,'B')i, =
j=1
if j = 1
if 2 < j < i 1
if j = i
if (i 1) < j < n
2(ni+1) ++ n+2il d,
+1 (di + d2 +" )+ dn1)
n +2i
2(i 2)d + 3d + 2(n i)d where d = max di = O(At)
= A
Bi = B'+ B"
0 0
d2 0 0
di d3 0 0
d2 d4 0 0
.j
(2) compute AiB" :
(a1)i,2d2
= (a1)i,(dj1 dj+)
0
2(ni+1) d
n d2
(j + 1)(n i + 1) (dj dj+)
i(n i k) (djl dj+)
n01
0
ifj = 1
if 1 < j < n
if j = n
if j = 1
if 1 < j < i
if j = i + k (k= 0,1, ,ni1)
if = n
where Id'I = max Idj1 dj+I = O(At h)
< O(At) + 3)(i O(At. h) + (n 1)(n i) O(At. h)
2 2
Therefore,
hll(Ax'BI)Il < hO (Ax'B') + II(Ai'B"),.)
= o(At).
Hence, for any i,
i(A. B")i,jl
j=l1
< 2(n i + 1) d n i+1(3 4 i)ld
S + d2+ I(3+4++id'
n+1 n+1
+ [(ni)+(ni1)+...+2+1] ld'I,
50
Lemma 5.3.2 Let A1 be the matrix described in Lemma 5.2.1 and A1 be the
diagonal matrix with Ai on the diagonal. Then
IAxI(A1A1 AixAi1)1 = O(At).
Proof: To prove this Lemma, we compute the maxnorm of Ai (A1A1 AIAI).
Let F = AxA1 A1A1. Then we can write F as
O fi
fi 0 f2
F = f2 0 f3 where f = e^+1(At) e(At)
0 f2 0 flf2
f 0 f3 0 0 f2 f3
= f2 0 f4 + 0 0 f3f4
= F'+F".
The same process in Lemma 5.3.1 is used to compute A1iF' and Ai'F"
respectively.
(1) compute A'F' :
(A1 F'),, = ((a1)i,jl (a)i,j+,)
2(ni+1l)f ifj
(n + 2i + 1) if
2i
n fj ifj>i
Therefore, for any i,
n
I(A= F'),, =
j=1
2(n l + f2l + + Ifi)+ +2i 1
n + 1 ,7+1
2i
+ (IAOi+l A + "fi2 + + If.)
n +1
< 2(i 1)If + 31f + 2(n i) f where f = max fI = O(At h)
= (2n + 1)O(At h)
= O(At).
(2) compute Ai'F" :
0
if j = 1
(Ai F")ij =
I
(a1)i,j_(fj1 fj) if j 0 1
0
(j 1)(n i + 1)(f f)
n+1 + k
i(n ni + 2 k) (f1 f3)
if j = 1
if 2 < j
if j=i + k (k= 2,3,. ,ni)
Hence, for any i,
SI(A 'F")i,j <  (1 +2++ i) f'
n+1
j= 1
+[(ni) + (ni 1) ++ 3 + 2] f',
where If'I = max fj_ fj = O(At h2)
< (i+ 1)io(At.h2)+ (n i+ 2)(n i 1)O(At.h2)
2 2
= O(At).
Therefore,
jIAi (AIA1 A1A1)II < II(AI'F')loo + I(A1 F")1io
< 0(At).
Theorem 5.3.3 Let A, and B1 be the matrices described in Lemma 5.2.1 and
Lemma 5.3.1, respectively. Then
IIC1DIl _< 1 + O(At),
where C = A1 Bx1 and D = A1Ax.
Proof: Note that
_lh AlA Bn llo ,Axi(AxAi AxAi)oo
IC1DI < A+ 1 h A1 + AA) (5.10)
21 llA Blloo 1 llAix1Bllo
whenever IIAi1B13II < 1, and where = max Ai < 1.
Therefore, from Lemma 5.3.1 and Lemma 5.3.2, we have
IIC'DII0 < A+ O(At) + O(At)
< 1+ O(At) .
5.3.2 MaxNorm for Matrix C1D in 2 and 3Dimensional Spaces
Now we extend the results from the 1dimensional case to 2 and 3dimensions.
First, we consider the 2dimensional case. Let A2 be a tridiagonal block matrix
with (A1 + 21)'s on the diagonal block, and I's on the superdiagonal and
subdiagonal block. Let B2 and A2 be diagonal block matrices with B1's and Al's
on the diagonal block respectively. Then
hA21B2 =
I+2A11 A,1
A'1 I+2Ai' A,1
h
A' I+ 2AI1 A,'
A1 I+ 2Ai1
AI1 B,
Ai1 B1
x
Ai1 B1
and we can rewrite it as
hA21B2 =
Al1 B1
AI1 Bx
hS2 (5.11)
AI1
where S2 is the matrix
I+2A11
A11
I+ 2A,
A1 I+ 2Af1 A,1
A,1 I+ 2A1
Thus,
hIA21B21. < hi AHA BIjlIjS2oo.
By Lemma 5.3.1, we have IAi1B1Io = O(At/h). Hence,
hIA2 B21loo < O(At)S2oo.
Note that S2 is difficult to analyze using the maxnorm. Hence, S2 is evaluated on
a computer to get a numerical value for JS21oo. For example, if A1 is 10 x 10 and
m = 10, so that S2 is 100 x 100, then S21oo 1.6687. If A1 is 20 x 20 and m = 20,
so that S2 is 400 x 400, then JS21K 2.0794. If A1 is 40 x 40 and m = 40, so that
S2 is 1600 x 1600, then IS2oo ; 2.5041. Therefore, we have
hIIA21B2[[o, = O(At).
A,1
(5.12)
Also, we have
A2'(A2A2 A2A2)
1
I + 2A11 A11
A11 I+2A1 A,1
A1' I+2A11
A1 A1A1 A1A1
x
Aj1 AIA1 AIA1
A,1 A1A1 AIA1
= S2
Ai~1 A1A1 A1A1
and then
IA2 (A2A2 A2A2)j IIAi'(A1A1 AiAi)oo S211,
= O(At), (5.13)
since IAi'(A1A1 AIA1)II = O(At) by Lemma 5.3.2. From Equation (5.10),
we have the estimate HIC1Dl, = 1 + O(At).
Three dimensions can be analyzed in the same manner. We write A as a
tridiagonal block matrix with (A2 + 21)'s on the diagonal block and with I's on
the superdiagonal and subdiagonal block. Write B and A as diagonal matrices with
B2's and A2's on the diagonal block respectively. Write
hA'B =
hSs
A2I
and
A'(AA AA)=
A21
where the matrix S3 is
I + 2A21
A1 I
A2
A2A2 A2A2
A2A2 A2A2
A~'2
+ 2A2 A2'1
A2' I+ 2A A21
A2' I+ 2A2'
To discuss the maxnorm, we still need to evaluate the factor S3. Using a computer
to evaluate the maxnorm of S3, we have: if As is 5 x 5, so that S3 is 125 x 125,
then IS3Ioo 5 1.9841. If A2 is 10 x 10, so that S3 is 1000 x 1000, then IS311 oo
2.5126. Similarly, we have
hIjA'B < h A21B2lootlSaloo
= O(At)
IA'(AA AA)II < IA21(A2A2 A2A2)IUISa3oo
= O(At)
by Equation (5.12) and Equation (5.13) from the 2dimensional space. Finally, we
have IIC'Doo = 1 + O(At) from Equation (5.10).
5.3.3 Numerical Results for Spectral Radius of C1D
Next we use the approximate value of y (Section 3.3) to calculate some
numerical results for spectral radius of C1D, (Table 51, Table 52 and Table
53.), in the domain Q1 when At = 0.1 s. Numerically, we observe that the spectral
radius of C'D is strictly less than 1.
Table 51: Spectral radius for C1D in one dimensional space
Matrix dimension (n x n) Spectral radius for C1D
n = 50 0.9990
n = 100 0.9991
n = 200 0.9992
n = 400 0.9992
n = 800 0.9992
n = 1600 0.9992
Table 52: Spectral radius for C D in two dimensional space
Matrix dimension (n x n) Spectral radius for C1D
n = 100 (m = 10) 0.9960
n = 400 (m = 20) 0.9984
n = 1600 (m = 40) 0.9989
n = 6400 (m = 80) 0.9991
Table 53: Spectral radius for C1D in three dimensional space
Matrix dimension (n x n) Spectral radius for C1D
n = 125 (m = 5) 0.9715
n = 1000 (m = 10) 0.9960
n = 8000 (m = 20) 0.9984
n = 27000 (m = 30) 0.9988
CHAPTER 6
NUMERICAL EXPERIMENTS
This chapter briefly reports some numerical experiments by using the finite
difference scheme described in Chapter 3. In these applications, the FORTRAN
package FISHPAK (see Appendix B) is utilized to deal with large sparse ma
trices in our model and to solve linear separable elliptic equations in three space
dimensions.
The conductivity function is computed as
a = 10p(z), (6.1)
where p(z) is the quartic polynomial (see Equation (3.21) and Table 31) and the
altitude z is measured in kilometers. Let us assume as in Hager et al. [7] that the
current per unit volume at position r is
V J = Apewirsp12 ANewrsN2,
where Sp and SN are the centers of the positively charged and negatively charged
current generators, respectively. The parameters Ap and AN are chosen so that
the total current generated by either the positive or the negative charge centers has
magnitude 1 A. Hence, we have
Ap h ew2dV = 1,
Jz>hi
where hi is the altitude for the positive charge center, R2 = Ir Sp 2, and z is the
altitude measured with respect to the positive charge center. The cylindrical shell
method is applied to the above integral, as shown in Figure 61, to obtain
(positive
zaxis (altitude)
P
Figure 61: Cylindrical shell of the volume integral
ew2 dV
ewh(p2+z2)2rpdzdp
J zht p20
= ( 27rpe,p2dp)
\Jp>0 /
(\z>hl
eWZdz)
7 (7) ) 1/2 1 7 1/2 )]
7 (r3/2 [l+ erf(v hl)]
where er is the error function 2
where erf is the error function defined as
erf (x) = 
V'77
f e` dt.
Jo
As a result, Ap is calculated as
S(w 3/2
Ap = 
7T
2
[1 + erf (v/hi)]'
and, similarly, AN can be obtained
A ( N3/2
AN= 
7r
2
[1 + erf(vwh2)]'
where h2 is the altitude for the negative charge center. Therefore,
V J = ( 312 2 Wl2 2 wSN 1
SW [1 + erf (vihl)] [1 + erf(v/( h2)]
is used in our model. In each experiment, the initial potential of the entire domain
is assumed to be zero. The centers of the current generators are positioned on the
zaxis, at altitudes 5 km () and 10 km (+).
Four experiments are conducted in this chapter. Section 6.1 uses different
values of w to find the breakdown time with the ambient conductivity and the
reduced conductivity. Section 6.2 explores a finer mesh for our model with an
acceptable accuracy for the potential 0. Section 6.3 and Section 6.4 check the
temporal error and the spatial error numerically respectively.
6.1 Experiment 1
This experiment considers the domain Q1, previously introduced:
i = {(x,y,z)\0 < z <100km, 50km < x,y< 50 km}
For conductivity, we use the ambient conductivity given in Equation (6.1). For
different values of w, we integrate the Equation (3.7) solving for potential 0. The
integration is stopped when the electric field reaches the breakdown threshold,
which is taken to be 250,000 V m1 [7]. We use the uniform mesh in all x, y,
and zdirections with mesh size 1000 m, and the time step is At = 1/16 s. As
shown in Table 61, the breakdown time is calculated based on different values of w.
Table 61 demonstrates that the breakdown time shortens rapidly as w
increases. For example, when w = 4, the breakdown occurs within 15.5625 seconds.
With w = 16, the breakdown time drops to 1.1850 seconds. This result is supported
by [7]: as w increases, the current becomes more concentrated at sp and SN.
Observe that when w = 1, there is no breakdown within one hour. According
to our numerical results, the maximum electric field approaches the limit 232429.88
Table 61: Breakdown time comparison for different values of w
breakdown time breakdown time
w value with the ambient conductivity with the reduced conductivity
(without cloud) (with cloud)
w = 1 00 64.7500 s
w = 4 15.5625 s 13.6875 s
w = 16 1.1850 s 1.75000 s
V/m which is less than the breakdown threshold 250,000 V/m, as t tends to
infinity. Therefore, the breakdown never happens when w = 1.
In a thundercloud, the conductivity is less than the ambient conductivity that
appears in Equation (6.1). Next, we consider the effect of reducing conductivity in
order to simulate a cloud. More precisely, the ambient conductivity is multiplied
by the factor 0.05 as described in [7]. The results are shown in the third column of
Table 61. Again, the breakdown time decreases dramatically as w increases. This
pattern is consistent with the previous computed breakdown times for the ambient
conductivity.
As can be seen in Table 61, reducing in conductivity to simulate a cloud leads
to the occurrence of breakdown when w = 1, while with the ambient conductivity,
charge is conducted into the surrounding atmosphere faster enough to prevent the
electric field from reaching the breakdown threshold. In the remaining experiments,
the conductivity is given by Equation (6.1) times 0.05 (to simulate a cloud).
6.2 Experiment 2
A finer mesh in our model is expected to obtain more accurate results.
However, such refinement is not only time consuming, but also increasing memory
usage. In order to maintain the accuracy of our results efficiently, a graded mesh
is considered. We place a large percentage of grid points near the charge centers
along the z direction where the potential changes rapidly. The grid along the
zaxis is chosen in the following manner [7]:
zk = kT for 1 < k < n (6.2)
and
zk = 2H + rs(k) for n < k < n, (6.3)
2
where H is the altitude of the positive charge center (10 km), T = 4H/n, and s is
chosen so that 2H + rs(n/ = 100 km is the height of the domain. According to
Hager et al. [7], the vicinity of the charge centers is critical to efficiency. Therefore,
the distribution of the grid points is uniform below 2H and geometric from 2H
to 100 km. We still use the uniform mesh for the x and ydirections as in
Experiment 1.
Let us consider two domains: Q1 as described in Section 5.1 and Q2 as follows,
Q2 = {(x, y, z)10 < z < 100 km, 25 km < x, y < 25 km}.
In this experiment, the potentials are found at the time 65 seconds with the time
step At = 1/2 s and the ambient conductivity reduced by the factor 0.05 when
w = 1. The graded mesh is applied with 200 grid points along zaxis (mesh size
= 200 m for 0 < z < 20 km). The mesh size is 1000 m for both in x and y
directions. Figure 62 shows the potentials in two domains, Q and Q2, for x = 0,
y = 0, and 0 < z < 100 km. The potential from Q1 is shown in a solid line and that
from Q2 is shown in a dashed line. The absolute difference in potential is about 105
and the relative difference is about 104 between these two domains. Therefore, we
cannot see any difference for the potentials in Figure 62. Figure 63 and Figure
64 are the potentials near the negative and positive charge centers respectively,
along the zaxis when x = 0 and y = 0.
Since the relative difference for potentials is quite small between these two
domains, we only focus on the smaller domain Q2 for the remaining experiments.
x 10
See Figure 6.4
4
2
0
2
4 S
See Figure 6.3
0 1 2 3 4 5 6
Altitude in meters
7 8 9 10
x 10
Figure 62: Potentials (i(0, 0, z), 0 < z < 100 km) in two domains, Q1 and Q2
) 5000 5100
Altitude in meters
Figure 63: Potentials (0(0, 0, z),
km) in two domains, Q, and Q2
4.6 km < z < 5.4 km) near the negative charge (5
 potential in the domain Qt
potential in the domain Q,
Altitude in meters x 10
Figure 64: Potentials (0(0, 0, z), 9.6 km < z < 10.4 km) near the positive charge
(10 km) in two domains, Q1 and Q2
xaxis in meters
2.5
x 10'
Figure 65: Potentials (O(x, 0, 5) and (x, 0, 10),
xaxis at the altitude 5 km and 10 km
25 km < x
< 25 km) along the
Further, we find the potentials (see Figure 65) along the xaxis at the altitude 5
km (the negative charge center) and 10 km (the positive charge center) when y = 0.
In Figure 65, the solid line represents the potential at the altitude 10 km and the
dashed line represents the potential at the altitude 5 km.
6.3 Experiment 3
According to the error estimate described in Theorem 4.2.4, when p = 1/2
(CrankNicholson scheme), the temporal error term should be smaller than for
other values of p. To verify the temporal error analysis numerically, the potential
is estimated at time t = 65 s, which is close to the breakdown time 64.75 seconds
when w = 1. The graded mesh (see Equation (6.2) and Equation(6.3)) is applied
with 200 grid points along zaxis (mesh size = 200 m for 0 < z < 20 km). The
mesh size is 1000 m for both in x and y axes. We compare the errors in different
time steps (1/2, 1/4, , 1/64) for p = 1/2 and p = 1. To analyze the error
numerically, we use the maximum norm:
I IAl = max
ij,k
where (OAt)i,j,k is the computed potential using the time step At at each mesh
point in the domain at t = 65 s, and bAt is the potential vector with elements
(At)ij,k. First, we consider the potential corresponding to the time step At =
1/128 s as the "exact" potential and subtract it from the potentials corresponding
to the compared time steps to obtain error estimates at all mesh points. The
maximum norm error at different time step is shown in the third column of Table
62 and Table 63. Second, the relative errors are estimated by the ratio
IIhAt al /2811oo
as shown in the fourth column of Table 62 and Table 63.
Table 62: Relative error in the potential at t = 65 seconds using the maximum
norm for CrankNicholson scheme (pt = 1/2).
Time step maximum norm maximum norm error CrankNicholson
(s) potential with respect to relative error
time step 1/128 s (p = 1/2)
1/2 564431185.91904 80.21113 1.42110 x 107
1/4 564431163.61478 20.30029 3.59659 x 108
1/8 564431157.94013 5.11413 9.06068 x 109
1/16 564431156.50102 1.27381 2.25681 x 109
1/32 564431156.13736 0.30471 5.39852 x 1010
1/64 564431156.04677 0.06052 1.07192 x 1010
1/128 564431156.04677 0.00000 0.00000
Table 63: Relative error in the potential at t
norm for Backward Euler scheme (p = 1).
= 65 seconds using the maximum
Time step maximum norm maximum norm error Backward Euler
(s) potential with respect to relative error
time step 1/128 s (0 = 1)
1/2 564422403.52880 35278.79275 6 2111I. x 10
1/4 564426772.69223 17375.01373 3.07832 x 10"
1/8 564428962.54654 8411.09363 1.49019 x 105
1/16 564430058.82113 3926.08120 6.95582 x 106
1/32 564430607.30153 1682.80248 2.98141 x 106
1/64 564430881.63001 560.96757 9.93864 x 107
1/128 564430881.63001 0.00000 0.00000
In Table 62 which corresponds to the CrankNichlinl,,Ti scheme, if the time
step At is multiplied by 1/2, the relative error is multiplied by 1/4. Hence, the
error for CrankNicholson scheme (pj = 1/2) behaves like O(At)2. In Table 63
which corresponds to the backward Euler difference scheme, if the time step At is
multiplied by 1/2, the relative error is also multiplied by 1/2. Therefore, the error
for backward Euler difference scheme (p = 1) behaves like O(At). This suggests
that the CrankNicholson scheme is far superior to the backward Euler difference
scheme.
Note that the forward difference scheme (pi = 0) does not work, since it is not
stable as we described in Section 4.1. Experimentally, taking p = 0, we find that
the maxnorm of the potential grows exponentially in t. The plot of loglo(IAtIloo)
versus t appears in Figure 66, with time step At = 1/8 s. Observe that the
maximum potential at t = 65 seconds is 1.04737 x 1089.
90
80
70
o
00
50
40
30
20
10 
0
0 10 20 30 40 50 60 70
Time in seconds
Figure 66: Maximum norm of potential versus time, where At = 1/8 s and t = 0
6.4 Experiment 4
The last study is to check the spatial error term. We take the time step At =
1/2 s and use CrankNicholson scheme p = 1/2. Four different graded meshes are
used in the domain Q2, where
Q2 = {(x, y,z)10 < z < 100 kni, 25 km < x, y < 25 kmn}
and the grid points are chosen according to Equation (6.2) and Equation (6.3)
along the zaxis:
Case 1: 30 uniform grid points in the x and ydirections, and 100 grid
points in the z direction (mesh size = 400 m for 0 < z < 20 km, and
exponentially increasing mesh for 20 < z < 100 km, from Az = 400 m near z
= 20 km to Az = 8043.61 m near z = 100 km).
Case 2: 60 uniform grid points in the x and ydirections, and 200 grid
points in the z direction (mesh size = 200 m for 0 < z < 20 km, and
exponentially increasing mesh for 20 < z < 100 km, from Az = 200 m near z
= 20 km to Az = 4652.41 m near z = 100 km).
Case 3: 120 uniform grid points in the x and ydirections, and 400 grid
points in the z direction (mesh size = 100 m for 0 < z < 20 km, and
exponentially increasing mesh for 20 < z < 100 km, from Az = 100 m near z
= 20 km to Az = 2629.65 m near z = 100 km).
Case 4: 240 uniform grid points in the x and ydirections, and 800 grid
points in the z direction (mesh size = 50 m for 0 < z < 20 km, and
exponentially increasing mesh for 20 < z < 100 km, from Az = 50 m near z
= 20 km to Az = 1462.03 m near z = 100 km).
Let "n, 02n, 4n,, and "8n denote the discrete potential evaluated along z
axis at z = 10 km, where n, 2n, 4n, and 8n are the numbers of mesh intervals
along zaxis in case 1, 2, 3, and 4 respectively. We focus on z = 10 km since the
the largest potential is achieved here. Suppose that 1 is the "exact" potential
at z = 10 km, and the sequence 02n 04,, .. is converging to the limit 0P.
Extrapolation techniques (see [3], Sectionl.8) are used to estimate various constants
in the expression of the error estimate.
First, based on Theorem 4.2.4, the spatial error is O(h'2). Hence,
n" = 4 + En, (6.4)
where En is the error term with the property IE, I < Ch2, and C is a constant, and
h = 1/n. In extrapolation techniques, we treat the inequality EnI < Ch2 as an
equality. Two elements in the sequence, say 'f and 2n, are used in order to solve
the constant C. Plugging "n and 02n into the error relation Equation (6.4), we
have
= E., (6.5)
02n 0 = E2n. (6.6)
Subtract Equation (6.6) from Equation (6.5) to obtain
Sn 02n = EnE2n
= Ch2 C()
= ,Ch2.
4
Therefore,
C 0 o2n
3 h2
4
where h = 1/n.
Any two potential vectors 0t and 2i can be employed in the same way as
above to find the constant C:
C = 3 (6.7)
4i2
Table 64: Estimates for the constant C generated by Equation (6.7)
n (0n _2n C
100
175594927.29868 2.34 x 108
200
25121434.49608 1.34 x 108
400
5574244.17420 1.19 x 108
800
In Table 64, only four different meshes are considered. If n < 100, the
potential vector is not accurate enough, and if n > 800, there is not enough
memory to do computations on our computer network. The maximum norms of
the difference of two successive potential vectors are shown in the second column
of Table 64. Observe that the maximum norm of the difference are getting
smaller. We can predict that it approaches 0 when n tends to infinity if more data
(potentials for larger n) provided. This shows that the sequence On is converging.
The estimates for the constant C are shown in the third column of Table 64. The
relative difference for the last two estimates of C is 0.13.
Next, let us assume that the spatial error is O(hP), that is,
n = + En,
where IE, I < ChP C and p are constants. Using the similar analysis as above, we
will compute the constants p and C. For three elements 0', 2", and 4", we have
o2n On = ( p I 1
(1= f 11
\2P 4P
Ci ( 1
= 1 (1 (6.9)
Equation (6.8) divides Equation(6.9) to get
2P = n 2n
2n 4n'
Therefore,
p = log2
P 1092 02n 04n '
and
n 2n
can be found, where h = 1/n.
Given any three elements 0", 3', and Ok with the property j2 = ik, the formula
for constants p and C are:
\(i j) i<(Oi 0j)
P = lg2 k I C = 3 O (6.10)
1
The computed constants p and C are shown in the third and fourth columns
of Table 65 respectively. The best estimate for the constant C is 1.45 x 108 from
Table 65, while the best estimate is 1.19 x 108 from Table 64. The relative
difference between them is 0.22. From the third column of Table 65, we can see
that the constant p is getting close to the theoretical value 2.
Table 65: Estimates for constants p and C generated by Equation (6.10)
n 2n p C
100
175594927.29868
200 2.81 2.04 x 108
25121434.49608
400 2.17 1.45 x 108
5574244.17420
800
CHAPTER 7
CONCLUSION
In this dissertation, Maxwell's equations are used to establish an electric field
model between the surface of the Earth and the ionosphere:
PDE: A+A+ z+ =0 (7.1)
at E ( E aZ E
BC: = 0, ifz = 0
= 250, 000, if z = Dz
f = c dz, if x = D or y = Dy
Domain: Q = {(x, y, z): \xl < Dx, \y < Dy, 0 < z < Dz}
The conductivity function a varies by many orders of magnitude between the Earth
and the ionosphere, which causes the partial differential equation (7.1) to be stiff.
We develop stable finite difference approximations to our model that take into
account the a variation. Error estimates are obtained by using the maxnorm. Four
numerical experiments are conducted to verify our theoretical analysis as described
in Chapter 4 and Chapter 5.
APPENDIX A
APPROXIMATION FOR CONDUCTIVITY FUNCTION
MATLAB programs are used to obtain the approximation for the conduc
tivity function. The detailed approximation for the conductivity function from
United Sates Air Force Cambridge Research Laboratories [27] is given in Appendix
A.1. The approximation from Lightning: physics and effects [21] is discussed in
Appendix A.2.
A.1 Approximation from Figure 9.6 [27]
Four points (see Table A1) on the curve B from Figure 9.6 [27] are chosen to
find an approximation for conductivity. The conductivity graph in Figure 9.6 [27]
contains the data up to the altitude 27.828 km (about 100 ft), while our model
domain is 100 km. All four points are chosen to have a cubic fit for the altitude
0 < z < 27.828 km, and the approximation can be obtained by MATLAB as
follows:
a (z) = 1013.40+0.1600zO.005175z2+O.00009077z3
For altitude 28.828 < z < 100 km, the last three points (P2, P3, and P4) are picked
to have a linear fit, and the approximation is
a(z) = 1012.88+0.06718z.
Figure A1 shows the graph of the approximation for the conductivity function.
A.2 Approximation from Figure 33
Thirteen points (see Table A2) in Figure 33 are picked to approximate the
conductivity function. The conductivity is chosen as the average of shaded region
in Figure 33 for the altitude above 30 km. The advantage for using this figure
Table Al1: Four points chosen from Figure 9.6 [27]
Point Height logo0 (conductivity)
(km) (1/ohmm)
P1 2.7432 13
P2 12.954 12
P3 20.7264 11.5
P4 27.82824 11
log Joconductivity), (Ohmmeter)
Figure A1: Electric conductivity as a function of altitude approximated from
curve B in Figure 9.6 [27]
to estimate the conductivity is that Figure 33 demonstrates a more completed
conductivity data than Figure 9.6 in [27].
Table A2: Thirteen points chosen from Figure 33
Point Height
(km)
P1 4.5
P2 12.5
P3 26
P4 40
P5 48
P6 52.5
P7 60
PA 70
P9 78.5
Pio 85
PI1 89
P12 95
P13 100
loglo (conductivity)
(S/m)
A quartic polynomial is found to
conductivity function is estimated by
approximate those points in Table 2. The
a = 10p(z)
where p(z) is the quartic polynomial
p(z) = 13.16 + 0.06351z + 0.0008682z2 0.000008674z3
+0.00000005781z4,
where the altitude z is measured in kilometers, and its graph is shown in Figure
A2.
78
100
90
80
70
S60
50
40
30
20
10
0 L'
14 12 10 8 6 4 2 0
log Pconductivity), S/m
Figure A2: Electric conductivity as a function of altitude approximated from
Figure 33
APPENDIX B
FISHPAK
The FORTRAN package FISHPAK developed by Swarztrauber and Sweet
is used for our numerical experiments. Documentation for the package is given
in ACM Algorithm 541 [24]. The original package is written in single precision.
In order to get more precise numerical results, a double precision version package
is obtained by a little adjustment. The subroutine POIS3D is used to solve for
potentials in our model. POIS3D solves the linear system of equations for unknown
X values,
C1*(X(I1,J,K) 2.*X(I,J,K) + X(I+1,J,K)) +
C2*(X(I,J1,K) 2.*X(I,J,K) + X(I,J+1,K)) +
A(K)*X(I,J,K1) + B(K)*X(I,J,K) + C(K)*X(I,J,K+1)
= F(I,J,K),
where I = 1, 2, ., L, J = 1, 2, M, and K = 1, 2,  ,N. The indices K1 and
K+1 are evaluated modulo N, i.e. X(I,J,0)=X(I,J,N) and X(I,J,N+1)=X(I,J,1).
The unknowns X(0,J,K), X(L+1,J,K), X(I,0,K), and X(I,M+1,K) are assumed
to take on certain prescribed values. The operation count for POIS3D is roughly
proportional to
L x M x N x (log2 M + log2 N + 5).
REFERENCES
[1] E. J. Adlerman and E. R. Williams, Seasonal variation of the global electrical
circuit, Journal of Geophysical Research, 101, 1996, pp.2967929688.
[2] M. B. Baker, H. J. Christian, and J. Latham A computional study of the
relationships linking lightning frequency and other thundercloud parameters,
Quarterly Journal of the Royal Meteorological Society, 121 1995, pp.1525
1548.
[3] W. W. Hager, Applied Numerical Linear Algebra, PrenticeHall, Englewood
Cliffs, NJ, 1988.
[4] W. W. Hager, Updating the inverse of a matrix, SIAM Rev., 31 1989,
pp.221239.
[5] W. W. Hager, A discrete Model for the Lightning Discharge, Journal of
Computational Physics, 144 1998, pp.137150.
[6] W. W. Hager, J. S. Nisbet, and J. R, Kasha, The evolution and discharge of
electric fields within a thunderstorm, Journal of Computational Physics, 82,
1989, pp.193217.
[7] W. W. Hager, J. S. Nisbet, J. R, Kasha, and W.C. Shann, Simulations of
electric fields within a thunderstorm, Journal of the Atmospheric Sciences, 46
1989, pp.35423558.
[8] L. C. Hale Middle atmospheric electrical structure, dynamics, and coupling.
Earths Middle Atmosphere Advances in Space Research. 4 1984. pp 175 1.,
[9] J. H. Helsdon, Jr. and R. D. Farley, A numerical modeling study of a
Montana thunderstorm: 1. Model results versus observations involving
nonelectrical aspects, Journal of Geophysical Research, 92 1987, pp.5645
5659.
[10] J. H. Helsdon, Jr. and R. D. Farley, A numerical modeling study of a
Montana thunderstorm: 1. Model results versus observations involving
electrical aspects, Journal of Geophysical Research, 92 1987, pp.56615675.
[11] J. H. Helsdon, Jr., R. D. Farley, and G. Wu, Lightning parameterization in
a storm electrification model, 1988 Proceedings International Conference on
Atmospheric Electricity, pp.849854
[12] J. H. Helsdon, Jr., G. Wu, and R. D. Farley, An intercloud lightning
parameterization scheme for a storm electrification model, Journal of
Geophysical Research, 97 1992, pp.58655884.
[13] H. W. Kasemir, A contribution to the electrostatic theory of a lightning
discharge, Journal of Geophysical Research, 65 1960, pp.18731878.
[14] D. R. MacGorman, A. A. Few, and T. L. Teer, Layered lightning activity,
Journal of Geophysical Research, 86 1981, pp.99009910.
[15] D. R. MacGorman, J. M. Straka, and C. L. Ziegler, A lightning parameteri
zation for numerical cloud models, Journal of Applied Meteorology, 40 2001,
pp.459478.
[16] E. R. Mansell, D. R. MacGorman, Ziegler, and J. M. Straka, Simulated
threedimensional branched lightning in a numerical thunderstorm model,
Journal of Geophysical Research, 107 2002, doi:10.1029/2000JD000244.
[17] L. Niemeyer, L. Pietronero, H. J. Wiesmann, Fractal dimension of dielectric
breakdown, Physical Review Letter, 52 1984, pp.10331036.
[18] F. Pockels, Uber das Magnetische Verhalten Einger Basaltischer Gesteien,
Ann. Phys. Chem., 63 1897, pp.195201.
[19] F. Pockels, Bestimmung Maximaler EntladungsStromStiirken aus Ihrer
Magnetisirenden Wirkung, Ann. Phys. Chem., 65 1898, pp.458475.
[20] F. Pockels, Uber die Blizentlandungen Erreicht Stromstdirke, Phys. Z., 2
1900, pp.306307
[21] V. A. Rakov and M. A. Uman, Lightning: physics and effects, Cambridge
University Press, Cambridge, New York, 2003.
[22] F. Rawlins, A numerical study of thunderstorm electrification using a three
dimensional model incorporating the ice phase, Quarterly Journal of the
Royal Meteorological Society, 108 1982, pp.779800.
[23] J. C. Strikwerda, Finite Difference Schemes and Partial Difference Equations,
Champman & Hall, New York, NY, 1989.
[24] P. N. Swarztrauber and R. A. Sweet, Algorithm 541: Efficient Foritran
subprograms for the solution of elliptic partial differential equations, A CM
Trans. Math. Software, 5 1979, pp.352364.
[25] T. Takahashi, Determination of lightning origins in a thunderstorm model,
Journal of the Meteorological Society of Japan, 65 1987, pp.777794.
[26] M. A. Uman, The Lightning Discharge, Academic Press, Orlando, Florida,
1987.
[27] United States Air Force Cambridge Research Laboratories, Handbook of
Geophysics, Macmillan, New York, 1960 (chapter on atmospheric electricity
revised in 1970).
[28] R. S. Varga, Matrix Iterative Analysis, PrenticeHall, Englewood Cliffs, NJ,
1962.
[29] H. J. Wiesmann and H. R. Zeller, A fractal model of dielectric breakdown
and prebreakdown in solid dielectrics, Jounral of Applied Physics, 60 1986,
pp.17701773.
[30] E. R. Williams, C. M. Cooke, and K. A. Wright, Electrical discharge
propagation in and around space charge clouds, Journal of Geophysical
Research, 90 1985, pp.60546070.
[31] C. T. R. Wilson, On Some Determinations of the Sign and Magnitude of
Electric Discharges in Lightning Flashes, Proceedings of the Royal Society of
London Series A, 92 1916, pp.555574.
[32] C. T. R. Wilson, Investigations on Lightning Discharges and on the Electric
Field of Thunderstorm, Philosophical Transactions of the Royal Society of
London Series A, 221 1916, pp.73115.
[33] C. L. Ziegler and D. R. MacGorman, Observed lightning morphology relative
to modeled space charge and electric field distributions in a tornadic storm,
Journal of the Atmospheric Sciences, 51 1994, pp.833851.
[34] C. L. Ziegler, D. R. MacGorman, J. E. Dye, and P. S. Ray, A model
evaluation of non inductive graupelice charging in the early electrification
of a mountain thunderstorm, Journal of Geophysical Research, 96 1991,
pp.833851.
[35] C. L. Ziegler, J. M. Straka, and D. R. MacGorman, Simulation of early
electrification in a threedimensional dynamic cloud model, 11th Interna
tional Conference on Atmospheric Electricity, Commission on Atmospheric
Electricity, 1999, pp.359362.
BIOGRAPHICAL SKETCH
ShuJen Huang was born in Taiwan in 1970. She was awarded a Bachelor of
Science degree in applied mathematics in 1992 from FuJen University, Taipei,
Taiwan. She worked as a teaching assistant at FuJen University for 3 years after
her graduation. In 1995, ShuJen started her graudate study in mathematics at the
University of Florida, from which she received her Ph.D. in mathematics in 2005.
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.
Villiam W. Hager, Chair
Professor of Mathematics
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doctor of Philosophy. A
Jay Gopalakrishnan
Assistant Professor of Mathematics
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.
Shari Moskow
Associate Professor of Mathematics
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.
iS. Pilyugin
Assistant Professor of Mathematics
I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.
Vladimir A. Rakov
Professor of Electrical and Computer
Engineering
2611111'1' Ar
J
This dissertation was submitted to the Graduate Faculty of the Department
of Mathematics in the College of Liberal Arts and Sciences and to the Graduate
School and was accepted as partial fulfillment of the requirements for the degree of
Doctor of Philosophy.
May 2005
Dean, Graduate School
LD
1780
20 DA
UNIVERSITY OF FLORIDA
3 1262 08554 258611111
3 1262 08554 2586

Full Text 
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EO5R0OPSG_127EA7 INGEST_TIME 20111102T16:58:32Z PACKAGE AA00004682_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES
PAGE 1
MULTISCA LE DISCRETIZATIO OF ELECTRICFIELD EQUATIO S B y SHUJEN H UANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE U I VERSITY OF FLORIDA I PART I AL FULFILLME T OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY U IVERSITY OF FLORID A 2005
PAGE 2
Copyright 2005 by ShuJen Hu ang
PAGE 3
To DengShan for now a nd for e ver.
PAGE 4
ACKNOWLEDGME TS First of a ll I would lik e to e xpr ess my gratitude to my supervisory committee chair, Profes so r Willi a m W. Hager. Without his assistance, t his dissertation c ould not have been completed. In addition, I appreciate his e ncouragem ent throughout this p e riod of the r esea rch. His confidence in me always strengthens my beliefs in my ability to complete this research. Second, I would a l so lik e to thank D r Jay Gopa l akrishn a n Dr. Shari Moskow. Dr. Sergei S. Pil y u g in a nd Dr. Vladimir A. Rakov for serving on my superv i sory com mittee. Their valuable suggestions have been very h e lpful to my r esea rch. Third, thanks go to m y officemates (Dr. Soonchul Park B eyza Hongchao, and Sukanya) and all colleagues and fri e nds in the Department of Mathematics at the U niv e rsit y of Florida. Their company alleviate d the stress a nd frustration of t hi s ti me. Last but not l east I wish to expr ess my spec ial thank s to my family: to my hu sband, D engS han for hi s lov e pati e n ce and e ncourage m ent ; to our d a u ghte r Marie, for bein g a glorious joy to us ; to my mother for h e r support and ass istan ce an d for h e lping me to take ca r e of my daught er throughout the dissertation writing process ; to my parents inlaw for t h e ir wholehearted und erstandi n g an d encouragement during my stu dies; an d to my si ste r for her un stopp in g support and e n couragem e nt. T h ey hav e always stoo d by me. Without t h e ir support, t hi s di ssertation could no t hav e be e n completed s u ccess fully lV
PAGE 5
TABLE OF CONTE TS ACK OWLEDGME TS LIST OF TABLES LIST OF FIGURES IV Vll Vlll KEY TO ABBREVIATIO S A D SYMBOLS ABSTRACT X Xlll CHAPTER 1 2 3 4 5 INTRODUCTION LITERATURE REVIEW 2 1 Dynamic Cloud Models 2. 2 Approaches for Lowering Charge due to Ligh t nin g 2. 2 1 Lightning Model with Bulk Flashes .... 2.2.2 Lig htning Model with Explicit Lig h tning Chann e l s 2. 2.3 H age r 's Model THE MODEL ....... 1 4 4 6 6 7 10 12 3.1 Formulation of the Equations 12 3.2 Discretized Equations . . 14 3 3 Important Rol e and Approximation for Co nductivity Function 22 3.4 Discussion: d = ,z (1 e r( 6.t)) / 1 26 STABILITY AND ERROR ANALYSIS. 4.1 Stability . . . . 4.2 Error Analysis ..... CONVERGENCE A ALYSIS. 28 28 3 1 38 5.1 MaxNorm Bound for Matrix c 1 = (A ~Br1 38 5.2 Estimates for A1 l . . . . . . . . . 4 0 5.3 MaxNorm Bound for Matrix c 1 D . . . . 4 6 5 3 1 Max No rm for Matrix c 1 D in IDimensional Space 47 5 3 2 MaxNorm for Matrix c 10 in 2and 3Dimensional Spaces 5 2 5 3 3 Numerical Results for Spectral Radius of c1 D . . . 5 7 V
PAGE 6
6 UMERJCAL EXPERIMENTS 6 1 Experiment 1 6 2 Exper im ent 2 6 3 Experiment 3 6.4 Experiment 4 7 CONCLUSION APPE DIX ... . A APPROXIMATION FOR CONDUCTIVITY FUNCTION A I Approximation from Figur e 9.6 [27) A.2 Approximation from Figure 3 3 B FISHPAK REFERENCES BIOGRAPHICAL SKETCH V l 5 9 61 62 66 69 74 75 75 75 75 79 80 83
PAGE 7
LIST OF TABLES Tab l e 3 1 Coefficients of the polynomial p( z) . . 25 5 1 Spectral radius for c 1 D in one dimensional space 57 5 2 Spectral radius for c 1 D in two dimensional space 58 5 3 Spectral radius for c 1 D in three dimensional space 58 6 1 Breakdown time comparison for diff e rent values of w 62 6 2 R e lativ e e rror in the potenti a l at t = 65 seconds usin g t h e maximum norm for Crankicholson schem e ( = 1 / 2). . . . . 67 6 3 Relativ e e rror in the potential at t = 65 seconds using the maximum norm for Backward Euler scheme ( = 1). . . . . 67 6 4 Estimates for the c onstant C generated by Equation (6.7) 71 6 5 Estimates for constants p and C generated by Equation (6.10) 73 A 1 Four points chosen from Figure 9 6 [27] 76 A 2 Thirteen points chosen from Figure 3 3 77 Vll
PAGE 8
LIST OF FIGURES Figure 3 1 Path P iPj connects the node Pi to the node Pj 19 3 2 Directed graph, G(M) for a tridiagonal matrix M with nonzero en tries on its diagonal superdiagonal and subdiagonal. . . . 21 3 3 Electrical conductivity CY and corresponding relaxation time T = ECY1 where c = 8.85 12 F m 1 versus altitude under a variety of geo physical conditions LL low latitude, wavy ; MLPS midlatitude typical night (highlatitude quite) ; AZTD auroral zone typical disturbed night ; MLD midlatitude day quite; MHL midhigh l atitude, typical of,....., 100 measurements ; REP, relativistic electron (energy from a few MeV to 10 M e V) precipitation event (unusual); PCA, polar cap absorption event ( an unusually lar ge flux of ener getic,,....., 100 MeV) so l ar protons within the polar cap). Adapted from Hale [8]. (So urce: reprinted with permission from Rakov and Uman [21], published in 2003, copyright Vladimir A Rakov and Martin A Uman publish e d by Cambridge University Pres ) 24 sin2 x 4 1 Graph of J( x) = ( )2 whe n a= 1.05. a+ cosx 30 sin2 x 4 2 Graph of f(x ) = ( )2 wh e n a= 1. a+ cosx 31 6 1 Cylindrica l she ll of the volume integral / dV . . . . . 60 6 2 Potentia l s ((0 0 z ) 0 z 100 km) in two domains f11 and f12 64 6 3 Potentials ( (0 0, z ) 4.6 km z 5.4 km) near the n egat iv e cha r ge (5 km) in two domains f11 and f12 . . . . . . . . 64 6 4 Potent i als ((0, 0 z ) 9 6 km ~ z 10.4 km) near the positiv e cha rge (10 km) in two domains f11 and f12 . . . . . . . . 65 6 5 Potentials (cp(x, 0 5) and cp(x, 0 10) 25 km x 25 km) along the xaxis at the altitude 5 km and 10 km . . . . . . . . 65 6 6 Maximum norm of pote ntial versus t i me where b.t = 1/8 s and = 0 68 Vlll
PAGE 9
A 1 Electric con du ct ivit y as a function of a ltitud e approx imated from curve B in Figure 9.6 [27] . . . . . . . . . . . . . 76 A 2 Electric conductivi ty as a function of altitude approximated from Fig ure 3 3 . . . . . . . . . . . . . . . . 78 l X
PAGE 10
KEY TO ABBREVIATIONS AND SYMBOLS The list s hown b e low gives a bri e f d esc ription of the major mathematical symbols d efi n e d in this diss e rtation. E e lectric fie ld J current d e nsit y CJ conductivit y c permittivity potential /Z d c c (l e7(Clt) ) CJz (J D {( x, Y z ) : 0 S z S Dz, lxl S Dx, IYI S D y } wh ere Dx, D y and D z are some prescribed va lu es. D1 {( x, y,z): lxl S 5 0 km, I Y I S 5 0 km, 0 S z S 100 km} D 2 {( x, y,z) : lxl S 25 km, IYI S 25 km, 0 S z S 100 km} X
PAGE 11
( 2 1 ) Ai = 1 2 1 . m x m m x m AJmxm ( A1 + I A I 2 I ) A1 + I I m2 x m2 ( B1 ) B2 = B1 m2 x m 2 ( A1 ) A2 = A1 m 2xm2 ( A2 + I A = _!_ I h 2 X I
PAGE 12
h C = AB 2 h D = AA+ (1) B 2 XU
PAGE 13
Abstract of Diss ertation Presente d to the Graduate School of the U niv e rsity of Florida in P a rtial Fulfi llm ent of the R e quir e m ents for the D egree of D octor of Philosoph y MULTISCALE DISCRETIZ AT I O OF ELECTRICF IELD EQUATIONS By ShuJen Hu ang Ma y 2005 Chair: William W. Hage r Major D epartme nt: Mathematics W e d eve lop stabl e fini te diff e r e nce approximations for equat i ons d escr ibi ng the e l ect ri c fie ld b etween the surface of the E arth an d t h e ionosph ere. In t hi s sc h eme w e introduc e a paramete r in the t im e discretization For e xample, = 0 and = l g ive us forward and backward difference methods respectively, and = 1/2 is the Crank Nicholson sc heme. The approximat i ons a r e uncondi t ion a ll y stabl e when 1/2 ::; ::; 1 and unstable whe n O ::; < 1/2. Maxnorm e rror esti m ates for the discr ete a ppro ximat ions are esta blish e d. We s how t hat t h e erro r i s O(~t)2 + O(h2 ) if = 1/2, and O(~t) + O(h2 ) othe rwi se W e conduct num e rical ex p eri m e nts to verify our analysis.
PAGE 14
CHAPTER 1 I TRODUCTIO T Most p eop l e are a fraid o f lig h t nin g because i t can destroy buildings and numerous systems critical to daily life. Moreover, it can kill peopl e Lo ss of U S property by li ghtning damage is a huge expense each year. The physical process es involved in li ghtning are the focus of inten sive r esea r c h As the particles, such as graup e l and i ce within a cloud grow in siz e and in crease in numb er some become cha r ged through colli sions. Small e r particles t e nd to ac quir e positiv e c h arge w hil e l a rger particles acq uir e more n e gativ e c h a r ge. T h ese particles tend to se p arate und e r the influ e nc e of updrafts a nd grav i ty until the upp er portion of the cloud obtains a net positive charge and the l ower portion of the cloud b ec om es n egat ively charged. Thi s separation of charge produces hug e e l ect rical pote ntial diff e r e n ce w i thin t h e cloud as well as between the cloud and t h e grou nd. As a result of these e l ect ri c field s a flas h occurs moving charge between po s i t iv e and n egat iv e r eg i ons of a thunderstorm. Uman [26] gav e a detaile d history of ea rly lig htning research. In t h e mid e ighte enth century, B e njamin Franklin perform e d the first scientific study of li ghtning. He d es ign e d an experiment that conclusive l y proved the e lectri c nature of lightning Little s i g nifican t progr ess was m a d e i n understanding t h e properties of lightning un t il the l ate nin etee nth century when photographic a nd spectrosco pi c tool s b ecame availabl e for li g htnin g researc h L i g htnin g cur rent measurements were first proposed b y Po c k e l s [18, 19 20]. H e a n a l yze d t h e magnetic fie ld indu ced b y lig htnin g currents to est imate t h e amount of curre nt. Wil so n [31 32] was t h e first researcher to u se the e l ectric fie ld m easure m e nts to est imate the structure of thunderstorm charges involv e d in li ghtn in g dischar ges H e won the obe l Prize 1
PAGE 15
2 for inventing the Cloud Chamber and made a major contribution to our present understandin g of lightning. Lightning r esearc h continu e d at a steady pace un t il the late 1960 s whe n the r esea rch b ecame part i c ularl y active. This in creased inter est was motivated by The danger of lightning to aircraft or spacecraft, Solidstate electronics used in c omput e rs and other devices Improved measurement and obs e rvation a l cap a bili t ies due to advanc in g tec hnology. Most li ghtnin g research i s done by physicists c h e mi sts, meteoro l ogists and e l ectr i cal engineers Hager [5, 6 7] was the first mathematician u s in g Maxwell' s eq uati o n s to develop a threedim ens i ona l mathematical e l ectr ical model to sim ulate a lightnin g di sc harge. His dis c harg e mod e l [6] was obtaine d by di sc retizing Maxwell s e quations to obtain a r e lation betwee n the pote ntial field a nd current d e nsit y du e to the mo t ion of c harg e d particl es moving und er the influ ence of w ind Spatial derivatives i n his equat ion were approximated by using volum e e l ements in space, whil e t h e tempora l derivatives were est imated b y a backward Euler scheme in time. Since co ndu ctivity is v e ry l arge in t h e r eg ion where the e l ectri c field reaches the breakdown threshold h e eva lu ate d the solutio n limi t as the cond u ctivity tends to infinity in the breakdown re g ion. In hi s mod e l [ 6], the output was the e l ectric field as a function of t im e and the inputs were currents generated by the flow of charged particles w ithin the t hund e r cloud und e r the influ ence of wind. This dissertation i s based on H ager s mathemati cal model. Some improvements are made compared to H age r s ea rli er work. For examp l e The equations are discretized in a diff erent way a ll owing us to so lve them by fast fourier transfo rm (FFT) rat h e r t han by C holesk y factorization. Co n seq uentl y finer meshes can be em plo yed and lar ger domain can be modeled. For examp l e, numerical results reported by Hager et al. [7], the top
PAGE 16
3 of the model domain was 48 kilometers, we measure the top of the domain at 100 kilometers (ionosphere). The huge increase in conductivity with altitude is taken in account in our scheme.
PAGE 17
CHAPTER 2 LITERATURE REVIEW Numerous studies in lightning from different aspects hav e been r eported in the past few decades This review is focused on numerical models of lightning discharges. 2.1 Dynamic Cl oud Models Atmospheric scientists Helsdon and Farley [10] and Ziegler et al. [33, 34, 35], considered the e l ectrom icroph ys ical proc esses for the charge exchange between interacting particles within the cloud domain to obtain dynamic cloud mod e ls. They divid ed water and i ce hydrom eteors into five categories: cloud wat e r cloud i ce, rain snow and highd ens ity precipitating ice (graupel/hail). Each of the five categories of hydrom eteo r was a ll owed to have a charge associated with it. This charge was then transported with the species of hydrometeor according to the transport equation =V'vQ+(UtQPa)+'vKm'vQ ) aQ 1 a (c5Q) 8t Pa 8 z c5t int er (2.1) where Q was the charge d ens ity (in cou l ombs per cubic meter) carried by hydrom eteor category (negative or positive) V was the velocity vector, Pa was the air d e nsity Ut was the massweighted mean termina l fall speed of the hydrometeor class Km was the nonlinear e ddy coefficient and (c5Q/c5t)inter was t h e charge ex cha n ged between categories of hydromet eo r due to interactions. The cloud model o f Helsdon and Farley [10] was twodimensiona l ; whi l e the models of Ziegler Straka, MacGorman, Dye, and Ray [33, 34, 35] were established in three dimensional space The first term on the righthand side of Equation (2. 1) was the advection term. The second term was the fallout term and it was equal to zero for cloudwater 4
PAGE 18
and cloudice charge The third term was the e ddymixing term, and the last term r e presented the charge exchanged between any two classes of hydrometeor during an interaction b e tw ee n those two classes. H e lsd on and Farl ey's mod e l [10] a l so considered the presence of small ion s. The growth and d ecay e quation govern in g the numb e r concentration of small ion s was 5 where n1 2 was t h e concentration of s mall ions with s ubs cript 1 r eprese n t ing positiv e ions and subscript 2 r eprese nting n egat iv e i o ns was the smalli o n mobil ity (in meters sq uar ed per volt p e r second pressure dependen t ) E was t h e two d imensional e l ect ric fie ld vector G was the ionge n e ration rate du e to cosm i c ray ionization (in ions per cubic meter per second, dep e ndent on height ) a was the ionion r eco mbina t ion coeffic ient (1.6 x 10 12 m3 /s), SRC was t h e so urc e term for small ions from processes other than cosmic r ay ionization and SI K was the si nk term for small i ons from pro cesses other t han recombination Zie g l e r e t al. [34 ] ex cluded t h e contr ibu tio n from free ions in their model because the y tho u g h t thi s co ntribution was very small within the cl o ud With the var iou s types of charges accounted for in the model at eac h grid point a total c h arge d e n s ity can be expresse d as where p was the total c h a r g e d ensity ( in co ul ombs p e r cub ic m e t e r ) a nd e w as the e l ectron i c charge (1.6 x 1019 C) The first term repres ente d t h e net c h a rg e du e to the diff ere n ce in s m a llion concentrations and it w as eq u a l to zero in the model of Zieg l e r et a l. [34] b eca use t h ey i g nor e d the contribution from free i o n s. T h e second term r epresented t h e net c h arge on the five categor ies of hydrometeor.
PAGE 19
With t h e total charge det e rmined they a pplied Poisson 's equat ion v'2 = [!__ c to obtain the scalar e l ect ri c pote nti a l in volts where c was the e l ect rical permit tivity of a ir (8. 8592 x 1012 C2 / ( m2)). Then the e l ectri cfie ld vector can b e so lv e d from the pote ntial using Gauss s law E = 'vcp. 2.2 Approac h es for Lowering Charge due to Lightning There are differ ent approaches for lower in g charge due to lightning. The n e utralization of charge b y li ghtn in g in t h e models with bulk flash es an d with exp li c i t lightnin g channels are discussed in Section 2.2.1 and Section 2 2 .2, resp ect i vely The approach u se d in the model of H ager e t al. i s discussed in S ect ion 2.2.3. 2.2.1 Ligh t ning Model wit h Bulk Flashes 6 Rawlin s [22] was t h e fir st scientist to add a s impl e li ghtnin g parameterization to a num er ic a l cloud model. When the e l ectr ic fie l d r eached the di sc harg e thres hold of 500 kV /m everywhe r e within the mod e l domain he s impl y reduc e d the net c harg e b y a constant fraction every wh e r e (arbitrari l y to 30% of t h e pre st r ik e va lu es ) This was the s implest r eprese ntation of the n e u tralization of charge by lightning On t h e other hand Takahas hi [25] s u ggeste d a mor e comp l ex parameterization of a li g htnin g flas h for hi s axisym metric, twodim e n s ion al mod el. His parameteri zation neutraliz e d 20 Co ul ombs of c h a r ge by removing that amount in eac h of t h e two oppositely c h a r ge d r eg ions whe n e ver t h e maximum e l ect ri c field had t h e va lu e of 340 kV /m. To neutralize the charge first h e found the maximum positive charge d e nsity and the m axi mum n egat iv e charge d e n s ity Iext, h e d etermi n ed the small est region abo u t each of t h e two maxima t hat conta i ned 20 Co ulombs F in ally h e
PAGE 20
subtracted positive charge equally from all charged particles in the positive region ; and similarly subtracted negative charge equally from all charged particles in the negative region. Ziegler and MacGorman [33] also applied a simple parameteriza t ion but used a threedimensional model to treat the bulk effect of lightning on storm charge in 7 a time step. They neutralized a fraction of the net c harge at all grid poin ts where the magnitude of the net charge density exceeded 0.5 nC/m3 to imitate the e ffect of lightning releasing the charge from its channels. The above bulkflash schemes are relatively easy to implement but lack realism. 2.2.2 Lightning Model with Explicit Lightning Channels Most lightning discharge models with explicit channel treat only a one dimensional or unbranched channel. The branching of lightning channels would be very difficult to treat analytically because of the lack of symmetry and the incomplete understanding of the processes. 2.2.2.1 The Model of Helsdon et al. Helsdon et al. [9, 10, 11, 12] estimated both the geometry and charge distri bution of an intercloud lightning flash in a twodimensional Storm Electrification Model (SEM). Since then, their parameterization has been extended to a three dimensional numeric a l cloud model. They applied t he concept of bidirectional leader [13]: The parameteriz e d lightning propagated bidirectionally ( initi ally parallel and antiparallel to the e l ect ric fie ld) from the point of initial breakdown. Their parameterization produced a s in g le, unbranched channel that traced the electricfield lin e from an initial point. Initiation propagation direction, and termination of the discharge were computed using the magnit ud e and dir ect ion of the electric fie ld vector as the dete rmining criteria. The charge redistribution associated with lig h tning was approximated by assum in g that the channe l re main ed electrically neutral over its total l ength.
PAGE 21
An initi a l criterion : a thresho ld of e lectri c field w as chosen with a va lu e o f 400 kV /m. The channe l was extende d in both directi ons a long the fie ld lin e un til the ambient electricfield magnitude fell b e low a certain threshold ( 1 5 0 kV /m) at the locations of the channeltermination points. They assumed that t h e linear charge density at a grid point, P along the channel was proportional to the differenc e between the potential at the point where the discharge originated, and the potential at P The linear charge density can b e presented b y Qp = k( p o), where Q p was the charge d ensity at P and p and 0 were the potentials at P and the initiation point of the discharg e r espect ively. The value of this proportionality co nstant k controlled the amount of charge transferred by t h e discharg e Charge distribution at eac h end of the channe l was adjusted in order to maintain net charge n eutrality of the flas h co n s istent with t h e theor y proposed by Kasemir [13]. They extende d the channel b y four grid points at each end, and adjusted t h e distribution of linear charge d ensity on t h e extensions to force the integral of lightning charge density for the whole channe l to b e zero. In this extended r eg ion they assumed that the charge d e nsity d ec r ease d like e0x2 where x i s the distance from the c hannel. 2.2.2 2 MacGorman' s Model 8 MacGorma n et a l. [15] s u ggested a li g htning parameterization that was conside r e d an extensio n of t h e H e lsdon et al. [12] parameterization in conjunctio n with some of the bulklightning parameterization methods presented b y Ziegler and MacGorman [33]. MacGorman et al. [15] d e veloped a para meterization to e n a bl e cloud mode l s to simulate t h e location and structure of individual li ghtning flas h es b y u sing the conceptual model of MacGorman et al. [14] and Williams et a l. [ 30]. T h eir parameterization [ 15] of li ghtning flas h structure procee d e d in two stages.
PAGE 22
9 First, H elsdon et al. [12] provided the initial development. From initiat ion a t a grid point, a flas h traced the e lectricfie ld lin e outward in both parallel and antiparallel directions until t h e magnitude of the ambient e l ectric fiel d at each end fell below some certain threshold value Secondly if one e nd of the channe l reached ground the parameterization terminated at that end but allowed the other end to continue developing. Charge estimati on and neutralization were paramete rized by appl y ing t h e technique proposed by Zieg l e r a nd MacGorman [33], except that Ziegler and MacGorman neutralized c harg e at all grid points h av in g l p(i,j, k)I p1 (w h ere p( i, j k) was the net charge density at the gr id point ( i, j k) and p1 was t h e minimum IP( i j k) I for a ll gr id points to be invol ved in lightning beyond initial propagation) throughout the storm, b u t their parameterization neutraliz e d charge only at suc h grid point within a sing l e localized flash. 2.2.2 3 Man se ll s Mode l Mansell et a l [16] proposed a li ghtning parameterization derived from the dielectric breakdown model t hat was developed by iemeyer et a l. [17] and Wiesmann and Zeller [29] to s imulate e l ectr i c discharges. They extended the dielectric breakdown model to a threedimensiona l domain to represent more realisti c e lectri c field and thunderstorm d y n amics In their work the stochast i c li g htnin g model (SLM) was an application of the WiesmannZeller model to sim ulate bidirectional discharges in the regions of varying net charge density (e.g., in an e l ectrified thund erstorm). Procedures for s imulating li ghtning flashes in t h e thunderstorm model were as follo ws. A flash occurred whe n t h e magnitude of the e l ectric field exceeded the initiation threshold E i n i t anywhe r e in the model domain. The lightning initiation point was chosen randomly from among a ll the points where the magnitude of the e lectric field is greater than 0.9Einit Both decisions for choosing the initiation thresho ld
PAGE 23
and the initiat ion point were suggested b y MacGorman et al. [15]. Positive and negative parts of the flash were prop agate d indep e ndentl y so that up to two new channel segments (positive and negativ e ) c ould b e a dded at eac h step. Bo t h ends had d efa ult initi a l propagation thresholds of 0 7 5Einit For flash n e utrality t h ey assume d that t h e chan n e l stru ct ur e would maintain overa ll charge neutrality as lon g as neither e nd reach e d the g r o und as suggested by Kasemir [13]. But, 10 for computational s impli city, their parameterization maintained nearneu tral ity (within 5 % ) by a technique of adjusting the referenc e pote n t ial to the growth of t h e lightning structure in s tead of adjusting the reference potenti a l of t h e c h annel. 2.2.3 Hag e r 's Model Hager et al. [5, 6, 7] proposed a t hr eed im ensiona l li ghtningdischarge model that produced bidirectional ax is y mmetri c I C an d CG flashes. The mod e l gen erated the di sc harg e region charge transfer a nd detailed cha r ge rearrangement associated with the flash Their approach to ligh t ning was quit e different from those in Section 2.2.1 and Section 2 2.2. Their breakdown model was based on Maxwell' s eq uati ons. They ass um ed t hat curre n t due to trans port o f c h arge und er the influence of wind was known T hey obtained an equat i o n gove rnin g the evo lu t ion of th e l ect ric potential und e r t h e ass ump t i on that t h e time d er i vative of t he magn e t i c fie l d can b e di s r egar d ed. After int eg rating this equation over boxes an d approximating derivatives by finite differences they obtained a n impli c i t system of diff erence equat ion s describing t h e evo lution of t h e e l ectr i c field. Their approac h to lightnin g was to let the conductivit y tend to infinit y wherever t h e e l ectr i c field reached the breakdown threshold. This ap pro ac h appeals to our basic conception of nature: Wh e n the e l ectr i c fie ld reaches breakdown threshold con du ctiv ity becomes very l arge as a p l asma forms.
PAGE 24
11 When the electric field reaches the breakdown threshold the electric potential changes instantaneously everywhere within the thundercloud The Inverse Matrix Modification Formula [4] was applied to evaluate this change: (2.2 ) where before was t h e electric potential b e fore discharge
PAGE 25
CHAPTER 3 THE MO D EL Hager et al. [6] were the first to use Maxwe ll 's equations to dev e lop a three dimens ional e lectri c mode l for a t h understorm. We first establish the mode l based on Maxwe ll s equations. 3 1 Formu lation of the Equations According to Hager et al. [6], the curl of the magnetic fie ld stre ngth H i s g iv e n bas e d on Maxw e ll 's equat ions: BE 'v x H = cat+ aE + J (3.1) where c is the p ermitti v ity, ai s the cond u ct i v ity, E is t h e e lectri c fie ld and J is the current densit y associated with the influenc e of the wind. In t hi s model we hav e the following assumptions: Assumption 1: The time d e rivativ e of the mag neti c fie ld i s zero i .e., 8 B = O => at 8 B 'v x E = = 0 at => E = 'v Therefore E is t h e gradient of a pote ntial Assump t ion 2: Let E8 b e the breakdown fie ld strength. T h e n t h e e l ectri c fie ld m ag ni t ud e i s alway s l ess t h a n or equa l to t h e break d own thrc hold ER. That i s I E I :S Ea. Assumption 3 : Wh e n the e l ectric field r eac h es the breakdown thresho ld E8 at some poin t th e conductivity tends to infinit y in a s m all n e i g hborhood at that point. That i s if IE(x, y,z)I = E a then o(x,y,z) = + 1 2
PAGE 26
13 In this di ssertat ion we do not impl e ment the lightning discharge described in Section 2.2 .3. Inste ad we focus on integratin g the equations up to the discharge initiation We take the div e rgenc e of Equation (3.1) to get 8E c'v at+ 'v oE + 'v J = 0 (3.2) Boundary conditions: In this model we consider a bound e d domain n b e twe e n the surface of the earth and ionosphere L e t us consider boundary conditions on an:
PAGE 28
The solution of Equation (3.6) is where sup e rscript n denot e s the time l e vel t n 1 = CJ / c, a nd 'Yz = CJz /E: i s t h e d e riva t ive of CJ/ E: with r espect to z Fir s t w e approxim ate t h e int eg r a l p a r t o f Equation (3.7), denoted as I: I= e'Y(Stn+1) ,z+ ds. lt n+l ( O'lp 'v. J) tn fJz E: The following steps are conducted. Step 1: Both fJ'ljJ/fJz and 'v J are evaluated at time (t n + tn+i) which is d enote d as time level n + H e nc e w e h a v e the followin g approximation Define 1 ey(.6.t) d= ,z 1 Then we can r e write the above approximation as follows : (fJ'ljJ)n+ 1 ey(.6.t) ('v. J r+ I~ d 0 + . z 'Y E: Step 2: Use (1 )(fJ'ljJ/oz r + (fJ'ljJ/oz r +1 to a pproximat e the t erm (fJ'ljJ/oz r+. H e nce, [ (fJ'ljJ) n (fJ'ljJ)n+ll 1ey(.6.t)('vJr+ I~d (1) + + , fJz fJz E: wh e r e O 1. 15
PAGE 29
16 So Equation (3. 7 ) i s t hen approximat e d b y the following e quation 6'lj;n+l = [ ( EJ'lj;) n ( EJ'lj;) n + 1 ] ).6'lj;n d ( 1 ) B z + oz 1 e'Y(L'.t) ( v' J)n+ 'Y E: ( 3.8) wh ere>.= e 'Y(L'.t), 1 = CJ/E:, and O::; u::; 1. ow w e con side r Equation (3.6) but ignor e the force t erm J (3.9) Similarly, t h e so lu t i o n of Equation (3.9) i s (3.1 0 ) a nd i t c a n b e approximat e d b y Next t h e s pati a l p a r t of E q uati o n ( 3 .11) ca n b e di sc retized by t h e finite diff ere nce m e thod. D efine wh e r e m i s a po sitiv e int ege r. The mesh p o in t (xi, Yi, zk) i s d efine d as wh ~ r e O ::; i ::; m + 1 0 ::; j ::; m + 1 a nd O ::; k ::; m + 1. L e t 'l/Ji,j k d e n ote t h e di sc r ete a ppr o xim a ti o n t o 'l/; at t h e p o int ( xi, Y i zk)For d o m a i n f21 : u s i ng t h e
PAGE 30
uniform m es h s ize for x y a nd z ( t hat i s 6x = 6y = 6 z = h ) we h ave t h e approxima t ion 8 . _ ~ i ,j, k + l ~i, j k 1 z'l' i J,k 2 h for the part ial d e rivativ e 8 ~/Bz at the point ( xi, Y J z k) T h e sec ond partial d e rivativ e 82~/8x2 a t the poin t ( xi, YJ, z k ) has an a pproximation d e noted b y O f~i,J,k : ,. 1 k 2 ''. k + ,,. l k 82.,,. _ 'I'+ J, '1',J, 'fl, ,1 1 'l' i J k h 2 Similarly, w e u se a pproximation s a J~i, J k a nd B J~i,J k to t h e sec o n d par t i a l d e riv a t ives 82~/8y2 an d 82~/8z2 : 82., .. _ ~i,j+l, k 2~i,j k + ~i,j1 k 2 'l' i J k h 2 and ,. k 1 2 k + k 1 82.1 _ '1'i,J, + 'l't J '1', 1 3 '1'i,J,k h 2 r es p ect iv e l y W e d efine t h e di sc r e t e L a pl ac i a n o p e rator 6,_h as foll ows : 3 6 h k = a2 1,. k '1'i 1 L 1 '1'i,J l = l Aft e r r e pl ac in g 6 b y 6 h and 8 ~/Bz b y the di sc r ete Bz~ in E qu a ti on ( 3 .11) w e ob t ain t h e followin g finit e diff e r e n ce sc h e me. 17 6 h.1 n+1 >. 6 h 1 n h [(l ) d k a 1 n dk O 1 n + 1 ] 'f'i,j k + k 'f'i, j k 2 h 2 Zl/'i,j,k + h2 Zl/'i,j,k (3.12) Equat i o n (3.12) y i e ld s a matri x v ector form (3.13) wher e \JI d enotes th e vector with c ompon ents ~ i J ,k T h e matri x A i s symmetri c and pos itive d e finit e and t h e structure of t h e m atrix A i s a bl oc k t ridi agona l
PAGE 31
matrix in the following form: 1 A = h 2 18 (3.14) where I m 2 i s t h e m2 x m2 id e nti ty matrix a nd A2 i s a n m2 x m2 b l ock t ridia g onal matrix in t h e form where I m i s a n m x m id e n t ity matrix a nd A1 i s a t ridiagonal matrix with 2 s on i ts diagonal and l' s on its superdiagonal and subdiagonal. The matrix A i s a diagonal matrix with m2 b l ocks on the diagonal and eac h blo c k h as ,\ on its di ago n al. The matrix B i s a di agona l blo c k matrix in the following form: 1 B = h 2 (3. 15)
PAGE 32
19 where the diagonal block B1 is an m x m tridiagonal matrix in the form Varga [28] gives the definition for an irreducible matrix and some theorems to determine the irredu c ibility of a matrix. Definition 3 .2.1, D e finition 3.2.2 The or e m 3 2.3 Definition 3 2 .4, and Theorem 3.2 5 ar e adopted from Varga [28]. Definition 3 .2.1 F o r n 2:: 2 an n x n matri x M is r e ducible i f th e r e exis t s an n x n permutation matrix P such that wh e r e M1 1 i s an r x r submatrix and M2 2 is an (n r) x (n r) s ubmatrix wh e r e 1 :s; r < n. If no s u c h p e rmutation matrix exists, th e n M is irre du cible. A directed graph is an alternative way to check irr e ducibilit y for a matrix. Let M be any n x n matrix, a finite directed graph G(M) can b e construct e d in the following way : Consider any n distinct points (nodes) P1 A , Pn in the plane. For ev e ry nonzero entry mi,j of the matrix, we connect a path direc t ed from the node P i to t he nod e P1 as shown in Figur e 3 1. p p I J Figur e 3 1: Path PJ11 conn e cts the node P i to the nod e P1
PAGE 33
20 Definition 3.2.2 A directed graph is strongly connected if, for any ordered pair of nodes P i and Pi there exists a directed path connecting Pi to Pi. Theorem 3.2.3 Ann x n matrix M is irreduc ible if and only if its directed graph G(M ) is strongly connec ted. Definition 3.2.4 Ann x n matrix M = (mi, j) is diagonally dominant if n l m I> '""lmI i,i 6 t J (3.16) j=l,#i for all 1 i n. Ann x n matrix M is irreducibly diagonally domina nt if M is irreducible and diagonally dominant, with strict inequa l ity in I nequality (3.16 ) for at l east on e i. Theorem 3.2.5 If M is an n x n irreducibly diagonally dominant matrix, then th e matrix M is nonsingular. Now, we rearrange Equat ion (3.13) to y i e ld Let C = A ~ B and D =AA+ ~(l ) B T h en we have w h ere t h e matrix C is a block tridi agona l matrix given b y 1 C = h 2 (3. 17) (3. 18)
PAGE 34
21 where C m 2 i s an m2 x m2 tridiagonal blo c k mat rix with m x m matrices 6 (3.19) on the diagon a l block and with n egat iv e m x m id e ntity matrices on both its supe r diagonal and s ubdi agonal blocks Lemma 3.2. 6 Let M b e any n x n tridiagonal matri x I f all entries o n t h e diag onal superdiagona l and s ubdiagonal of M are nonzero, th e n M is i rr educible 000 00 ... ... > > " " ........ ,.__ii_, _______ ___, ____ p n Figure 3 2: Directed graph, G(M), for a tridi agona l matri x M with nonzero entrie s on its diagonal superd i agona l and subdiagonal. Proof: Since the t ridi ago nal matrix M has nonzero e ntri es on its diagonal supe rdiagon a l and subdiago n a l the directed graph of t h e matrix M G(M), i s s hown in Figure 3 2. From Figure 3 2 it i s easy to see t hat G(M) i s stron g l y connecte d and t h e r efore M i s irreduc ibl e by Theorem 3.2 3 Lemma 3.2.7 The matrix C in (3.18) is irr e ducibl e Proof: First, w e want to show that each d i agona l b l ock Cm2 of the m atrix C i s irreducible By Lemma 3 2.6 eac h di agonal block (3. 19) o f C m 2 is irr educ ible; h ence t h ere exists a path for any two vertices associate d with a diagonal blo c k (3.19). Suppose we are g i ven vert i ces i < j in two distinct diagona l blocks ; the subdi agona l blocks of C m 2 provides a path between i and i + m and between i + m and i + 2m etc. Event u a ll y we reach a vertex k in the block containing j. And
PAGE 35
22 since the vertices k and j are in the same irreducible diagonal block there ex ists a path from k to j. H ence, there exists a path from i to j. Now the matrix C has irreducible diagonal blo cks C m2. The structure of C is similar as the structure of Cm2 Similarly we can show that the matrix C is irreducible. Therefore C = A ~B is irreducibly diagonally dominant by Defini tio n 3.2.4 and Lemma 3.2. 7. Based on Theorem 3 2.5 C is nonsingul ar. Therefore, we hav e the iteration where C = A ~B, D =AA+ ~(1)B, and c10 is the amplification matrix. 3.3 Important Role and Approximation for Conductivity Function The conductivity function varies by many orders of magnitude between the surface of the Earth and the ionosphere Figure ?? and Figure 3 3 show that grows more than the factor 101 0 S/m from the surface of the Earth to the ionosphere which l eads the differential Equation (3.6) to be a very stiff quation. For a stiff equation dy dt = ky, (3.20 ) where k is a very larg e positive numb er, its solution is y(t) = ekty(O) and y(t) tends to zero quickly. Let us consider the following three schemes: 1. Euler e xplicit scheme for Equation (3.20) : Therefore yn = (1 k6tty(O). Whe n k is very larg e, (1 k6t) < 1 and yn + +oo as n incr eases. H e nc e Euler ex plicit scheme do es not work. 2 Euler implicit scheme for Equation (3.20) Yn+l yn ____ = kyn+ l 6t
PAGE 36
Therefore, yn = ( 1+!.6t ) n y(0) and O < 1+!.6t < 1 whe n k i s large H ence yn + 0 as n + oo. The scheme works. 3. The exact formul a: Therefore, yn+ I = ( eMtr y(0). The factor ( ek.6tr tends to zero much faster than ( 1+!.6t ) n when k i s large. This scheme works better than Euler implicit scheme. 23 Therefore, the above explanation demonstrates that the exact formul a is the best method to solve a stiff equation. The conductivity function, a plays an important role in the Equation (3.6). That i s the reason that the exact solution (3.7) is chosen in Section 3 2. At the beginning of the research, we use Figure 9 6 in [27] to find an approximation for the conductivity function a (se e Appendix A .l). For O z 27.828 km: a( z ) = 1O13.40+0.1600z0.005175z2+0.00009077z3 For z > 27.828 km: a( z ) = 1012. s s+ 0 .06718z and Then we have { a z = (0.1600 O O1035 z + O.OOO2723z 2)(ln 10 ) ; = 0 1546 for O :S z '.S 27 828 for z > 2 7 .82 8 { : z ( ~ ) = ( 0.0103 5 + O .0OO5446z)(ln 10) for O :S z '.S 27. 828 (az ) = O for z > 27.828 o z a Figure 9 6 [27] only provide s conductivity d a t a up to 27 .8 km (ab o u t 100 0 00 ft). In order to produce more accurate results, we use Figure 1.3 in the book
PAGE 37
120 110 100 90 BO E 70 J,t'. .g 60 ::, ~50 40 30 20 10 MLPS I I I I I I I Night .. "Nom inal"_.,'ilDay I / /// " ..::= ___ .,,,,.,,./z/'I I = / Nov '69 PCA, day ~ ,,1' ,.,,.~MLD ,., Nov '69 PCA. night 101 1013 10 12 1011 10 1 0 109 108 101 108 10"'5 10~ 103 102 101 10 Conductivity, S m 1 103 102 101 1oO 10 1 10 2 103 10_. 10 5 10 6 101 10a 1cr9 10 10 10 11 Relaxation time, s 21 Figure 3 3 : Electrical conduc t ivity CJ and corresponding relaxation time T = c
PAGE 38
2 5 Lightning: phy sics a nd effect s [21]. In the book [21], Rako v and U m a n g iv e a figur e adapted from Hale [8] that es timates conducti vity at altitude up to 120 km. Figure 1.3 in [21] i s r eproduce d in Figure 3 3 of t hi s t h es i s w i t h p e rmi ss i on o f t h e authors In Appendix A.2 w e g iv e t h e foll owi n g a p prox imat i o n b ase d o n Fi gure 3 3: CJ= 10 r(z) where p(z ) is a quart i c polynomial 4 p(z ) = L aiz i, ( 3.2 1 ) i = O wher e z is in kilometers and a i are co effici e nts of the pol y n omial p(z) as s h ow n i n Tabl e 3 1: Co effic i ents of t h e p o l y n o mi a l p( z ) i a i 0 1.316 X 101 1 6.3 5 1 X 10 2 2 8 682 X 10 4 3 8 6 7 4 X 10 6 4 5.781 X 10 8 T a bl e 31. The n w e hav e and CJz = ( 0 063 5 1 + 0.001736 z 0.00002602 z2 + 0 0000002312 z3)( ln 10 ) (J : z ( : z ) = (0.001 7 36 0.0000 5 20 4 z + 0.0000006936 z2)(l n 10 ) For the a ppro x imation s estimated from F i gure 3 3 w e c a n s how t hat a n d % z ( ~ ) are b ounde d in the d omain D1 ( d e fin e d in ( 3 .5)) as foll o w s : I ~ I ~ o 4793
PAGE 39
and I : z (:z ) I 0.007987 3.4 Disc u ssion: d = ]z(l e 1 ( 6. t ))h In this sect i on we discuss t h e parameter d, where 26 d = (l e1(6.t))l, / Z (3 .22) 1 =
PAGE 40
27 Lemma 3.4.2 Ford as d efine d in (3.22) and E C2 we have Proof: By the mean value t heor e m we h ave where %z ( 7 ) i i s the derivative of ( 7 ) with resp ect to z at som e zi between Z i a nd Zi+l, ~i and 1i are t h e value a nd the derivative o f at some Z i between zi and zi+l Sin ce D i s bound e d a nd a E C2 both 7 and % z ( 7 ) a r e bound e d and we h ave ldi di+1I < 1(1 e r;(~tl ) : z ( ~ )ihl + 11i(6.t) ( ~ )i+i h i < ,,i(6.t):z ( ~ ) /I+ 11i(6.t) ( : z )i+l hi 0(6.t h)
PAGE 41
CHAPTER 4 STABILITY A D ERROR A ALYSIS 4.1 Stability This section uses Fourier analysis in our study of stability. We use Fourier analysis on the grid of integers, Z or hZ where hZ is defin ed by hZ = { hm : m E Z}. Strikwerda [23] gives the Fourier transform as follows: Definition 4.1.1 If v is a grid function in 2 d efined for all integers m its Fouri er transform is given by 00 v(~) = _1_ L eim{Vm V27r m=oo for~ E [1r,1r], andv(1r) = v(1r). If th e spacing between the grid points is h, by changing variables w e hav e for E [1r / h 1r / h]. ext, we consider 'l/; a gri d function whose spacing between the grid points is h Then we have 1. 28
PAGE 42
2. 1 v'21r L eimh0('1/Jm+l 'l/Jm1)h m=oo 00 Theorem 4.1.2 Th e finit e diff e r ence sc h e me, Equation (3 .1 2), is unconditionally stab l e wh e n 1/2 S S 1 and is unstab l e wh e n O s < 1/2 Proof: Take Fourier t r ansform of Eq uati o n (3.12) to obtain 3 L (2cos0j 2),J;n+l j=l 3 >.k L(2cos0 j 2),J;n dkh(l )isin03,J;n dkhisin 0 3,J;n+l_ j = l 2>.k E(cos0 j 1)dkh(l )isin03,J;n_ 2I::(cos0j 1) + d k hisin03 Define t h e a mplification factor as 2>.k L ( cos 0 j 1) d kh(l )i s in 0 3 p= 2I::(cos0jl)+ dkhisin03 T h e n the magnitud e of p i s g iv en b y I 12 = 4>.% ( E (cos01 1))2 + d%h2(1)2sin203 p 4(I::(cos011))2+d%h22sin203 29 For 1/2 S S 1 we h ave IPl2 S 1 since >.k = er,.,(~t) s 1. H e n ce t h e sche m e i s un con dition a ll y stable. For O S < 1/2, w e want to s how I P l2 S 1 + J<(f:).t ) where K i s a constant, does not h o ld for a ny 0 j E [ 7!", 7!"]. If IPl2 s 1 + J<(f:).t), then we h av e
PAGE 43
30 4>.% (I)cos 0 i 1)) 2 + d%h2 ( 1 )2 sin 2 03 < 4 ( I)cos 0j 1) r + 4K(6.t) (I)cos 0 j 1 ) r +d%h2 2 s in 2 03 + K(6.t)d % h2 s in2 0 3 4(>.% 1 K(6.t)) (I)cos 0i 1)) 2 :S d % h2( K ( 6.t)2 + 2 1) sin2 03 4(>.% 1 K(6.t)) :S d%h2(K(6.t) + 2 1) sin 2 0 3 2 ( L (cos0i 1)) s i n 2 03 4(1 + K(6.t)) >.%) 2:: dkh2 ( 1 2 K(6.t) 2 ) 2 (L(cos0i 1 )) 3 2 1 (\ \ \ 1 2 3 X sin2 x Figure 4 1 : Graph of J(x) = ( )2' wh e n a = 1.05. a + cos x ow consider the function where a = cos 01 + cos 02 3 a n d 5 :S a :S 1. Let f x = sin 2 x ( ) ( a + COS X) 2
PAGE 44
31 where 1r s x ::; 7f, a nd 5 s a ::; 1. If a =fa 1 t h e graph of f (x) i s shown in Figure 4 1. We can find the maximum of f(x ) i s a/(a2 1)312 When a te nds to 1, the maximum of J(x) approaches to infinity. If a= 1, the graph of J(x) is shown in Figure 4 2. We can see the supremum of J(x) is infinity. Henc e, 300250 20 Y15 10 3 2 1 2 3 X sin2 x Figure 4 2: Graph of J(x) = ( )2, when a= 1. a+ cosx sup (""'( ))2 = +oo. 1r/h'.5.03'.5.1r/h 6 cos 0 j 1 ( 4.1) Therefore it i s impossible to hav e the inequality 4(1 + K(t::. t)) A~) d~h2(12 K(t::.t)2 ) sin2 03 2 (~(cos 01 1)) In other words the sc heme is unstable when O s < 1/2. 4.2 Error Analysis In this section t h e error generated from the discretizing process i s anal yzed. Lemma 3.4.1 and Lemma 4.2.1 are applied in the error anal ys is. L emma 3.4.1 g i ves
PAGE 45
the estimate for d and Hager [6] ment i ons the e rror estimat e for int e gration as described in Lemma 4 2.1. 32 Lemma 4.2.1 L e t Wk 00([a, bl) d e note the space of functions d e fin e d on the int erval [a, b] with d e rivat i ves through orde r k e s s e ntially bounde d If g E W1 00 and f E W2 00, the n we hav e the e stimat e 1 1 b f(z)g(z)dz f(c) 1 b g(z)dzl (b ~t)3 llf'llllg'I I + ( b ; 4a)3 llf" l lllgll, whe r e II II d e notes the essential s upr e m 'um and c is t h e midpoint of t h e interva l [a, b]. Proof: wh e r e 1 1 b f(z)g(z)dz f(c ) 1b g(z)dzl 11b (J( z ) f(c)) g(z)dxl 1 1 b (lz J'( s)ds ) g(z)dz l 1 1 b (lz [J'(c) + 1 s J"(r)d r ] d s ) g(z)dzl < 1 1 b (1z f ( c )ds ) g(z)dz l + 11b (1z 1 s J"(r)d r d s ) g(z) dzl I 1 1 1 b (1z J' ( c)ds ) (g(c ) + l z g'(s)ds ) dzl < 1 1b(z c )J'( c)g(c)dz l + 1 1b(z c)J'(c) (1z g'(s)d s ) d z l < l l f'll llg'll 11b ( z c)2dz l ( b ~t)3 l l f'll llg'II,
PAGE 46
and T hi s proves the lemma. ote that i s considered small in Lemma 4.2.2 Lemma 4 2.3 and Theo rem 4.2 .4. 33 Lemma 4.2.2 If w e use th e finit e diffe renc e schem e (3 .12} to di scretize th e int e gral e quation (3.10 ) th e n the truncation e rror is wh e r e m = 3 if= 1/2 and m = 2 if =I= 1/2, Proof: Defin e R1 as follows: l tn+ I a'lj; ( a'lp)n+ l tn+I e'Y(Stn+1),zdS /Z e'Y(Stn+1lds t,. a z a z tn l tn+ I B'lj; 1 e1'(6.t) ( B'lj;) n+ e'Y(Stn+l)fzdS ,z t n a z 'Y a z l tn+I e'Y(Stn+il,z B'lj; ds d ( B'lj;) n+ t n a z a z wher e d = 'Yz(l e 'Y{b.tl)/,Apply Lemma 4.2 1 to o btain t h e follow i ng estimate: < C(Llt)3 (4.2) D efine R2 as follows : (a'lj;)n+ [ ( a'lj;) n ( a'lj;) n +1 ] R 2 = d a z d ( 1 ) a z + a z
PAGE 47
34 where O ::; ::; 1. By Tay l or s expansion, w e hav e 8'lj; f)z (tn) = 8'lj; ( t::.t) (f)'lj;) 1 ( t::.t)2 (f)'lj;) (t .! ) + (t .! ) + (t .! ) + ... f)z n+2 2 f)z n+2 2 2 f)z n+2 t tt and f)'lj; f)z (tn+ 1 ) = f)'lj; ( t::.t) ( f)'lj;) 1 ( t::.t) 2 ( f)'lj;) (t .! ) + (t .! ) + ( t .! ) + .... f)z n+2 2 f)z n+2 2 2 f)z n+2 t u Hence This implies that if = 1/2, (f)'lj;) n (f)'lj;) n+l (f)'lj;)n+ 2 (1 ) f)z + f)z f)z = O(t::.t) and if i= 1/2, ( f)'lj; ) n ( f)'lj;) n + 1 ( f)'lj; ) n + (1 ) f)z + f)z f)z = O(t::.t). T h erefore (4.3) whe r e m = 3 if= 1 / 2 a nd m = 2 if i= 1 /2. Let R3 b e d efined as follows: R = ,6 J n+l )..,6oJ n (6h01 n + I ).. 6h01,n. ) 3 'I' 'I' '1't,J,k k '1'i,J, k
PAGE 48
By Tay l or s e xpansion w e have Combining (4.2) (4. 3) and (4.4) w e s ee that when t h e so lution t o ( 3.10 ) is substituted in (3.12) we obtain where m = 3 if= 1/2 and m = 2 if i= 1/2. 35 ( 4.4) ow if w e add t h e for ce t erm J b ac k to Equatio n ( 3.10 ) t h e n the corre spo nd ing finite d i ff e rence sc h e m e is A k 6h'l/J~j, k [(1 ) ~ ~ B z'l/J~j, k + ~ ~ Bz'l/Jftf] 1 e ,k(Ll.t) (v' J ) n + 1 k E ( 4 5) Lemma 4.2.3 If w e use th e fin i t e diff e r ence sc h e m e (4. 5 ) to d is cr e t ize Equ ation (3. 1}, th e n i t has th e t runcation e rror whe r e m = 3 if= 1/2 and m = 2 if:/= 1/2. Proof: W e do n o t includ e t h e for ce t e rm J in L emma 4.2 .2. o w w e on l y need to foc u s on t h e e rror caus e d by estimating the J term wh e n t h e fini te diff e r e nce sch e m e ( 4 5) i s appl i e d t o E qu atio n ( 3.7) In t h e s c h e m e (4.5) w e ap proxi m a t<' t h<' v' J t erm in (3.7 ) by t h e v' J term in ( 4.5 ) B y L e mm a 4 2 1, w e h a ve 1tn+ 1 y' J (v' J )n+ 1tn + 1 e,
PAGE 49
36 Combin i ng thi s est imate with L emma 4 .2.2 w e find that when the o luti o n to (3. 7) i s substituted in ( 4 .5), we obtain where m = 3 if = 1/2, and m = 2 if =/. 1/2. The finite difference scheme ( 4 .5) y i e ld s a matrixvector re lation of t h e form h 1 e1'(t,.t) 1 Awn+i AAwn = [( 1 )Bwn + Bwn +1 ] I n (4.6) 2 'Y E where I n correspo nd s to the term V J of ( 4. 5) e v a lu ate d at t im e l eve l n + The next theorem gives u s the maxnorm error estimate wh e n the solutio n to (3. 7) is substituted into Equation ( 4.6). D efine t h e maxnorm for a vector x as where x = ( x1 x2 . Xnf, a nd d efine the maxnorm for a matri x M as n IIMlloo = l !\~ L l mi,jl j = l wh e r e M = ( mi,j) i s a n n x n matrix Theorem 4.2.4 L et Er denot e th e continuous solution of (3. 7) evalu ated at t = ni0..t and at the m esh points of th e grid { (xi, y1 zk)}. L et q,n b e the vector with components '1/Ji,j, k in Equation (4. 6). Then we have w h e r e m = 2 if= 1 / 2 and m = 1 if=/. 1/2. Proof: In sert e n into Eq uation ( 4.6) and w e get
PAGE 50
37 where Tn is the truncation e rror and where m = 3 if= 1/2, and m = 2 if =I1/2. Then by subtracting Equation (4.7) from Equation ( 4.6) we have (4.8) where C = A ~ B and D =AA+ ~(1 )B. Following from Equation (4.8), it is found n l w n e n = M n ( w0 E>0 ) L M 1c 1rnll, l = O (4.9) where M = c 10 ote that we prove JJMJ J00 = 1 + O(.6.t) in Sec t ion 5.3 for O ::; ::; 1 if I is small. Therefore JJMJl00 can be wr itten as IIMJl00 ::; 1 + c.6.t ::; ec6.t_ Taki ng t h e maxnorm of Equat i on (4.9) then we have ote that IJC 11100 is bounded by a constant which is independent of h (see L emma 5 1.2). Consequent l y if l l w0 800 = O(h2 ) then where m = 2 if= 1/2, and m = 1 if =I1/2.
PAGE 51
CHAPTER 5 CO VERGENC E A ALYSIS 5 1 Maxorm B ound for Matri x c 1 = ( A ~ B r1 In this sect i o n we obtain a m ax n orm boun d for c 1 = ( A / L B ) 1 whe r e A and B a r e d e fin e d i n Equat ion ( 3 1 4 ) and E quati o n (3. 1 5) r especti v e ly. L emma 5 1.1 from Varga s Matrix I t e rat i v e Analys i s [28] is us e d to s h ow that t h e matrix C i s inv ertibl e and the maxnorm of t h e mat rix c 1 i s bounde d b y a con stant whi c h i s ind e p endent o f h Lemma 5 .1.1 If M = (mi,i) is a r e a l irr e d u c i bly diagonally domi nant n x n matr ix wi th mi, i 0 f or i =f. j an d mi,i > 0 f or all I i n, t h en M 1 > 0 i e all e l e m ents of M 1 ar e p o si t ive Lemma 5.1.2 If h is s uffi c i e ntly s mall, th e n th e m at r ix C is invertible Mor e over 11c 11100 is boun d e d b y a c o ns ta n t whic h is i ndepend en t of h. Proof: The proof o f t hi s r es u l t is the sam e proof u se d in [6] for a s li g h t l y diff e r e n t matrix. ote that a ll t h e e l e m ents in the matri x C sati sfy t h e condit i o n s in L emma 5.1.1 except the e l e m ents on t h e subdi agona l of C The e lements o n t h e subdiagonal hav e the form 1 + ~ d w h e r e dis a b ounde d p os i t i ve nu m b e r fro m S ection 3 .4, and O l. In orde r to u se L emma 5 .1.1 w e mus t sati s f y the condition h 1 + d < 0 2 ( 5 .1) According l y i f h i s s uffici entl y small c 1 0 i .e ., t h e e l e m e nts of c 1 a r e p os i t ive b y Lemma 5. 1.1. H e n ce I I C 11 100 = I I C 111100, wh e r e 1 i s the vecto r wi t h eve r y component equa l s to 1. Furthermore IIC11100 11c1F l loo if F 1. 3 8
PAGE 52
39 ote that the matrix C i s the coefficient matrix corresponding to the expres sion D, h.J k dkO o / k 'f/i,J, Z'f/t,J , (5.2) where 6 h is the discrete Laplacian operator O z is the discr ete 8'1/; / o z and d k = Next we show that 11c11100 :::; 2a2 where a = hf, and I i s the maximum va lu e of the ind ex i associated with points in the set Nh. Let us defin e 0i,j, k = (h i)2 Observe that Moreover defin e w i,j,k = (hi)2 for (i j k) outsid e h_ If (i, j, k) E we c h oose wi,j,k to satisfy 6hw k dk w k = 0. i,J, z t,J (5.3) Then from Equation (5.3) we obtain Wi,j k = i [ W i+ l j k + Wi1, j k + Wi, j + I k + Wi,j1, k + ( 1 + %dk) W i ,j k+l + (1 ~dk) W k 1 ] (5.4) 2 i,J, for (i j, k) E Nh_ If his sufficie ntl y sma ll the condit i on (5.1) is satisfied that is 1 ~dk 2: 0. Hence, by Equation (5.4) we h ave where 8Nh is the set of boundary points of N11. H ence,
PAGE 53
for all ( i, j k) E Nh. If 'lj;i,j,k = 0i,j,k + wi,j,k, t h en we have 1Pi,j ,k = 0 t0,.h.1 k dko n l k > l f.fli,J, z l.fli,J, if (i,j, k) E fJN\ if (i, j k) E Nh. 40 Let '11, 8 and n denote the vectors formed from 'lj;i,j,k, 0i,j,k, and wi,j,k for (i,j, k) E Nh. Take F = Cw 2". 1 to ob t ain I I C 1 lloo < IIC1Flloo = 1 1'1111 118 + fllloo < IIBl loo + llfllloo < a2 + a2 = 2a2 where a is a constant which is ind ependent of h. Remark: As noted in t h e co ndition (5. 1 ), t he matrix C will satisfy the conditions of Lemma 5.1.1 if h is suffic i e ntl y small. In fact, using t h e est im ate for d from Section 3. 4 we obtain the condit ion of h '.S 4 .173 kilometers wh e n = l for the domain D1 5 2 Est i mates for A1 l In this section we w ill get some pre li minary resu lts for finding the maxnorm of the amplification matrix c 1 D Lemma 5.2.1 L e t A1 b e th e n x n symme tric tridiagonal matr ix wit h 2's on th e diagonal and with 1 's on th e supe rdiagonal and on th e subdiagonal Th en the ith row of th e inve r se matrix of A1 has the form ( A1 1)i = n: 1 [ni+l ,2(ni+l), ,i(n i+l),i(ni),i(nil), 2 i,i].
PAGE 54
Proof: Note that we compute the (i jlh entry of the matrix A1 i A1 b y t h e formu l a n :~::)A1 1)i, k(A1)k,J k = l (a1)i,j1Ujl, j + (a1 ) i,jaj, j + (a1)i, j+laJ+l,j, where (a 1 )i,j denot e s the ( i,j)th entry of the matrix A1 1 For j = i, we h ave 41 1 ( A1 l A1)i,i = n + 1 [ ( i 1) ( n i + 1) ( 1) + i( n i + 1 ) (2) + i( n i) ( 1)] = 1 For j # i we consid e r two cases: 1 ( A1 l A1)i1 = [(j l)(ni+ 1)( 1) + j(ni+ 1)(2)+(j + l)(ni+ 1)(1)] = 0 n+l if j = 1 2, i1, and 1 ( A1 l A1) i1 = [ i(ni(ml))(l)+i(nim)(2)+i( n i (m+ l))(l ) ] = 0 n+ 1 if j = i + m w h ere m = 1 2 , n i Therefore A 1 A= I Remark: Dr. J. Gopala kri s hn a n suggests to use the finite element method to find A1 1 Consider the two point value problem: { U" =Jon (0, 1), U(O) = U(l) = 0 Then the Green s fun c t i on G( x, s ) i s g iv e n by { (1 s ) x G( x, s ) = ( 1 x ) s if X < S, if X > S. L e t O = X o < X1 < < Xn+ l = 1 with Xi+i Xi = 1/(n + 1) for all i = 0 1 ... n. Defin e a(u v) = 11 u ( x ) v '( x)dx,
PAGE 55
and S = { v : v ( x ) is continuous on [0, 1], linear on e ach interval [xi,xi+1 ] for all i = 0 1 ... n and v(0) = v( l ) = 0}. The di sc r e t e G ree n s func t i o n Gi(x ) E S i s defined by 42 a(Gi, v) = v (xj) for all v E S ( 5.:) for any grid v erte x XjThen and { (1 Xj) X Gi(x) : = G(x Xj) = (1 x ) x1 if X X j , if X 2: Xji(n+lj) ifi<. ( n + 1)2 J G( ) j ( n + l i ) f . Xi' X j = ( n + l) 2 1 'I, > J Let { c/>i} fi l b e a basis of S where ( x ) = { I if X = Xi, 0 if X = Xj (j =/ i) Then t h e di sc r e t e Gree n function c a n b e written as G1(x) = (I x ;) t,,x.+ (~/' , ~ / 'x,) x; = (I x ;) t ,x, + Ct.1 ,t.1 ,x,) x;, s inc e X o = 0 and Xn+l = 1. Thus
PAGE 56
a ((1 Xj) t >kxk + X j ( t
PAGE 57
44 Lemma 5.2.2 L e t A1 b e th e matrix d esc r i b e d in Lemma 5.2 1 L et T b e then x n tri diagonal matrix with O s on th e diagonal with 1 's on th e s up e rd i agonal a n d with 1 'son th e su bdiagon a l The n IIA1 1Tll00 = n l. Proof: Since w e h ave n 2)A1 i)i,k(T)k,J k = l (a1)i,J1tJ1 J + (a1)i,J+1tJ+1, 1 2(n + 1 i) , if j = 1 2 , i1 n+l n + 2i 1 n+l 2i n+ l 2 i n+ l if j = i if j = i + 1 i + 2 , n 1 i f j = n The n t h e infinity norm for the matri x A1 1T i s n max L /(A1 1T)i j / t j = l max i1 ++n i [2(ni+l)( ) ln+2il l 2i ( )] i n+l n+l n+l [ 4i2 + 4i(n + 1 ) 2(n + 1) l 2i ( n + 1) 1 ] max+. i n+l n+l Now focus on the term in s id e brackets of (5.7): (5.7 )
PAGE 58
45 4 i2 + 4i(n + 1) 2(n + 1) l 2 i (n + l)j + n + l n+l [ 4i2 4i(n + 1 ) + (n + 1)2] + (n + 1)2 2(n + 1) l 2 i (n + 1) 1 '+n+l n+l [2i(n+ 1)]2 + l2i(n+ 1) 1 ( ) + n l n+l (5.8 ) < n l for a ll 1 i n where i and n are both integer s. The maximum value of Equation (5 8) i s n l when [2i (n + 1)]2 + l2i (n + 1)1 = 0 ow w e consider two cases: Case 1: n i s an odd number. The condition (5 9) is satisfied i f we choose i = (n + 1)/2. Putting t hi s i into Equation (5 7) w e obtain n 1. Case 2 : n i s a n eve n numbe r. The condition (5 9) i s satisfie d if we c h oost' i = n/2. Putting this i into Equation (5.7 ) w e obta in n 1. (5.9 ) Lemma 5.2.3 L e t J i = e ,;+i(~t) er;(~t), and I E C2 Th e n l f il = O(i6.t h) and Proof: Both equalities are prove d by the mean value theore m : (1) where 1i and ry~ are the value and t h e derivativ e of I at some zi between Zi and Z i + l Hence, I /i i = O(i6.t h).
PAGE 59
(2) ( 6.t)h ( er;+,(6t)'Y~+1 er;(6tly:) (6.t)h ( e7 ; (6t ) (1i(6.t) )1ih + et(6t)1i h) where 'Yi and "y~ are the value and the derivativ e of I at some zi between Zi and zi+1 'Yi+i and 'Y~+1 are the va lu e and the derivative of I at some zi+1 between =I =II Zi+l and Zi+2 and 'Yi, Ii, and Ii are the value the first d er ivative, a nd t h e seco nd d erivat ive of I at som e Zi b e tween zi a nd Zi+ 2 H e nce, l fi /i+1I < (6.t)h ((1i)2(6.t )h+fh) 0(6.t. h2 ) 5 3 MaxNorm Bound for Mat rix c1 D The maxnorm bound for the matrix c 10 wh e re C = A ~B, and 46 D =AA+ ~(l )B are discussed in this section. Only t h e case= l i s dis c uss e d for the convergence analysis. Note that
PAGE 60
H ence 11c1011 < + ~!!. IIA1 Blloo IIA1 (AA AA)lloo oo_ 21~IIA1 Blloo + l~IIA 1 B l loo whenever ~IIA 1 B1100 < 1 and where = max Ai. Our goa l in this sec tion i s to show that 11c1Dlloo = 1 + O(flt). 5.3 1 Maxarm for Matrix c1 D in 1Dimensional Spac e The preliminary results from Section 5.2 are used to show that IIC1D1100 = 1 + O(flt) in one dimensional space. 47 Lemma 5.3.1 L e t Ai b e the matrix d e scrib e d in L emma 5.2 .1. L e t Bi b e then x n tridiagonal matrix with O s on th e diagonal with d i 's on t h e s up e rd ia g on al a n d with di 's on th e s ubdiagonal Th e n Proof: The matrix Bi has the form 0 di d 2 0 d 2 d 3 0 d3 Bi= d 4 0 d 4 where each d i = O(flt). We can write the matrix Bi as
PAGE 61
B1 B' + B11 0 d1 0 0 0 0 d2 d2 0 0 di 0 d3 d1 d3 0 0 + d2 0 d4 d2 d4 0 0 We will compute A1 1B and A1 1B11 respective ly. (1) compute A1 1B': where (a1)i,j is the (i, j)th entry of A1 i H ence, Therefore, for any i n L /(A11B')i,j/ j=l 0 if _j = 1 2(n + i l) d if 2 _< J _< ,,1 n + 1 jl n + 2i 1 dj1 if j = i n+l 2i d 1 n + 1 1 if ( i 1) ::; j ::; n 2(n i + 1) / n + 2 i 1 / (d1 + d 2 + + d 2) + d:_1 n + l i n+l 1 2i + n + 1 (di+ di+l + + dn1) < 2( i 2)d + 3d + 2(n i)d, where d = max di= 0(.6.t) 0 ( ~t) 48
PAGE 62
49 1 II (2) comp ute A1 B : if j = l ( A 1B11 ) = 1 t,J 0 if j = n 2(n i + 1) n+l rlj=l (j + 1 ) ( n i + 1 ) (d d ) if 1 < 1 < i n + l J1 1+1 = i(ni k)(d d ) n + l J1 J + l if j = i + k (k = 0 1 , n i l ) 0 if j = n H ence, for any i, n L l (A1 lB'\,jl < 2(ni+l)d ni+l(3 4 .)ld' I 2 ++ + +i n+l n+l j=l i + n + 1 [(n i) + (n i 1) + + 2 + 1] J d' I where J d' I = max J dj1 dj+1l = 0(6.t h) < 0(6.t) + (i + 3)(i 2) 0(6.t h) + (n i + l )(n i ) O(l:!..t h) 2 2 = 0 ( ~t) Therefore hjj(A1 lB1) l loo < h (IJ(A1 1B')lloo + IJ(A1 l B ")JI ) = 0(6.t)
PAGE 63
Lemma 5 .3.2 L e t A1 b e t h e matrix d esc r i b e d in L emma 5 2 1 a n d A 1 b e th e diagonal matrix with >.i on th e diago n al. Th en 5 0 Proof: To prov e this L emma w e c omp u t e t h e m a xnorm of A1 1( A1A1 A1A1 ) L e t F = A1A1 A1A1 The n w e can write F as 0 Ji Ji O h F h O h 0 h !1 0 h h O f4 0 + 0 0 hh 0 T h e sam e pro cess in L emma 5.3 1 i s u se d t o co mpu t e A1 1F1 a n d A1 1F r espec t i v e ly. (1) compute A1 1F' : ( A 1F' ) 1 i,J 1 ((a 1 k j 1 (a1k1+1) 2 (n i + 1 ) ! i f J < i n + 1 1 (n + 2 i + 1 ) J if j = i n+l 1 if j > i
PAGE 64
51 Ther efo re, for any i, n L l (A1 lF'kjl = 2(n i + I ) I n + 2i I I n + 1 ( lfil + l h l + + l.fi11 ) + n + 1 l.fil j=l 2 + n: l (lh+1I + l.fi+2I +'''+Ifni) < 2(i I ) l f l + 3 lfl + 2(n i)lfl, where I l l = max llil = O(t::.t h) = (2n + I )O(t::.t h) = O(t::.t) 1 II (2) compute A1 F : if j = I if j =/= I 0 if j = I = (j I)(n i + I) (ff) if 2 :s; j :s; i + 1 n+I 1 1 1 i(n i + 2 k) (f ) f k (k ) n + 1 1 1 1 i j = i + = 2 3 , n i Henc e, for any i, n L l(A1 lF ")i,j l < ni+I n+1( I + 2 + + i) I J I j = l i + n + 1 [(n i) + (n i I)++ 3 + 2] lf'I where lf'I = max l fj1 fjl = O(t::.t h2 ) < ( i + I) i O(t::.t. h2 ) + (n i + 2)(n i I ) O(t::.t h2) 2 2 = O(t::.t).
PAGE 65
Therefore IIA1 l(A 1A1 A1A1) l loo < ll(A1 l F )lloo + ll( A 1 lF ") l loo < O( ~ t) 52 Theorem 5 .3.3 L e t A1 and B1 b e th e matrices d escri b e d in L emma 5.2.1 and Lemma 5. 3.1, respectively. Th e n Proof: ate t hat w h enever ~IIA1 1B1'100 < 1 a nd w h ere X = max>. i 1. Therefore from Lemma 5.3.1 and Lemma 5.3 2 we have 11c 1D1ioo < X + O(~t) + O(~ t) < 1 + O( ~ t) 5.3 2 MaxNorm for Matrix c 1 D in 2and 3Dimens i onal Spaces ow we extend the resu lts from the 1 dimensional case to 2a nd 3 dimensions First, we consider the 2dimens i ona l case Let A2 be a tridiagona l block matrix with ( A1 + 2 I ) 's on the d i agona l b l ock and I's on the sup e rdiagona l and s u bdiagona l block. Let B 2 and A 2 be d i agona l b l ock matrices with B1 s and A1 s
PAGE 66
on the di ag on a l blo c k resp ect ively. Then hA2 1B2 = h X I + 2 A1 1 A 1 1 A 1 1 I + 2 A1 1 and we can rewrite i t as A 1 1 Ai l I + 2A1 l A 1 1 A 1 1 5 3 I (5.11)
PAGE 67
where S2 i s t h e matrix Thus, I + 2A1 l A 1 1 A 1 1 I + 2 A1 1 A 1 1 Ai l I + 2 A1 1 A 1 1 By Lemma 5.3.1, w e hav e IIA1 1B11100 = 0(6.t/h). H e n ce A 1 1 I + 2 A1 l 54 1 Note that S2 is difficult to anal yze using the maxnorm. H e n ce, S2 is eva lu a ted on a computer to get a numerical valu e for 11S21100. For e xampl e if A1 is 10 x 10 and m = 10 so that S2 is 100 x 100 then IIS2ll00 1.6687. If A1 i s 20 x 20 a nd m = 20, so that S2 is 400 x 400, then I I S21100 2 0794. If A1 is 40 x 4 0 and m = 40, so t hat S2 is 1600 x 1600 then IIS2ll00 2.5041. Therefore we have (5.12)
PAGE 68
Also, we h ave X and then I + 2A1 l A 1 1 A 1 1 I + 2 A1 1 A 1 1 55 1 JJA2 l(A2A2 A2A2)lloo < JJA1 l(A1A1 A1A1)lloollS2lloo O(.6.t) (5.13 ) since IIA1 1( A1A1 A1A1)l l00 = O(.6.t) by Lemma 5.3.2. From Equation (5.10) we have the est imate 11c1 D1100 = 1 + O(.6. t). Three dim e nsions can be anal yzed in the same manner. We writ e A as a tridi ago n a l block matrix with (A2 + 2I)'s on t h e diagonal block and w ith I's on the superdiagona l and subdiagona l block. Write B and A as diagona l matrices with
PAGE 69
B2 's and A2 s on the diagonal block respectiv e ly. Write hSa and where the matrix Sa i s I + 2 A21 A1 2 A 1 2 A21 I + 2 A21 A 1 2 56 1 A1 2 I + 2 A21 To discu ss the maxnorm we still n ee d to eva lu ate t h e factor Sa Using a computer to e valuate the maxnorm of Sa, we have: if A2 is 5 x 5, so that Sa i s 12 5 x 12 5 the n IISalloo 1.9 841. If A2 is 10 x 10 s o that Sa is 100 0 x 1000 t h e n I ISalloo
PAGE 70
2.5126. Similarly w e hav e and hllA lBlloo < hllA2 lB2lloollS3lloo O(b.t) IIA l(AA AA)lloo < IIA2 l (A2A2 A2A2)lloollS sll O(b.t) 5 7 by Equation (5.12) and Equation (5.13) from t h e 2dimensiona l space. Finally we h ave 11c 1Dll00 = 1 + O(b.t) from Equati on (5.10). 5.3.3 um er ical R esu l ts for Spect r a l Radiu s of c 1 D ext we use the approximate va lu e of 1 ( Sect i on 3.3) to c alcu late so m e num er i cal results for spectra l radius of c 1 D (Tab l e 51, Tab l e 52 and Tab l e 53.), in the domain 01 when b.t = 0 1 s. umerically, we observe t hat the spectral radius of c 1 D is strict l y l ess than 1. Tabl e 5 1: Spect r a l rad iu s for c 1 D in one dim ens ion a l space Matrix dim e nsion ( n x n) n = 50 n = 100 n = 200 n = 400 n = 800 n = 1600 Spectral radius for c 1 D 0 .9990 0.999 1 0 9992 0.9992 0 9992 0 .9992
PAGE 71
Tabl e 5 2: S p ect r a l r a diu s for c 1 D in two dimensional pace Matri x dimension ( n x n) n = 100 (m = 10) n = 400 ( m = 20) n = 1600 (m = 40) n = 6400 (m = 80) Spectra l radius for c 1 D 0.9960 0.99 84 0 9989 0 9991 Table 5 3: Spectral r a diu s for c 1 D in three dimens ion al space Matrix dimens i o n ( n x n ) n = 125 (m = 5) n = 1000 (m = 10) n = 8000 (m = 2 0 ) n = 27000 (m = 30) Spectral radiu s for c 1 D 0.9715 0.9960 0.998 4 0.9988 58
PAGE 72
CHAPTER 6 UMERICAL EXPERIME TS This chapter bri efly reports some numerical experiments by u s in g t h e finite difference scheme des c ribed in Chapter 3 In t h ese a pplications t h e FORTRA package FISHPAK (see Appendix B ) is u t ilized to d ea l with l arge spar. e ma trices in our model a nd to so lv e lin ear separable e llip t i c eq uati ons in t hr ee space dimension s. The conductivity function is compute d as a = 10p(z) (6.1 ) where p(z) is the quartic pol y nomial (see Equation (3.21) and Tabl e 31) and the altitude z i s m easured in kilom ete r s. Let u s ass um e as in Hager et al. [7] that the curre n t p er unit volume at position r is where Sp and SN are t he centers of the posi t ively c h arge d a nd n egative l y charge d current generators respectively. The parameters Ap a nd AN are ch ose n so that the total current generate d b y either t h e positive or t he nega t iv e charge centers has magnitude 1 A. H e n ce we h ave Ap { e wR2 dV = l }z?_h1 wher e h1 i s the altitude for t h e posi t iv e charge cente r R2 = I r Sp 12 and z is the altitude m easure d with re spect to the posi t iv e charge ce n ter. The cy lindri ca l shell method is applied to the above integ r a l as sho wn in F igur e 6 1 to obtain 59
PAGE 73
60 z axis (altitude) p dp Figure 6 1: Cy lindrical s h e ll of t h e volum e in tegral/ dV wher e e r f is the e rror function d e fin e d as 2 t 2 e r f ( x ) = .,fa J o et dt. As a r esult, Ap i s calc ul ated as ( w ) 3 / 2 2 Ap = ; [ 1 + e r J( jwh1)]' a nd si milarly, AN can be obtained
PAGE 74
61 where h2 is the altitude for t he negative charge center. The refore is used in our model. In each experiment, the initial potential of the entire domain is assumed to be zero. The centers of the current generators are posi t ioned on the zaxis at altitudes 5 km () and 10 km ( + ). Four experiments are conducted in this chapter. Section 6.1 us es different values of w to find the breakdown time with the ambient conductivity and the r educe d conductivity. Section 6 2 explores a fin e r mesh for our mode l with an acceptable accuracy for the potential Section 6.3 and Section 6.4 c heck t h e temporal error and the spatial error nume rically respectively. 6.1 Experiment 1 This experiment considers the domain 01 previously introduced: 01 = {( x, y z ) I O :S z '.S 100 km, 5 0 km :S x y '.S 50 km} For conductivi ty, we use the ambient conductivity given in Equation ( 6 1). For different values of w we integrate the Equation (3 7) solving for potential The integration is stopped when the electric field reaches t he breakdown t hreshold which is taken to b e 250 000 V m 1 [7]. W e u se the uniform m es h in a ll xy and zdirection s with mesh size 1000 m and the time step is 6.t = 1 / 16 s s s hown in Tabl e 6 1 the breakdown time i s calc ulated base d on diff e r ent values of w Table 61 demonstrates that the breakdown time shortens rapidly as w increases. For example, when w = 4 the breakdown occurs within 15 .5 625 seconds With w = 16 t h e breakdown time drops to 1.1850 seconds. T hi s r esult i s supported by [7]: as w in creases t h e current b ec omes more co n centrated at Sp and SN. Observe that wh e n w = l there i s no breakd own within o n e hour. A cco rd i n g to our numerical r esults the maximum e lectri c fie ld approac h e s t h e limit 232429 88
PAGE 75
Table 6 1: Breakdown time comparison for differ ent values of w w value w=l w=4 w = 16 breakdown time with the ambient conductivity (without cloud) 00 1 5.5625 s 1.1850 s breakdown time with the reduced conductivity (with cloud) 64. 7500 s 13 687 5 s 1.75000 s V /m which is less than the breakdown threshold 250 000 V / m as t tends to infinity. Therefo re t h e breakdown n eve r happe ns when w = 1. 62 In a thundercloud the conductivity is l ess than the ambient cond uctivity that appears in Equation (6. 1) ext, w e consider the e ffect of reducing conducti vity in order to simulate a cloud. More pre c i se ly the ambient conductivity is multiplied b y the fa cto r 0 0 5 as described in [7]. The results a r e shown in t h e third column of Table 61. Again the breakdown time d ec r eases dramatically as w increases. This pattern is consistent with the previous computed breakdown t ime s for t h e ambient conductivity. As can b e see n in Table 61 r educing in conductivity to s imulate a cloud leads to the occurrence of breakdown when w = 1 w hil e wit h the ambi e n t conductivity charge i s conducted into the surrounding atmosph e r e faster eno u g h to prevent t h e e l ect ri c fie ld from reac hing t h e breakdown thr es hold. In t h e remaining expe rim ents, the conductivity is give n b y Equation ( 6 1 ) t im es 0 .05 (to s imulate a cloud) 6.2 Experiment 2 A fin er m es h in our mode l is expecte d to obtain more accurate res ul ts. How eve r s u c h r e fin ement i s not onl y t im e consuming, but a l so incr eas in g m emory u sage In ord e r to m a intain t h e accuracy of our r esults effic i entl y a g rad ed me h i s considered. W e pl ace a l arge percen t age of g rid points near th e c h a rgf' rc'nt c'r s a l o n g t h e zdirect i on w h ere t h e potential chan ges rapidl y T h e g rid a long the
PAGE 76
63 zaxis is c hos e n in t h e following mann e r [7] : n for 1 < k < 2 ( 6 .2) and for'!.!.< k < n 2 ( 6.3 ) whe r e His t h e a l t i t u de o f t h e p os i t ive c h a r ge center (10 km ), T = 4H/n, a nd sis c hos e n so t hat 2H + TS(n/2 ) = 100 km is the h e ig h t of the d o main. Acc ordin g to Hag e r e t al. [ 7], t h e v icinity of the charg e c ente rs i s c ritic a l to efficiency. T h e r e for e the dis t ribu t ion o f t h e g rid points i s unifo r m b e low 2 H a nd geo metric fr om 2 H to 100 km W e still u se the u n iform mes h for t h e x and y d ir ect i ons as i n Exp e rim e n t 1 L e t u s co n s id e r two d o main s : 01 as d esc rib e d in Sect i on 5.1 a n d 02 as follows, 02 = {( x y z ) I 0::; z::; 100 km 2 5 km::; x y::; 2 5 k m } In this e xp e rim e n t t h e pote n t ials are found at the t im e 6 5 sec ond s with t h e time ste p t::i..t = 1 / 2 s an d t h e ambi e n t c ondu ctiv ity r e du ce d b y t h e factor 0 0 5 w h e n w = l. The g r a d e d mes h i s appli e d with 200 g rid p oints a l ong zaxis ( m e sh s i z e = 200 m for 0 ::; z ::; 20 km ). The mesh s ize i s 1000 m for both in x a n d y dir ect ions. Fi gure 6 2 s how s t h e pote nti a l s in tw o d o m a in s 01 a n d 02 for x = 0 y = 0 a nd 0 ::; z ::; 100 km The p ote ntial fro m 01 i s s hown in a so lid lin e an d t hat from 02 i s s hown in a d as h e d lin e T h e a b solut e diff e r e n ce i n poten t i a l i s abo u t 105 a nd t h e r e lative diff ere n ce i s a bou t 10 4 b etwee n t h ese two do m a in s T h e r efor e w e c anno t see a n y diff e r ence for t h e p ote ntial s in Fi g ur e 6 2 F igure 6 3 a nd Fi g ur e 6 4 are the p ote n t i a l s n ea r t h e n ega tive a nd pos i t i v e c h a r ge cente r s r espect iv e l y a l o n g t h e zaxis w h e n x = 0 a nd y = 0 Since t h e r e lative diff e r e nce for po t e n t i a l s i s quite s m a ll between t h ese two d o m a i ns w e on l y foc u s o n t h e s m a ll e r dom a i n 02 for t h e r emai n in g e xp e rim ents.
PAGE 77
X 108 6.,~,,~~~~~~ 2 0 > .!: 4 2 .; o :g Q) 0 a. 2 4 See F i g ure 6 .4 l See Figu re 6 3 p o t e nti a l i n th e domain D.1 po t e nt ia l in the domain 6~~~~~~~~~~~~ 0 2 3 4 5 6 Altitude in meters 7 8 9 1 0 X 104 F i gu r e 6 2: Pote nti a l s ((0, 0 z ) 0 z 100 km ) in two d o m a in s, n1 and n2 X 108 5. 15 5. 2 5. 25 5.3 .'!J 5.35 0 > 5.4 i: 0 a. 5.45 5. 5 5. 55 5.6 5 65 4600 4 700 p ote nti a l in th e doma in n pot e n t ial in the d o main 4 8 00 4 900 5000 Altitude i n m e ters 5100 5200 5300 5400 6 4 Fig ur e 6 3: P ote ntial s (( 0 0 z), 4 6 km z 5.4 km ) n ea r t h e n egat iv e c h a r ge (5 km) in two dom a ins, n1 and D 2
PAGE 78
!I? 0 > .!: i: ., a a. X 109 5 65 ~ ~~5 6 5 55 pote nti a l in th e domain .n1 potential in the domain 5 15 .__ __ _._ __ _,._ __ __._ __ ___. ___ ..__ _ ...J._ __ '' 0 96 0 97 0 98 0 .99 Altitud e i n m eters 1 .01 1 02 1 03 1 04 x10' 65 Figure 6 4: Potentials (cp(0, 0 ,z), 9.6 km ~ z 10.4 km) near the pos i t iv e c h a r ge (10 km) in two domains, 01 and 02 X 107 2 5~!I? .!: i: 2 poten ti a l a t th e a ltitud e 10km 1 5 potential a t the a ltitud e 5km a o 5 a. ' \ \ 0. 5 \ I I I / I 1 l...''L''''''_J 2 5 2 1. 5 1 0 5 0 0 5 xaxi s i n m e t ers 1 5 2 2 5 X 10' Figure 6 5: Potenti a l s (cp(x, 0 5) and cp(x, 0 10) 25 km x 25 km) a l ong t h e xaxis at t h e a l i t ud e 5 km and 10 km
PAGE 79
66 Further, we find t he potentials (see Figure 6 5) a long the xaxis at the altitude 5 km (the negative charge center) and 10 km (the positive charge center) when y = 0. In Figure 6 5, the so lid lin e represents th potential at the altitude 10 km and the dashed lin e represents the potential at the altitude 5 km. 6.3 Experim nt 3 According to the error estimate described in Theorem 4.2.4, when = 1/2 (Crankicholson scheme) the tempora l error term should be smaller than for other values of. To verify the tempo r a l error anal ysis numeri ca lly the potent i a l is estimated at time t = 65 s, which i s close to the breakdown tim 64. 7 seconds when w = 1. The graded mesh (see Equation (6.2) and Equation(6.3)) is applied with 200 grid points a long zaxis (mesh size = 200 m for 0 ::; z ::; 20 km). The mesh size is 1000 m for both in x and yaxes. We compare the erro rs in different time steps (1/2 1/4, 1/64) for = 1/2 and = 1. To analyze the e rror numerically we use the maximum norm: where (t::.t)i,i,k is the computed potential using the time step 6.t at each mesh point in the domain at t = 65 s, and t::.t is the potential vector with elements ( t::.t)i,i,k First, we consider the potential cor r esponding to the time step 6.t = 1 / 128 s as the "exact potential and subtract i t from the potentia l s corresponding to the compared tim steps to obtain e rror estimates at all mesh point T h e maximum norm e rror at diff e ren t time step i s hown in the third column of Table 62 and Table 63. Second the relative errors are estimated by th ratio ll4tt::.t 1;12slloo ll1;12slloo as shown in the fourth column of Table 62 and Tabl e 63
PAGE 80
67 Table 6 2: Relative error in the potential at t norm for Crankicholson scheme ( = 1/2). 65 seconds using th maximum Time step (s) 1/2 1/4 1/8 1/16 1/32 1/64 1/128 maximum norm potential 564 431185.91904 564431163 61478 564431157 94013 564431156 50102 564431156.13736 564431156. 04677 564 431156.04677 maximum norm error Cranki cholson with respect to relative error time step 1/128 s ( = 1/2) 80 21113 1.42110 X 107 20.30029 3 59659 X 10s 5 .11 413 9 06068 X 109 1.27381 2.25681 X 109 0.30471 5.39852 X 10lO 0.060 52 1.07192 X 10lO 0.00000 0 00000 Table 6 3: Relative error in the potential at t norm for Backward Eul er scheme ( = 1) 65 seconds using the maximum Time ste p maximum norm maximum norm error Backw a rd E ul e r ( s) potential with r es pect to r e lativ e error time step 1 / 128 s ( = 1 ) 1 / 2 5 6 4422403. 52 8 80 3 5 27 .792 T 6.2 5 033 X 105 1/4 564426772. 69223 173 75 .013 7 3 3 0 7 32 X lQ1/8 5 6 4 428962. 54654 8 4 11.09363 1.4 9019 X 10 5 1/16 564430058.82113 3926 08120 6.95582 X 106 1/32 564430607 .30153 1682 80248 2.98141 X 10 6 1/64 564430881.63001 560 96757 9.93864 X 10 7 1/128 564 430881.63001 0.00000 0 00000
PAGE 81
6 In Table 62 which corresponds to the Cranki c holson ch emr, if t h time' step b..t is multiplied by 1/2, the relative error is multiplied by 1/4. H ne e the error for CrankNicholson scheme ( = 1/2) behaves lik e O(b..t)2. In Table 63 which corresponds to t he backward Euler difference scheme, if the time step b..t is multiplied by 1/2, the relative error is a l so multiplied by 1/2. Therefore the error for backward Euler difference scheme ( = 1) behaves like O(b..t) Thi sugge ts that the Crankichol son scheme i s far superior to the backward Euler difference scheme. ote that the forward difference scheme ( = 0) does not work since it is not stable as we de crib ed in Section 4.1. Experimentally taking = 0 we find that the maxnorm of the potential grows exponentia lly int. The plot of log10(llt.tlloo) versus t appears in Figure 6 6, with time step b..t = 1/8 s. Observe that the maximum potential at t = 65 seconds i s 1.04737 x 1089 90 80 ,e 70 ... 0 0 60 C: E ::, ., 50 E 00 2 40 .9 30 20 10 0 0 10 20 30 4 0 50 6 0 70 Time in sec ond s Figure 6 6 : Maximum norm of potential v e rsu s tim e wh e r e 6. t = 1/8 s and 1 1 = 0
PAGE 82
69 6 .4 Exp e rim e n t 4 The last study is t o ch ec k the s p a tial e r ro r t e rm. W e take the ti m e step 6.t = 1/2 s and us e Crankichol son sch e m e = 1 / 2 Four diff e r ent g r a d e d m es hes a r e us e d in the dom a in fh whe r e 02 = {(x, y ,z ) I 0 z 100 km, 25 km~ x,y 2 5 k m }, and t h e g rid p o ints are c h ose n a cco rdin g to E quati on ( 6.2 ) an d Eq uati o n ( 6 3) along the zaxis : Cas e 1: 30 uniform g rid points in the x and ydirect ion s and 100 gr id points in the z direc t ion (m e sh siz e = 4 00 m for O z 20 km and exp on e n t iall y in c r eas in g m es h for 20 z 100 km fro m 6.z = 4 00 m near z = 20 km t o 6.z = 80 4 3.61 m near z = 100 km ) Cas e 2 : 60 uniform grid points in the x a nd ydirect i on s an d 2 00 gr id point s in t h e zdir ect i on (mesh s ize = 200 m for O z 20 km and e xpon e nti a lly in c r easing mesh for 20 z 100 km from 6. z = 200 m near z = 20 km t o 6.z = 4 6 52.4 1 m near z = 100 km ) C ase 3: 120 uniform g rid points in the x and ydirect ion s an d 4 00 gr id points in t h e z dir ect i on ( mesh size = 100 m for O z 20 k m and e xpon e n t i ally in c r eas in g mesh for 20 z 10 0 km fr o m 6.z = 100 m nea r z = 20 km to 6. z = 2629 .65 m near z = 100 km ) Case 4: 2 4 0 uniform gr id points in t h e xa nd y di rect ion s an d 8 00 gr id points in the zdir ect i on (mesh s ize = 5 0 m for O z 20 k m and e xpon e n t i a ll y in c r eas in g m e sh for 20 z 100 km fro m 6.z = 5 0 m near z = 20 km t o 6.z = 1 4 62.03 m near z = 100 km ) Let
PAGE 83
70 the l argest potential is achieved here Suppose that /'0 is the exac t potential at z = 10 km and the sequence
PAGE 84
Any two potential vectors qi and 2i can be employed in the same way as above to find the constant C: qi_ i C= 3 4i2 Table 64: E t im ates for t h e co nstant C gen e rated by Equation (6. 7) n cpn cp2n C 100 17 5594927.29868 2.34 X 108 200 25121434.49608 1.34 X 108 400 5574244 17420 1.19 X 108 800 71 (6.7 ) In Table 64, only four different mes hes ar e considered. If n < 100 t h e potential vector i s not accurate enough, and if n > 800 there is not enough memory to do computation s on our co mpu ter network. The maximum norms of t h e diff erence of two succe sive potential v e ctors are shown i n the sec ond c olumn of Table 64. Obs e rve that t he maximum norm of the diff e rence are getting smaller. We can predict that i t approac h es O whe n n t e nds to infini ty if mor e data (potent ial s for l arger n) provided. This s hows that the s e quence n is c onver g ing The estimates for the co nstant C a r e s hown in t he t hird colum n of Tab l e 64. The relative diff erence for the l ast two estimates of C i s 0 .13. ext, let us as ume that the spatia l e rror is O ( h P ) that i s
PAGE 85
72 where IEnl ::; ChP C and pare constants. Using the similar anal ys is as above we will compute the constants p and C. For three e lem ents P, 2n, and 4n, we have (6.8)
PAGE 86
73 Table 6 5: Estimates for constants p a nd C generat e d by Equation ( 6.10 ) n qP n p C 100 175594927.29868 200 2.81 2.04 X 108 25121434.49608 400 2.17 1.45 X 108 5574244 17420 800
PAGE 87
CHAPTER 7 CONCLUSION In this dissertation Maxwell's eq uations are used to establish an electric field model between the surface of the Earth and the ionospher e : PDE: ~ 6 +
PAGE 88
APPENDIX A APPROXIMATIO FOR CO DUCTIVITY F U CTIO MATLAB programs ar e used to o btain t h e a pproximati on for t h e con du c tivity function. The d e tailed approximation for the c ondu c tivit y fun c t i o n fr o m Unit e d Sates Air Force Cambridge R esear c h Laborator ies [27 ] i s g i v e n in App e ndi x A .1. The approxima t ion from Lightn ing: phy s i c s a nd effects [21] i s disc u ssed in Appendix A.2 A 1 Approximation from Figur e 9.6 [27 ] Four point s (se e Table A 1) on t h e curv e B from Figur e 9.6 [ 27 ] a r e c h ose n to find an approximation for conductivit y The c onductivit y graph in F i g ur e 9 6 [ 2 7 ] contains the d ata up t o the a ltitud e 27. 828 km (ab o u t 100 ft), whil e o ur m ode l domain is 100 km. All four points are chos e n to h ave a cubi c fit for t h e al t itud e 0 :S z :S 27. 828 km and the approximation c an b e obtaine d b y MAT LAB as follows: O" ( z ) = 1o13 .4 0 + 0.1 600z0 0051 7 5 z2+ 0 0000 9077 z3 For altitude 28.828 :S z :S 100 km the last three points (A, P3 and P4 ) are pick e d to have a linear fit and the approximation is a( z ) = 1012 88 + 0 0671 8 z Figur e A 1 s hows the graph of the approxim a tion for the condu c tivi ty fun c tion A.2 Approximation from Figur e 3 3 Thirte en points ( see Tabl e A2) in Figur e 3 3 a r e pick e d to a pproxim ate t h e condu c tivit y fun c tion The c ondu c tivity i s c h o s e n as t h e av e r age of s h a d ed r e g i o n in Figur e 3 3 for the a lti t ud e a bov e 30 km. The a d va ntage for u s i n g t hi s figure 7 5
PAGE 89
Table A 1: Four points chosen from Figure 9.6 [27] Point Heigh t log10 (conductivity ) ( km ) (1/oh mm ) Pi 2.7432 13 A 12.954 12 ?3 20.7264 11.5 ?4 27.82824 11 100 90 80 ] 70 ., 0 60 ~ :;: 50 40 30 20 10 0 14 13 12 11 10 9 8 7 6 lo g 1Jconducti vity), (Oh m meter) I Figure A 1: Electric conductivity as a function o f altitude approximated from curve B in Figur e 9.6 [27] 76
PAGE 90
to estimate t h e conductivity i s that Figure 3 3 d emonstrate s a more completed conductivity data than Figur e 9.6 in [27]. Table A 2: Thir teen po ints c ho s e n from F i g u r e 3 3 Point Pi A P 3 P 4 Ps p 6 P 7 P s P g P10 P n Pi2 P13 H e igh t (km) 4.5 12 .5 26 40 48 52.5 60 70 78.5 85 89 95 100 l og10 (conduct i vity) (S/m) 13 12 11 10 9 8 7 6 5 4 3 2 1 A quartic pol ynom ial is found to approximate those points in Table 2. T h e conductivity function is est imated b y (J = lQP (z) w h ere p(z) i s the qua r tic polynomial p( z ) = 13 16 + 0 063 5 1 z + 0 0008682 z2 0.000008674z3 +o. 00000005 781 z4 where the altitude z is m easure d in kilometer s and its graph i s s h own in Figure A 2 77
PAGE 91
100 90 80 ] 70 QJ "Cl 60 2 50 40 30 20 10 0 14 12 10 8 6 4 2 0 log 1Jconductivity) S / m Figure A 2: Electric conductivity as a function of altitude approximated from Figure 3 3 78
PAGE 92
APPENDIX B FISHPAK The FORTRA packag e FISHPAK d e v e lop e d b y Swarztrauber and Sweet is used for our nume ri c al ex periments. Documentati o n for the pack age i s g i v e n in ACM Algori thm 54 1 [ 24]. The original p ac k age i s w ritten i n s in g l e preci s ion. In orde r to get more pre cis e nume ric a l r esults, a double prec i s i o n vers i o n p ac k age is obtained by a little adjust m ent. The subroutine POIS3D is us e d t o s ol ve for potentials in our mode l. POIS3D solv e s the linear s ystem of equations for unknown X values Cl *(X(I1 J K) 2.*X(I, J K ) + X(I+l, J K )) + C2 (X(I ,Jl, K) 2 X ( I J K ) + X(I,J+ l K )) + A(K)* X ( I J K1 ) + B(K)* X ( I J K) + C( K )* X ( I J K + l ) = F(I, J K) whe r e I= 1 2 , L J = 1 2 , M and K = 1 2 , N The indice s K1 and K+l are evaluate d modulo N i.e X(I, J ,0)=X(I, J ,N) and X ( I J + l )=X(I J ,1 ). The unknowns X ( 0 J K) X (L+l, J K) X(I, 0 K ) and X ( I ,M+l K ) a r e assu m ed to take on c ertain presc rib e d v a lu e s T h e o p e rati o n count for POIS 3D i s ro u gh l y proportiona l to L x M x N x (log2 M + log2 + 5) 79
PAGE 93
REFERE CES [1] E. J. Adl e rman and E. R. Williams S eas onal var i ation of the g lob a l e l ect ri cal circuit Journal of G eophysic al R esea r ch, 101, 1996 pp.29679 29688. [2] M. B. Baker H. J. Christ i an and J Latham A comput ional s t udy of t he relationships linking li ghtning fre quenc y and other t hund e rcl oud parameters Quart e rly Journal of the Royal M e t eoro logi ca l S ociety, 121 1995, pp. 1 525 1548 [3] W. W. Hager Applie d Numerica l Linear Alge bra, Prent iceHall Englewoo d Cliffs J 1988. [4] W. W H ager Updating the inv e r se of a matrix SIAM Rev. 31 1989 pp.221239. [5] W. W Hager A di sc r ete Mode l for t h e Lig htnin g Disch a r ge J ournal o f Computatio nal Physics 144 1998 pp .137 1 50. [6] W W Hag e r J S Nisbet, and J R K asha, The evo lu t i on and di sc h arge of electric fie ld s within a thund ersto rm J o urnal of Computational Physics, 82, 1989 pp. 193 217 [7] W. W Hager J. S isb e t J. R Kash a and W.C. Shann Simu lati ons of e lectri c fie ld s within a thund e rs to rm J ournal of t h e Atm osph eric S ciences 46 1989 pp. 35423558 [8] L. C. H a l e Midd l e atmospheric e l ectrical stru c t ure d y n a mics, a n d co u pling. Earths Middl e Atmosphere Advances in Space R e s earch. 4 1 984. p. 1 75 I 6 [9] J H H e l sdon Jr. and R. D F arley A num e r i ca l mod e ling tudy o f a Montan a thunderstorm: 1. Mod e l r es ults v ers u s observations i nvol v in g non e l ect rical aspects, Journal of G eophysica l R e search 92 198 7 pp.56455659. [10] J H Hel s don Jr. a nd R. D. Farley A num er ical m ode lin g s tudy of a Mont a n a thunderstorm: 1. Mod e l results v ersus observatio n s invo l v i n g e l ec trical as p ects J ourna l of G e ophysi c a l R e sear c h 92 1987 pp 56615 67 5 [11] J H H e l s don Jr. R. D. F a rl ey a nd G. Wu, Lig h t nin g param eteriza t i o n in a stor m e l ect rification model 19 88 P rocee ding s International Con f e r ence on Atmosp h e ric Ele ctricity, pp.849 854 80
PAGE 94
81 [12] J H. Helsdon Jr. G. Wu, and R. D. Farley An intercloud lightning parameterization sc h eme for a storm electrification mod e l Journal of Geophysic al R ese arch 97 1992 pp.5865 5884. [13] H. W. Kasemir A contr ibuti on to the e l ectrostat i c theory of a li ghtning discharge, J ournal of G e ophysical R esearc h 65 1960 pp 1873 1878. [14] D. R. MacGorman, A. A. Few and T. L. Teer, Layered lightning activity, Journal of Geophysical Research 86 1981 pp.9900 9910. [15] D. R. MacGorman J. M. Straka, and C. L. Ziegler, A lightning parameteri zation for numerical cloud models Journal of Applied Meteorology 40 2001 pp.459 478. [16] E. R. Mansell D R. MacGorman, Ziegler and J. M. Straka Simulated threedimensiona l branched li ghtning in a num er ical thund e rstorm mod e l Journal of G eophysica l R esearch, 107 2002 doi:10.1029 / 2000JD000244 [17] L. iemeyer L. Pietronero, H. J Wiesmann Fractal dimension of dielectric breakdown Physical R eview L etter, 52 1984, pp.1033 1036. [18] F Pockels Uber das Magnetische Verhalten Einger Basaltischer Gesteien Ann. Phys. Chem., 63 1897 pp.195 201. [19] F Pockels B esti mmung Maximaler EntladungsStromStark e n aus Ihr e r Magnetisirenden Wirkung Ann. Phys Ch e m ., 65 1898 pp. 458 475. [20] F. Pock e ls Uber die Blizentlandun gen Erreicht Stromstark e Phys. Z., 2 1900 pp.306307 [21] V. A. Rakov and M. A. Uman Lightning: physics and effects Cambridge University Press Cambridge, New York, 2003 [22] F Rawlins A numerical study of thunde rstorm e l ectrification using a t hre e dimensional mode l incorporating the ic e phas e Quarter ly Journal of the Royal M e t e orological Socie ty, 108 1982 pp.779 800 [23] J C Strikwerda Finite D ifference S c h e m es and P artial D ifj' e r ence Equations, Champman & Hall ew York NY 1989. [24] P . Swarztrauber and R. A Sweet Algorithm 541: Effic ient Fort r a11 subprograms for the solution of e lli ptic partial differentia l quati o n s A C M Trans. Math. Software, 5 1979 pp.352 364. [25] T. Takahashi D ete rmination of lig h tning origins in a thunderstorm model Journal of the Meteoro l ogica l Society of Japan, 65 1987 pp.777 794 [26] M. A. Uman The Lightning D isc harge, Academic Press Orl a ndo Florida 1987.
PAGE 95
[27] United States Air Force Cambri dg e Research Laboratories H andbook of Geophysics, Macmillan e w York 1960 ( chapte r on atmospheric e l e ctri city revised in 1970). [28] R. S. Varga Matrix I terative Analysis PrenticeHall, Englewood C liff s J, 1962. [29] H. J. Wiesmann and H R. Zell e r A fractal model of dielectric breakdown and prebreakdown in solid dielectrics Jounral of Applied Physics 60 1986, pp.1770 1773. [30] E. R. Williams C. M. Cooke, and K. A. Wright, Electrical discharge propagation in and around space charge clouds, Journal of G e ophysical R esearc h 90 1985 pp.6054 6070. 82 [31] C. T R. Wilson On Some D eterm ination s of t h e Sign and Magnitude of Electric Di scharges in Lightning Flashes, Proceedings of the Royal Soci e ty of London S e ries A 92 1916, pp. 555 574. [32] C. T. R. Wilson, Investigations on Lightning Discharges and on the Electric Field of Thunderstorm, Philo sop h ical Transactions of th e Royal Society of London S e ries A 221 1916, pp.73 115. [33] C. L. Ziegler and D R MacGorman, Observed li ghtning morpholo g y relativ e to modeled space c h arge and e l ect ri c field distributions in a tornadic storm Jou rnal of th e Atmosph e ric Sci e nc e s 51 1994 pp.833 8 51. [34] C L. Ziegler D R. MacGorman J E Dye, and P S Ray A mode l eva lu ation of non inductive graupe l ice charging in the early electrification of a mountain thunder storm, Journal of Geophysical R e s e arch 96 1991 pp.833 851. [35] C. L. Ziegler, J. M. Straka, and D. R. MacGorman Simulation of e arly e l ectrification in a threedimensional dynamic cloud model 11th Int e rn a tional Conf e r e nc e on Atmospheric El e ctricity Comm i ssion on Atmosphe r i c Electricity 1999 pp.359 362.
PAGE 96
BIOGRAPHICAL SKETC H ShuJen Huang was born in Taiwan in 1970. She w as awa rd ed a Bachelor of Scie nc e d eg ree in a ppli e d math e matics in 1992 from F uJ en Un i versity, Taipei, Taiwan. Sh e worked as a teaching ass istan t at FuJ en Uni versity for 3 years after h e r graduation In 1995 ShuJen started h e r gra u date study in mathe matics at the University of Florid a, from whic h she receiv e d h e r Ph. D in mathematics in 2005.
PAGE 97
I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scope and quality as a dissertation for the degree of Doctor of Philosophy. w illiam \V. Hager Chair Professor of 1 fathematics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scope and quality as a dissertation for the degree of Doctor of Pru ~ Jay Gopalakrishnan Assistant Professor of Mathematics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate. in scope and quality as a dissertation for the degree of Doctor of Philosophy. ~tilt /tfJ4itrJv,..__ Shari Moskow Associate Professor of 1Iathematics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scope and quality, as a dissertation for the degree of Doctor of ,r,,'~~i S. Pilyugin Assistant Professor of 1Iathematics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scop e and quality as a dissertation for the degree of Doctor of Philosophy. Vladimir A. RakoY Professor of Electrical and Computer Engineering
PAGE 98
This dissertation was submitted to the Graduate Faculty of the D epartment of Mathemat ics in t h e College of Liberal Arts and Sci ences and to the Grad u ate School and was accepted as partial fu l fillm ent of the requirements for the degree of Do ctor of Philosophy. May 2005 D ean, Grad uate School
PAGE 99
' LD 1780 2005 H B 1'f UNIVERSITY OF FLORIDA II I II IIIIII Ill I l l lllll l llll I I IIIIII I III III I Ill lllll II IIII IIIII I I 3 1262 08554 2586

