8:40-9:40

Element-Free Galerkin (EFG) Methods
Ted Belytschko, Department of Mechanical Engineering, Northwestern
University, Evanston, IL
tedbelytschko@nwu.edu

In many problems of computational mechanics, such as the simulation of
failure by crack growth or of manufacturing processes such as casting,
one is forced to deal with moving discontinuities. Conventional
computational methods such as finite element, finite volume or finite
difference methods are not well suited to these problems since the
most viable strategy for dealing with moving discontinuites - to
remesh in each step of the evolution so that the mesh lines remain
coincident with the discontinuities - introduces many difficulties.
The objective of meshless methods is to eliminate at least part of
this structure by constructing the approximation entirely in terms of
nodes.

The element-free Galerkin (EFG) method is a meshless method for
solving partial differential equations which uses only a set of nodal
points and a CAD-like description of the body to formulate the
discrete model. The EFG method employs moving least-square (MLS)
approximants to construct an approximate solution for the PDE. The
main three components are: a compactly supported weight function
associated with each node, a basis (usually consisting of
polynomials), and a set of coefficients that depend on position. One
attractive property of MLS interpolants is that their smoothness can
be prescribed by the smoothness of the weights. Some new variants of
the element-free galerkin method which are based on kernel
approximations with corrections on the functions or their derivatives
are described.  Then examples are given of solutions of crack growth
in two and three dimensions.

 back to top

---

---

Contibuted Papers Ia (9:45-10:45)

A New Formulation of the Modified Nodal Integral Method
Bradley Wescott and Rizwan-uddin, Department of Nuclear Engineering,
University of Illinois at Urbana-Champaign
wescott@students.uiuc.edu

Nodal integral methods are coarse-mesh methods that rely on transverse integration and
analytical solutions of the averaged one-dimensional equations to yield solutions in
less CPU time than conventional schemes.  An efficient formulation of the recently proposed
modified nodal integral method has been developed to further reduce computation time when solving
nonlinear partial differential equations with a nonlinear convection term such as Burgers'
equation and the Navier-Stokes equation.  In this formulation, by adding and subtracting a
linearized convection term, in which the node-averaged velocity at the previous time step
multiplies the spatial derivative, the node-interior approximate analytical solution is
developed in terms of this previous time step node-averaged velocity.  This leads to a set of
discrete equations with coefficients which only need to be evaluated once each time step for each
node, resulting in a significant reduction in computing time when compared with the original
MNIM formulation.  A numerical scheme using these node-averaged velocities at the previous time
step--to be referred to as M2NIM--has been developed for the two-dimensional, ime-dependent
Burgers' equation.  The method is shown to be second order and to posses inherent upwinding.
When compared with MNIM, numerical results show a significant reduction in the computation time
without sacrificing accuracy.  In addition to the standard implicit formulation, an explicit M2NIM
based on the evaluation of the space-averaged u and v nodal velocities at the next time step
(u_xy and v_xy) using the previous time step solution has also been developed.  For the next
time step, given u_xy and v_xy, all other discrete variables can be evaluated.  Results for
the implicit and explicit M2NIM will also becompared.

---

A Numerical Scheme for the Convection-Diffusion Equation in Arbitrary Geometries
Allen Toreja and Rizwan-uddin, Department of Nuclear Engineering,
University of Illinois at Urbana-Champaign
atoreja@uiuc.edu

The convection-diffusion equation (CDE) can be used to determine contaminant transport in
groundwater flow.  A previous application of a hybrid nodal-integral/finite-element method to
the heat conduction problem has demonstratede the feasibility of combining the nodal integral
method (NIM) with a finite element type method (FEM) [1].  In that approach, the computational
domain was divided into rectangular and triangular nodes.  The rectangular nodes were treated using
the standard NIM [2], and a finite element type approach was applied to the triangular nodes.  In
this work, we have developed a new hybrid scheme for the CDE using a set of trial functions with
exponential terms for the triangular nodes.  The exponential terms of these trial functions are
necessary to adequately treat the convective terms of the CDE.  In the limit of pure diffusion,
these trial functions recover the linear trial functions used in the previous application.
Numerical results for benchmark problems and comparisons will be presented.

1.  Rizwan-uddin.  Hybrid Nodal-Integral/Finite-Element Method for Arbitrary Geometries.
    Trans. Am. Nucl. Soc., 1997, 76, 170-172

2.  Michael, E.P.E., Dorning, J.J., and  Rizwan-uddin.  Nodal Integral Method for the Convection-Diffusion Heat Equation.Trans. Am. Nucl. Soc., 1993, 69 239-241

---

Modified Nodal Method for Navier-Stokes and Energy Equations
Fei Wang and Rizwan-uddin, Department of Nuclear Engineering,
University of Illinois at Urbana-Champaign
feiwang1@uiuc.edu

Two-dimensional, time-dependent Navier-Stokes and energy equations
coupled via the Boussinesq approximation are solved using the newly
developed advanced coarse-mesh nodal method. The transverse-integrated
energy equation was first solved for temperature, and solution for
temperature was plugged into the continuity and momentum equations to
solve for the velocities and pressure. After eliminating the "pseudo"
source terms using the constraint equations, the problem was closed
with 11 equations for 11 unknowns. Instead of moving all the
non-linear terms to the right-hand side of the nodal integral
equations as pseudo source terms, which is done in the conventional
nodal method, the method features a linearized convection term kept on
the left-hand side. Fortran code was developed using the developed
scheme and successfully applied to lid-driven problem to produce the
exact linear distribution velocity profile, and Hagen-Poiseuille flow
to produce the exact parabolic velocity profile. Compared with the
conventional nodal integral method, this modification is expected to
greatly improve the accuracy and the convergence speed of the
scheme. We are currently applying the method to two-dimensional,
time-dependent problems.

back to top

---

---

Contributed Papers Ib (9:45-10:45)

Inexact approximate simplified Newton iterations for implicit Runge-Kutta methods
Laurent Jay, Department of Mathematics, University of Iowa, Iowa City, IA
ljay@math.uiowa.edu

The major difficulty in the implementation of fully implicit Runge-Kutta (IRK) methods is to
solve the arising systems of nonlinear equations. In this talk we show that the use of inexact
approximate Newton methods can be very efficient for this purpose.

---

Design and Implementation of DIRK Integrators for Stiff Systems
Roger Alexander, Department of Mathematics, Iowa State University, Ames, IA
alex@iastate.edu

In this talk I will describe the development of some new DIRK formula
pairs for the integration of stiff systems.

The basic integrator in each pair has four stages, with the first
(formally) explicit, and is stiffly accurate and L-stable.

It turns out that an A-stable integrator of order four of this
form must evaluate derivatives at points outside the set. So we take
the basic integrator to have order three and consider alternatives
allowed by the remaining free parameters:

1) Minimize the truncation error of the formula, and pair it with
   a second-order error estimator ("DIRK(3,2)").

2) Pair the formula with a fourth-order error estimator
   ("DIRK(3,4)").

I will report experiments with these formulas comparing different
predictors for the stages, and comparing their performance with the
well-known Hairer-Wanner code SDIRK4, a DIRK(4,3) formula pair.

---

A New Approach for solving Thomas-Fermi Equation
A.M.Wazwaz, St. Xavier University, Chicago, IL
WAZWAZ@vaxd.sxu.edu

The Thomas-Fermi equation which plays an important role in mathematical
physics will be investigated.  We propose a new approach to develop
a nonperturbative approximate solution for this equation.  This approach
rests on the modified decomposition method combined with the diagonal
Pade' approximants.  The approximation obtained for the initial slope
y'(0) shows significant improvements over existing approximations.

back to top

---

---

Contributed Papers IIa (11:00-12:00)

Computational Methods for Kadomtsev-Petviashvili Equations
E.H. TWIZELL, Department of Mathematics and Statistics, Brunel
University Uxbridge, Middlesex, England
E.H.Twizell@brunel.ac.uk

Computational techniques based on a linearized, implicit,
finite-difference scheme and on a predictor-corrector approach
are presented for the solution of the Kadomtsev-Petviashvili (KP)
equation and its generalized form (GKP).

An important advantage to be gained from the use of the former
method over the latter, which is conditionally stable, is the
ability to vary the mesh length, thereby reducing computational
costs.

Numerical results portraying a single line-soliton solution, the
interaction of two line-solitons, a lump-like soliton, and the
interaction of two lump solitons will be reported for the KP
equation.

---

Solving Spherical PDEs using Radial Basis Functions
Tanya Morton, Mike Neamtu, and Larry L. Schumaker, Department of
Mathematics, Vanderbilt University, Nashville, TN
morton@math.vanderbilt.edu

Consider the three topics:  PDEs, radial basis functions, and
spheres. Each combination of two of these three topics gives
a research field of considerable interest; solving PDEs using
radial basis functions, solving PDEs defined on spheres, and
radial basis functions on spheres. Our aim is to consider
the intersection of all three topics, that is, the use of
radial basis functions to solve spherical PDEs. We will present
convergence order results for this problem, and illustrate
these results numerically.

---

Wavelet-Galerkin Method for the Vibration of Cable Carrying a Moving Mass
M. AL-Qassab and S. Nair, Department of Mechanical, Materials and
Aerospace Engineering, Illinois Institute of Technology, Chicago, IL

The dynamic behavior of an elastic cable supported between two points
on the same level and vibrating in a plane due to a moving mass along
its length is investigated here using multilevel representational
wavelets. The equations of motion are derived using the Hamilton's
principle which resulted in nonlinear partial differential
equations. The equations of motion are linearized and the solution is
obtained by discretization following the procedure of Galerkin method
using anti-derivative Wavelets to represent the spatial variable
function. The anti-derivative wavelets guarantee satisfaction of the
boundary conditions at the supports. After discretization, equations
of motion are reduced to a system of ordinary differential equations
which are integrated over time using implicit method. examples are
given to examine the effect of different wavelets on the
representation of the dynamic shape of the cable.

back to top

---

---

Contributed Papers IIb (11:00-12:00)

Generating Orthogonal Polynomials for Exponential Weights on a Finite Interval
Raymond C. Y. Chin, Department of Computer and Information Science,
Indiana University Purdue University Indianapolis, Indianapolis, IN
rchin@cs.iupui.edu

Generating orthogonal polynomials given a positive weight is simple in
theory but difficult in practice.  The actual construction of the
polynomials is fraught with numerical instabilities.

In theory, the process of constructing orthogonal polynomials is the
Gram-Schmidt orthogonalization process applied to a subspace of
polynomials, called the Stieltjes procedure.  The sequence of
orthogonal polynomials satisfies a three term recurrence whose
coefficients are the weighted centroids and the ratios of the weighted
norms of the successive polynomials. From another perspective, it is
also a multiple scales problem.  The scales are the inherent scale of
the weight and the polynomials and the interval length.

In practice, the calculations involve integrating a $2n$th degree
polynomial for each $n$th-order coefficient and marching forward via
the three term recurrence relation.  This process is extremely
delicate and is susceptible to numerical perturbations. Numerical
stability is a major concern.

In this presentation, we will discuss an alternative formulation that
focuses on a set of nonlinear recurrence equations governing the
coefficients of the three-term recurrence and the values of the
orthogonal polynomials at the end points. This set of algebraic
equations is redundent with a core of linearly independent
equations. We will subject these equations to an asymptotic analysis
to extract their local behaviors and to select the final set of
equations.  We will discuss resulting asymptotics induced domain
decomposition method to the orthogonal polynomial generation problem.

---

An Efficient Block-Matrix Algorithm for Evaluating Orthogonal
Approximations to EXP(K) Arising in Reduced Hartree-Fock (RHF)
Calculations
Mike Minkoff, MCS, Argonne National Laboratory, Argonne, IL, and Paul
Sutton, North Central College, Naperville, IL
minkoff@mcs.anl.gov

The orthogonal transformation EXP(K) incorporates non-redundant SCF
coefficients in a skew-symmetric matrix K.  In the RHF case at least
half of the elements are members of two diagonal null blocks.  By
using a block-matrix approach to evaluating K and its powers, an
orthogonal approximation is obtained without using these null blocks.
We compare this approach with the alternative of approxming
EXP(K) by applying Taylor's series expansions with recursive
doubling and reorthogonalization procedures. This latter approach is
used in production large-scale parallel RHF software.  Issues of
relative size and matrix blocks in a distributed software
implementation are examined.

---

Estimating the Hausdorff Dimension
David Stewart, Department of Mathematics, University of Iowa, Iowa City, IA
dstewart@math.uiowa.edu

The Hausdorff dimension of a set is often confused
with the fractal (or capacity, or box-counting)
dimension, but the two can be different even
for quite simple sets.  This difference will
be explained, along with a numerical algorithm
which can distinguish between them, and some
numerical results.

back to top

---

---

1:30-2:30

Boundary Value Problems for ODEs
Lawrence F. Shampine, Department of Mathematics, Southern Methodist
University, Dallas, TX
lshampin@mail.smu.edu

Solving boundary value problems (BVPs) for ODEs is not the routine
matter that solving initial value problems has become.  We'll look
at some simple examples that show a BVP might have a finite number
of solutions, an infinite number, or none at all.  Parameters that
are to be determined as part of the problem and singularities are
not unusual.  bvp4c is a collocation code for BVPs that controls
the size of residuals.  We'll talk about what collocation is and
about some new theoretical results that make the unusual control
of residuals very attractive. In the course of solving a couple of
numerical examples, we'll talk about how facilities of the MATLAB
problem solving environment can be exploited to obtain a strikingly
simple user interface.  bvp4c and a tutorial with many examples are
available gratis by anonymous ftp.

back to top

---

---

Contributed Papers IIIa (2:30-3:10)

PVODE, An ODE Solver for Parallel Computers
George D. Byrne, Illinois Institute of Technology, Chicago, IL, and
Alan C. Hindmarsh, Lawrence Livermore National Laboratory, Livermore, CA
gdbyrne@mediaone.net

PVODE is a general purpose solver for stiff and non stiff ordinary
differential equations (ODEs) and is based on CVODE, which is written in
ANSI standard C.  PVODE uses MPI (message passing interface) and a
variation of the vector module in CVODE to achieve parallelism.  The
idea is to solve large-scale ODEs in the SPMD  (Single Program Multiple
Data) environment with distributed memory.  Often such large scale ODEs
arise in the solution of partial differential equations.  To date, PVODE
has been run on an IBM SP-2, a Cray T3D, a Cray T3E, and a cluster of
Sun workstations.  It is currently being used to simulate tokamak edge
plasmas at Lawrence Livermore National Laboratory.  Some numerical
results will be given.

This work was done in collaboration with Alan C. Hindmarsh, Lawrence
Livermore National Laboratory.

---

On The Impicit Integration Of Differential-Algebraic Equations of Multibody Dynamics
Dan Negrut, Mechanical Dynamics, Inc., Ann Arbor, MI, and Edward J. Haug, Department of Mechanical Engineering, University of Iowa, Iowa City, IA
dnegr@adams.com, haug@nads-sc.uiowa.edu

Three methods for the state-space based implicit integration of
differential-algebraic equations of multibody dynamics are summarized and
numerically compared.  In the state-space approach, the time evolution of
a mechanical system is characterized using a number of generalized
coordinates equal with the number of degrees of freedom of the system.
In this paper these independent generalized coordinates are a subset of
the Cartesian position coordinates and orientation Euler parameters of body
centroidal reference frames.  Depending on the method, the independent
generalized coordinates are implicitly integrated and dependent quantities
(including Lagrange multipliers) are determined to satisfy constraint
equations at position, velocity, and acceleration levels.  Five computational
algorithms based on the proposed methods are used to simulate the motion of
a stiff 14-body vehicle model.  Results show that the proposed methods deal
effectively with challenges posed by stiff mechanical system simulation.
A comparison with a state-space based explicit algorithm for the simulation
of the same model indicates a speed-up of approximately two orders of
magnitude.

back to top

---

---

Contributed Papers IIIb (2:30-3:10)

The QR Algorithm for Orthogonal Hessenberg Matrices
Greg Ammar, Department of Mathematical Sciences, Northern Illinois
University, De Kalb, IL
ammar@math.niu.edu

Unitary Hessenberg matrices are closely related to Szego polynomials,
and have a structure that is amenable to exploitation in eigenvalue
computations.
We will give an overview of these points, and pay particular attention to
Gragg's unitary Hessenberg QR algorithm.  We'll outline a derivation of
the UHQR algorithm using a device for describing the efficient
implementation of single-bulge chasing procedures on unitary Hessenberg
matrices.

We will then turn our attention to the efficient implementation of the QR
algorithm on (real) orthogonal Hessenberg matrices.  Using our device
for deriving the UHQR algorithm, we will show that the double-bulge
chasing sweep that arises from the Francis shift strategy can be
implemented by interleaving three single-bulge chasing sweeps.
The resulting OHQR algorithm avoids the additional storage and
computation associated with the complex arithmetic that is required
when the single-shift UHQR algorithm is applied to a real matrix.

---

Efficient Evaluation of Scattered Data Interpolants Obtained Using a Thinning Algorithm
Greg Fasshauer, Department of Applied Mathematics, and Jesus Miranda, Department of Computer Science, Illinois Institute ofTechnology, Chicago, IL
mirajes@charlie.iit.edu

In [1] the following algorithm for scattered data interpolation was proposed:
Given a set of scattered centers, ${\cal X} \subset \RR^d$, and a function
$f : \RR^d \to \RR $ to be matched on $\cX$, a nested sequence of subsets of $\cX$,
$\cX_1 \subset \cX_2 \subset \cdots \subset \cX_M = \cX$, is built by using Delaunay triangulations.
The subsets are chosen according to the density of the current subset. The density is computed by using Delaunay triangulations.  At the $k$-th step one matches the error function
$$
 f - \sum_{i=1}^{k-1} s_i, \qquad k = 1, \ldots, M,
$$
on $\cX_k$  by computing the coefficients of the $k$-th interpolant
$$
s_k(x) = \sum_{j=1}^{N_k} c_j^k \phi_{\alpha_k}(\|x-x_j^k\|),
$$
where $\phi_{\alpha_k}$  is a smooth compactly supported radialfunction and $\alpha_k$ is related to the density of the set $\cX_k$. The Conjugate Gradient method is used in order to obtain the $c_j^k$'s.  Our presentation focusses on the fast evaluation of the error function. We describe various algorithms and give numerical comparisons.

1. Floater, M.~S., and A.~Iske. Multistep scattered data interpolation using compactly supported
radial basis functions. J.\ Comput.\ Applied Math.\ 73 (1996), 65--78.

back to top

---

---

Contributed Papers IVa (3:30-4:30)

Recent Developments in Post-Processing Solutions to Partial Differential Equations
Mike Gosz, Department of Mechanical, Materials, and Aerospace
Engineering, Illinois Institute of Technology, Chicago, IL
gosz@mmae.iit.edu

Because of the tremendous advances in computer speed and memory
capacity thathave occurred over the past ten years, we now have the
capability to solve large-scale,three-dimensional partial differential
equations efficiently and at low cost.  Because three-dimensional
solutions to partial differential equations are inherently difficult
to visualize, the need for post-processing and visualization tools is
readily apparent.  In the presentation, a shareware program called the
Persistence of Vision Ray-Tracer (POV-RAY) is described for
visualizing solutions of partial differential equations.  The program
creates three-dimensional, photo-realistic images using a rendering
technique called ray-tracing. Ray-tracing is a rendering technique
that calculates an image of a scene by shooting rays into the
scene. The scene is built from shapes, light sources, a camera,
textures, and other special features.  It produces high quality images
with realistic reflections, shading, perspective and other effects.
To demonstrate the usefulness of the method, rendered animations of
solutions to problems solved by the finite element method are
presented.  The problems range from flexible structures subjected to
ocean waves, to crack propagation in three-dimensional solids.  Ideas
for utilizing the program for post-processing solutions obtained by
other numerical techniques, such as meshless methods, are also
discussed.

---

Finite Elements for Elliptic Interface Problems with Strongly Discontinuous Coefficients and Solutions
Daoqi Yang, Department of Mathematics, Wayne State University, Detroit, MI,
yang@math.wayne.edu

An iterative finite element algorithm is proposed for numerically
solving two-phase steady-state generalized Stefan interface problems
with discontinuous solutions, conormal derivatives, and coefficients.
This algorithm employs finite element methods and
iteratively solves smaller subregion problems for each
phase with good accuracy, and exchanges information at the interface
to advance the iteration.
The finite element grids in different phases do not have to match
each other at the interface. Numerical experiments are performed to
show the accuracy and efficiency of the algorithm for capturing
discontinuities in the solutions and coefficients.
One surprising property of the algorithm is that its accuracy does not
deteriorate as the discontinuity in the coefficients becomes worse.
That is, the accuracy remains the same for continuous problems and
strongly discontinuous problems. Another surprising property is that
its conditioning becomes better as the discontinuity gets worse. In
other words, the stronger the discontinuity, the faster convergence.
Numerical examples on matching and non-matching unstructured grids
are given with coefficient discontinuity jumps in the order of
$10^3, 10^5, 10^{10}, 10^{50},$ and $10^{100}.$

---

The linear rational collocation method for PDEs
Richard Baltensperger and Jean-Paul Berrut, Universite de Fribourg,
Mathematiques, Fribourg/Perolles, Switzerland
Jean-Paul.Berrut@unifr.ch

The global pseudospectral method of lines for solving hyperbolic or parabolic differential equations for a function $u(x,t)$ of the spatial variable $x$ and the time variable $t$ can be described as follows: One chooses interpolation points (nodes) $x_i$ in the spatial interval, replaces $u$ for every $t$
by a polynomial $\widetilde u$ interpolating between the $x_i$'s, inserts $\widetilde u$ into the equation and collocates at the $x_i$'s. This leads to a system of ODEs in $t$ that is then solved by a time marching technique. In the latter, the vectors of values of $\widetilde u$ are multiplied by differentiation matrices $D$ whose entries are the derivatives of the fundamental Lagrange polynomials at the nodes.

For precise and well--conditioned interpolation, the $x_i$'s must be chosen in such a way that their preimages on the circle by the application $\arccos$ are almost equidistant (e.g., \v Ceby\v sev or Legendre points), which make them concentrate near the extremities of the interval. This results in
ill--conditioned differentiation and highly nonnormal $D$'s which for stability require smaller time
steps than difference or element methods using equidistant nodes.

For the \v Ceby\v sev case, Kosloff and Tal-Ezer in 1993 suggested to shift the points toward a more regular distribution by means of a conformal map. In view of the chain rule, the corresponding change of variable results in differentiation matrices that are sums of products of the $D$'s with diagonal matrices of values of derivatives of the map. These new matrices are closer to normality than the $D$'s and therefore increase the timestep necessary for the stability of the time marching technique.

In our talk we will show how the same can be achieved by replacing the interpolating polynomial by some linear rational interpolant, i.e., a rational interpolant whose denominator is independent of the interpolated function.

back to top

---

---

Contributed Papers IVb (3:30-4:30)

A Domain Based Multilevel Preconditioning Technique For Solving General Sparse Linear Systems
Jun Zhang, Department of Computer Science, University of Kentucky, Lexington, KY
jzhang@cs.uky.edu

We introduce a domain based multilevel block incomplete LU
preconditioner (BILUTM) for solving general sparse linear
systems. This preconditioner combines a high accuracy
incomplete LU factorization with an algebraic multilevel
recursive reduction. In the first level the coefficient matrix is
permuted into a two by two block form using (block) independent set
ordering and an ILUT factorization for the reordered matrix is performed.
The reduced system is the approximate Schur complement associated
with the partitioning and it is obtained implicitly as a byproduct
of the partial ILUT factorization with respect to the complement of
the independent set. The incomplete factorization process is
repeated with the reduced systems recursively. The last reduced system
is factored approximately using ILUT again. The successive reduced
systems are not stored. This implementation is efficient in controlling
the fill-in elements during the multilevel block ILU factorization,
especially when large size blocks are used in domain decomposition type
implementations. Numerical experiments will be presented to show the
robustness and efficiency of the proposed technique for solving some
difficult sparse linear systems.

---

Multigrid Method and High Order Compact Scheme for Solving Boundary
Layer Problems on Nonuniform Grids
Lixin Ge (speaker) and Jun Zhang, Department of Computer Science,
University of Kentucky, Lexington, KY
lixin@csr.uky.edu  jzhang@cs.uky.edu

A fourth order compact finite difference scheme (nine-point compact
difference scheme) and a multigrid cycling algorithm are employed
to solve the two dimensional convection diffusion equations with
boundary layers. The computational domain is first discretized on a
nonuniform (graded) grid to resolve the boundary layers. A grid
transformation technique is used to map the nonuniform grid to a
uniform one. The fourth order compact scheme is applied to the
transformed uniform grid. A multigrid method is then used to solve
the resulting linear system. We conduct analyses to show how the
grid stretching may affect the computed accuracy of the solutions
from the fourth order compact scheme. We also demonstrate that the
stretching may also affect the convergence rate of the multigrid
method. Numerical experiments are used to show that a graded mesh
and a grid transform are necessary to compute high accuracy solution
with the fourth order compact scheme for convection diffusion problems
with boundary layers. Accuracy comparisons between the standard central
difference scheme and the present fourth order compact scheme are given.

---

A Sparse Approximate Inverse Technique for Solving General Sparse Matrices
Jiqing Zhang (speaker) and Jun Zhang, Department of Computer Science,
University of Kentucky, Lexington, KY
zhangj@cslab.uky.edu  jzhang@cs.uky.edu

Krylov subspace methods are popular choices of iterative solution
techniques due to their ability to handle large general sparse
matrices. It is found that the performance of these Krylov subspace
methods can be significantly enhanced by coupling them with a
suitable preconditioner, and that the convergence rate of a
preconditioned iterative method is usually dictated by the quality
of the preconditioner. Therefore, designing and identifying efficient
preconditioning techniques have motivated many researches.

In this paper, we propose a forward factored sparse approximate
inverse for inverting nonsymmetric and indefinite matrices. This method
is extracted from a dense matrix inverse technique that is based on
computing two sets of biconjugate vectors. The forward sparse
approximate inverse is computed as a factored form and is used as a
preconditioner to work with some Krylov subspace methods. Sparsity
patterns of the factored inverse are exploited to reduce computational
cost.

We tested the forward factored sparse approximate inverse techniques
for solving a few linear systems with general sparse matrices arising
from realistic applications. Numerical experimental results are presented
to make a comparison between this method and a backward factored sparse
approximate inverse techniques to show the effectiveness and efficiency
of the new sparse approximate inverse preconditioner. The results show
although both methods are very much cheap in terms of the cost of
construction on sequential machines, the computation of the forward
decomposition method is relatively cheaper and runs faster than that of
the backward method.

Among several other advantages, the main advantage of the new preconditioner
is its robustness. We found that, in general, the construction cost and the
solution cost are comparable, unless the numbers of iterations are very
small, or the matrices are very large.

back to top