Citation
Multiscale discretization of electric-field equations

Material Information

Title:
Multiscale discretization of electric-field equations
Creator:
Huang, Shu-Jen
Publication Date:
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 )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 2005.
Bibliography:
Includes bibliographical references.
General Note:
Printout.
General Note:
Vita.
Statement of Responsibility:
by Shu-Jen 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 non-profit 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
ELECTRIC-FIELD EQUATIONS















By

SHU-JEN 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

Shu-Jen Huang
















To Deng-Shan,

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, Deng-Shan, 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 parents-in-law 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 Max-Norm Bound for Matrix C-1 = (A B) .. 38
5.2 Estimates for A -1 ........................... 40
5.3 Max-Norm Bound for Matrix C-1D ................ 46
5.3.1 Max-Norm for Matrix C-1D in 1-Dimensional Space 47
5.3.2 Max-Norm for Matrix C-1D in 2- and 3-Dimensional Spaces 52
5.3.3 Numerical Results for Spectral Radius of C-1D ...... 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 3-3 ..

B FISHPAK ...................

REFERENCES ...................

BIOGRAPHICAL SKETCH ............


FUNCTION















LIST OF TABLES
Table page

3-1 Coefficients of the polynomial p(z) .... 25

5-1 Spectral radius for C-1D in one dimensional space ... 57

5-2 Spectral radius for C-1D in two dimensional space ... 58

5-3 Spectral radius for C-1D in three dimensional space ... 58

6-1 Breakdown time comparison for different values of w ... 62

6-2 Relative error in the potential at t = 65 seconds using the maximum
norm for Crank-Nicholson scheme (pA = 1/2). ...... 67

6-3 Relative error in the potential at t = 65 seconds using the maximum
norm for Backward Euler scheme (pA = 1). ...... 67

6-4 Estimates for the constant 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















LIST OF FIGURES
Figure page

3-1 Path PXPj 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 a and corresponding relaxation time T = Ea-1,
where E = 8.85-12 F m-1, versus altitude under a variety of geo-
physical conditions. LL, low latitude, "wavy"; MLPS, mid-latitude
typical night (high-latitude, quite); AZTDN, auroral zone, typical
disturbed night; MLD, mid-latitude day, quite; MHL, mid-high-
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
4-1 Graph of f(x) = sn+ )2, when a = -1.05. ........... 30
(a + cos X)2
sin x
4-2 Graph of f(x) = cosx)2, when a = -1. .............. 31
(a + cos X)2

6-1 Cylindrical shell of the volume integral f dV 60

6-2 Potentials (0(0, 0, z), 0 < z < 100 km) in two domains, Q1 and Q2 64

6-3 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

6-4 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

6-5 Potentials (O(x, 0, 5) and O(x, 0, 10), -25 km < x < 25 km) along the
x-axis at the altitude 5 km and 10 km .. ..... 65

6-6 Maximum norm of potential versus time, where At = 1/8 s and p = 0 68










A-1 Electric conductivity as a function of altitude approximated from curve
B in Figure 9.6 [27] . 76

A-2 Electric conductivity as a function of altitude approximated from Fig-
ure 3-3 . . 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 = e-0(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 = h-I



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
ELECTRIC-FIELD EQUATIONS

By

Shu-Jen Huang

May 2005

Chair: William W. Hager
Major Department: Mathematics

We develop stable finite-difference 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 Crank-Nicholson scheme. The approximations are unconditionally stable when

1/2 < p < 1 and unstable when 0 < p < 1/2. Max-norm 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,

Solid-state 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 three-dimensional 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 high-density 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 mass-weighted 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 two-dimensional; while the models of Ziegler, Straka,

MacGorman, Dye, and Ray [33, 34, 35] were established in three-dimensional space.

The first term on the right-hand side of Equation (2.1) was the advection term.
The second term was the fallout term, and it was equal to zero for cloud-water










and cloud-ice charge. The third term was the eddy-mixing 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 small-ion

mobility (in meters squared per volt per second, pressure dependent), E was the

two-dimensional electric field vector, G was the ion-generation rate due to cosmic

ray ionization (in ions per cubic meter per second, dependent on height), a was

the ion-ion recombination coefficient (1.6 x 10-12 m3/s), SRC was the source term

for small ions from processes other than cosmic-ray 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 10-19 C). The first term represented the net charge due to

the difference in small-ion 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 10-12 C2 / ( N m2)). Then the electric-field 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 pre-strike

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, two-dimensional 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 three-dimensional 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 bulk-flash 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 two-dimensional 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 electric-field 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 electric-field magnitude fell below a certain threshold (150 kV/m)

at the locations of the channel-termination 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 e-ax2, 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 bulk-lightning 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 electric-field 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 three-dimensional domain to represent more
realistic electric field and thunderstorm dynamics.

In their work, the stochastic lightning model (SLM) was an application of
the Wiesmann-Zeller 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 near-neutrality

(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 three-dimensional lightning-discharge 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 A-iU(UTA-1U)-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*- + V-aE + 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 fair-weather 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 10-9.
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(S-t+) 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 (s-t+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 e-Y(At) / \n+i (V. J)n+)
e ( ( ) ) + ( /J)

Define
d 1 e-y(At)7
-y
Then we can rewrite the above approximation as follows :
(/-n+ d e-Y(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_(s-t+1)y^z ds, (3.10)


and it can be approximated by

Ao++-=AAA -d (1-A) 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,k-1
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 + Vi-l,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 + Oij-lk
2' Oi,j,k h- 2

and
i,,k i,j,k+1 20i,j,k + i,j,k-l
)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 matrix-vector 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 3-1.



P P


Figure 3-1: 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 Pr-1 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

super-diagonal and sub-diagonal 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 n-1 n

Figure 3-2: 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 3-2. From Figure 3-2, 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 C-1D 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 3-3 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) = e-kty(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 = (e-kIt)n y(O). The factor (e-kAt)" 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) = 10-13.40+0.1600z-0.005175z2 +0.00009077z3

For z > 27.828 km:

a(z) = 10-12.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 -






10-14 10-"1 10-12 10-11 10-so 10. 104 10-7 10-8 10" 10-,4 103 10"2 10- 10
Conductivity, S m-1

1 I I i i t -l l1 lI
10s 102 101 111 10-1 10-2 103 10- 10 5 10-6 10-' 10- 104 10-10 10-11
Relaxation time, s

Figure 3-3: Electrical conductivity or and corresponding relaxation time T = Eo-i,
where e = 8.85-12 F m-1, versus altitude under a variety of geophysical conditions.
LL, low latitude, "wavy"; MLPS, mid-latitude typical night (high-latitude, quite);
AZTDN, auroral zone, typical disturbed night; MLD, mid-latitude day, quite;
MHL, mid-high-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 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 3-3 of this thesis, with permission of the

authors. In Appendix A.2, we give the following approximation based on Figure

3-3:

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 3-1: Coefficients of the polynomial p(z)

i a,

0 -1.316 x 101
1 6.351 x 10-2
2 8.682 x 10-4
3 8.674 x 10-6
4 5.781 x 10-8


Table 3-1. 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 3-3, 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

max-norm 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 e--Y(At) z
d =
= 1 e--y^A L s


= e-7 At,


where 0 < < 1. Since e-7'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 -e-1 -+(At)) ( A
v0 i 0 i+1

+ (1- e--i(At) -) (1 e- 7i(At)) --
( i ri+1
= (1 e-'Y(At)) --) -- + (e--At) e-A(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

) = e-imvm
Sm==-oo
for ( E [-7r,7r], and (-7r) = (7r). If the spacing between the grid points is h, by

changing variables we have



()= m-0oo

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 e-imhO(O.+l 20m + Om-1)h
V2=s e

000

S m=(--oos
= (2 cos 0 2) 1 e-imhmh
v m= -oo
= (2 cos 0 2) m,










2.



1 -oo
1 e-imh (m+ Om-1
M=00



= 2isin0 1 e-i 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 |p|2 < 1, since Ak = e-k(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 p|l2 < 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 4-1: 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 4-1. 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 4-2. 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 4-2: 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)d-r 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= e-Y(s-+1)-zds -z e-(s-t+')ds
ft. 8z 09z
tn l n7z

=I e7 -tn+1) -yds d Pl
et -?S-tn )l ) rhA"ds


= 7z(1 e-7(at))/7. Apply Lemma 4.2.1 to obtain the following estimate:
= **eT+ \ i) ] et+l
R1 = e s-t y)7P ds -7z- f e-(s-tl)ds


< C(At)3.


(4.2)


Define R2 as follows:

R2=d -\z


I


d (1-p) +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 + I-Uijk]
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 matrix-vector 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 max-norm error estimate when the solution to

(3.7) is substituted into Equation (4.6). Define the max-norm for a vector x as


|ix|oo = max |xil,
1
where x = (X1, x2,.. ., x)T, and define the max-norm for a matrix M as
n
in
l||oo = 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 e-p(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 = C-1D ('n -_ a) C-1Tn, (4.8)


where C = A AB and D = AA + h(1 tp)B. Following from Equation (4.8), it
is found
n-1
r" = M" (o0 _0) ZMIC-1 n-_-1, (4.9)
1=0
where M = C-1D. Note that we prove IMIoo = 1 + O(At) in Section 5.3 for 0 <
pi < 1 if 7 is small. Therefore, |IM|oo can be written as I|M|Ioo < 1 + cAt < ecat.
Taking the max-norm of Equation (4.9), then we have

11 e"no < e"^ct |I1 6o11o + IIC-1lloo 11.r-1 -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 Max-Norm Bound for Matrix C-1 = (A ,/B)

In this section we obtain a max-norm bound for C-1 = (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 max-norm of the matrix C-1 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 M-1 > 0, i.e., all

elements of M-1 are positive.

Lemma 5.1.2 If h is sufficiently small, then the matrix C is invertible. Moreover,

IIC-11 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, C-1 > 0, i.e., the elements of C-1 are positive
by Lemma 5.1.1. Hence, IC-'Ioo = FIC-11Hoo, where 1 is the vector with every

component equals to 1. Furthermore, IIC-I[ _< 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 e-k ) (A7/7)). Let

Nh = {(i,j,k) : (xi, yj, Zk) E }.

Next, we show that IC-11loo < 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 + i-l,j,k + Wi,j+l,k + 2i,j-l,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 _< llC-1Fllo = 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 max-norm
of the amplification matrix C-1D.

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-) [n-i+1,2(n-i+l),-- ,i(n-i+l),i(n-i),i(n-i-1), ,2i,i].
n +1









Proof: Note that we compute the (i, j)th entry of the matrix A1- 'A by the
formula

(Ai-A1)i,,j = -(Aj-1)i,k(Ai)kj
k=1
= (a-1),jjajj, + (a-1)i,jaj,j + (a-1)i,j+laj+,j,

where (a-1)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
(A1-1A)i,j = 1 [(j-1)(n-i+1)(-1)+j(n-i+1)(2)+(j+1)(n-i+1)(-1)] = 0,

ifj= 1,2,.---. ,i-1, and

1
(Ail-1A)i', = -- [i(n-i-(m-1))(-1)+i(n-i-m)(2)+i(n-i-(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 two-point 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+ 1-i) 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=i-lorj=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

(AI-IT)i,j = E(Aj-1)i,(W)kj
k=1
= (a-1)i,j-tj-l,j + (a-1)ij+ltj+l,j

= (a-)i-l (a -)i,j+1


we have


(Ai-T)i,j <





Then, the infinity norm for


IIAr1'TIK


-2(n+l-i) ifj=l,2,.i -1

-n + 2i 1i)
if 3 = i

-2, if j =i+1,i + 2,-- ,n-1
n+1
2i
2i if j = n
n+1'

the matrix Ai-'T is


n
= maxE (A-1'T)i,3
j=l
[2(n- i + ~1) -n+2i-1 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)

< n-l,

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) e-7(), and 7 e C2. Then I|l = O(At h) and

I +1i = O(At. -h2).

Proof: Both qualities are proved by the mean value theorem:

(1)

fi g-7.+l(At) _g-7.(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+ = e--7i+(At) e-7i(At) (e-i+2(At) e-,+(At
= e-i(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 Max-Norm Bound for Matrix C-1D
The max-norm bound for the matrix C-1D, 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

C-1D = (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- A-1B A-BA+ I A-'B A-AA-AA)
2 ( 2 /\ 2









Hence,
<-1 + Ah II A-BI +, IIA-I(AA AA)||oo
21 IIA-IBII 1 IIA-IB|lo

whenever /,lA-1'Boo < 1, and where A = max A,. Our goal in this section is to
show that IC-Dlloo = 1 + O(At).
5.3.1 Max-Norm for Matrix C-1D in 1-Dimensional Space
The preliminary results from Section 5.2 are used to show that IIC-1DIk =
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

h|1Al-'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':

(AI-B'),,j = dj-1 ((a-1)i,j-1 (a-1)i,+l)

where (a-1)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(n-i+1) ++ |-n+2i-l d,

+-1 (di + d2 +" )+ dn-1)
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 A-iB" :


-(a-1)i,2d2

= (a-1)i,(dj-1 dj+)

0


-2(n-i+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,- ,n-i-1)

if = n


where Id'I = max Idj-1 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

+--- [(n-i)+(n-i-1)+...+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

IAx-I(A1A1 AixAi1)1 = O(At).

Proof: To prove this Lemma, we compute the max-norm 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 fl-f2
-f 0 f3 0 0 f2 -f3
= -f2 0 f4 + 0 0 f3-f4




= F'+F".

The same process in Lemma 5.3.1 is used to compute A1-iF' and Ai-'F"
respectively.
(1) compute A-'F' :

(A-1 F'),, = ((a-1)i,j-l (a-)i,j+,)
-2(n-i+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


(a-1)i,j_-(fj1- fj) if j 0 1


0

(j 1)(n- i + 1)(f f)
n+1 + k
i(n -ni + 2 k) (f-1 f3)


if j = 1


if 2 < j
if j=i + k (k= 2,3,.- ,n-i)


Hence, for any i,


SI(A 'F")i,j < --- (1 +2+---+ i) f'
-n+1
j= 1
+---[(n-i) + (n-i- 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(A-1 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

IIC-1DIl _< 1 + O(At),

where C = A1 Bx1 and D = A1Ax.

Proof: Note that

_lh AlA- Bn llo ,|Ax-i(AxAi AxAi)|oo
IC1DI < A+ 1 h A-1 + AA-) (5.10)
21- llA Blloo 1- llAix-1Bllo

whenever |IIAi-1B13II < 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 Max-Norm for Matrix C-1D in 2- and 3-Dimensional Spaces

Now we extend the results from the 1-dimensional case to 2- and 3-dimensions.
First, we consider the 2-dimensional 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


hA2-1B2 =

I+2A1-1 -A,-1

-A-'1 I+2Ai-' -A,-1


h


-A-' I+ 2AI-1 -A,-'

-A-1 I+ 2Ai-1

AI-1 B,

Ai-1 B1


x




Ai-1 B1

and we can rewrite it as


hA2-1B2 =

Al-1 B1

AI-1 Bx


hS2 (5.11)


AI-1









where S2 is the matrix


I+2A1-1

-A1-1


I+ 2A,


-A-1 I+ 2Af-1 -A,-1
-A,-1 I+ 2A-1


Thus,


h|IA2-1B21. -< hi AH-A BIjlIjS2|oo.

By Lemma 5.3.1, we have IAi-1B1Io = O(At/h). Hence,

hIA2- B21loo < O(At)||S2||oo.

Note that S2 is difficult to analyze using the max-norm. Hence, S2 is evaluated on
a computer to get a numerical value for J|S21oo. For example, if A1 is 10 x 10 and
m = 10, so that S2 is 100 x 100, then ||S21|oo 1.6687. If A1 is 20 x 20 and m = 20,

so that S2 is 400 x 400, then J|S2|1K| 2.0794. If A1 is 40 x 40 and m = 40, so that
S2 is 1600 x 1600, then |IS2||oo ; 2.5041. Therefore, we have


hIIA21-B2[[o, = O(At).


-A,-1


(5.12)










Also, we have


A2-'(A2A2 A2A2)
-1
I + 2A1-1 -A1-1

-A1-1 I+2A1- -A,-1




-A-1' I+2A1-1

A1- A1A1 A1A1

x


Aj-1 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 HIC-1Dl, = 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


A2-I


and


A-'(AA AA)=

A2-1


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 max-norm, we still need to evaluate the factor S3. Using a computer

to evaluate the max-norm of S3, we have: if As is 5 x 5, so that S3 is 125 x 125,

then |IS3I|oo 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| A2-1B2|lootlSaloo

= O(At)


IA-'(AA- AA)|II < IA2-1(A2A2- A2A2)IUISa|3oo

= O(At)


by Equation (5.12) and Equation (5.13) from the 2-dimensional space. Finally, we
have IIC-'Doo = 1 + O(At) from Equation (5.10).

5.3.3 Numerical Results for Spectral Radius of C-1D

Next we use the approximate value of -y (Section 3.3) to calculate some

numerical results for spectral radius of C-1D, (Table 5-1, Table 5-2 and Table

5-3.), 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 5-1: Spectral radius for C-1D in one dimensional space

Matrix dimension (n x n) Spectral radius for C-1D

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 5-2: Spectral radius for C- D in two dimensional space

Matrix dimension (n x n) Spectral radius for C-1D

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 5-3: Spectral radius for C-1D in three dimensional space

Matrix dimension (n x n) Spectral radius for C-1D

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 3-1) 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 = Ape-wir-sp12 ANe-wr-sN2,


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 e-w2dV = 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 6-1, to obtain


















(positive


z-axis (altitude)


P


Figure 6-1: Cylindrical shell of the volume integral


e-w2 dV


ew-h(p2+z2)2rpdzdp
J z-ht p20


= ( 27rpe-,p2dp)
\Jp>0 /


(\z>-hl


e-WZdz)


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 -w-SN 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

z-axis, 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 m-1 [7]. We use the uniform mesh in all x-, y-,

and z-directions with mesh size 1000 m, and the time step is At = 1/16 s. As

shown in Table 6-1, the breakdown time is calculated based on different values of w.

Table 6-1 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 6-1: 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 6-1. 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 6-1, 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










z-axis 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 y-directions 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 z-axis (mesh size

= 200 m for 0 < z < 20 km). The mesh size is 1000 m for both in x- and y-

directions. Figure 6-2 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 10-4 between these two domains. Therefore, we

cannot see any difference for the potentials in Figure 6-2. Figure 6-3 and Figure

6-4 are the potentials near the negative and positive charge centers respectively,

along the z-axis 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 6-2: Potentials (i(0, 0, z), 0 < z < 100 km) in two domains, Q1 and Q2


) 5000 5100
Altitude in meters


Figure 6-3: 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 6-4: Potentials (0(0, 0, z), 9.6 km < z < 10.4 km) near the positive charge
(10 km) in two domains, Q1 and Q2


x-axis in meters


2.5
x 10'


Figure 6-5: Potentials (O(x, 0, 5) and (x, 0, 10),
x-axis at the altitude 5 km and 10 km


-25 km < x


< 25 km) along the










Further, we find the potentials (see Figure 6-5) along the x-axis at the altitude 5

km (the negative charge center) and 10 km (the positive charge center) when y = 0.

In Figure 6-5, 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

(Crank-Nicholson 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 z-axis (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

6-2 and Table 6-3. Second, the relative errors are estimated by the ratio

IIhAt al /2811oo


as shown in the fourth column of Table 6-2 and Table 6-3.













Table 6-2: Relative error in the potential at t = 65 seconds using the maximum
norm for Crank-Nicholson scheme (pt = 1/2).

Time step maximum norm maximum norm error Crank-Nicholson
(s) potential with respect to relative error
time step 1/128 s (p = 1/2)

1/2 564431185.91904 80.21113 1.42110 x 10-7
1/4 564431163.61478 20.30029 3.59659 x 10-8
1/8 564431157.94013 5.11413 9.06068 x 10-9
1/16 564431156.50102 1.27381 2.25681 x 10-9
1/32 564431156.13736 0.30471 5.39852 x 10-10
1/64 564431156.04677 0.06052 1.07192 x 10-10
1/128 564431156.04677 0.00000 0.00000


Table 6-3: 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 2-111I. x 10-
1/4 564426772.69223 17375.01373 3.07832 x 10-"
1/8 564428962.54654 8411.09363 1.49019 x 10-5
1/16 564430058.82113 3926.08120 6.95582 x 10-6
1/32 564430607.30153 1682.80248 2.98141 x 10-6
1/64 564430881.63001 560.96757 9.93864 x 10-7
1/128 564430881.63001 0.00000 0.00000












In Table 6-2 which corresponds to the Crank-Nichlinl,,Ti scheme, if the time

step At is multiplied by 1/2, the relative error is multiplied by 1/4. Hence, the

error for Crank-Nicholson scheme (pj = 1/2) behaves like O(At)2. In Table 6-3

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 Crank-Nicholson 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 max-norm of the potential grows exponentially in t. The plot of loglo(|I|AtIloo)

versus t appears in Figure 6-6, 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 6-6: 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 Crank-Nicholson 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 z-axis:

Case 1: 30 uniform grid points in the x- and y-directions, 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 y-directions, 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 y-directions, 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 y-directions, 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 z-axis 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 = En-E2n

= 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 6-4: 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 6-4, 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 6-4. 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 6-4. 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 6-5 respectively. The best estimate for the constant C is 1.45 x 108 from
Table 6-5, while the best estimate is 1.19 x 108 from Table 6-4. The relative
difference between them is 0.22. From the third column of Table 6-5, we can see
that the constant p is getting close to the theoretical value 2.






























Table 6-5: 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 max-norm. 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 A-1) 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) = 10-13.40+0.1600z-O.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) = 10-12.88+0.06718z.

Figure A-1 shows the graph of the approximation for the conductivity function.

A.2 Approximation from Figure 3-3

Thirteen points (see Table A-2) in Figure 3-3 are picked to approximate the

conductivity function. The conductivity is chosen as the average of shaded region

in Figure 3-3 for the altitude above 30 km. The advantage for using this figure













Table A-l1: Four points chosen from Figure 9.6 [27]

Point Height logo0 (conductivity)
(km) (1/ohm-m)

P1 2.7432 -13
P2 12.954 -12
P3 20.7264 -11.5
P4 27.82824 -11


log Joconductivity), (Ohm-meter)


Figure A-1: Electric conductivity as a function of altitude approximated from
curve B in Figure 9.6 [27]










to estimate the conductivity is that Figure 3-3 demonstrates a more completed

conductivity data than Figure 9.6 in [27].

Table A-2: Thirteen points chosen from Figure 3-3


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

A-2.









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 A-2: Electric conductivity as a function of altitude approximated from
Figure 3-3















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(I-1,J,K) 2.*X(I,J,K) + X(I+1,J,K)) +

C2*(X(I,J-1,K) 2.*X(I,J,K) + X(I,J+1,K)) +

A(K)*X(I,J,K-1) + 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 K-1 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.29679-29688.

[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, Prentice-Hall, Englewood
Cliffs, NJ, 1988.

[4] W. W. Hager, Updating the inverse of a matrix, SIAM Rev., 31 1989,
pp.221-239.

[5] W. W. Hager, A discrete Model for the Lightning Discharge, Journal of
Computational Physics, 144 1998, pp.137-150.

[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.193-217.

[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.3542-3558.

[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.5661-5675.

[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.849-854










[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.5865-5884.

[13] H. W. Kasemir, A contribution to the electrostatic theory of a lightning
discharge, Journal of Geophysical Research, 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
three-dimensional 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.1033-1036.

[18] F. Pockels, Uber das Magnetische Verhalten Einger Basaltischer Gesteien,
Ann. Phys. Chem., 63 1897, pp.195-201.

[19] F. Pockels, Bestimmung Maximaler Entladungs-Strom-Stiirken aus Ihrer
Magnetisirenden Wirkung, Ann. Phys. Chem., 65 1898, pp.458-475.

[20] F. Pockels, Uber die Blizentlandungen Erreicht Stromstdirke, Phys. Z., 2
1900, pp.306-307

[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.779-800.

[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.352-364.

[25] T. Takahashi, Determination of lightning origins in a thunderstorm model,
Journal of the Meteorological Society of Japan, 65 1987, pp.777-794.

[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, Prentice-Hall, 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.1770-1773.

[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.6054-6070.

[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.555-574.

[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.73-115.

[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.833-851.

[34] C. L. Ziegler, D. R. MacGorman, J. E. Dye, and P. S. Ray, A model
evaluation of non inductive graupel-ice charging in the early electrification
of a mountain thunderstorm, Journal of Geophysical Research, 96 1991,
pp.833-851.

[35] C. L. Ziegler, J. M. Straka, and D. R. MacGorman, Simulation of early
electrification in a three-dimensional dynamic cloud model, 11th Interna-
tional Conference on Atmospheric Electricity, Commission on Atmospheric
Electricity, 1999, pp.359-362.















BIOGRAPHICAL SKETCH

Shu-Jen Huang was born in Taiwan in 1970. She was awarded a Bachelor of

Science degree in applied mathematics in 1992 from Fu-Jen University, Taipei,

Taiwan. She worked as a teaching assistant at Fu-Jen University for 3 years after

her graduation. In 1995, Shu-Jen 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 Go-palakrishnan
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


2611-111'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 UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EO5R0OPSG_127EA7 INGEST_TIME 2011-11-02T16:58:32Z PACKAGE AA00004682_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES



PAGE 1

MULTISCA LE DISCRETIZATIO OF ELECTRIC-FIELD EQUATIO S B y SHU-JEN 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 Shu-Jen Hu ang

PAGE 3

To Deng-Shan 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 eng-S 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 -in-law 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 Max-Norm Bound for Matrix c -1 = (A ~Br1 38 5.2 Estimates for A1 l . . . . . . . . . 4 0 5.3 Max-Norm Bound for Matrix c -1 D . . . . 4 6 5 3 1 Max -No rm for Matrix c -1 D in I-Dimensional Space 47 5 3 2 Max-Norm for Matrix c -10 in 2-and 3-Dimensional Spaces 5 2 5 3 3 Numerical Results for Spectral Radius of c-1 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 = ECY-1 where c = 8.85 -12 F m -1 versus altitude under a variety of geo physical conditions LL low latitude, wavy ; MLPS mid-latitude typical night (high-latitude quite) ; AZTD auroral zone typical disturbed night ; MLD mid-latitude day quite; MHL mid-high 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 x-axis 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 e-7(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 = A--B 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 ELECTRIC-F IELD EQUATIONS By Shu-Jen 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. Max-norm 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, Solid-state 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 three-dim 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 high-d 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 mass-weighted 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 two-dimensiona 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 right-hand side of Equation (2. 1) was the advection term. The second term was the fallout term and it was equal to zero for cloud-water 4

PAGE 18

and cloud-ice charge The third term was the e ddy-mixing 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 small-i 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 ion-ion 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 10-19 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 ll-ion 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 10-12 C2 / ( m2)). Then the e l ectri c-fie 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 three-dimensional 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 bulk-flash 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 two-dimensional 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 electric-field 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 electric-field magnitude fell b e low a certain threshold ( 1 5 0 kV /m) at the locations of the channel-termination 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 e-0x2 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 bulk-lightning 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 lectric-fie 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 three-dimensiona 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 Wiesmann-Zeller 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 near-neu 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 ee-d im ensiona l li ghtning-discharge 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+ a-E + 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 o-E + '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(S-tn+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 -e--y(.6.t) d= ----,z -1 Then we can r e write the above approximation as follows : (fJ'ljJ)n+ 1 e--y(.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 1-e--y(.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 ub-di 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 c-10 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) = e-kty(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 = ( e-Mtr y(0). The factor ( e-k.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 ) = 1O-13.40+0.1600z-0.005175z2+0.00009077z3 For z > 27.828 km: a( z ) = 10-12. 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 10-1 10-13 10 -12 10-11 10 -1 0 109 10-8 101 10-8 10"'5 10~ 103 102 10-1 10 Conductivity, S m -1 103 102 101 1oO 10 -1 10 -2 103 10_. 10 -5 10 -6 10-1 10--a 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 3-1. 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 e-1(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 e-im{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 e-imh0('1/Jm+l -'l/Jm-1)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::(cos0j-l)+ 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 = e--r,.,(~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, 300-250 -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(1-2 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(S-tn+1),z-dS /Z e'Y(S-tn+1lds t,. a z a z tn l tn+ I B'lj; 1 e-1'(6.t) ( B'lj;) n+ e'Y(S-tn+l)fz-dS -----,z t n a z 'Y a z l tn+I e'Y(S-tn+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 matrix-vector re lation of t h e form h 1 -e-1'(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 max-norm error estimate wh e n the solutio n to (3. 7) is substituted into Equation ( 4.6). D efine t h e max-norm for a vector x as where x = ( x1 x2 . Xnf, a nd d efine the max-norm 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 -1rn-l-l, 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 max-norm 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 max-norm 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 IIC-11100 11c-1F 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 11c-11100 :::; 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 + Wi-1, j k + Wi, j + I k + Wi,j-1, 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 < IIC-1Flloo = 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 max-norm 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 [n-i+l ,2(n-i+l), ,i(n -i+l),i(n-i),i(n-i-l), 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 (a-1)i,j-1Uj-l, j + (a1 ) i,jaj, j + (a-1)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)(n-i+ 1)( -1) + j(n-i+ 1)(2)+(j + l)(n-i+ 1)(-1)] = 0 n+l if j = 1 2, i-1, and 1 ( A1 -l A1) i1 = --[ i(n-i-(m-l))(-l)+i(n-i-m)(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+l-j) 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 (a-1)i,J-1tJ1 J + (a-1)i,J+1tJ+1, 1 -2(n + 1 i) ----, if j = 1 2 , i-1 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 ---i-1 +-----+-n i [2(n-i+l)( ) l-n+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) e--r;(~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 ( e--r;+,(6t)'Y~+1 e--r;(6tl-y:) (6.t)h ( e-7 ; (6t ) (-1i(6.t) )1ih + e-t(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 Max-Norm Bound for Mat rix c-1 D The max-norm 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 11c-1011 < + ~!!. IIA1 Blloo IIA1 (AA AA)lloo oo_ 21-~IIA-1 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 11c-1Dlloo = 1 + O(flt). 5.3 1 Max-arm for Matrix c-1 D in 1-Dimensional 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 (a-1)i,j is the (i, j)th entry of A1 i H ence, Therefore, for any i n L /(A1-1B')i,j/ j=l 0 if _j = 1 2(n + i -l) d if 2 _< J _< ,,-1 n + 1 j-l -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 + + dn-1) < 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 J-1 1+1 = i(ni -k)(d -d ) n + l J-1 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(n-i+l)d n-i+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 dj-1 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 h-h 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 -(a-1k1+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.fi-11 ) + 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) (f-f-) 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 < n-i+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 fj-1 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 Max-Norm for Matrix c -1 D in 2-and 3-Dimens i onal Spaces ow we extend the resu lts from the 1 dimensional case to 2-a nd 3 dimensions First, we consider the 2-dimens 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 max-norm. 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 11c-1 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 -A-1 2 A 1 2 A21 I + 2 A21 A 1 2 56 -1 -A1 2 I + 2 A21 To discu ss the max-norm we still n ee d to eva lu ate t h e factor Sa Using a computer to e valuate the max-norm 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 2-dimensiona 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 5-1, Tab l e 5-2 and Tab l e 5-3.), 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 3-1) 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 e-t 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 z-axis 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 z-direction 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 6-1 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 6-1. 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 6-1 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 z-axis 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 z-axis ( 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 z-axis 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 x-axi 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 x-axis 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 x-axis 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 z-axis (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 6-2 and Table 6-3. Second the relative errors are estimated by th ratio ll4tt::.t -1;12slloo ll1;12slloo as shown in the fourth column of Table 6-2 and Tabl e 6-3

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 10-7 20.30029 3 59659 X 10-s 5 .11 413 9 06068 X 10-9 1.27381 2.25681 X 10-9 0.30471 5.39852 X 10-lO 0.060 52 1.07192 X 10-lO 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 10-5 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 10-6 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 6-2 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 Crank-Nicholson scheme ( = 1/2) behaves lik e O(b..t)2. In Table 6-3 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 max-norm 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 z-axis : Cas e 1: 30 uniform g rid points in the x and y-direct 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 y-direct 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 y-direct 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 6-4: 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 6-4, 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 6-4. 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 6-4. 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 600z-0 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 A-2) 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 m-m ) 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(I-1 J K) 2.*X(I, J K ) + X(I+l, J K )) + C2 (X(I ,J-l, K) 2 X ( I J K ) + X(I,J+ l K )) + A(K)* X ( I J K-1 ) + 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 K-1 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 ice-Hall Englewoo d Cliffs J 1988. [4] W. W H ager Updating the inv e r se of a matrix SIAM Rev. 31 1989 pp.221-239. [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 three-dimensiona 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 Entladungs-Strom-Stark 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 Prentice-Hall, 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 three-dimensional 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 Shu-Jen 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 u-J en Un i versity, Taipei, Taiwan. Sh e worked as a teaching ass istan t at Fu-J en Uni versity for 3 years after h e r graduation In 1995 Shu-Jen 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