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