NASA Contractor Report 178074

lCASE REPORT NO. 86-12

leASE

NASA-C((-l (~) Of"!

NASA-CR-178074

19860013879

A COMBINED INTEGRATING AND DIFFERENTIATING MATRIX FORMULATION

FOR BOUNDARY VALUE PROBLEMS ON RECTANGULAR DOMAINS

w. D. Lakin

Contract Nos. NASl-17070, NASl-18l07

March 1986

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NI\SI\

National Aeronautics and

Space Administration

Langley Research center

Hampton. Virginia 23665

U,r,GL[:Y RESEARCH GENTER

UeRARY, NASA

E:,~.:~TON, VIRGINIA

A COMBINED INTEGRATING AND DIFFERENTIATING MATRIX FORMULATION

FOR BOUNDARY VALUE PROBLEMS ON RECTANGULAR DOMAINS

W. D. Lakin

Old Dominion University

and

Institute for Computer Applications in Science and Engineering

Abstract

Integrating and differentiating matrices allow the numerical integration

and differentiation of functions whose values are known at points of a

discrete grid. Previous derivations of these matrices have been restricted to

one-dimensional grids or to rectangular grids with uniform spacing in at least

one direction. The present work develops integrating and differentiating

matrices for grids with non-uniform spacing in both directions. The use of

these matrices as operators to reformulate boundary value problems on

rectangular domains as matrix problems for a finite-dimensional solution

vector is considered. The method requires non-uniform grids which include

"near-boundary" points. An eigenvalue problem for the transverse vibrations

of a simply supported rectangular plate is solved to illustrate the method.

Research was supported by the National Aeronautics and Space Administration

under NASA Contract Nos. NASI-17070 and NASI-18107 while the author was in

residence at ICASE, NASA Langley Research Center, Hampton, VA 23665-5225.

nl @3BdOS#-

i

I° Introduction

Rotating beam configurations have traditionally been used to study the

vibrations and aeroelastic stability of rotating structures such as helicopter

rotor blades and propeller blades. More recently, models involving elastic

plates have been proposed to include the effects of spanwise variations in

material properties. The fourth-order boundary value problems associated with

both the beam and plate models do not, in general, have useful closed form

solutions. Consequently, most theoretical work on these problems has been

asymptotic or numerical in nature.

In one approach to the numerical solution of these problems, harmonic time

dependence is assumed to reduce the governing partial differential equation to

a differential equation in space variables which includes an eigenvalue. For

beam models, this is an ordinary differential equation. The fundamental

derivative which represents beam curvature may now be taken as a new dependent

variable, and the eigenvalue problem for the beam can be reformulated as an

integro-differential equation (White & Malatino, 1975; Kvaternik, White, &

Kaza, 1978; Lakin, 1982). This equation may be conveniently expressed using

integral, differential, and boundary evaluation operators. The operator

equation for the continuous solution may further be converted to a matrix

operator equation for a finite-dimensional solution vector by evaluating the

continuous equation at a finite set of discrete grid points which span the

interval of interest. A key question is now the manner in which the matrix

operators are approximated.

For beam models, one method for approximating the integral and

differential operators involves the use of integrating matrices (Vakhitov,

1966; Hunter, 1970; Lakin, 1979) and differentiating matrices (Hunter &
-2-

Jainchell, 1969; Lakin, 1985). In the simplest terms, these matrices provide,

respectively, a means of numerically integrating and differentiating a

function whose values are known at a finite set of discrete grid points. A

key property of both integrating and di[ferentiating matrices is that their

derivation requires only knowledge of the grid points, and no information is

needed about the function to be numerically integrated or differentiated. In

the case of a beam model with its single space variable, this property allows

the integrating and differentiating matrices based on one-dimensional grids to

be used directly as approximations for the integral and differential operators

in the matrix operator form of the eigenvalue problem. The result of this

approximation is a striaghtforward matrix eigenvalue problem which can be

solved by standard methods. This approach has proved capable of efficiently

handling a wide variety of beam problems including beams with concentrated

masses, follower forces, and point loadlngs (Lakin, 1982).

For vibration and buckling problems which involve two-dimensional elastic

plates, removal of the time dependence from the original boundary value

problem yields an eigenvalue problem which continues to be governed by a

partial differential equation. By analogy with the one-dimensional case, it

would seem desirable to reformulate this eigenvalue problem as a matrix

integro-partial differential equation for a finite-dimensional solution vector

on a two-dimensional grid of discrete points. Integrating and differentiating

matrices based on two-dimensional grids could then he used to approximate the

respective operators resulting, again, in a standard matrix eigenvalue

problem.

The present work will explore the potential of this approach by

considering an eigenvalue problem associated with the transverse vibration of
-3-

a simply supported rectangular plate. This problem consists of the biharmonic

eigenvalue equation

_2 ___2_2 u(x,y) = 12 u(x,y) (I.I)

(_x2 + _y2 j

for 0 J x J A, 0 J y J B, and the boundary conditions

u(0,y)= Uxx(0,y)= 0,

u(A,y) = Uxx(A,y) = 0,

(1.2)

u(x,0) = Uyy(X,0) = 0,

u(x,B)= u (x,B)= O.

YY

In section 2, equation (1.1) will be reformulated as an integro-partial

differential equation consistent with the form of the boundary conditions

(1.2). Because (1.2) involves conditions on u(x,y) itself at all four

boundaries, the present approach retains u itself as the dependent

variable. Conversion to a matrix eigenvalue problem will now require the

derivation of appropriate integrating and differentiating matrices based on a

two-dimensional rectangular grid of discrete points.

One type of integrating matrix for a function of two variables has been

previously derived by Lakin (1981). This matrix may be used in twodimensional problems whose reformulation is possible using an integrating

matrix alone, e.g., the plate analogue of a beam with cantilevered boundary
-4-

conditions. Unfortunately, this matrix is not suitable for the present

purposes as its derivation requires that the spacing of the grid points be

uniform in at least one direction. The boundary conditions (1.2) lead to a

reformulation which will require the use of differentiating matrices to

approximate partial derivatives with respect to both x and y. To preserve

accuracy, the computational grid must now include points "close" to all four

boundaries (Lakin, 1985) giving non-uniform grid spacing in both directions.

Generalized integrating matrices for functions of two variables whose

values are known on non-uniform rectangular grids are derived in the

Appendix. Differentiating matrices which approximate partial derivatives on

non-uniform rectangular grids are also derived, as are matrices which evaluate

quantities at boundary grid points. The only restriction in this derivation

is that the grid sub-units [xj ! x ! xj+ 1, Yk ! Y ! Yk+l ] should be

rectangles.

In section 3, appropriate integrating and differentiating matrices are

used to approximate the reformulated eigenvalue problem by a matrix eigenvalue

problem involving a "stacked" column vector. Numerical calculations are

presented for a rectangle with A = 2 and B = I on a 7-by-7 grid

(including near boundary points). Despite the coarseness of this grid, good

agreement with the exact eigenvalues of (1.1) and (1.2) is obtained.

2. Reformulation of the Eigenvalue Problem

Before considering the two-dimensional problem (I.i) and (1.2), it is

useful to briefly consider the reformulation process for the one-dimensional

problem of a cantilevered beam. For the beam, integrating the fourth-order
-5-

ordinary differential equation twice with respect to the space variable x

(say) gives an integral equation for the curvature w"(x) and allows the

boundary conditions at the free end to be explicitly invoked. Boundary

conditions at the free end now enter explicitly through relation of the

fundamental derivative w" to w" and the original dependent variable w.

Thus, the appropriate reformulation makes use of all boundary conditions.

The initial steps required to obtain an appropriate reformulation of the

two-dlmenslonal problem (I.I) and (1.2) are a straightforward generalization

of the one-dimensional procedure. Equation (I.I) may be integrated twice with

respect to x from x to A making use of the boundary condition Uxx(A,y)

= 0. This result may then be integrated two additional times with respect

to y from y to B making use of three other boundary conditions at x = A

and y = B, i.e.,

u(A,y) = u(x,B) = Uyy(X,B) = 0. (2.1)

These steps give the equation

B

f (_ - y)[Uxx(X,_) + (A - X)Uxxx(A,_)]d_

Y

A

+ f (n - x)[Uyy(_,y) + (B - X)Uyyy(n,B)]dn x

(2.2)

+ 2(A - x)[ux(A,y) - Ux(A,B) ] + 2(B - y)[Uy(X,B) - Uy(A,B)]

B A

+ 2(A - x)(B - Y)Uxy(A,B) + 2u(x,y) = _2 f (_ _ y) f (n - x)u(n,_)d_ d_.

y x
-6-

The corner consistency condition u(A,B) = 0 has also been explicitly used in

obtaining (2.2).

Equation (2.2) contains three different types of second partial

derivatives. A question to be answered is therefore whether, as in the single

variable beam case, the boundary conditions at x = 0 and y = 0 can be

accounted for by designating one type of second derivative as a fundamental

derivative and new dependent variable. The form of the boundary conditions is

critical in deciding this question.

The second partial derivatives Uxx and Uyy are related to the

variable u through the relations

x

/ (q - X)Uxx(_,y)dq = u(x,y) - u(0,y) - xu (0,y) (2.3) 0 x

and

Y

f (_ - (x,$)d_ = u(x,y) - u(x,0) (2.4)

0 y)Uyy - yUy(X,0).

Thus, for either Uxx or Uyy alone to be a viable candidate for the

fundamental derivative role, the boundary conditions (1.2) would have to

specify that the respective first partial derivative vanish for x or y

equal zero.

In contrast to (2.3) and (2.4) the boundary conditions in (1.2) at

x = 0 and y = 0 insure that integrated terms vanish when the mixed second

derivative u is related to u through the expression

y x

f f Uxy(q,_)dn d$ = u(x,y) - u(x,O) - u(O,y) +

u(0,0). (2.5)

0 0
-7-

This is remeniscent of the single variable case. Unfortunately, the higher

derivative conditions at x = 0 and y = 0 are not on Uxy itself but on

Uxx and Uyy. The mixed derivative can be directly related to Uxx and

Uyy through an appropriate quadrature and partial differentiation. However,

the relationship is such that neither of the two remaining higher derivative

boundary conditions can be satisfied directly. Consequently, the mixed

derivative alone is also not suitable as a new dependent variable.

Two possibilities now remain. First, the boundary value problem (i.I) and

(1.2) can be reformulated in terms of three fundamental derivatives, i.e., all

three types of second partial derivatives can be designated as separate

dependent variables. The original variable u could then be related to the

new variables through (2.5) satisfying two of the four conditions at x = 0

and y = 0. However, while attractive in principle, this option will be

impractical in practice as it leads to matrix eigenvalue problems involving

large matrices. For example, in the case of the relatively coarse 7-by-7

rectangular grid used in the present calculations for the sample problem, 147-

by-147 matrices would be required if the three derivatives are retained as

distinct dependent variables.

A second possibility is to retain the original variable u(x,y) in the

reformulation of (1.1) and (1.2). This choice requires that two additional

partial derivative operations be approximated. It also requires that equation

(2.2) be modified to explicitly take account of the four boundary conditions

at x = 0 and y = 0. However, when the resulting reformulation is

approximated by a matrix eigenvalue problem for a finite-dimensional solution

vector, considerably smaller matrices are needed than if fundamental

derivatives were used. For example, in calculations on the 7-by-7
-8-

rectangular grid, 49-by-49 matrices are sufficient for the basic development,

and the final eigenvalue problem involves a 25-by-25 matrix.

To obtain the desired reformulation, equation (2.2) is first evaluated

at x = 0, at y = 0, and at the corner point x = 0 and y = 0. The three

relations which result now explicitly use all conditions at the boundaries

(0,y), (x,O), and satisfy the corner consistency conditions at (0,0). In

particular:

A B

f n[Uyy(_,y) + (B - y)Uyyy(B,B)]d_ + f A(_ - Y)Uxxx(A,_)d _

0 y

+ 2A[Ux(A,y) - Ux(A,B)] -- 2(B - y)[Uy(0,B) - u Y (A,B)] (2.6)

B A

+ 2A(B - Y)Uxy(A,B ) = X2 f (_ - y) f _u(_,$)d_ d$,

y 0

B A

0 _[Uxx(X,$) -- (A - X)Uxxx(A,_)]d_ + fx B(B - X)Uyyy(n,B)d_

+ 2(A - x)[ux(A,0 ) - Ux(A,B)] + 2B[Uy(X,B) - Uy(A,B)] (2.7)

B A

+ 2B(A - X)Uxy(A,B) = X2 f _ _ (_ - x)u(_,_) dn d_, 0 x

B A

A_Uxxx(A,$)d_ + f B_Uyyy(_,B)dn 0 0

+ 2A[Ux(A,O) - Ux(A,B) ] + 2B[Uy(O,B) - Uy(A,B)] (2.8)

B A

+ 2ABUxy(A,B) = X2 f _ _ Bu(B,_)dn dE. 0 0
-9-

Subtracting (2.6) and (2.7) from (2.2) and adding (2.8) now gives

Ify_ + y [XUxxx(A,$) - Uxxx(X,_)]d _

0

+ xn + xf [yUyyy(n,B) - Uyyy(n,y)]dn x

+ 2X[Ux(A,0) - Ux(A,y)] + 2y[Uy(0,B) - Uy(X,B)]

(2.9)

+ 2XyUxy(A,B) + 2u(x,y)

%21 y A B x

-- x f _u(n,_)d_ d_ + yf f nu(_,_)d_ d$

x y 0

B A y x

+ xyfy fx u(n,_)dn d_ + 0f f0 _u(_,_)d_ d_[ .

This equation is the reformulation required for appoximation of the eigenvalue

problem (I.I) and (1.2) using integrating and differentiating matrices on the

rectangular grid.

3. Approximation By A Matrlx Eigenvalue Problem

The first step in approximating equation (2.9) by a matrix eigenvalue

problem for a finite-dimensional solution vector is to discretlze this

equation on a rectangular grid G of discrete points. This will allow the

integrating, differentiating, and boundary matrices on this two-dimensional

grid to be used as operators to approximate the corresponding operators in the

continuous equation.
-i0-

The two-dimensional grid G may be formed from the cross-product of

appropriate one-dimensional grids in the x- and y-directions. In

particular, let Gx be the one-dimensional grid of N discrete points

0 = x I < x 2 ... < x N = A (3.1)

which discretizes the interval 0 J x _ A, and let Gy be the one-dimensional

grid of M points

0 = Yl < Y2 < "'" < YM = B (3.2)

which discretizes the interval 0 _ y J B. Neither Gx nor Gy need to have

uniform spacing. Indeed, in actual implementation, to preserve accuracy it

will be necessary to choose the spacings x2 - Xl' XN - XN-l' Y2 - YI' and

YM - YM-I relatively small as the formulation will involve differentiating

matrices (Lakin, 1985). The two-dimensional grid G for the continuous

region of the boundary value problem (I.i) and (1.2) may now be taken as the

set of NM discrete points

G = {(xi,Yj):xiCGx,YjEGy }. (3.3)

Thus, the subunits of the grid G are rectangles which need not have equal

areas.

The straightforward format for displaying values of the solution u(x,y)

at the discrete grid points of G is an N-by-M matrix [U] with elements

Uij = u(xi,Yj). Unfortunately, this format is unsuitable for the desired

reduction to a matrix eigenvalue problem. Rather, it is convenient to arrange
-11-

the NM values of u on the grid as a "stacked" NM-by-i column vector. In

particular, the flnite-dimensional solution vector {u} is taken to have

elements

uk = u(xi,Y j) with k = N(j - I) + i (3.4)

where i = 1,...,N and j = I,...,M. It should be noted that this format

for {u} is x-oriented, i.e., the M groups of N consecutive elements in

the stacked vector give values of u for a fixed value of y in Gy while

x varies in Gx. Because of this orientation, as indicated in Appendix A,

matrices which approximate integrals and derivatives with respect to y will

require extra operations in their construction. A general flow chart for the

construction of both x- and y-operation matrices is given in Figure I.

Once equation (2.9) has been discretized, integrating, differentiating,

and boundary evaluation matrices on the grid G may be used to obtain a

matrix eigenvalue problem which provides the required approximations. In

particular, the eigenvalue problem for (2.9) may be written as

[G]{u} = %2[H]{u} (3.5)

where [G] = [GI] + [G2] + 2 [G3],

[GI] = ([JRY][Y] + [Y][JTY])([X][BA][DX3] - [DX2]), (3.6)

[G2] = ([JRX][X] + [X][JTX])([Y][BB][DY3] - [DY2]), (3.7)

and
-12-

[G31 = [XI([BA0] - [BAI)[DX] + [YI([BB0] - [BB])[DY]

(3.8)

+ [X][Y][BA][BB][DX][DY] + [I].

In (3.6) through (3.8), [JRX] and [JTX] approximate x-integrals from 0

to x and from x to A, repectively, on the two-dimensional grid G

while [DX], [DX2], and [DX3] approximate first, second, and third partial

derivatives with respect to x on G. Similarly, [JRY] and [JTY]

approximate y-integrals from 0 to y and from y to B, respectively,

on G while [DY], [DY2], and [DY3] approximate first, second, and third

partial derivatives with respect to y. The matrices [BA] and [BB]

evaluate quantities at boundary points with x = A and y = B, respectively,

while [BA0] and [BB0] give values at the boundary points (A,O) and

(O,B). The derivation of these integrating, partial differentiating, and

boundary matrices on the grid G is discussed in Appendix A. [I] is the

NM-by-NM identity matrix. The NM-by-NM matrices IX] and [Y] are

diagonal matrices such that if position k in the stacked solution vector

corresponds to the point (xi,Yj) of G, then Xkk = ×i and Ykk = Yj"

The matrix [H] in (3.4) may be written as the sum of four matrices

[H] = [H1] + [H2] + [H3] + [H4] with

[HI] = [X][JRY][JTX][Y], (3.9)

[H2] = [Y][JTY][JRX][X], (3.10)

[H3] = [X][Y][JTY][JTX], (3.11)
-13-

and

[H4] = [JRY][JRX][X][Y]. (3.12)

To enhance accuracy in calculation of the lower eigenvalues, it is convenient

to further re-write (3.5) as

[A]{u} = (.o{u} (3.13)

with

[A] = [G]-I[H] and _ = I/X2. (3.14)

To test the accuracy of this matrix eigenvalue approximation to the

continuous problem (1.1) and (1.2), equation (3.13) was solved on a rectangle

with A = 2 and B = I. G x was taken to be the seven point grid consisting

of the two end points x I = 0 and x7 = 2, two near-boundary points x2 =

0.0001 and x 6 =1.9999, and three interior points x 3 = 0.5, x 4 = 1.0, and

x 5 = 1.5. Gy was taken to be the seven point grid consisting of the two

endpoints Yl = 0 and Y7 = 1.0, two near-boundary points Y2 = 0.00005

and Y6 = 0.99995, and three interior points Y3 = 0.25, Y4 = 0.5, and

Y5 = 0.75. The implementation of (3.6) through (3.12) on this grid thus

involves 49-by-49 matrices. However, the size of the matrix [A] in the

matrix eigenvalue problem (3.13) can be reduced by noting that u = 0 at all

boundary points. Rows and columns of [A] corresponding to boundary points

may thus be deleted leading to a 25-by-25 matrix.

Exact solution of the boundary value problem (I.I) and (1.2) are of the

form

u pq(X, y) = Cpq sin --_-- p_x sin q_yB (3.15)
-14-

where p and q are integers. The corresponding eigenvalues are

Ipq = (_)2 + (_)2 . (3.16)

For the present test case with A = 2 and B = i, the exact values of the two

smallest eigenvalues are

Iii = 12.337 and 121 = 19.739. (3.17)

Approximate values for these eigenvalues were obtained by solving the matrix

eigenvalue problem (3.13) on the grid G. Differentiating matrices were based

on fourth degree polynomials while integrating marices were based on fifth

degree polynomials. This even/odd degree scheme allows grid points to be

centered as much as possible within the sliding subgrids on which the

differentiating and integrating matrices are based. The computations give the

values

111 = 12.553 and 121 = 17.635. (3.18)

It must be remarked that G in this test problem is a relatively coarse grid

which, when boundary and near-boundary points are ommitted, has only a total

of nine points in the interior of the two-dimensional region. Values of these

approximations could be improved through the inclusion of additional interior

points.
-15-

4. Concluding Remarks

The present work has examined an extension of integrating and

differentiating matrix methodology to partial differential equations involving

two space variables. Matrices which approximate integrals and derivatives on

one-dlmensional grids are used as a starting point to develop matrices which

approximate integrals and partial derivatives on two-dimenslonal rectangular

grids. The method requires that the original boundary value problem be

reformulated to take account of all boundary conditions. Integrating,

differentiating, and boundary matrices may then be used as operators to

approximate the boundary value problem by a standard matrix problem for a

stacked, finite-dimensional solution vector. The inclusion of near boundary

points in the grid helps to prevent the degradation of accuracy at boundaries

associated with differentiating matrices.

While only two-dimensional rectangular domains have been explicitly

considered in the present work, a further generalization of the method to

three-dimensional domains is relatively straightforward. This is due to the

use of a stacked column vector format for the solution vector which allows

matrices on the higher dimensional grid to be obtained from matrices on the

underlying one-dimensional grids. The primary requirement in going to three

space dimensions is the use of appropriate change matrices analogous to the

matrices [CXy] and [CYX] of Appendix A. These matrices will shuffle the

order of the stack to orient it with respect to a given variable and then

restore the original orientation after an operation with respect to that

variable has been approximated.
-16-

Appendix A: Partial Di[ferentlatlng, Integratlng, and Boundary Matrice8

For Rectangular Grids

Let G be the rectangular grid in (3.3) formed from the two onedimensional grids Gx and Gy in (3.1) and (3.2). Let {u} be the

NM-by-i "stacked" column vector defined in (3.4) which gives values of

u(x,y) at the points of G. Then, for the present boundary value problem on

a rectangle, it is necessary to derive two types of integrating matrices and

three partial differentiating matrices for each of the two space variables

x and y.

Consider first the matrices which approximate operations with respect to

x. The "0-to-x" integrating matrix [JRX] on G is an NM-by-NM matrix

such that the NM-by-I column vector [JRX]{u} contains approximate values

x

of the integral f u(B,y)dB at the points of G. In particular,

0

x2 A x2

[JRX]{u} _ (0,f u(_,Yl)dn,...,f u(n,Yl)dn,0,_ u(n,Y2)dn,... ,

0 0 0

(AI)

A x2 A

u(_,Y2)dn,''',0,f u(_,YM)d_,''', f u(_,YM)d_) r.

0 0 0

Similarly, the "x-to-A" integrating matrix [JTX] leads to approximations

A

at points of G of the integral f u(n,y)dn so that

x

A A A

[JrX]{u} _ (f u(n,Yl)d_,...,f u(_,Yl)d_,0,f u(_,Ym)dn,...,0,... ,

0 XN_ 1 0

(A2)

A A

f u(n,YM)dn,''',f u(n,YM)dn,0) r.

0 XN_ 1
-17-

The x-partial derivative matrix [DX] is an NM-by-NM matrix such that

[DX]{u} contains approximate values of _u/_x at points of G. Thus,

[DX]Iu } = I_u _x (xI'Yl)'''''_-x- _u (xN'Yl)'''''-__u (xl,YM),. "",_-_u (xN,YM)) T • (A3)

The matrices [DX2] and [DX3] which lead to approximations of second and

third partial derivatives of u(x,y) with respect to x at points of G

have similar definitions.

The stacked column vector {u} has been constructed so as to be xoriented. It thus consists of M segments which contain N elements

apiece. In each of these segments, x varies through Gx for a constant

value of y in Gy. The matrices in (AI) to (A3) approximate x-operations

for fixed values of y. Consequently, use of the present stacked column

vector format will allow x-operation matrices on the rectangular grid G to

be constructed from the corresponding N-by-N matrices on the one-dimensional

grid Gx. This construction is most easily accomplished through definition of

a "diagonalizing" mapping from the set of N-by-N matrices to the set of NMby-NM matrices.

Let p and q be integers, and let [A] be a p-by-p matrix. The

diagonalizing mapping Diag(p,q,[A]) then assigns to [A] the pq-by-pq

matrix [B] obtained by placing q matrices [A] along the diagonal of

[B] and taking all other elements of [B] to be zero, i.e.,

[B] = Diag(p,q,[A]) = [A]. • (A4)

"[A
-18-

The x-operation matrices on G may be formed by applying this mapping with

p = N and q = M to appropriate matrices on the grid G.

Consider first the construction of [JR]. Let f(x) be a function whose

values are known at the points of G, and let {f} be the N-by-I column

vector which contains these values. Further, let [jrx] be an N-by-N

integrating matrix which approximates integrals of f(x) from 0 to x on

Gx so that

[jrx]{f} ~ f f(n)dn • (A5)

0

Comparing equation (At), segment-by-segment, with (A5) now shows that

[JRX] = Diag(N,M,[jrx]). (A6)

The matrix [jrx] (and hence [JRX]) is not unique, but depends on both the

number of points included in the sliding subgrids of Gx and the manner in

which f(x) is approximated on these subgrids. The present work uses the

integrating matrix [jrx] for one-dimensional non-unlform grids developed by

Lakin (1979).

To obtain the second required integrating matrix [JTX] on the

rectangular grid G, let [jtx] be an integrating matrix on Gx such that

[jtx]{f} : f f(n)dn . (A7)

x i

Then, a segment-by-segment comparison of (A2) and (A7) shows that
-19-

[JTX] = Diag(N,M,[jtx]). (A8)

Differentiating matrices on one-dimensional non-uniform grids have been

derived by Lakin (1985). Let the matrices which approximate first, second,

and third derivatives of f(x) at points of Gx be denoted, respectively,

by [dx], [dx2], and [dx3]. As is the case with integrating matrices, these

differentiating matrices are not unique. Further, in the usual case where the

sliding subgrids contain fewer than all N points of G x, [dx2] and [dx3]

cannot be obtained by simply squaring or cubing the matrix [dx]. Rather,

these matrices must be obtained directly from approximations to the second and

third derivatives of f(x) on the sliding subgrids.

The three matrices [DX], [DX2], and [DX3] which approximate partial

derivatives with respect to x of u(x,y) at points of G may now be

constructed from [dx], [dx2], and [dx3] using the diagonalizing mapping.

In particular,

[DX] = Diag(N,M,[dx]),

[DX2] = Diag(N,M,[dx2]),

and (A9)

[DX3] = Diag(N,M,[dx3]).

For consistency, the differentiating matrices used in the present work were

based on the same subgrid approximation scheme as was used for the integrating

matrices on Gx.

Consider next derivation of the matrices which approximate operations with

respect to y on the rectangular grid G. The required integrating matrices
-20-

for this variable are the "0-to-y" integrating matrix [JRY] which

Y

approximates the integral f u(x,_)d_ at points of G, and the "y-to-B"

0

B

integrating matrix [JTY] which approximates the integral f u(x,_)d_. The

Y

product of these NM-by-NM matrices with {u} are the NM-by-1 column

vectors

Y2 Y2

[JRY]{u} _ (0,...,0,f U(Xl,_)d_,...,f U(XN,_)d_,... '

0 0

(AI0)

B B

U(Xl,_ld$,''', I U(XN,_IH$)T

o 0

and

B B

[JTY]{u} _ (f U(Xl,_)d_,...,f U(XN,_)d_,... '

0 0

(All)

B B

f U(Xl,_)d$,''', f U(XN,_)d_,0,..-,0)r.

YM- I YM- i

Let g(y) be a function whose values are known at the points of the onedimensional grid Gy, and let {g} be the M-by-I column vector which

contains these values. Further, let [jry] and [jty] be integrating

matrices on Gy such that

[jry]{g} ~ lyl f g(_)d_ iand [jty]{g} _ I_ g(_)d$ 1. (AI2)

0 Yi

Because {u}, as defined, is x-oriented, the integrating matrices [JRY]

and [JTY] on G cannot be formed using the single mapping (A4) on the

corresponding matrices for the onemdimensional grid Gy. As [JRY] and
-21-

[JTY] approximate integrals with respect to y for fixed values of x, two

additional mappings which convert from x- to y-orientation and from yto x-orientatlon will also be required.

If the values of u(x,y) at the NM points of the rectangular grid G

are written as an NM-by-I stacked column vector {v} whose format is yoriented, then the elements of {v} are

vk = u(xi,Y j) where k = M(± - i) + j. (AI3)

Thus, the vector {v_ consists of N segments which contain M elements

apiece. In each segment, y varies in Gy for a fixed value of x in G x.

If {v} and not {u} had been chosen as the format for the solution vector,

then [JRY] and [JTY] could be formed from [jry] and [jty] directly

using the diagonalizing mapping (A4) with p = M and q = N.

An x-oriented vector {u} may be associated with its corresponding

y-oriented vector {v} through a mapping Cxy from the set of NM-by-I

column vectors into itself so that Cxy({u}) = {v}. In symbolic terms, if the

values of u(x,y) on the rectangular grid G are arranged in an N-by-M

array, then the effect of applying the mapping Cxy is to produce an

M-by-N array which is the transpose of the original. For the present

purposes, the mapping Cxy may be carried out by multiplying {u} by an NMby-NM matrix [CXY] so that

[CXY]{u} = {v}. (AI4)

The matrix [CXY] may be written as a stack of N, M-by-NM matrices
-22-

I1

[cxY(I)]

[CXY] = : . (AI5)

CXY (M) ]

If [e.JJ is the j-th unit vector in N-dimensional real space, i.e., a row

vector with a one in the j-th position and zeros in the other N-I

positions, then each matrix in (A15) can be written in the form

[ejj . [ej 0 ]I

[CXY(J)] = " ] • • (AI6)

0 "[ej

The matrices in (AI5) thus have the row vector [ej] along the diagonal and

zeros elsewhere.

The mapping Cxy is one-to-one and hence invertable. Let the inverse

mapping be denoted by Cyx so that Cyx({v}) = {u}. Thus, when applied to

a y-oriented vector, Cyx restores the standard x-orientation. This mapping

may be carried out by multiplying {v} by the NM-by-NM matrix [CYX] so

that

[CYX]{v} = {u}. (AI7)

Because of the inverse relationship of Cxy and Cyx, [CYX] = [CXY]-I.

Having defined the mappings which change the format of a stacked column

vector from x- to y-orientation and back again, the procedure which uses

the M-by-M integrating matrix [jry] on Gy to produce the the vector

[JRY]{u} defined in (AI0) may now be given. The vector {u} is first
-23-

multiplied by [CXY] to produce a y-oriented format. It is then multiplied

by the matrix Diag(M,N,[jry]) to give a vector which consists of M

segments, each of which approximates the integral of u(x,y) from 0 to y

for a fixed value of x. Finally, the original x-orientation is restored

through multiplication by [CYX]. This process implies

[JRY] = [CYX] Diag[M,N,[jry])[CXY]. (A18)

The matrix [JTY] may likewise be formed in this manner from the integrating

matrix [jty] on Gy. Hence,

[JTY] = [CYX] Diag(M,N,[jty])[CXY]. (A19)

The NM-by-NM matrix [DY], which approximates partial derivatives with

respect to y on the rectangular grid G, is defined by

_u 8u _u 8u T.

[DY]{u} ~ (_ (xl,Yl),...,_ (xN,Yl),...,_ (x1,yM),...,_ (xN,YM)) (Am0)

Let [dy] be an M-by-M differentiating matrix on the one-dimensional grid

Gy such that

[dy]{g} ~ {g'}. (A21)

Then, [DY] can be formed from [dy] through the relation

[DY] = [CYX] Diag(M,N,[dy])[CXY]. (A22)
-24-

Similarly, let [dy2] and [dy3] be matrices which approximate second and

third derivatives of g(y) on Gy. Then, replacing [dy] in (A22) by

[dy2] or [dy3] leads to the matrices [DY2] and [DY3], respectively,

which approximate second and third partial derivatives of u(x,y) with

respect to y on G.

Matrices which evaluate quantities at the boundaries x = A and y = B

and at the corner points (A,0) and (0,B) are the final items needed to

construct the matrices [G] and [HI in the matrix reformulation of (I.I)

and (1.2). The NM-by-NM boundary matrix [BA] is such that [BA]{u} gives

values of u(A,y) at points of the rectangular grid G. In particular,

[BA]{u} = (u(A,Yl),...,u(A,Yl),U(A,Y2),... ,

(A23)

u(A,Y2),''.,u(A,B),...,u(A,B)) r.

[BA] may be written as a stack of M, N-by-NM matrices. If [bail is the

j-th matrix in this stack (j = I,...,M), then the element in the i-th

row (i =I,...,N) and k-th column of [baj] is unity if k = jN and zero

otherwise. Similarly, the NM-by-NM matrix [BB] is such that [BB]{u}

gives values of u(x,B) at points of G, i.e.,

[BB]{u} = (U(Xl,B),...,U(XN,B),...,U(Xl,B),...,U(XN,B)) T. (A24)

[BB] may also be written as a stack of M, N-by-NM matrices.

[BB] = [b_] . (A25)

[bb]
-25-

However, each of the M matrices in the stack (A25) is identical. The righthand block of N columns of [bb] is simply the N-by-N identity matrix.

All other elements of [bb] are zero.

The matrices [BA0] and [BB0] evaluate quantities at the corner

points (A,O) and (O,B), respectively. Both of these NM-by-NM matrices

consist of a single non-zero column which contains all ones. For [BA0],

the N-th column is non-zero. The (NM-N+I)-st column is non-zero in the

case of [BB0].
-26-

References

[i] Hunter, W. F., "Integrating matrix method for determining the natural

vibration characteristics of propeller blades," NASA TN D-6040, December

1970.

[2] Hunter, W. F. & Jainchell, R. A., "Integrating and differentiating

matrices," NASA Langley Working Paper 733, March 1969.

[3] Kvaternik, R. G., White, W.F., & Kaza, K. R. V., "Nonlinear flap-lagaxial vibrations of a rotating beam with arbitrary precone angle," Paper

78-491, AIAA/ASME 19th Structures, Structural Dymanics, and Materials

Conference, Bethesda, MD, April 1978.

[4] Lakin, W.D., "Integrating matrices for arbitrarily spaced grid points,"

NASA CR-159172, November 1979.

[5] Lakin, W. D., "Two-dimensional integrating matrices on rectangular

grids," NASA CR-165668, February 1981.

[6] Lakin, W. D., "Integrating matrix formulation for vibrations of rotating

beams including the effects of concentrated masses," NASA CR-165954, June

1982.

[7] Lakin, W. D., "Differentiating matrices for arbitrarily spaced grids,"

Int. J. Numer. Meth. Engrg., 23 (1986), pp. 209-218.
-27-

[8] Vakhitov, M. B., "Integrating matrices as a means of numerical solution

of differential equations in structural dynamics," Izv. VUZ.

Aviatsionnaya Tekn. 3 (1966), pp. 50-61.

[9] White, W. F. & Malatino, R. E., "A numerical method for determining the

natural vibration characteristics of rotating nonuniform cantilever

blades," NASA TM X-72751, October 1975.
-28-

m

One-dimensional grid

Gs:Sl<S2< ... < sK

K= N (s= x) or K= M (s= y)

FormmatricesonGs

[jrs] [jts] [ds] [ds2][ds3]

Orientation

changematrices

[CXY] [CYX]

[Diag(N,M; [])I[[CYX]Diag (M,N;E])[CXY]

x-operation

I

matrices

.

y-operation

I

matrices

[JRX][JTX] [JRY] [JTY]

[DX][DX2][DX3] [DY][DY2][DY3]

Figure I. Flowchart illustrating the formation of matrices which

approximate integrals and partial derivatives with respect to

x and y on the rectangular grid G.

1. Report No. NASA CR-178074 12. Government Accession No. 3. Recipient's Catalog No.

ICASE Report No. 86-12

4. Title and Subtitle 5. Report Date

A COMBINED INTEGRATING AND DIFFERENTIATING March 1986

MATRIX FORMULATION FOR BOUNDARY VALUE PROBLEMS ON 6. Performing Organization Code

RECTANGULAR DOMAINS

7. Author(s) 8. Performing Organization Report No. W. D. Lakin 86-12

9. Ifllgfi~~€~afio¥i~~~ ~IflS:i.cations 10. Work Unit No.

in Science

and Engineering

Mail Stop 132C, NASA Langley Research Center 11. qqA~'lC.ifIJUm; ~ASl-18107

Hampton, VA 23665-5225

13. Type of Report and Period Covered 12. Sponsoring Agency Name and Address Contractor Report

National Aeronautics and Space Administration 14. Sponsoring Agency Code

Washington, D.C. 20546 505-31-83-01

15. Supplementary Notes

Langley Technical Monitor: Submitted to J. Engineering

J. C. South Mathematics

Final Report

16. Abstract

Integrating and differentiating matrices allow the numerical integration

and differentiation of functions whose values are known at points of a

discrete grid. Previous derivations of these matrices have been restricted to

one-dimensional grids or to rectangular grids with uniform spacing in at least

one direction. The present work develops integrating and differentiating

matrices for grids with non-uniform spacing in both directions. The use of

these matrices as operators to reformulate boundary value problems on

rectangular domains as matrix problems for a finite-dimensional solution

vector is considered. The method requires non-uniform grids which include

"near-boundary" points. An eigenvalue problem for the transverse vibrations

of a simply supported rectangular plate is solved to illustrate the method.

17. Key Words (Suggested by Authors(s)) 18. Distribution Statement

integrating matrix, differentiating 64 - Numerical Analysis

matrix, boundary value problem, twodimensional

Unclassified - unlimited

19. rf~~rit~sC~arfr~a this report) 120tffie~1rAt~sCifr~H>f this page) 21. ~'O of Pages 122'lo~e

For sale by the National Technical Information Service, Springfield, Virginia 22161

NASA-Langley, 1986