Yu Wang

DYNAMICS OF UNSYMMETRIC PIECEWISE-LINEAR/NON-LINEAR SYSTEMS USING FINITE ELEMENTS IN TIME

The dynamic response and stability of a single-degree-of-freedom system with unsymmetric piecewise-linear/non-linear stiffness are analyzed using the finite element method in the time domain. Based on a Hamilton's weak principle, this method provides a simple and efficient approach for predicting all possible fundamental and sub-periodic responses. The stability of the steady state response is determined by using Floquet's theory without any special effort for calculating transition matrices. This method is applied to a number of examples, demonstrating its effectiveness even for a strongly non-linear problem involving both clearance and continuous stiffness non-linearities. Close agreement is found between available published findings and the predictions of the finite element in time approach, which appears to be an efficient and reliable alternative technique for non-linear dynamic response and stability analysis of periodic systems.

INTRODUCTION

In a large class of vibrating systems there exist piecewise-linear or piecewise-non-linear stiffnesses, due to the presence of gaps, clearances, backlash and impacting components. Bearings, gears, impact print hammers and offshore structures are a few examples that fall into this category. In the design of these systems, it is essential to obtain a comprehensive understanding of their dynamic behavior in order to prevent excessive levels of noise and vibration. For a simple piecewise-linear system, its dynamic response may be solved exactly because the exact solution form is known in the intervals which the system possesses constant stiffness characteristics [START_REF] Shaw | The dynamics of a harmonically excited system having rigid amplitude constraints[END_REF]2]. However, the characteristics of a general system are inherently non-linear and the prediction of its response is not an easy task, particularly if the stiffness is piecewise-non-linear.

In this paper a numerical study is presented of unsymmetric piecewise-linear/non-linear systems subject to general periodic excitations. Such systems are often used as models of compliant offshore structures and articulated mooring towers [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF] and their dynamics have been previously investigated with well-illustrated complex characteristics in their response, including chaotic vibrations. A common solution procedure for the strongly non-linear system is the numerical time integration, which can give both transient and steady state solutions. However, there are some difficulties in this approach. For example, it is usually very time-consuming to obtain the system response with respect to a system parameter or parameters, especially in parameter regions where multi-valued responses exist. Furthermore, numerical integration with digital computers may encounter numerical difficulties for some types of integration algorithms, particularly when the non-linearity becomes very strong [START_REF] Singh | 22nd Biennial Mechanisms Conference[END_REF]5].

An alternative approach is the harmonic balance method. With assumed waveforms of vibration, approximate solutions for a single degree-of-freedom oscillator are determined in reference [6]. A more general harmonic balance procedure is developed in references [7][8][9] and has been applied to the two different unsymmetric piecewise-linear systems given in references [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF] and [6]. The approach is known to be limited by its ability of solving the resulting non-linear algebraic equations. For a number of different implementations of the procedure [START_REF] Ling | [END_REF][START_REF] Cameron | [END_REF][12][13], it is found that the strength of non-linearity, the number of harmonics sought in the solution and the relative magnitudes of higher harmonics with respect to the fundamental are, among others, the major factors that are instrumental for the success of the method [START_REF] Cameron | [END_REF]13]. In particular, if an impulse-like force, for example, generated by a rigid stop, is applied to the system, a large number of higher harmonics are inevitably required, giving rise to a large number of system equations. In this case, it is usually difficult to achieve convergence in the numerical solution.

Another possible approach is the use of the finite element in time method (FET) based on Hamilton's principle [14][15][16][17], in which the solution for all the spatial degrees of freedom at all time steps within a given time interval of interest is sought via a set of algebraic equations. Some recent development has shown that this approach offers several potential advantages over numerical solution of ordinary differential equations in the time domain, including the flexibility in formulating the problem directly from system energy expressions, the greater accuracy at specific time points of interest, and the use of adaptive finite elements to improve computational efficiency and accuracy [14,[17][18][19].

For obtaining the periodic response of the piecewise-linear/non-linear system, moreover, two specific features of the finite elements in time make the approach an ideal tool. First, the periodic solution can be made readily available by assembling a number of time elements and imposing the appropriate periodic boundary constraint relations, in contrast with time-stepping methods, which require initial conditions to start the integration procedure but cannot impose periodic boundary conditions. For a harmonic balance procedure, the periodicity condition of a solution is attained by the periodic functions in Fourier series expansions. For the time finite element method, there is no such restriction imposed on the interpolation functions and the conventional shape functions can be used. In fact, it is easy to show that the Fourier series approach is a special case of finite elements in time in which harmonics are used as trial and shape functions for a single time element [19].

The second advantageous feature is the straightforward determination of stability of the periodic solution. Since the corresponding transition matrix for the analysis of stability of small perturbations about the periodic solution is a by-product of the FET procedure, Floquet's stability theory can be readily applied without any special effort.

These unique characteristics of the FET method have created a wide range of interest in the field of dynamics and control, where it has been successfully used for analysis and design of analytic-non-linear rotorcraft systems [14][15][16]19], for dynamics of multi-body systems [17,[START_REF] Mello | Time finite element methods for large rotational dynamics of multibody systems[END_REF] and for optimal control [START_REF] Hodges | [END_REF]. It appears that the method is well-suited as an efficient and consistent approach to predicting the dynamic response of piecewise-linear/non-linear systems. The focus of this paper is to apply the FET method to the unsymmetric piecewise-linear (and -non-linear) systems in order to demonstrate its effectiveness. The first objective of this investigation is to develop the solution of the equation of motion for piecewise-linear/non-linear systems with finite elements in time. Next, a number of techniques for efficient numerical calculations and parametric studies are presented. Finally, the predictions of this method, when applied to several piecewise linear and non-linear problems, are discussed. The stability of the periodic solution is analyzed using Floquet's theory.

EQUATION OF MOTION

The response of the system analyzed in the present investigation is described by a non-linear equation of motion, expressed in the form

mz¨+ cz˙+ k 1 z + k 2 g(z) = F(t), (1)
where m, c and k 1 denote the mass, viscous damping coefficient and linear spring constant, respectively. The second spring gives rise to an unsymmetric restoring force k 2 g(z) for the system. If the second spring is linear with a spring constant of k 2 , then the overall system stiffness is defined to be piecewise-linear, as shown in Figure 1, where the clearance value is denoted to be d and is assumed to be equal to or greater than zero. For a more general case, the second spring can have a non-constant stiffness, for example, a Hertz contact spring, representing an unsymmetric piecewise-non-linear system. The system is assumed to be subject to a general periodic excitation F(t) with a circular frequency of v.

The equation of motion can be written into a non-dimensionalized form by introducing the following non-dimensional parameter:

u = vt, v 2 n = k 1 /m, z = c/(2mv n ), b = k 2 /k 1 , P(u) = F(u)/k 1 , r = v/v n . The equation of motion becomes r 2 z0 + 2zrz' + z + bg(z) = P(u), (2)
where ( • )' denotes 1( • )/1u, and r and z are the frequency ratio and damping ratio, respectively. Without loss of generality, the amplitude of the excitation P(u) is assumed to be unity in what follows.

FINITE ELEMENT METHOD IN THE TIME DOMAIN

The finite element method in the time domain for obtaining solutions to dynamics problems is based on a weak form of Hamilton's law of varying action [14,18]. The commonly known as Hamilton's weak principle can be expressed in its general displacement form as

g t F t I (dL + dW) dt = dz • p= t F t I , (3)
which describes the real motion of the system between the two known times t I and t F . Here, L and W denote the Lagrangian and the work done by the generalized forces of the system, respectively, and p = 1L/1z˙is the set of the generalized momenta. For the piecewise-linear/non-linear system described by equation ( 2), the variational equation becomes

g u F u I {dz' • (r 2 z') + dz • [-2zrz' -z -bg(z)] + dz • P(u)} du = dz • p= u F u I . (4)
Here, p = r 2 z' and the time interval of interest, T = u F -u I , equals the period of steady state response of the system. For example, T = 2p for the fundamental periodic response of a harmonic excitation. In this variational equation, z and dz are simply required to be continuous, and p appears only as a discrete quantity at the ends of the time interval. These boundary terms are strictly enforced as natural boundary conditions, allowing p to be determined without differentiation of z [14,17,18]. Because of these special features, the use of the finite element method is the most powerful tool for an easy approximation of equation ( 4).

In a finite element approximation to the dynamic equation ( 4), the time interval T is divided into a finite number (N e ) of time elements, usually equally spaced for convenience. Similar to the standard finite element technique in spatial discretization, a set of nodal displacements Z i for each time element i are defined, and the displacement z and velocity z' within the element are interpolated among the nodal displacements Z i by using standard shape functions, for example, the Lagrange's polynomials, as

z(u) = N i (j)Z i and z'(u) = N' i (j)Z i , j = (u -u i )/(u i + 1 -u i ), i = 1, 2, . . . , N e , (5, 6)
where the vector N for the elemental shape function has (k + 1) components of polynomial function if the order of interpolation is chosen to be k. Correspondingly, the nodal displacement vector Z i has (k + 1) temporal components for the time element. The Lagrange family of interpolation functions is used in this study. Thus the equation of motion ( 4) can be put in the following variational form for time element i:

g u i + 1 u i {dZ T i [N' T i r 2 N' i Z i -2zN T i rN' i N i Z i -bN T i g(z)] + dZ T i N T i P(u)} du = dZ T i B i , (7)
where B i is a vector of (k + 1) elements containing linear momenta of the dynamic system at the two ends of the time element, B i = {-p(u i ), 0, . . . , 0, p(u i + 1 )} T . Since dZ i is completely arbitrary, the finite element approximation at the element level is given by

A i Z i + G i (Z i ) + P i = B i , (8)
where , 11) The variational equation ( 4) can now be represented by finite element approximation equations that correspond to the nodal point displacements of the assemblage of finite elements, as a sum of integrations over the time intervals of all time elements. Furthermore, for periodic steady state response, the momentum vector vanishes since the periodic boundary condition requires both displacements and velocities to be identical at u 1 and u N e + 1 , or z(u 1 ) = z(u N e + 1 ) and p(u 1 ) = p(u N e + 1 ). Therefore, the dynamic response is described at the system level by

A i = g u i + 1 u i N' T i r 2 N' i du -2z g u i + 1 u i N T i rN' i du -g u i + 1 u i N T i N i du, (9)
G i = -b g u i + 1 u i N T i g(z) du, P i = g u i + 1 u i N T i P(u) du. ( 10
AZ + G(Z) + P = 0, (12)
where, with a summation notation representing the element assemblage operation,

Z = s N e i = 1 Z i , (13)
A = s N e i = 1 A i = s N e i = 1 6g u i + 1 u i N' T i r 2 N' i du -2z g u i + 1 u i N T i rN' i du -g u i + 1 u i N T i N i du 7 , (14)
G = s N e i = 1 G i = -b s N e i = 1 g u i + 1 u i N T i g(z) du, (15)
O = s N e i = 1 P i = s N e i = 1 g u i + 1 u i N T i P(u) du. (16)
Thus, the finite element approximation results in a set of non-linear algebraic equations in N e unknown nodal displacement vectors. In view of a numerical solution method-for example, the Newton-Raphson method-it is often necessary to derive a linear approximation of equation (12). For a given state of vibration Z (n) in an iteration solution procedure, a neighboring state can be expressed by adding an increment,

Z (n + 1) = Z (n) + DZ (n) . ( 17
)
For a small increment, the non-linear equation ( 12) can be expanded into a set of Taylor series, and a set of linear incremental equations are obtained by retaining only the linear terms in the expansions, in the form

H (n) + K (n) DZ (n) = 0, (18)
where (n) . (19,[START_REF] Mello | Time finite element methods for large rotational dynamics of multibody systems[END_REF] The matrix K is often referred to as the tangent matrix in finite element analysis. Equations ( 17) and ( 18) represent an iteration process of the Newton-Raphson method, from which the steady state response Z can be calculated. For a converged solution, DZ (n) becomes nearly zero.

H (n) = A (n) Z (n) + G (n) + P, K (n) = A (n) + (1G/1Z)
It should be pointed out that the linearization procedure can be performed either on the variational equation ( 4) of Hamilton's weak principle or after the finite element approximation. Conceptually, these linear approximations to the non-linear algebraic system are equivalent, although they may yield different implications in the formulation of analytical equations and in numerical computation. This issue is similar to that involved in different implementations of the harmonic balance method [22]. For the one-degree-of-freedom problem considered in this investigation, all formulae are worked out in explicit forms, regardless of the stage of linearization. Furthermore, for unsymmetric piecewise-linear systems, the formulation can be greatly simplified by explicit expressions for the integrals involving the restoring force of the second spring. In addition, for a general parametric study, one can use various numerical procedures to obtain the solution of the nonlinear algebraic equations ( 12) for a variable parameter. For example, a bifurcation diagram of the unsymmetric piecewise-linear system with variation of stiffness ratio b can be easily traced by first seeking the solution for a small value of b of weak non-linearity and then increasing the parameter and obtaining the corresponding solution with a predictor-corrector method. These are the subjects of the following section on numerical calculations.

NUMERICAL CALCULATIONS

evaluation of integrals

Calculation of the solution of the non-linear equation (12) or the linearized equation ( 18) requires evaluation of the integrals of each time element in matrix A and vectors G and P in equations ( 14)- (16). The integral components of matrix A and vector P are constants and, therefore, they need to be calculated only once. In addition, if identical interpolation functions are used for different time elements, the corresponding integrals are identical. With polynomial interpolation functions (e.g., Lagrangian polynomials), each integral component of the matrix A can be made exact analytically without quadrature approximation.

The evaluation of the elements of vector G and matrix 1G/1Z,

G i = -b g u i + 1 u i N T i g(z) du, ( 21
)
1G i /1Z j = -b g u i + 1 u i N T i 1g(z)/1zN j du (i, j = 1, 2, . . . , N e ), (22)
requires a knowledge of the zeros of the function (z(u) -d) with the interval [u i , u i + 1 ]. Since the displacement z within the interval is interpolated by polynomials among the temporal nodes Z i , these zeros are the roots of a polynomial and therefore can be easily determined in numerical computation of a computer program. Given the m roots, j 1 , j 2 , . . . , j m , and the sign of function (z -d), s 1 , s 2 , . . . , s m + 1 , corresponding to the intervals [j 0 , j 1 ], [j 1 , j 2 ], . . . , [j m , j m + 1 ] with j 0 = u i and j m + 1 = u i + 1 , the integrals ( 21) and ( 22) are expressed, respectively, as a sum of integrals over the (m + 1) intervals in the form

G i = -b s m + 1 k = 1 h k g j k j k -1 N T i g(z) du, ( 23
)
1G i /1Z j = -b s m + 1 k = 1 h k g j k j k -1 N T i 1g(z)/1zN j du. ( 24
)
The corresponding step functions are given as

h k = 6
1 for s k q 0 (z e d),

0 for s k Q 0 (z Q d), (25)
For a piecewise-linear system, g(z) = h(zd) and 1g(z)/1z = h, with h to be the step function defined above. The integrals ( 23) and ( 24) become definite integrals of polynomials and, therefore, can be made exact analytically without numerical quadrature. Consequently, this will increase efficiency of numerical computation for the solution. In addition, the restoring force in the time interval of the ith element is determined only by the nodal displacement Z i of the corresponding element. Therefore, 1G i /1Z j = 0 for j $ i.

static condensation

The finite element method in the time domain has the crucial property that the dynamic response of a system is sought through propagation forward in time. As a result, the resulting tangent matrix K (20) is in a block-diagonal form and the set of linearized equations (18) can be solved efficiently. The solution of equation ( 18) can be obtained without actually assembling the global tangent matrix K and the residual vector H. This property is similar to that of finite elements in space used in solid mechanics. There exists a standard procedure of static condensation in the space finite element method, by which the elemental degrees of freedom that are not involved in satisfying the inter-element compatibility are eliminated before the element is assembled into the overall system equations, yielding fewer unknowns in the system equations [START_REF] Bathe | Finite Element Procedures in Engineering Analysis[END_REF]. A similar procedure for finite elements in time is described in reference [15] for obtaining the solution of linear equations (18) efficiently. It is briefly discussed as follows.

Consider the tangent matrix K, the displacement and force vectors, DZ and P, and the boundary conditions B of linearized equation (18), to be partitioned as

$ K bb K ib K bi K ii %6 DZ b DZ i 7 + 6 H b H i7 = 6 B b 0 7 , ( 26
)
where b and i correspond the temporal nodes at the ends and in the interior of current time element separately, and B b = {-p(u i ), p(u i + 1 )} T denotes the boundary conditions. The lower part of equation ( 26) can be solved for DZ i as

DZ i = -K -1 ii H i -K -1 ii K ib DZ b . ( 27
)
When this equation is substituted into the upper part of equation (26) for DZ i and the terms are collected, the remaining equation for DZ b can be written as

K bb DZ b + H b = B b , (28)
where the reduced tangent matrix and force vector are given by

K bb = K bb -K bi K -1 ii K ib , H b = H b -K bi K -1 ii H i . (29, 30)
When the next element is generated, the condensation procedure can be applied identically, resulting in a similar set of reduced equations. These two sets of reduced equations are further condensed by imposing the compatibility conditions of the displacement vector and the boundary conditions at the common node of the two elements, yielding the elimination of the common node. The above procedure is repeated for each remaining element, and the final set of equations is written as

$ K II K FI K IF K FF% 6 DZ I DZ F7 + 6
H I H F7 = 6 B I B F7 , (31)
where I and F denote the initial and final time nodes, respectively, and B I = -p I = -p(u I ) and B F = p F = p(u F ). Finally, for the periodic steady state response it is required that

DZ F = DZ I and p I = p F . ( 32
)
Therefore, the initial and final nodal displacements can be solved from equation (31) as

DZ I = DZ F = -(K II + K IF + K FI + K FF ) -1 (H I + H F ). ( 33
)
This procedure is the same as Gaussian elimination. Once the boundary nodes have been computed, then the condensed unknown nodes can be recovered by a back substitution process for each corresponding condensation. The procedure is continued for all nodes in the system, allowing a complete solution of the equations to be obtained.

parametric solution

The finite element in time method is ideally suited to parametric studies for obtaining solution diagrams, since it can step from a known state of vibration to a neighboring state corresponding to an incremental change in one of the governing parameters of the system. The solution of the non-linear algebraic equations (12) with a variable parameter defines a path in the corresponding state-parameter space, or the solution diagram. A method for following the path is discussed next, which employs a predictor to move along a path direction and a corrector to allow accurate path following.

Consider the solution diagram, defined by equation (12) as

F(Z, l) = A(l)Z + G(Z, l) + P = 0 (34)
for a variable parameter l, which can be the frequency ratio r, the stiffness ratio b, the clearance value d, or any parameter of interest in the system. Define w to be the arc-length along the solution path. Then, the derivative of equation ( 34) with respect to arc-length is given as

1F 1Z dZ dw + 1F 1l dl dw = 0, ( 35
)
where 1F/1Z is the tangent matrix K defined in equation ( 18) by linearization. The above derivative indicates where the path is going, allowing a step to be moved in that direction. The arc-length is defined to be

B dZ dw B 2 + 0 dl dw 1 2 = 1. ( 36
)
Suppose that the current point on the path corresponds to a solution Z (n) and a parameter value l (n) : then, the next point farther along the path of step size h is defined as (37,38) where DZ (n) and Dl (n) satisfy the derivative and unit arc-length equations ( 35) and (36); i.e.,

Z (n + 1) = Z (n) + hDZ (n) , l (n + 1) = l (n) + hDl (n) ,
DZ (n) = -00 1F 1Z 1 -1 1F 1l 1 (n) Dl (n) , Dl (n) = 21 >0B00 1F 1Z 1 -1 1F 1l 1 (n) B 2 + 1 1 1/2
. (39, 40)

This procedure is known as the Euler predictor, which is of first order [START_REF] Zangwill | Pathways to Solutions, Fixed Points, and Equilibria[END_REF]. Higher order predictors (e.g., Adams-Bashforth) are also available using higher derivatives of the solution equation (34).

A predictor is essentially an approximation method. The predicted solution can significantly drift off the solution path in one step or a sequence of predictor steps, and it must be corrected to return more closely to the solution path [START_REF] Zangwill | Pathways to Solutions, Fixed Points, and Equilibria[END_REF]. The classical method is a Newton-Raphson corrector with a sequence of iterations constrained to lie perpendicular to the predictor direction and passing through predicted solution Z (n + 1) and l (n + 1) . At convergence, the corrector gives rise to an accurate solution.

Let DZ (n + 1) and Dl (n + 1) be the corrector step from the current solution Z (n + 1) and l (n + 1) , satisfying the Newton-Raphson method such that

$ 1F 1Z 1F 1l % (n + 1) 6 DZ Dl 7 (n + 1) = -F (n + 1) . ( 41
)
The corrector direction is defined by

{DZ T Dl} (n) 6 DZ Dl 7 (n + 1) = 0. ( 42
)
On combining the above equations, the corrector step is the solution of the set of linear equations

$ (1F/1Z) (n + 1) DZ (n) (1F/1l) (n + 1) Dl (n) %6 DZ Dl 7 (n + 1) = 6 -F (n + 1) 0 7 . ( 43
)
The new value of the corrected point DZ (n + 2) and Dl (n + 2) is given by 1) . (44, 45)

Z (n + 2) = Z (n + 1) + DZ (n + 1) , l (n + 2) = l (n + 1) + Dl (n +
A finite number of corrector steps will be taken so that the norm of F is smaller than some selected small number. Repetition of the above predictor and corrector process will follow the solution path quite well, as shown in the numerical results presented in section 6.

STABILITY OF STEADY STATE SOLUTION

The stability of any periodic steady state solution obtained above is investigated by considering small perturbations of the solution and the boundary condition in equation ( 31) for the initial and final time nodes [15,17]; hence, equations (31) becomes

$ K II K FI K IF K FF% 6 dZ I dZ F7 = 6 -d p I d p F 7 . (46)
According to Floquet's theory, the stability of the system is determined by the eigenvalues of the transition matrix that relates the initial and final perturbations. This transition matrix happens to be a by-product of the FET procedure, and is readily found by rearranging equation (46):

6 dZ F d p F 7 = $ -K -1 IF K II K FI -K FF K -1 IF K II -K -1 IF -K -1 IF K FF% 6 dZ I d p I 7 . (47)
Therefore, the same solution algorithm can be used for both the periodic response and stability calculations. There is no need for a special effort in calculating the transition matrix, as in, for example, references [START_REF] Hsu | [END_REF]26]. This is useful in assessing the characteristics of the non-linear system. The stability of the steady state response is determined by the eigenvalues of the transition matrix defined in equation ( 47). If the eigenvalue l with the larger modulus satisfies =l= Q 1, then the periodic solution examined is asymptotically stable. When =l= q 1, the solution is unstable. Bifurcations occur if =l= = 1, resulting in qualitative changes in system dynamics [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF]. There are three types of bifurcations. The first type is a saddle-node bifurcation corresponding to the merging of two periodic solutions. It occurs when l = 1 and results in the classical jump phenomena. The second type is a flip bifurcation associated with period doubling of the response. It occurs when l = -1. Finally, if a pair of complex eigenvalues have modulus one, the so-called Hopf bifurcation happens, resulting in an amplitude modulated response of the system. For the system considered in this investigation, the first two types of bifurcations are observed as discussed in the results presented in the next section. It was previously shown that the Hopf bifurcation does not happen in a single-degree-of-freedom system with piecewise-linear stiffness [2].

RESULTS

A simple unsymmetric bilinear system is considered first, in which d = 0 and b = 4. The system is under a harmonic excitation P(u) = cos u. This system has been studied in reference [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF] using direct numerical integration and in references [7,9] using the harmonic balance method. The numerical results obtained by the finite element method are shown in Figure 2, where the steady state amplitude is defined as half of the peak-to-peak of the response. The frequency response is plotted as a function of the dimensionless frequency ratio r, with damping ratio z = 0•1. These results are found to agree very well with those obtained by using direct integration and harmonic balance method.

The FET method yields all possible steady state solutions, including the sub-periodic and the unstable solutions, for a wide range of frequency ratios. The order of the sub-periodic vibration is denoted by m in the figure . A bifurcation from a stable solution on the fundamental response to a stable solution on the sub-periodic of order 2 is found to occur at r = 2•27. The stability analysis reveals that it is a flip bifurcation. This bifurcation was investigated in reference [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF] by using a Poincare´map. Numerical calculation of the transition matrix in a harmonic balance analysis also showed the period doubling phenomena [9]. Resonant peaks of sub-periodic response appear regularly to the right of the fundamental resonance with increasing order of sub-periodic vibration. Furthermore, to the left of the fundamental resonance, the FET method shows a number of peaks occurring irregularly. This behavior was also investigated in reference [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF], where the correlation with the cyclic moment of the Poincare´point was shown.

The mth sub-periodic response is defined to have a period of oscillation equal to m times the fundamental period, i.e., T = 2mp, for the harmonic excitation. To obtain the sub-periodic response without modifying the FET formulation, one simply takes the excitation to be P(u) = cos mu and shifts the obtained solution with the frequency divided by m.

To show the convergence of the obtained solutions and to find an adequate number of time elements and an adequate order of interpolation functions, this problem has also been computed with a different number of time elements and different order of Lagrangian polynomials for the shape function. The results of the average amplitude obtained are listed in Table 1. It is observed that the use of six time elements and quartic polynomials is sufficient to obtain a well converged solution with three to four digit accuracy. Therefore, six time elements and quartic polynomials are used for the numerical results shown hereafter. The convergence development with subsequent iteration of the FET method for t = 1•2 is shown in Figure 3. When the stiffness ratio b increases, the strength of the non-linearity increases. The frequency response of the case of b = 9, which also agrees very well with previous results [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF]9], is shown in Figure 4.

As discussed in section 4.3, the effect of increasing the stiffness ratio can be fully investigated by seeking the solution diagram for the variable parameter b with the predictor-corrector method presented in section 4. unstable. This phenomenon seems to be a characteristic feature of this type of non-linear systems, and it is one of the possible mechanisms of occurrence of chaotic vibration of the systems [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF]9].

Cases involving a finite clearance value d are also examined. The bilinear system considered has a clearance d = 0•5, with the damping ratio z = 0•1. Two cases are studied for b = 9 and 99. In the first case, with the modest level of non-linearity, the fundamental resonance is shown (Figure 6) to occur in the neighborhood of frequency ratio of 1•519, which corresponds to the natural frequency of the undamped, zero clearance system [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF]. In the range of frequency ratios above 1•723, the system behavior is linear, as the maximum amplitude of vibration is less than the clearance value. It is well known that the piecewise non-linearity introduces a hardening effect in the system dynamics if b q 0. For the second case with strong non-linearity (b = 99), it is shown in Figure 6 that the hardening effect becomes highly pronounced. In the neighborhood of the 1•723 frequency ratio, a complex hardening phenomenon with a number of consecutive jumps occurs, as clearly shown in a regional plot (Figure 7). Multiple steady state solutions exist in a number of regions of frequency ratio. Stability analysis reveals that one of the multiple solutions is unstable and that a saddle-node bifurcation occurs at the point at which the unstable response merges with a stable response. In a lower frequency range (r = 1•221-1•285), the fundamental response becomes unstable. One of the real eigenvalues of the transition matrix is less than -1, indicating flip bifurcations at the transitions between the stable and unstable solutions.

In the last example, a piecewise-non-linear system with both clearance and continuous stiffness non-linearities is considered, in which d = 0•5, z = 0•1 and b = 9. The restoring force of the second spring is given by g(z) = (z -d) n , where n denotes the type of non-linearity of the second spring. Two types are examined here: n = 3, cubic; and n = 1•5, referring to the Hertzian spring. The frequency responses for these cases are shown in Figure 8. As expected, the hardening effect of the non-linear spring becomes stronger as the stiffness non-linearity represented by n increases. The overall characteristics of the response seem to be dictated by the clearance non-linearity. However, multiple solutions with jump transition exist in the case of n = 3. It is previously known that the combination of the clearance and continuous stiffness non-linearities can yield highly complex dynamics, including chaotic motions [START_REF] Kahraman | [END_REF]. The detailed study of the dynamic behavior of such systems will be deferred to another paper.

CONCLUSIONS

A method of finite elements in the time domain for determining the steady state vibration of unsymmetric systems with piecewise-linear/non-linear stiffness characteristics has been described in this paper. The method is based on a Hamilton's weak principle that allows for a simple choice of interpolation function and an easy problem formulation. The method offers a number of special features that can be utilized for accurate numerical results with only a few elements and for minimal computation effort. The stability of the steady state solution is directly determined as the transition matrix of Floquet's theory is readily available.

This method is applied to unsymmetric piecewise-linear and piecewise-non-linear systems and the numerical results are seen to provide excellent agreement with previous findings of a direct numerical integration method and a harmonic balance method. Both stable and unstable steady state solutions are determined. A parametric solution of the system dynamics is performed with a predictor-corrector method. This possibility permits one to investigate strongly non-linear systems with both clearance and continuous stiffness non-linearities.

The finite element in time method has the flexibility for adaptive change of the order of interpolation functions and the number of time elements. Therefore, the number of temporal nodes can be easily changed from one iteration step to the next, as only the linearized equations are assembled from each time element in each step. This fact may lead to a systematic and effective impelementation for obtaining steady state solutions over a wide range of system parameters and for obtaining solutions with greater accuracy at certain specified time locations. This method is currently under development by the author in applying the method to problems with discontinuities in the states and the state equations (such as systems with pre-loaded springs and dry friction).

Figure 1 .

 1 Figure 1. (a) The physical model of an unsymmetric system. (b) Unsymmetric piecewise-linear restoring force.

Figure 2 .

 2 Figure 2. A frequency response plot for d = 0, z = 0•1 and b = 4. ----, Stable; ----, unstable.

Figure 3 .

 3 Figure 3. The convergence of the finite element method; d = 0, z = 0•1, b = 4 and r = 1•2. The integers 1, 2, . . . , 7 correspond to the following iteration steps: ----, 1; ••••, 2; • -• -, 3; --, 4; • • •, 5; ^, 6; W, 7.

3 .

 3 The case with z = 0•1, d = 0 and r = 2•89 is analyzed as a parametric study of the effect of the stiffness ratio on the steady state response. As shown in Figure 5, period doubling occurs at b = 1•9858, with the appearance of a stable sub-periodic response of order m = 2. When the system becomes more non-linear as b increases, there are regions of b, e.g., at b = 22 or b = 50, in which this sub-periodic response is unstable. It can be shown in some regions that all responses are

Figure 4 .

 4 Figure 4. A frequency response plot for d = 0, z = 0•1 and b = 9: ----, Stable; ----, unstable.

Figure 5 .

 5 Figure 5. The steady state response with increasing stiffness ratio b, where z = 0•1, d = 0 and r = 2•89: ----, Stable; ----, unstable.

Figure 6 .

 6 Figure 6. The steady state solution with z = 0•1 and d = 0•5: ----, Stable; ----, unstable.

Figure 7 .

 7 Figure 7. Multiple jumps in a strongly non-linear system: z = 0•1, d = 0•5 and b = 99: ----, Stable; ----, unstable.

Figure 8 .

 8 Figure 8. The response for both clearance and continuous stiffness non-linearities; z = 0•1, d = 0•5 and b = 9: ----, Stable; ----, unstable.

Table 1

 1 Average amplitudes obtained with different numbers of time elements and different orders of polynomials

			Order
		ZxxxxxxxxxxxxxxCxxxxxxxxxxxxxxV
	Element	3	4	5	6
	4	2•53389	2•52567	2•52664	2•52612
	6	2•52333	2•52598	2•52601	2•52594
	8	2•52479	2•52608	2•52594	2•52596

ACKNOWLEDGMENTS

This work is being carried out as part of a project supported by the National Science Foundation under grant No. MSS-9308663, which is gratefully acknowledged. The author would also like to than Professor Fred Tasker for his input as to the finite element in time method.


M A Christon

D W Roach

The numerical performance of wavelets for PDEs: the multi-scale finite element

The research summarized in this paper is part of a multi-year effort focused on evaluating the viability of wavelet bases for the solution of partial differential equations. The primary objective for this work has been to establish a foundation for hierarchical/wavelet simulation methods based upon numerical performance, computational ef®ciency, and the ability to exploit the hierarchical adaptive nature of wavelets. This work has demonstrated that hierarchical bases can be effective for problems with a dominant elliptic character. However, the strict enforcement of orthogonality in the usual L 2 sense is less desirable than orthogonality in the energy norm. This conclusion has led to the development of a multi-scale linear ®nite element based on a hierarchical change-of-basis. This work considers the numerical and computational performance of the hierarchical Schauder basis in a Galerkin context. A unique row-column lumping procedure is developed with multi-scale solution strategies for 1-D and 2-D elliptic partial differential equations.

Introduction

Wavelets are a relatively new mathematical tool that dissect data, functions, and differential operators into components with an associated resolution that is matched to the scale of each component. The name ``wavelet'' or ``ondelette'' was coined in the early 1980s by French researchers Morlet, Arens, Fourgeua, Giard and Grossman [START_REF] Morlet | Wave propagation and sampling theory[END_REF][START_REF] Morlet | NATO ASI Series[END_REF][START_REF] Grossmann | Decomposition of hardy functions into square integrable wavelets of constant shape[END_REF]. However, functions with the attributes of wavelets have been known for almost 100 years. Meyer [START_REF] Meyer | Wavelets Algorithms and Applications[END_REF] points out that there are seven primary origins for wavelets that date from around 1930 with the Haar wavelet dating back to 1909. However, the literature from this era does not explicitly use the term ``wavelet''.

In recent years, the application of wavelet bases to the solution of partial differential equations (PDEs) has evolved to the point where there are a number of com-peting formulations that include, but are not limited to wavelet-Galerkin, wavelet-collocation, and reproducing kernel methods. Despite the growing number of waveletbased formulations and solution algorithms for PDEs, the ®eld is still very new, and many technical issues remain. That is to say, there is promise in the approach, but there is a clear need for fundamental research that characterizes the numerical and computational performance of wavelets for the solution of partial differential equations.

A complete historical perspective of wavelet based approaches to solving PDEs is beyond the scope of this paper. However, a survey of the literature on wavelet based methods for solving PDEs may be found in Christon et al. [START_REF] Christon | An investigation of wavelet bases for grid-based multi-scale simulations ± Final Report[END_REF]. Over the past 8 years, there has been a great deal of attention paid to the wavelet-Galerkin method because of the generality that the Galerkin method provides, and the ease with which alternative bases may be implemented and tested. In addition, the trend in recent years has been away from traditional Daubachies wavelets with L 2 orthogonality and towards wavelets that are constructed with weaker orthogonality properties yielding biorthogonal and semi-orthogonal wavelets (see for example Dahmen and Micchelli [START_REF] Dahman | Biorthogonal wavelet expansions[END_REF]). This approach has led to the idea that wavelets may be most useful for solving PDEs if they are used to simply ``complete'' a subspace with the concomitant relationship between grid scales. A related approach that combines a spectral method with biorthogonal wavelet bases adapted to non-square domains and higher spatial dimensions may be found in Canuto et al. [START_REF] Canuto | The wavelet element method part ii: Realization and additional features in 2d and 3d[END_REF].

In this work, orthogonality with respect to the energy norm is used with a nodal change-of-basis to develop multi-scale ®nite elements that make use of a hierarchical Schauder basis. In 1-D, the uniform stability of the Schauder basis in H 1 yields uniformly bounded condition numbers, i.e., independent of re®nement level (grid resolution), for elliptic operators in a Galerkin form. This property was previously described by Zienkiewicz [START_REF] Zienkiewicz | The Finite Element Method[END_REF] as the ``diagonality'' feature for one dimensional hierarchical shape functions. However, the relationship between uniform stability, a nodal change-of-basis and the wavelet transform was not considered. In the ensuing discussion, the theoretical issues surrounding the hierarchical Schauder basis are presented, and the details of a uniformly stable basis are made concrete. The use of the wavelet transform in conjunction with the assembly of the mass and stiffness operators for a hierarchical basis is presented in the framework of 1-D and 2-D multi-scale ®nite elements. In addition, a comparison between the hierarchical Schauder basis and the linear ®nite element basis is presented. The application of an ad-hoc rowcolumn lumping procedure for the hierarchical basis is outlined, and its effectiveness for solving elliptic boundary value problems in one and two-dimensions is demonstrated.

The hierarchical Schauder basis

In this section, the formulation issues associated with the hierarchical Schauder basis are outlined. For this development, the boundary value problem under consideration is Àer 2 u cu f on X ; 1 with homogeneous Dirichlet boundary conditions u 0 on C : 2

The weak form of the problem, after introducing the ®nite dimensional subspace, V k , is

X eru k Á rU k dC X cu k U k dC C f U k dC : 3
Here, u k P H 1 0 , and U k is a ®nite dimensional basis for the subspace V k at scale k. In matrix form, this may be written as

M k K k fu k g fF k g 4
where M k is the mass matrix, and K k is the stiffness matrix associated with scale k. Alternatively, in a bilinear form, aU k ; U k u k FU k . Now, consider an alternative basis for V k ,

W k W T k U k ; 5
where W k is a nonsingular N Â N matrix. From this transformation, u k Du T W may also be found by solving

aW k ; W k Du FW k 6
for the multi-scale representation of the ®eld, Du. The linear system, A W k , that derives from the change of basis may also be viewed as the linear system obtained by preconditioning

A U k M k K k with W k , i.e., A W k W T k A U k W k : 7
Note that the superscript W (U) indicates that W k (U k ) is used as both the test and trial function in obtaining the weak form.

The matrix W k is a wavelet transform and provides the mechanism for decomposing a ®eld into multiple scales. In order to understand the role of the wavelet transform, consider that is a one-sided sequence of nested ®nite-dimensional subspaces of H, a Hilbert space, such that V k H. Now, de®ne the ``coarse-grid'' wavelet subspace, W 0 V 0 and, for k ! 1, choose W k so that completion of subspace V kÀ1 may be achieved via a direct sum,

V k V kÀ1 È W k : 9
Here È denotes a direct sum, but not an orthogonal direct sum.

Beginning with Eq. ( 5), the wavelet transform may be written recursively in terms of a two-scale transform as

W k T k T kÀ1 0 0 I kÀ1 ! Á Á Á T 1 0 0 I 1 ! 10
where I k is the N W Â N W identity matrix with N W dimW k .
Here, T k is a two-scale transform such that

U kÀ1 W k & ' T T k U k : 11
As an aside, one very important aspect of this formulation is that it does not require the decomposition matrix W À1 k . This is signi®cant because it permits the relaxation of strict orthogonality in the selection of the wavelet bases.

One choice for W k that satis®es Eq. ( 9) is the hierarchical Schauder basis [3, 10±13]. This basis consists of the elements of the linear ®nite element basis at multiple scales where the relationship between W k and U k is simply

w j k / 2jÀ1 k ; j 1; . . . ; Nnp k ; 12
where U k is the linear ®nite element basis. Here the superscript j indicates the node number, k indicates the scale, and Nnp k is the number of node points at scale k.

The Schauder basis is illustrated in Fig. 1.

With the relationship between basis elements in the Schauder basis given by Eq. ( 12), the construction of the two-scale transformation may be written in terms of the grid spacing. In order to de®ne the wavelet transform, let H k be the Nnp k À 1ÂNnp kÀ1 À1 matrix H k h j;j H k j H ;j and let G k be the Nnp k À 1 Â Nnp kÀ1 matrix

h j;j H k D j H k D j kÀ1 ; j H 2j À 1 1; j H 2j D j H 1 k D j kÀ1 ; j H 2j 1 0; otherwise V b b b b b b b b b b X 14 g j;j H k d j H ;2jÀ1 : 15
In one dimension, the element length from

x jÀ1 k to x j k is de®ned as D j k x j k À x jÀ1 k .
The action of h j;j H k applied at scale k is to average, while g j;j H k acts to inject at scale k. This can be seen in Fig. 1 where the relationship between basis elements is given by

/ j kÀ1 j H h j;j H k / j H k 16 w j k j H g j;j H k / j H k :
17

With the Schauder basis de®ned, attention is turned to the question of stability. To begin the discussion on stability, several de®nitions are required.

De®nition 1 Given a basis U / 1 ; /

V 0 & V 1 & Á Á Á & V k Á Á Á & H ; 20 such that U k is a basis for V k .
De®nition 3 Given a nested sequence of subspaces of a Hilbert space, H, and their respective bases fU k g kPZ (normalized in H), then fU k g kPZ is uniformly stable in H if there exists constants a; b > 0 independent of k such that

aku k k 2 `2 ku T k U k k 2 H bku k k 2 `2 k 0; 1; . . . 21
With these de®nitions, the relationship between the stability of a basis and the condition number of an oper-ator discretized in the basis is outlined. In speci®c, it will be demonstrated that a uniformly stable set of bases yields uniformly bounded condition numbers for the discrete operators independent of the re®nement level.

Lemma 1 Let U k be a basis for a subspace, V k , of a Hilbert space H such that

A U k is positive de®nite. Let a kA U À1 k k À1 2 and b kA U k k 2 . Then aku k k 2 `2 ku T k U k k 2 H bku k k 2 `2 ; 22 and condA U k b a : 23
Proof The condition number of a symmetric matrix is given by condA

U k kA U k k 2 kA U k À1 k 2 .
For the stability condition, note that

ku T k U k k 2 H ku k k 2 `2 u T k U k ; u T k U k À Á H u T k u k u T k A U k u k u T k u k : 24 Moreover, for symmetric matrices b kA U k k 2 qA U k max u k u T k A U k u k u T k u k

25

which gives us the upper bound. For the lower bound, observe that for symmetric positive de®nite matrices

a kA U À1 k k À1 2 max u k u T k A U À1 k u k u T k u k 3 À1 min u k u T k A U k u k u T k u k : 26
The proof of Lemma 1 states that the best stability bounds are provided by the smallest and largest eigenvalues of the discretized matrix ± a well known result. Now consider the stability of a linear ®nite element basis, shown in Fig. 1a for k 2, where H L 2 0; 2 and V k spanf/ j k ; j 1; . . . ; 2 k1 À 1g; k 0; 1; . . . ; 27

U k / 1 k ; / 2 k ; . . . ; / 2 k1 À1 k T : 28 Lemma 2
The sequence of bases U k for the linear ®nite element basis is uniformly stable in L 2 0; 2.

Proof Let the discretized matrix A U k be a diagonally scaled ®nite element mass matrix. In 1-D this is a symmetric tridiagonal matrix consisting of 1 on the diagonal and 1=4 on the non-zero off diagonals. Note that for the ensuing proof, U k has been normalized in a manner that is equivalent to diagonal preconditioning applied to the ®nite element mass matrix. Using Gershgorin's theorem, the eigenvalues for A U k lie in the interval 1=2; 3=2 giving

1 2 ku k k 2 `2 ku T k U k k 2 L 2 3
2 ku k k 2 `2 ; 29 and condA U k 3. Therefore, the sequence of bases, U k , is uniformly stable in L 2 .

The stability associated with the ®nite element basis in L 2 is consistent with the empirical observation that the consistent mass matrix is well behaved in terms of its condition number. In practice, this is re¯ected in the ability to easily solve mass matrix dominated problems with simple iterative techniques.

Theorem 1 (Poincare Â-Friedrich) For bounded domains, the H m norm and semi-norm are equivalent in the sense that there exists constants a; b > 0 such that ajuj H m kuk H m bjuj H m for u P H m 0 : 30

Recall that the norm and semi-norm for the Sobolev space H 1 for one space dimension are

kuk H 1 ku H k L 2 kuk L 2 juj H 1 ku H k L 2 :
Lemma 3 After re-normalization with respect to the H 1 semi-norm, the sequence of bases U k is not uniformly stable in H 1 0 0; 2.

Proof Using the equivalence of the norm and semi-norm, a sequence of vectors u k is constructed for which uniform bounds are not possible. If

u k 1; 1; 1; . . . ; 1; 1 z}|{ 2 k1 À1
for k 0; 1; 2; . . . ; 31

then ku k k 2 `2 2 k1 À 1 and ju T k U k j 2
H 1 1 for k 0; 1; 2; . . . ; 32

and stability for the basis can not be demonstrated. Therefore, using the equivalence of the H 1 norm and semi-norm, there are no uniform bounds for the sequence of bases U k in H 1 0 . Next, the hierarchical Schauder Basis is considered where U 0 corresponds to the linear ®nite element basis described above. Let W k span / j k ; j 0; 2; 4; . . . ;

2 k1 À 2 n o 33 V k V 0 a k j1 W j 34 W k / 0 0 ; / 0 1 ; / 2 1 ; . . . ; / 0 k ; / 2 k ; . . . ; / 2 k1 À2 k T : 35
The hierarchical basis for k 2 is shown in Fig. 1.

Lemma 4 After re-normalization with respect to the H 1 semi-norm, the sequence of bases W k for the Schauder basis is uniformly stable in H 1 0 0; 2.

Proof The H 1 semi-norm of W k is equivalent to the L 2 norm of the orthogonal Haar wavelet basis. The Haar basis is an orthogonal basis for L 2 0; 2 and thus A W k is the identity with respect to the semi-norm for k 0; 1; 2; . . . ; I.

Lemma 5 After re-normalization with respect to the L 2 norm, the sequence of bases W k for the Schauder basis is not uniformly stable in L 2 0; 2.

Proof De®ne the sequence u k which gives the hat function centered at one with support 1 À 2 k ; 1 2 k as u k 1; À1=2; À1=2; 0; À1=4; À1=4; 0; 0; 0; 0; À1=8; À1=8; 0; 0; 0; . . . : 36

Thus,

ku k k 2 `2 ! 1 ; 37 and ku T k W k k 2 L 2 2 3 1 2 k : 38 Because ku T k W k k 2 L 2 3
0 as k 3 I the stability condition in Eq. ( 3) can not be met.

The multi-scale finite element

This section outlines the 1-D and 2-D multi-scale ®nite elements. The multi-scale elements are based upon the hierarchical Schauder basis with a change-of-basis incorporated at the element-level in order to make use of the well-known ®nite element assembly procedure [START_REF] Hughes | The Finite Element Method[END_REF].

Before embarking on a description of the 1-D multiscale ®nite element, a brief review and interpretation of the multi-scale transformation is presented. Recall from Eq. ( 5) that W k W T k U k where W k is the multi-scale transformation matrix. In order to make this transformation concrete, consider the following example. Beginning with a 1-D grid consisting of 5 grid points and 4 linear ®nite elements, the nodal basis will be decomposed into a coarse-grid consisting of two elements and the associated ``psuedo-wavelets''. This decomposition is shown schematically in Fig. 2.

Remark 1

The term pseudo-wavelet is used here to indicate that the elements of W k in the hierarchical Schauder basis do not possess the property that their zeroth moment is zero. However, the elements of W k used to complete the subspace V kÀ1 are uniformly stable in H 1 . Thus, the term pseudo-wavelet seems appropriate. In this example, the wavelet transform,

W 1 W T 1 U 1 , is / 3 0 w 2 1 w 4 1 V X W a Y 1=2 1 1=2 1 0 0 0 0 1 P R Q S / 2 1 / 3 1 / 4 1 V X W a Y ; 39
where the subscript indicates the scale, and the superscript indicates the node number. From this, it is clear that the wavelet transform performs an averaging procedure to obtain the coarse-grid basis elements, U 0 , and an injection to obtain the coarse-grid pseudo-wavelets, W 1 .

The decomposition of nodal variables, u, may be accomplished with the wavelet transform, but its inverse is required to obtain the coarse-grid coef®cients, i.e., Du W À1 u, where Du is the multi-scale component of the ®eld. Relying on the inverse wavelet transform is impractical because the orthogonality constraint between the wavelets and scaling functions has been relaxed in the hierarchical Schauder basis. However, incorporating the hierarchical Schauder basis at the element level yields solution algorithms that compute the multi-scale solution directly and rely only on the reconstruction algorithm, i.e., u WDu. Thus, given a multi-scale representation of the ®eld that corresponds to the multi-scale basis elements in Fig. 2, the reconstruction algorithm is simply

u 2 1 u 3 1 u 4 1 V X W a Y 1=2 1 0 1 0 0 1=2 0 1 P R Q S u 3 0 Du 2 1 Du 4 1 V X W a Y : 40
As an aside, the wavelet transform is comprised of two components and both may be viewed in terms of a discrete convolution. Using the nomenclature introduced above for the wavelet transform, T 1 H 1 jG 1 , where

H 1 1=2 1 1=2 V X W a Y ; G 1 1 0 0 0 0 1 P R Q S : 41
In this simple example of a two-scale decomposition, W 1 T 1 , but in general, the wavelet transform is computed recursively according to Eq. ( 10).

The 1-D multi-scale element

The description of the multi-scale element begins with the linear ®nite element for which the shape functions are

N i 1 2 1 À n i n : 42
Here, n is the natural coordinate, n i is the nodal value of the natural coordinate, i 1; 2, and À1 n 1.

The concept of scale is introduced at the element level by injecting degrees-of-freedom (DOF) that are supported by the pseudo-wavelets of the hierarchical Schauder basis. The 1-D multi-scale element is shown in Fig. 3 where a single ``internal'' degree-of-freedom located at n 0 in the element is introduced at Scale-1, two DOF are introduced at Scale-2, and four DOF at Scale-3.

At Scale-1, the pseudo-wavelet is

w 1 n 1 n if À 1 n 0 1 À n if 0 n 1 .

& 43

More generally, the psuedo-wavelets for the multi-scale DOF may be written in terms of the translates and dilates of wn as

w k n w ñ ; 44 where ñ 2 kÀ1 1 n À 2j À 1 ; 45 and 2 2Àk j À 1 n 2 2Àk j 2 1Àk À 1 2 2Àk j 2 1Àk À 1 n 2 2Àk j 2 2Àk À 1 : 46
Here k indicates the scale, and j indicates the translates in the element parametric space (À1 n 1). The derivatives of the shape functions yield constant functions that are orthogonal to the derivatives of the pseudo-wavelets at all scales. The derivatives of the pseudo-wavelets yield Haar wavelets as illustrated in Fig. 3, and at any given scale they are orthogonal with the derivatives at all other scales in the 1-D multi-scale element.

The reconstruction algorithm may be viewed as an element-by-element procedure that relies only on the multiscale information in each element. The reconstruction is shown schematically in Fig. 3 where the DOF located at n À3=4 is computed as a linear combination of the detail, Du, at scales 0À3. Here, a canonical node-numbering scheme is used where the node numbers are n 1 ; n 2 at Scale-0, n 3 at Scale-1, n 4 ; n 5 at Scale-2, etc. as shown in Fig. 3.

With this numbering scheme, the reconstruction for the DOF located at n 6 may be written as

u 6 N 1 u 1 0 N 2 u 2 0 w 3 1 Du 3 1 w 4 2 Du 4 2 w 6 3 Du 6 3 ; 47
where each of the basis elements is evaluated at n À3=4. Thus, the reconstruction algorithm begins with the interpolant of the coarse-grid solution and injects re®nements, or detail, up to the desired scale. The reconstruction may be written more generally as

u 22 kÀ1 j k Nnpe i1 N i u i 0 Nscale k1 w 22 kÀ1 j k Du 22 kÀ1 j k 48
where the basis elements (N i , and w k ) are evaluated at locations corresponding to the DOF location in the parent element.

The 1-D multi-scale operators

The computation of the stiffness for the model problem, Eq. ( 1), is a straight-forward procedure that begins with the coarse-grid stiffness. At the element level, i.e., Scale-0, this is simply

K e 0 e h 1 À1 À1 1 ! ; 49
where h is the node-spacing for the coarse-grid. Making use of orthogonality, the diagonal stiffness entries for the scale DOF, i.e., for k > 0, are

K 22 kÀ1 j k 2 k1 e h : 50
The resulting multi-scale element-level stiffness matrix is

K e e h 1 À1 À1 1 4 8 8
. . .

2 k1 P T T T T T T T T R Q U U U U U U U U S : 51
With this form of the element stiffness, only the coarsescale terms contribute to the element assembly procedure since all internal DOF simply require a point evaluation and do not rely on information outside the element. Here, a uniform mesh re®nement was applied to inject the scale DOF, but this is not a necessity for the multi-scale element.

The mass matrix computation is somewhat more involved because it requires the calculation of inner-products that involve the basis elements across scale. The computation of the multi-scale mass matrix consists of the following

52

This series of element-level integrals leads to the ®ngerdiagonal matrix structure shown in Fig. 4 with the inset coarse-grid element mass matrix.

A row-column lumping procedure is shown schematically in Fig. 4 for the multi-scale DOF at n 7 in the element. Note that the traditional row-sum lumping for the coarsegrid mass matrix can be used for the entries corresponding to the coarse-grid, but this type of mass lumping can not be used for the multi-scale DOF in the element. This concept is discussed further in the comparison between the multi-scale and linear ®nite elements.

A multi-scale algorithm

As an example of a multi-scale algorithm, consider the model problem, Eq. ( 2), with e 1, c 0, f x 1, and essential boundary conditions, u0 0; u1 0. In Fig. 3. Basis elements and their piecewise derivatives for the one-dimensional multi-scale element with three re®nement scales. The derivatives have been scaled by 1=k for scales k 1; 2; 3 this example, three scale solutions are computed. Scale-0 corresponds to the coarse-grid solution using two elements. Scale-1 corresponds to the injection of one multiscale DOF per element, while Scale-2 corresponds to the injection of two multi-scale DOF per element. This may be seen in Fig. 5 where the multi-scale DOF are shown relative to the elements (e 1 ; e 2 ) of the one-dimensional grid. At the ®nest resolution, the multi-scale solution corresponds to the solution associated with 8 linear ®nite elements. The multi-scale algorithm proceeds as follows. 4. Repeat 1±3 until e is smaller than some user-speci®ed criteria. 5. Perform the element-by-element multi-scale reconstruction using Eq. (48).

For the example problem, the multi-scale solution for k 1 corresponds to Du 4 1 Du 5 1 1=32, and for k 2, the Du 6 2 Du 9 2 1=128, and Du 7 2 Du 8 2 1=128. After the reconstruction procedure, the scale DOF yield solution values that interpolate the exact solution ± a result that is expected for linear problems. Similar results have been obtained for problems with non-linear source terms, and for problems with inhomogeneous essential and natural boundary conditions.

The algorithm presented for the 1-D multi-scale element possesses the property that all scale injection relies only on element-local data and does not require a re-solve of the coarse-grid problem to improve the solution. For problems that involve a signi®cant mass-matrix contribution, the correction procedure outlined in the following sections is required when the row-column lumped mass is used.

The 2-D multi-scale element

Attention is now turned to the 2-D multi-scale element. As in the 1-D case, the 2-D bilinear element provides the element-level components of the global basis functions. To begin, Fig. 6 shows a four-patch of bilinear ®nite elements with the injected multi-scale DOF corresponding to k 1. The con®guration of the multi-scale DOF in the parent element is shown in Fig. 7. Like the shape functions, the components of the pseudo-wavelets at the element level take on a value of 1 at the DOF location, and they are zero at all other node locations. where i 1; 2; 3; 4, and À1 n; g 1. At the ®rst scale, k 1, the psuedo-wavelets are w 5 1 1 À jnj1 À jgj; w 6 1 1 2 1 n1 À jgj; w 7 1 1 2 1 À jnj1 g; w 8 1 1 2 1 À n1 À jgj; 56 w 9 1 1 2 1 À jnj1 À g :

In a more general way, the pseudo-wavelets for the multi-scale DOF may be written in terms of the translates and dilates of the basis functions at scale k 1. The pseudo-wavelets in two-dimensions are w m k n; g w m 1 ñ; g ; 57 where ñ 2 kÀ1 1 n À 2j À 1 g 2 kÀ1 1 g À 2j À 1 ; 58

j k À 1 2 2Àk j À 1 ñ 2 2Àk j 2 2Àk À 1 2 2Àk j À 1 g 2 2Àk j 2 2Àk À 1 :

59

Here k indicates the scale, j indicates the translates in the element parametric space (À1 n; g 1), and m 5; 6; 7; 8; 9 for the element-local numbering of the pseudo-wavelets. With the basis elements de®ned this way, the use of recursion at the element level as indicated in Fig. 6 can be used to automate the computation of the mass and stiffness operators with a given scale of resolution, Nscale.

There are several key points regarding the 2-D multiscale element. First, the orthogonality of the derivatives of the pseudo-wavelets is not preserved in two dimensions ± even on an orthogonal grid. The ®nger-diagonal matrices that arise from this discretization can lead to excessive storage costs if the matrices are used without threshholding, row-column lumping, or element-by-element procedures. Surprisingly, the row-column lumping can be applied to both the multi-scale mass and stiffness operators, albeit only for the rows and columns corresponding to the multi-scale DOF. The 2-D multi-scale element is compatible with many h-adaptivity strategies being implemented in ®nite element codes today, however, a detailed discussion of this aspect is beyond the scope of this paper. The ability to use this element as a change-of-basis preconditioner is just beginning to be explored. Finally, the numerical performance of the 1-D and 2-D multi-scale element in terms of dispersive behavior, rate of convergence, etc. is identical to the linear (bilinear) element since any multi-scale solution can be cast in terms of the reconstructed solution in the ®nite element basis at the ®nest grid scale.

Multi-scale vs. linear elements

In this section, a comparison of the storage, matrix conditioning, and computational complexity for the hierarchical Schauder and linear ®nite element bases is presented in the context of the model elliptic problem developed in Sect. 2. The approximations inherent in using a lumped multi-scale mass matrix are discussed, and an ef®cient row-column lumping procedure with a lumped-mass correction algorithm is presented to illustrate that only the largest wavelet coef®cients need be corrected in the approximate row-column lumped-mass solution. The row-column lumping procedure is applied to both the multi-scale mass and stiffness matrices in 2-D, and the lumped-mass correction algorithm is extended to account for the approximations introduced by the rowcolumn lumped mass and stiffness matrices. Although a uniform discretization is considered in the ensuing discussion, it should be noted that in 1-D the hierarchical Schauder basis retains all of its properties for a non-uniform grid. However, in two-dimensions, orthogonality in the H 1 sense is lost even for uniform grids.

1-D comparison

The non-zero ®ll pattern for the mass and stiffness matrices are shown in Fig. 8 for both the hierarchical Schauder and the linear ®nite element bases. The formulae for the number of non-zeros associated with the mass and stiffness operators for both the hierarchical Schauder and ®nite element bases are shown in Table 1. Because the hierarchical Schauder basis diagonalizes the stiffness, Fig. 8. The non-zero entries in the mass and stiffness matrices for the hierarchical Schauder and linear ®nite element bases for a mesh with 64 nodes only Ndof (number of DOF) words of storage are required for the diagonal, although in practice a coarse-grid ®nite element matrix is required for the multi-scale element. The consistent mass matrix, although sparse, requires increased storage relative to the tri-diagonal ®nite element mass matrix due to the ®nger-diagonal structure. In comparison to the linear ®nite element, the storage for the multi-scale mass matrix scales with both the number of DOF, Ndof, and the scale, k.

As demonstrated in Sect. 2, the ®nite element mass and the hierarchical Schauder stiffness matrices are both well conditioned. In fact, the hierarchical Schauder stiffness is ideal in 1-D since the nodal change-of-basis results in a diagonal stiffness operator when the full wavelet transformation is applied as a preconditioner. In contrast, the ®nite element stiffness and the Schauder mass matrices are both poorly conditioned. To be more precise, the ®nite element basis is uniformly stable in L 2 , but not in H 1 . However, the hierarchical Schauder basis is uniformly stable in H 1 , but not in L 2 .

In order to illustrate the differences in stability between the ®nite element and the hierarchical Schauder bases, consider the condition numbers associated with the mass and stiffness matrices as shown in Table 2. Here, multiple levels of mesh re®nement are considered with k 1 corresponding to a mesh with 3 nodes (2 elements). The growth of the condition numbers for K k is seen to be proportional to Oh À2 ) for the linear ®nite element bases, while the condition numbers for the mass matrix are bounded asymptotically at 3. The condition number associated with the combined ®nite element mass and stiffness (M k K k ) is dominated by the stiffness since e c 1 in this case. In contrast, the mass matrix for the hierarchical Schauder basis yields a condition number that grows approximately as Oh À3=2 while the condition number for the stiffness is uniformly bounded at 1. The combined mass and stiffness operator (M k K k ) for the Schauder basis is also uniformly bounded indicating the dominance of the stiffness.

Although the multi-scale mass matrix associated with the hierarchical Schauder basis is not as well conditioned as the ®nite element mass matrix, its primary drawback is the storage associated with the ®nger-diagonal ®ll pattern. However, the entries in the ®nger diagonals decrease in amplitude with increasing scale suggesting that the relative importance of the matrix entries become negligible with increasing scale. This effect is re¯ected in the integrals for the multi-scale mass matrix in Eq. ( 52). Therefore, either a mass lumping procedure or threshholding based on the relative size of the entries can be used to limit the matrix ®ll-in. In order to test this idea, an ad-hoc procedure for mass lumping is considered. In a ®nite element setting, mass lumping refers to the row-sum lumping procedure used to obtain a diagonal mass matrix. However, in the context of a multi-scale basis, the physical interpretation of mass lumping is not simple, and the idea is perhaps less well founded for the multi-scale element than for the linear ®nite element.

Regardless of these issues, experimentation with mass lumping strategies has suggested that a viable lumping procedure for the multi-scale mass matrix is to sum the values to the left and above the diagonal entry to the diagonal as shown in Fig. 4. The use of the row-column lumping procedure in Algorithm 1 permits the point-wise evaluation of the multi-scale DOF for problems where the mass matrix is required, but the subsequent reconstruction procedure can only yield an approximate solution. However, inspection of the resulting approximate solution, Du, obtained with the lumped multi-scale mass indicates that a small correction to the largest wavelet coef®cients can dramatically improve the solution quality. In other words, the row-column lumping procedure cannot be used over the entire range of scales without degrading the solution accuracy, i.e., the in¯uence of mass lumping is ``felt'' more at the coarser grid scales. A lumped-mass correction algorithm has been developed to account for the approximate nature of the solution and proceeds as follows.

Algorithm 2 Lumped-mass correction algorithm 1. Row-column sum lump the multi-scale mass matrix resulting in a diagonal matrix, Ã M L k K k . Here, the superscript L indicates a lumped operator. For the multiscale element, only those entries in the mass matrix associated with the multi-scale DOF are lumped. 4. Permute the rows and columns of A to match the permutation of Dũ using a consistent multi-scale mass only for the Dũ N unknowns. Here, the permuted matrix, A M k K k is formed with the consistent mass matrix, M k . 5. Solve the smaller system using the consistent mass matrix for Du N , i.e., the coarse-grid correction,

(k) Ndof M k K k M k K k M k K k M k K k 1 3
A N Du N f N À A R Dũ R : 61
Here, A N is the operator for the coarsest grid, and A R corresponds to the rows of A associated with Dũ R but formed with the consistent mass matrix. An ef®cient means of computing the matrix-vector product, A R Dũ R is through the use of the wavelet transform in Eq:7.

Multiplication of Du by the matrix W is equivalent to the application of the wavelet transform, i.e., the reconstruction algorithm, which can be implemented in an O(Nnp) operation algorithm. It should be noted that multiplication by W T is not equivalent to the inverse wavelet transform, i.e., W T W T I. However, the properties of W T permit implementation of a matrix-vector product using convolution techniques similar to the reconstruction algorithm.

In practice, the DOF associated with the coarse-grid tend to be the largest and are the ones that need to be corrected. Thus, the permutation step can be replaced with a solve for the coarse-grid DOF followed by the injection of the multi-scale DOF. This algorithmic approach is exactly what the multi-scale ®nite element does since the coarsegrid solution is computed with only the injection of scale relying on the hierarchical Schauder basis. Figure 9 illustrates the effectiveness of the lumped-mass correction algorithm for the model problem with e 1, c 1, and f Àx 8 x À 2 3 . In Fig. 9a, the exact solution is shown with the approximate solution computed with the lumped multi-scale mass matrix. In Fig. 9b, the approximate solution has been improved by correcting only the three largest wavelet coef®cients.

2-D comparison

In two dimensions, the behavior of the hierarchical Schauder basis changes rather signi®cantly relative to the 1-D basis. The most apparent changes are the additional ®ll-in for the multi-scale mass matrix and the structure of the stiffness matrix which has a ®nger diagonal ®ll pattern rather than being diagonal. This is shown in Fig. 10 and is due to the multi-scale interaction of the tensor product basis in 2-D. The storage for the ®nger diagonal mass matrix grows as 4k 2 Ndof while the storage for the stiffness is proportional to 4kNdof where Ndof 2 k1 À 1 2 . In comparison, the storage for the ®nite element mass and stiffness matrices scales linearly in the number of DOF as $9Ndof. As in the 1-D case, the amplitude of the entries in the ®nger-diagonals of the multi-scale mass matrix decrease with increasing scale.

In addition to the change in the ®ll pattern, the O1 condition number for the 1-D Schauder stiffness becomes Ologh À2 in 2-D. Again, the condition number for the ®nite element stiffness grows as Oh À2 regardless of the dimensionality. Table 3 illustrates how the condition number varies with mesh resolution for the model problem in Eq. ( 1). Here, the L 2 stability of the ®nite element basis is re¯ected in the uniformly bounded condition number for the mass matrix and the Oh À2 ) growth of the condition number for the stiffness. In contrast, the condition number for the multi-scale mass matrix grows faster than Oh À2 although ®rm estimates have not been made. For the stiffness, the condition number for K k grows as logh À2 . Although the condition number for the multi-scale mass matrix degrades rapidly, the condition number for M k K k is dominated by the stiffness and is much better behaved than the ®nite element basis as shown in Table 3.

Following the 1-D mass lumping, the concept of rowcolumn lumping is applied in 2-D. However, here the lumping procedure is applied to both the multi-scale mass and the stiffness matrices. Although the row-sum lumping procedure is not valid for the ®nite element stiffness, the row-column lumping is possible for the multi-scale Fig. 9a,b. Solutions, u, on 0 x 2 obtained using the hierarchical Schauder basis with mass-lumping and a no corrections and b 3 wavelet coef®cient corrections representation of the stiffness ± an artifact of the nodal change-of-basis. In other words, the multi-scale representation of the stiffness does not retain the ``row-sum to zero'' property of the nodal stiffness matrix. In order to illustrate this property, consider the application of Algorithm 2 in two-dimensions to the model problem with e c 1 and f xy2 À x2 À y. However, in this case, the stiffness is replaced by the ``lumped'' stiffness, i.e., Ã M L k K L k . Figure 11b andc show the lumped-approximate solution with no wavelet correction and the associated error relative to the exact solution. The approximate solution was obtained by performing one vector divide and a multiscale reconstruction resulting in a very reasonable initial solution. In Fig. 11c, the peak in the error is associated with the largest approximate wavelet coef®cient. Figure 11d illustrates the lumped-approximation after the largest wavelet coef®cient has been corrected following Algorithm 2. Nearly an order of magnitude reduction in the error has been achieved by correcting only the largest wavelet co-ef®cient. In Fig. 11e, the peaks in the error are associated with the next 8 largest wavelet coef®cients. Although the row-column lumping procedure yields an approximate solution, it is thought that this approach may generate a good initial solution for an iterative procedure or alternatively provide a multi-scale preconditioner.

Although the hierarchical Schauder basis is promising, the question of computational complexity remains. In 2-D, the change in ®ll pattern and matrix conditioning make computational complexity a more important issue in the comparison of the Schauder and linear ®nite element bases. In order to evaluate the computational complexity associated with the Schauder basis, a brief comparison of the cost of solving the model problem using a conjugate gradient method was performed. In this comparison, complete ®ll-in of the multi-scale mass and stiffness was accounted for ignoring the advantages of row-column lumping.

In order to perform the comparison, both the number of non-zero entries in the matrix, and the number of iterations required to solve the problem are required. Table 4 shows the number of non-zero entries and associated iteration count for the Schauder and linear ®nite element bases for the model problem with 0 e 1000 and c 1. For problems where the mass matrix dominates, e 0, the uniform stability of the ®nite element basis in L 2 results in a computational cost that is bounded independent of the mesh resolution. In contrast, the hierarchical Schauder basis results in a lower-order computational cost for problems where the stiffness dominates due to the multi-scale preconditioning of the stiffness ± despite the additional ®ll-in associated with the ®nger-diagonal matrices. This result is further emphasized by the data in Table 5 for a ``purely'' elliptic problem with no mass matrix contribution. In summary, Table 6 provides order of magnitude estimates for the computational complexity for the hierarchical Schauder and linear ®nite element bases. Here, the number of ¯oating point operations required to solve the model problem grows as 9Ndof Ndof p for the ®nite element basis. In comparison, the number of ¯oating point operations required to solve M k K k using the Schauder basis scales as 48k 3 Ndof and as 48k 2 Ndof for the purely elliptic problem. From this data, the Schauder basis does eventually have a lower computational cost, but the bene®t is only apparent for relatively large problems. For the ``purely'' elliptic operator, K k , the hierarchical Schauder basis wins, i.e., has lower computational cost, for mesh resolution exceeding Ndof 10 5 . However, for M k K k , the Schauder basis does not win until the mesh resolution exceeds Ndof ! 10 8 .

(k) Ndof M k K k M k K k M k K k M k K k 1 9 4

Summary and conclusions

In the search for a wavelet basis with good numerical and computational performance for multi-scale simulations, the following ``shopping list'' of characteristics was developed as the goal for the ideal wavelet basis. The attributes include: (1) Low order, e.g., linear, for computational ef®ciency, (2) Consistent reproduction of polynomials, e.g., reproduce f1; x; y; xyg in two-dimensions, (3) Nodal, i.e., possesses the Kronecker delta property, (4) Hierarchical: V 1 V 0 È W 1 , (5) Element based ± compatible with isoparametric ®nite elements, (6) Analytic expressions for the basis elements / and w, (7) Easy treatment of boundary conditions, (8) Good numerical performance, e.g., dispersion characteristics, truncation error, etc., (9) Appropriate for both Eulerian and Lagrangian computations, [START_REF] Yserentant | On the multi-level splitting of ®nite element spaces[END_REF] Computationally ef®cient decomposition and reconstruction of ®elds, and (11) Extensible to multiple spatial dimensions. Based on these characteristics and the results of this exploratory effort the following conclusions are drawn. The hierarchical Schauder basis, and the multi-scale elements, are prototypical of what the ideal multi-scale basis should be although this basis does not satisfy all 11 characteristics of an ideal basis. Unfortunately, the storage and computational cost associated with the ®nger-diagonal matrices from this type of basis is a signi®cant penalty. However, the use of ad-hoc lumping procedures ameliorates this problem and offers the potential for the development of fast, simple preconditioners. Currently, the real value of the multi-scale elements lies in the application to elliptic problems. A preconditioner based on the multi-scale element would be of great value in applications where a dominant elliptic component is present, e.g., in time-dependent incompressible ¯ow and quasi-static electro-magnetics.

As demonstrated in the discussion of the hierarchical Schauder basis, it is dif®cult to construct a wavelet basis that is stable in both L 2 and in H 1 , i.e., for all possible combinations of mass and stiffness operators. The design of wavelet bases that are customized for a speci®c partial differential equation remains an open topic of active research, and the use of wavelet bases for the solution of partial differential equations remains a research topic that is currently centered in the mathematics community. However, hierarchical solution procedures that use wavelets tailored to the physical problem appear to be the most viable candidates for using wavelet bases.

Fig. 1 . 13 The

 113 Fig. 1. The hierarchical Schauder basis for scale k 2 with 8 elements and 9 nodes showing a the equivalent ®nite element basis, b the composite basis, c the basis functions for Scale-2, d the basis functions for Scale-1, and e the basis functions for Scale-0

Fig. 2 .

 2 Fig. 2. One-dimensional two-scale decomposition of a ®nite element nodal basis at scale-1

10 X e cw 1 N j dX; M 11 X e cw 1 w 1 dX M 12 X e cw 1 w

 10111121 i N j dX; M 01 X e cN i w 1 dX; M 2 dX; . . .

Algorithm 1 3 .

 13 Multi-scale solution algorithm 1. Form the coarse-grid operator K 0 , and solve the coarsegrid problem, K 0 U 0 F 0 . 2. For each element, inject one scale DOF and solve for the wavelet coef®cient, Du. fDu k g K k À1 fF k g 53 Compute the termination measure for the scale DOF injection. One possibility for the termination measure relies on stopping when the scale DOF are small relative to the overall solution. (Other stopping criteria have not been investigated.) e Du 22 kÀ1 j k ku k k 54

Fig. 4 .Fig. 5 . 6 .

 456 Fig. 4. Finger diagonal structure of multi-scale element mass matrix

Fig. 7a±f . 1 , b w 5 1 , c w 6 1 , d w 7 1 , e w 8 1 , and f w 9

 7a±f111119 Fig. 7a±f. Parent element and multi-scale DOF for k 1. a reference element and basis elements for multi-scale DOF with k 1, b w 5 1 , c w 6 1 , d w 7 1 , e w 8 1 , and f w 9 1

Fig. 10 .

 10 Fig. 10. The non-zero entries in the mass and stiffness matrices associated with the hierarchical Schauder and linear ®nite element bases in 2-D for a 32 Â 32 mesh

Fig. 11a±e .

 11a±e Fig. 11a±e. The exact and approximate solutions using the lumping procedure with and without wavelet correction for a

Table 1 .

 1 Formulae for the number of non-zeros in the 1-D mass and stiffness matrices where Ndof is the number of degrees-offreedom, Ndof 2 k1 À 1, and k indicates the scale

	Basis ± matrix	No. of non zeros
	Linear ®nite element ± Mass matrix (M k )	3Ndof)2
	Linear ®nite element ± Stiffness matrix (K k ) 3Ndof)2
	Schauder ± Mass matrix (M k )	( 2 k)1)
		(Ndof + 1) + 3
	Schauder ± Stiffness matrix (K k )	Ndof

  Calculate the approximate wavelet coef®cients, Dũ, by dividing the right-hand-side by the diagonal entries of Ã.3. Permute the components of Dũ so that the N largest coef®cients are at the top. Let Dũ N be the largest N lumped-approximate wavelet coef®cients and Dũ R be the remaining coef®cients such that

	2. permDũ	&	Dũ N Dũ R	'	:		60
							2.09	6	4	7	1	1.23
					2	7	2.71	25	18	18	1	1.27
					3	15	2.92	103	74	49	1	1.29
					4	31	2.98	414	295	121	1	1.29
					5	63	3.00	1659	1181	288	1	1.29
					6	127	3.00	6640	4724	667	1	1.29
					7	255	3.00	26560	18900	1517	1	1.29

Table 4 .

 4 Number of non-zeros and iteration count for the model problem using the a

		Mesh parameters	No. of non-zeros Iteration count
	Schauder and b linear ®nite	Scale (k)	Ndof		e 0	e 0:1	e 1	e 10 e 1000
	element bases in 2-D. The iteration count is based on using Jacobi preconditioned conjugate gradient with 0 e 1000, c 1	a Schauder basis 1 2 3 225 9 49 4 961	4 9 729 6889 51529	3 11 50 187	3 10 27 48	3 10 25 42	3 9 22 37	3 9 22 37
		5	3969	335241	488	73	61	52	51
		6	16129	1990921	979	107	87	70	59
		b Linear ®nite element basis
		1	9	4 9	3	3	3	3	3
		2	4 9	3 6 1	6	6	7	7	7
		3	225	1849	7	13	15	16	16
		4	961	8281	6	27	31	31	31
		5	3969	34969	5	54	62	63	63
		6	16129	143641	3	110	125	126	126

Table 5 .

 5 Storage requirements, matrix ®ll, and iteration count using Jacobi preconditioned conjugate for the Schauder basis with the ``purely'' elliptic model problem with e 1, c 0

	Mesh parameters		No. of	Iteration
			non-zeros	count
	Scale (k)	Ndof
	1	9	33	3
	2	49	329	9
	3	225	2265	22
	4	961	13113	37
	5	3969	68985	51
	6	16129	342265	59

Table 6 .

 6 Computational complexity for the linear ®nite element and the hierarchical Schauder basis in 2-D where Ndof = (2 k+1 )1) 2 , and k indicates the scale

	Basis	No. of	Iteration	Floating point
		non-zeros	count	operations
	Linear ®nite element	9Ndof	p	Ndof	9Ndof	p	Ndof
	(M k K k ) Schauder basis (M k K k ) 4k 2 Ndof Schauder basis (K k ) 4 kNdof	12k 12k	48k 3 Ndof 48k 2 Ndof

Sandia is a multiprogram laboratory operated by Sandia

Corporation for the United States Department of Energy under Contract DE-AC04-94AL85000.