When will a stiffness matrix become non-symmetric?

Hi,

A case in which the stiffness
matrix becomes non-symmetric is when the stiffness characteristic is highly
nonlocal or when the nonlocal effects become significant at a reduced scale of
study. According to the nonlocal theory, the stress at any material point is a
function of not only the strain at that point but also the strains at all
material points in the neighborhood. Theoretically, This could be a source of
non-symmetry in the stiffness matrix when you discretize the problem for a body of
finite dimensions using FEM or any other method. In
practice, however, there would be a cut-off radius which limits the
neighborhood of a point causing the stiffness matrix to be symmetric at the interior material points of the finite-dimension body. However,
even in this case, the non-symmetry is still expected at the material points close to the boundaries.

Regards,

Amir Shakouri 

Can't cite a reference off the cuff (and will have to think hard as to precisely where I read it), but as far as I know, if a variational formulation is available for a problem, or, to put it more specifically: if a potential function can be defined for a problem, then the resulting stiffness matrix is symmetrical.

It is the existence of a potential function which, I think, provides both the necessary and sufficient conditions for Maxwell's reciprocity theorem to hold. (Someone please correct me if I was wrong there.) It is Maxwell's reciprocity theorem which is usually cited as the basis for the symmetry of the stiffness matrix, in the more engineering-oriented literature---and, I don't know about the more mathematically oriented literature (nor do I care much(!!)).

Caution: Sometimes, people take "variational" and "weak" as synonymous terms, equivalent in all respects, which, I think, is erroneous. A weak formulation can always be had via application of the method of weighted residuals: by lowering the differential order and thereby weakening the demand being made on continuity (using integration by parts to lower that order). However, the possibility of a weak formulation via weighted residuals does not necessarily mean that a corresponding potential function can be defined.

As far as my own understanding goes, the symmetry is guarunteed only if a potential function exists.

A potential function cannot be defined for nonlinear problems.

Ajit

- - - - -
[E&OE]

Dear Khong:

I do get a sense of what you have been trying to get at, but no, I can't off-hand think of a demonstration that's simple enough. (Would like to know if someone can offer one.) For examples of non-symmetric stiffness matrices, search the literature on plates and shells. Hope this helps. Bye for now.

Update a few minutes later: Also see CFD (though, here, I was trying to think of examples from mechanics of solids). Bye.

--Ajit
- - - - -
[E&OE]

I'll try. One can illustrate this for linear problems (no need for nonlocal or nonlinear effects), and I will pick boundary-value problems in one dimension (extension to multi-dimensions will be apparent) for it'll be easier to understand. Have used LaTeX notation below.

One must distinguish between the continuous formulation (weak or variational) and the choices made to construct the discrete (FEM) equations.  As a mathematical model, consider 1-D diffusion (Laplace operator) with a source function f, which is governed by -u'' = f in the domain (0,1) (Lu = f). Here, L is a self-adjoint operator, which points to the associated matrix being symmetric (Hermitian in general:  see this for its many properties). In contrast, the non-selfadjoint stationary convection/advection-diffusion equation is: -u'' + ku' = 0 in the domain (0,1) (Lu = 0), where L is now non-selfadjoint. The latter (non-selfadjoint PDEs) is typical of what is realized in fluid mechanics (Navier-Stokes equations in its entire generality).

One can use the method of weighted residuals to derive the weak form in each case. In addition, for the self-adjoint (diffusion) problem, there also exists a variational formulation that is equivalent to the weak form (more on this later). For the self-adjoint problem with essential boundary conditions prescribed one has: a(u,w) = \ell(w) \forall w \in V as the weak statement, where a(u,w) = \int_0^1 u' w' dx, and u and w are the trial and test functions, respectively. For the convection-diffusion problem (non-selfadjoint problem), the bilinear form a(.,.) is defined as: a(u,w) = \int_0^1 [u'w' - ku'w]dx. On using appropriate discretizations (say FEM) for u and w, the specific discrete representation of the stiffness matrix is obtained.

With C^0 FEM, one expands u and w with the same nodal basis functions \phi_b, i.e., u(x) = \phi_b (x) u_b (sum on b) and the weak statement is `tested' by w(x) = \phi_a (x) (a = 1,2,,...N) so that we obtain N equations with N unknowns. The stiffness matrix (prior to imposing the essential BCs) is: K_ab = \int_0^1 \phi_a^' \phi_b^' dx which is clearly symmetric. This is known as the Bubnov-Galerkin approach.  The discrete equations that stem from using an ansatz in the variational formulation (minimizing the potential energy functional) for this self-adjoint problem will be identical to what is obtained with the Bubnov-Galerkin method.  However, note that if we chose different basis functions as test functions in the weak form, say w(x) = \psi_a (x) then the stiffness matrix would be K_ab = \int_0^1 \psi_a^' \phi_b^' dx, which is clearly now non-symmetric. This is the Petrov-Galerkin approach : a few meshfree methods also adopt this choice of expanding the trial and test functions with different basis functions.

Now, for the convection-diffusion problem, even with the Bubnov-Galerkin approach, one obtains K_ab = \int_0^1 [phi_a^' phi_b^' - k phi_a phi_b^'] dx which is non-symmetric. Clearly, if a Petrov-Galerkin method is used (which is the preferred choice), the stiffness matrix will also be non-symmetric. To summarize, the symmetry/non-symmetry in the FEM stiffness matrix depends, both, on the underyling weak form and the selection (linear combinantion of basis functions) of the trial and test functions in the FE approach.

The question: "... under what situation will a stiffness matrix in a FEM problem become non-symmetric?"

The answer: "...the symmetry/non-symmetry in the FEM stiffness matrix depends, both, on the underyling weak form and the selection (linear combinantion of basis functions) of the trial and test functions in the FE approach."

Ok. A few points:

(i) Thanks for a neat clarification re. u and w. For pointing out that it's not just the governing equation but also the basis functions.

(ii) Is there anything significant that gets added in saying that the differential operator be self-adjoint instead of saying that the relevant F and d variables (vectors, really) obey the required kind of a reciprocity relationship? From a physical viewpoint, the answer to the immediately last question, to my limited knowledge, is: nothing. Speaking in terms of function spaces and operators may be a mathematically more abstract way of putting it. However, if one asks for an example of a problem wherein a self-adjoint operator leads to a symmetrical matrix while also violating reciprocity relationship, will it be possible to come up with one? Of course not. Vice versa (a problem wherein the reciprocity relationship isn't violated but the operator isn't self-adjoint)? Again not... If wrong (on either count), I would like to get corrected.

(iii) As the Wiki entry for Petrov-Galerkin [^] indicates, it seems to be rather about the odd-/even-ness of the order of differential terms. Simple mindedly: Even-ordered terms decompose nicely and lead to symmetry. Potential functions involve even-orders because two first-order terms are being related together [update on Jul. 1, 2012: replaced the more stupid "multiplied" by the more intelligent "related", in this statement]. The higher order potentials can be decomposed into second-order ones (e.g. the biharmonic eqn). Availability of a potential function implies that of a corresponding functional. (Or, is it the other way around? Should the existence of a functional be the more basic idea and that of a potential function, the secondary one? Don't know about it; can't even guess!) 

(iv) Now, a question (because it, I think, leads us closer back to the original question): If the trial and test functions are chosen to be different even if the operator is self-adjoint, then what is the physical significance of such a choice? ... BTW, note, the physical significance being sought pertains to the use of the trial and test functions, i.e. to the problem situation. On the side of the solution situation, the physical significance is fairly obvious. You couldn't, e.g., expect to "interchange" voltages and currents across any two resistors in a passive network (simplest case) if the matrix weren't to be symmetric. The point is, what's the physical significance of the differing functions on the side of the problem or at the stage of conducting the procedure.

Ajit

- - - - -

[E&OE]

Dear Khong,

Defining the functional of a selected variational principle can help you to solve your problem.

Finite Element matrices are symmetric. This is due to the equation [K]=[Tt][K*][T] where [K] : elementary matrix in Global Coordinate System (C.S), [K*] : elementary matrix in Local C.S, [T]: transformation matrix (Orthogonal matrix!). For example in structural mechanics the total potential energy gives U=(1/2)*integral(sigmat*epsilon)dv = (1/2)*{qt}[K*]{q} and [K*]=integral([Bt][D][B]dv) where U: strain energy, {q}:d.o.f vector, [B]=[derivatives of N], N:interpolation functions, [D]: material properties matrix. Your Boundary Conditions have to be included by the mean of a modified formulation of the functional as Constaint Equations with Lagrange Multipliers added to your d.o.f. This can generate an extended stiffness matrix (including your constraints) which is also symmetric.

Mohammed lamine