I 

NASA Contractor Repo+44&l.711 I # 

IWE REPORT NO. 88-50 

ICASE 

MULTIGRID METHOD FOR STABILITY PROBLEMS 

Shlomo Ta'asan 

Contract Nos. NAS1-18107 and NAS1-18605 

September 1988 

(HASA-CR-185825) MULTIGRID MBTHOD FOB 

98 e CSCL 128 

'Bi 8 9- 2 7 4 20 

STABILITY PROELKCIS Final Report (ICASE) 

Unclas 

G3/64 0199021 

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 

NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 

National Aeronauticsand 

Space Administration 

Ln0kY-w 

Hamptm.Virginia 23665 
Multigrid Method for Stability Problems 

Shlomo Ta’asan * 

Department of Applied Mathematics and Computer Science 

The Weizmann Institute of Science 

and 

Institute for Computer Applications in Science and Engineering 

Abstract 

The problem of calculating the stability of steady state solutions of differential equations is 

treated. Leading eigenvalues (i.e., having maximal real part) of large matrices that arise from 

discretization are to be calculated. An efficient multigrid method for solving these problems is 

presented. The method begins by obtaining an initial approximation for the dominant subspace 

on a coarse level using a damped Jacobi relaxation. This proceeds until enough accuracy for 

the dominant subspace has been obtained. The resulting grid functions are then used a8 an 

initial approximation for appropriate eigenvalue problems. These problems are being solved 

first on coarse levels, followed by refinement until a desired accuracy for the eigenvalues has 

been achieved. The method employs local relaxation on all levels together with a global change 

on the coarsest level only, which is designed to separate the different eigenfunctions as well as 

to update their corresponding eigenvalues. Coarsening is done using the FAS formulation in a 

non-standard way in which the right hand side of the coarse grid equations involves unknown 

parameters to be solved for on the coarse grid. This in particular leads to a new multigrid 

method for calculating the eigenvalues of symmetric problems. Numerical experiments with a 

model problem demonstrate the effectiveness of the method proposed. Using an FMG algorithm 

a solution to the level of discretization errors is obtained in just a few work units (less than lo), 

where a work unit is the work involved in one Jacobi relaxation on the finest level. 

*Supported in part by the National Aeronautics and Space Administration under NASA 

Contract Nos. NAS1-18107 and NAS1-18605 while the author was in residence at ICASE, NASA 

Langley Research Center, Hampton, VA 23665, and in part by the Air Force Office of Scientific 

Research, United States Air Force under Grant AFOSR-86-0127. 

i 
1 Introduction 

Stability analysis plays an important part in hydrodynamics (41 and in other physical sciences [6]. 

Basically, one is interested in the growth (decay) rate of infinitesimal perturbations to steady state 

solutions. In rectangular geometries with simple enough steady state solutions one can use Fourier 

transform techniques to reduce the resulting stability problem to an eigenvalue problem for an 

ODE with parameter(s) [4]. General geometries or complicated solution structure do not permit 

this and eigenvalue problems in two and three space dimensions need to be solved. 

A problem related to the stability of a given steady state solution is the problem of transition, 

when a nonlinear problem with a parameter is given and the value of that parameter at which 

steady state solutions loose their stability is required. In fact, in cases of marginal stability this is 

the relevant problem instead of finding the stability of a given steady state solution, since in such 

cases the discretization (in space) might have shifted the eigenvalue from, say, a stable region to 

an unstable one. 

Steady state solutions are often obtained by time marching algorithms. In such cases the 

stability to small perturbations can be studied during the solution process. However, efficient 

multigrid solvers for steady state problems do not involve usually any time marching. The solver 

is designed for solving the steady state problem directly. This allows obtaining solutions which are 

physically unstable. A way to recover the information about the stability of a given steady state 

solution is therefore needed. 

One way of obtaining that information is to use time marching for the time dependent problem 

obtained by linearizing the original problem, around the steady state solution, using a random 

initial condition. The discretization has to be such that it preserves stability. That is, if a certain 

component is decaying oscillating or growing, the same should happen for the discretization. CrankNickolson time stepping has this property. 

The effectiveness of such a method depends on the closeness of the dominant eigenvalues to 

the imaginary axis. The closer they are the more time steps are needed to determine the growth 

rate of small perturbations to the steady state solution under study. For symmetric problems the 

situation is much better than the non-symmetric ones, since the time steps can be taken as larger 

and larger for larger times. In the non-symmetric case the time stepping depends also on the 

imaginary part of the dominant eigenvalues, and this can be much larger than the real part. In 
2 

other words, the time stepping can depend on the oscillation rate rather than the decay rate which 

may be much smaller. An improvement of the above method is to use time marching until the 

subspace of dominant eigenvalues has been determined to a desired accuracy, then obtaining from 

it the set of eigenfunctions and eigenvalues as in ['I]. 

The above approach is inappropriate when the steady state solution has been obtained by an 

efficient multigrid solver since it is more expensive than getting the steady state solution using a 

multigrid method. Moreover for transition problems it can be very time consuming and a different 

approach should be used. 

In this paper an alternative way for calculating the stability of steady state solutions by solving 

directly, using a multigrid algorithm, an appropriate eigenvalue problem. This approach could be 

modified for transition problems. The difference is only in a global step (explained later) that is 

performed on the coarsest level. The exact treatment of transition problems by multigrid, together 

with some physical examples will be described elsewhere. 

Since the eigenvalue with maximal real part is required, the multigrid processes must not 

contradict this requirement. The relaxation should be such that it will relatively damp all the nondominant eigenfunctions. Jacobi or a damped Jacobi. relaxation achieves this purpose. Kaczmarz 

on the other hand, should not be used since it damps more the eigenfunctions whose eigenvalues are 

away from zero. As in [2] a single grid method which involves a relaxation method and a step for 

updating the eigenvalue is constructed first. In non-symmetric problems the relaxation (Jacobi. ) 

couples the eigenfunctions that correspond to complex conjugate eigenvalues. Therefore, the step 

for updating the eigenvalue, (the global step) performs also a projection onto one of the subspaces 

of the two conjugate eigenvalues, say, with non-negative imaginary part. 

When the dominant eigenvalue is not separated enough from the others, in the sense that their 

real parts are close, one has to include more eigenfunctions and to solve for the dominant eigenspace 

which is of higher dimension. 

The coarsening strategy for the multigrid method is done using an FAS formulation since the 

eigenvalue problem is a non-linear one. However, since the global change done on the fine level 

affects the coarse grid right hand side significantly in this problem, the standard FAS should be 

modified to include this change of the right hand side. The resulting FAS is such that the right 

hand side of the coarse grid equations involves some parameters to be solved on that level. The 
3 

resulting algorithm has fast algebraic convergence rates. Moreover, it leads to a new multigrid 

method for symmetric eigenvalue problems. 

The algorithm is performed in an FMG version, where the problem is solved first on coarse levels 

using the approximation obtained there as an initial guess for the solution of finer grid problems. 

On each level a fixed number of basic multigrid cycles are performed, yielding a solution to the 

level of discretization errors in just a few work units (about lo), where a work unit is the work 

involved in one Jacobi relaxation on the finest level. 

Numerical experiments with a two dimensional model problem are given. The results demonstrate the effectiveness of the proposed method. Solution to the level of discretization errors is 

obtained with the 1-FMG algorithm (which uses one basic multigrid cycle per refinement and cost 

about 10 WU). 

2 The Problem 

Let A be a linear real valued operator. We wish to compute its eigenvalues with maximal real part. 

That is 

A!# = A@ with Real(A) mazimized. (2.1) 

In general, for asymmetric operators, the eigenvalues and eigenfunctions are complex. The above 

eigenvalue problem can then be written in terms of real quantities as 

or as 

depending on whether the dominant eigenvalue is real or complex. In this notation X = Real(A), 

p = Imag(A), + = Real(!#), 3 = Imag(9). 

Since the operator A is real its eigenvalues come in pairs of complex conjugates. Hence, if (A, p) 

is a solution to (2.3), also (A,-p) is a solution. Moreover, each of these eigenvalues is double. If 
4 

(4, t,b)T corresponds to (A, p) then so does (9, -+)T . Similarly, ($J, 4)T and (4, -$J)~ correspond 

to (A, -4. 

3 Single Grid Algorithm 

One way to get functions in the subspace spanned by the eigenfunctions with eigenvalue of maximal 

real part is to simulate a time dependent process, that is, solving approximately the initial value 

problem 

ut = AU 

U(0) = random 

See (I. Goldhirsch et a1 1987). The dominant part of the solution at large time is a linear combination of the desired eigenfunctions. This leads us to the following. 

3.1 Initial Approximation 

The initial approximation for the eigenvalue problem is obtained from the following process applied 

repeatedly, 

(a) Relax using a damped Jacobi. relaxation AW = 0. 

(b) W + W/llWll* 

The stopping criterion for the process is as follows. Let 41 = W, tjj = A4j-1 (j > l), and let 

dn) be an nxn matrix whose (ij) element is < 4i,4j >. Define 1 to be the minimal positive integer 

for which det + 0. The dominant subspace U is of dimension 1, and the process can be 

stopped when det D('+l) is small enough. The dimension of U determines how many eigenfunctions 

are to be simultaneously solved for in order to obtain fast convergence. 

After an initial approximation for the dominant subspace has been obtained a solution of the 

eigenvalue problem can be started. The algorithm consists of two steps which are performed successively. The first is a local relaxation process on equation (2.2) or (2.3) depending on the eigenfunction. Its role is to damp the error in the components corresponding to non-dominant eigenvalues. 

The analogy of damped Jacobi relaxation with explicit time stepping suggests that using it with 

an appropriate under-relaxation parameter may serve our purpose. The local relaxation is then 

followed by a step referred to as the global step which is explained in detail next. 
The Global Step 

(a) dim U = 1. In this simple situation X can be updated according to 

(b) dim U = 2. Assume for simplicity of exposition that the dominant subspace is spanned by a 

pair of complex conjugate functions. In that case the approximation we get at an intermediate step 

will be a linear combination of the four eigenfunctions mentioned in section 2, and the eigenfunctions 

corresponding to other A’s. The latter converge to zero in the relaxation process (relative to these 

four), while these four corresponding to the maximal X are coupled by the relaxation process. If our 

approximate solution is in the subspace spanned by the two eigenfunctions corresponding to one of 

the p’s, say, the non-negative one, which belong to maximal A, then X and p can be calculated by 

certain inner products. Since in general this is not the case, we have to perform a projection onto 

the two dimensional subspace that corresponds to one of the p’s, say, the nonnegative one. 

Let (4, $)’ be an approximate solution of (2.3) which lies mainly in the relevant four dimensional 

subspace. A projection into a two dimensional subspace which correspond to (A, p) is to be done: 

Each of the two components of an eigenfunction is approximately a linear combination of the current 

6 and $. That is, 

Substituting this in (2.3) and taking inner products of each of the two equations 4 and $, we get 

a 4x4 system given below 
6 

where C and D are given by 

and 2, = (a11,(~12)~, % = (021,~~22)~. This is equivalent to the following (complex valued) 

generalized 2x2 eigenvalue problem 

(C - pD) p = 0 (3.10) 

where p = X f ip. The solution of the above eigenvalue problem gives X = Real(p), p = Imag(p) 

and 3 = ReaZ(aJ,& = Imag(g). Under our assumption the matrix D is nonsingular since 4 and 

+ are linearly independent. Hence, the above 2x2 eigenvalue problem is a standard one, and has 

two eigenvalues. 

(c) dim U = n. 

Let U be spanned by the set of real valued functions (41,. . . ,#,,). A projection into the subspace 

corresponding to non-negative imaginary part (41,. . . , &)' t E(41,. . . , q5JT is described. Let C 

be an nxn matrix whose (ij) element is < A4i,4j > , and D is an nxn matrix whose elements are 

< +i,4j >. D is non-singular since the set of 4's is linearly independent. The projection step is 

done by solving the nxn eigenvalue problem (3.10) with the new definition for C and D together 

with the constraints ll(Er~5)i11~ - = 1, (i = 1,. . . ,n) . It gives n eigenvalues where each complex 

valued one appears together with its conjugate. With each real eigenvalue we can associate a real 

eigenfunction that solves a real eigenvalue problem like (2.2), while each complex one, together with 

its complex eigenvector, defines a pair of functions, corresponding to the real and the imaginary 

part of the complex eigenfunction associated with it, and a problem like (2.3). In an analogous 
7 

way we define a set of n real vectors, constructed from the eigenvectors of (3.10), which define an 

nxn real valued matrix whose rows are the above real n valued vectors and whose action on the set 

(41,. . . ,4,J defines a new set of functions which is a better approximation to the corresponding 

set of eigenvalue problems. 

4 Multigrid Algorithm 

From the above single grid method we construct a multigrid algorithm. The relaxation defined 

for the single grid algorithm involved a local step and a global step. The local step which uses 

the damped Jacobi relaxation has the property of damping errors that belong to non-dominant 

eigenfunctions on that grid. These functions are mainly high frequency functions which cannot be 

represented on a coarser grid. The local step serves therefore as a smoother for our problem. The 

global step updates the global quantities (the eigenvalues), as well as performing a projection of 

the solution onto an appropriate subspace that corresponds to the eigenvalues under consideration. 

This projection step is designed to control the important smooth eigenfunction to be calculated. 

In trying to construct a multigrid algorithm we have a few options in treating the global step. 

The first is to employ it on all levels of discretization. Such an approach is in contrast to a basic 

principle in multigrid; that is, global changes should be performed on coarse grids only [l]. Note 

that these global changes involve in most problems the convergence of the smooth components and 

therefore it makes sense to perform them on coarse grids. The effect of the global step on the 

non-smooth components is not important since these components should converge to zero as the 

algorithm proceeds. However, one has to make sure that the size of the non-smooth components is 

small relative to the relevant smooth ones before performing the global step, since small errors in 

the highly oscillatory components can cause large errors in the global quantities. This is another 

reason why it is preferred to perform the global step on coarse grids, where small errors in the high 

frequency functions have smaller effect on the outcome. 

4.1 Coarse Grid Correction 

For clarity assume first that U is spanned by a pair of complex conjugate eigenfunctions. In that 

case the relevant eigenvalue problem is a discretization of (2.3). The standard coarse grid equations 
8 

in a multigrid algorithm for that case using the FAS formulation are 

where 

(4.3) 

(4.4) 

Th Hh (w ) = AHIfwh - IfAhWh 

‘f = r: + llrh IIH - Ilwhlli H h2 

where wh is any fine grid function. 

The right hand side of the coarse grid equation depends heavily on the current approximation 

on the fine level, hence it would be different if a global step was carried out on the fine grid. In 

other words, the fine grid change 

would change the right hand side of the coarse grid problem to 

Since the global step may change the solution on the fine grid significantly, and as a result also 

the right hand side of the coarse grid equations, it is better if the coarse grid right hand side 

takes this global change into account. The role of the equations l14hll& = rr , and Ilt)hll& = rf 

is to prevent us from getting a zero solution to our problem. Any other equation that will have 

the same effect is equally good. Since the projection step leads to complicated equations for these 

constraints, we deviate here for these constraints, from the standard FAS. We basically ask that the 

coarse grid solution is of norm one for both the components under study. This leads to a significant 

simplification. This together with taking into account the global step as if done on the fine level 

leads us to the following coarse grid equations. 
In this equation the unknowns are +H,~H,X,~ and ~~11,~~12,~~21,~~22. The iterative method for 

solving this problem will involve Jacobi relaxation for 4H,~H keeping the rest of the unknowns 

fixed, followed by a global step. The global step requires the solution of an eigenvalue problem like 

(3.10) with the following C and D 

(4.10) 

A corresponding eigenvector defines the four CY'S as in the single grid method. Having solved this 

problem approximately for dH, @, A, p and 0111, ~~12, ~~21~~22 the interpolation to be used is 

(4.11) 

We now come to the general case. 

dim U = n. Assume that there are r complex eigenvalues and n-r real eigenvalues corresponding 

to this subspace. In this case the fine grid problem is of the form 

(Ah - r")f = 0 (4.12) 

Ah nz n 

(4.13) 

(4.14) 
10 

and 

h - r= 

Following the previous arguments, the coarse grid problem is then 

Hh (AH - r)EP - = Ezh (& ) 

II(E&H)ill2 = 1, (i = 1,. . . , n) 

where 

(4.15) 

nz n 

(4.16) 

(4.17) 

(4.18) 

and E is an nxn real valued matrix, whose rows are obtained from a solution of the following 

eigenvalue problem 

det(C - pD) = 0, (4.19) 

HH where D is an nxn matrix whose elements are < di , q5j > and C is an nxn matrix whose elements 

are < AH# - rf(43),+$ >. Since D is invertible there are n eigenvalues where complex ones 

appear in pairs of conjugates. A complex eigenvalue p correspond to a pair (Xj, pj) for some j 5 1 

while a real one corresponds to a real Xj for a j > 1. The real and the imaginary part of a complex 

eigenvector will define two rows in E that correspond to the corresponding pair (A,,pj), while a 

real valued eigenvector will define one row in E in a place corresponding to that of the eigenvalue 

associated with it. 

Having solved the coarse grid problem approximately the following interpolation is to be used: 

(4.20) 
11 

c 

4.2 Multigrid Algorithm 

Consider a sequence of grids nk(k 5 m) with mesh sizes hk satisfying 2hk+l = hk . Suppose on 

each grid operators Ak are given in such a way that Ak, (k < m) is an approximation to Ak+'. 

Assume also that interpolation operators from coarse to fine grids, and a restriction operator 

I:-', from fine to coarse grids, are given. The problem is to solve the following equation given on 

the finest level 

(A" - ")+" - = 0 , IlI;(&")ill2 = 1 (i = 1,. . .,fa). 

On levels (k < m) the equations to be solved are of the form 

1k (Ak - - = Erk -m (+k+l), - )ill2 = 1 (i = 1,. . . , fa). 

where 

(4.21) 

(4.22) 

zE(g") = 0 (4.24) 

where gk+l is any fine grid n-vector function. 

Note that the finest grid problem, Le., for k=m, can be written as the general k-grid one with 

the choice E=I, the nxn identity matrix. 

Given an approximate solution {(+", yY""T, A, p, E} to the problem (4.22) the multigrid cycle 

for improving it is denoted by 

and is defined recursively as follows: 

(4.25) 

if k=l then solve (4.22) by enough local relaxations together with global steps 

else 
12 

0 Perform vl local relaxation sweeps on (4.22), starting with (+k, $k, A, p, E) and resulting 

in a new approximation ($k, $k, A, p, E) 

-k-l- k 

0 Starting with +k-l = Ii-l$k,+k-l = Ik + make 7 successive cycles of the type 

(C)k-',$k-l,A,p,E) c MG(k - l,+k-l,+k-l,A,p,E) 

0 Interpolate correction according to the formula $k = $k + - Ij-lp) . 

0 If k = m make the change 6" e E6", 

0 Perform v2 relaxation sweeps on (4.22) starting with ($', $, A, p, E) and yielding 

then set E = I. 

(cjk, $", p, E), the final result of (4.25). 

FMG Algorithm 

In order to obtain full efficiency, the first approximation on a given level is obtained from a 

solution of the same problem on the next coarser level, which itself has been calculated in a similar 

way. The resulting algorithm is called (FMG) and is described next. Let n%-, be an interpolation 

operator (usually of higher order than I,"-,). Given the problem (4.21) , the N-FMG solution of 

that problem is: 

N-FMG Algorithm 

Set rn = 1. 

Obtain initial approximation as in section (3.1). 

Calculate (+', +', A, p) the solution of (4.21) for m=l by several relaxations. 

for m = 2,. . . , M do: 

0 Calculate +m n~-,+~-l, 4" + nz-l+m-l 

0 Perform the cycle (+m, +", A, p, E) + MG(m, 4", $m, A, p, E) N times. 

5 Numerical Examples 

Numerical experiments were carried out with the following operator 
13 

in (O,l)x(O,l), with a0 x bo < 0. 

The eigenvalues of this operator are 

x f 1p 

where X is an eigenvalue of A and p = dm. 

Uniform discretization was used, where A was discretized using the standard 5-point Laplacian. 

was the 9-point full 

weighting operator. The following parameters were used in the multigrid cycling: u1 = 2,uz = 

1,7 = 1. Initial approximation on the coarsest level was obtained as in section (3.1). Results for 

the 2-FMG algorithm are given in Tables 1-4. The L2 norm of the residuals as well as the values 

of Ah and ph are given at the end of each cycle. Although a 2-FMG algorithm was used the results 

clearly show that 1-FMG algorithm solves the problem to the level of discretization errors. Observe 

that the imaginary part p of the eigenvalue is exact on coarse grids since the truncation errors in 

approximating the off diagonal elements of A are zero, and only those elements contribute to the 

non symmetry of the problem. Also note the insensitivity of the results to the closeness of the 

eigenvalue to the imaginary axis. 

was bi-cubic interpolation, I:" was bi-linear interpolation and 

References 

[l] Brandt A. (1984): Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics. 

Monograph, 191 pages. Available as GMD studien 85 from GMD-AIW, Postfach 1240, D-5205, 

St. Augustin 1, W. Germany. 

[2] Brandt A., McCormick S. and Ruge J. (1983): Multigrid algorithms for differential eigenproblems. SIAM J. Sci. Stat. Components. 4, pp. 244-260. 

[3] Brandt A. and Ta'asan S. (1985): Multigrid method for nearly singular and slightly indefinite 

problems. in Multigrid Methods I1 ( W. Hackbusch and U. Trottenberg, Eds, ) Lecture Notes 

in Mathematics 1228, Springer-Verlag. pp. 99-121. 

[4] Drazin P. G. and Reid W. H. (1981): The Theory of Hydrodynamic Stability, Cambridge University Press, New York. 
14 

2 

I [5] Wilkinson J. H. (1965): The Algebraic Eigenvalue Problem, Clarendon, Oxford University Press. 

[6] Stanly H. E.(1971): Introduction to Phase- Transition and Critical Phenomena, Oxford University Press, New York. 

, [7] Goldhirsch I., Orszag S. A., Maulik B. K. (1987): An Efficient Method for Computing Leading 

Eigenvalues and Eigenvectors of Large asymmetric Matrices, J. of Sci. Components., vol. 2, 

NO. 1, pp. 99-58 

1 .279e+00 -19.4809 1.00000 

2 .150*01 -19.4863 1.00000 

I 

3 

4 

1 .298e-01 -19.6737 1.00000 

2 .151*02 -19.6757 1.00000 

1 .312*02 -19.7228 1.00000 

2 .203*03 -19.7233 1.00000 

2 

3 

4 

Table 1: 00 = 1. bo = -1. c2 = 0. 

1 .279e+00 -19.4809 3.16227 

2 .150*01 -19.4863 3.16227 

1 .298*01 -19.6737 3.16227 

2 .151*02 -19.6757 3.16227 

1 .312*02 -19.7228 3.16227 

2 .203+03 -19.7233 3.16227 

L I 

Table 2: a0 = 10. bo = -1. c2 = 0. 
15 

level cycle ))residucrls))2 Ah 

2 1 .279e+00 0.51912 

2 .15Oe-O1 0.51370 

3 1 .298e-O1 0.32634 

2 .151e-O2 0.32421 

Ph 

3.16227 

3.16227 

3.16227 

3.16227 

4 I 1 I .312e-02 I 0.27714-1 3.16227 

2 .203e-03 0.27766 3.16227 

Table 3: a0 = 10. bo = -1. c2 = 20. 

3 

4 

1 

1 .298e-01 0.32634 1.OoooO 

2 .151e-O2 0.32421 1.OOOOO 

1 .312e-O2 0.27714 1.00000 

2 .203e-03 0.27766 1.00000 

Table 4: 00 = 1. bo = -1. c2 = 20. 
Report Documentation Page 

1. Report No. 

NASA CR-181711 

ICASE Report No. 88-50 

2. Government Accession No. 3. Recipient’s Catalog No. 

I 

7. Author(s) I 8. Performing Organization Report No. 

4. Title and Subtitle 

MULTIGRID METHOD FOR STABILITY PROBLEMS 

Shlomo Ta’asan 

5. Report Date 

September 1988 

6. Performing Organization Code 

I 88-50 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 

Mail Stop 132C, NASA Langley Research Center 

Hampton, VA 23665-5225 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, VA 23665-5225 

and Engineering 

11. Contract or Grant No. 

NAS1-18107, NASl-18605 

13. Type of Report and Period Covered 

Contractor Report 

14. Sponsoring 4gency Code 

Final Report 

_- 

17 Key Words (Sugyested by Authorls)) 

n,*: rig, z;tattlitv, eigenvalues 

18 Distribution Statement 

64 - Numerical Analysis 

70 - Physics (General) 

Unclassified - unlimited 

I I I 

NASA-Langley, 1988 NASA FORM 1626 OCT 86 

19 Security Classif (of this report) 

Unclassified 

20 Security Classif (of this page) 21 No of pages 22 Price 

Unclassified 17 A0 2