## Domain boundary and its regularity

Until now we have worked with open bounded sets without paying special attention to their boundaries. This will change in this section, since we will need to use the unit outer normal vector to the boundary and calculate surface integrals. Let us begin with introducing the notion of a domain Definition A.52 (Domain in Rfi) A subset il C R'' is said to be a domain if it is nonempty, open and connected. Set Q C Rrf is said to be connected if every two points in Q can be connected by a continuous...

## Dyi

Theorem A.6 (Conditions for C(U, V) to be a Banach space) Let U be a normed space and V be a Banach space. Then C(U. V) is a Banach space. Proof Let L,, , be a Cauchy sequence in C(U. V) and u e U arbitrary. We show that L,, is a Cauchy sequence in V This is clear if u 0. If . 0 consider an arbitrary c > 0. There exists an index n() so that L, L < f M r for all ni.n > )(). Then Lln)i - L u - < L , L m ( < f for all (. ) > o. and indeed L,,u is a Cauchy sequence. Since V is a Banach...

## Embeddings of Sobolev spaces

Sometimes we need to decide whether all functions that belong to a Banach space U also lie in another Banach space V. Thus we are asking if the following implication holds This is equivalent to the question whether the identity operator 1 U > V is continuous. The reader already knows the answer in some situations. For example, when ft C Kd is a bounded domain, U Lp(ft) and V LV(Q). In this case the answer is positive if q < p (see Lemma A.29, Paragraph A.2.10). We also know the answer when...

## Inner Product Spaces

Some linear spaces can be endowed with inner product, which is a binary operation similar to the dot-product u.v) R u v (A.56) of vectors in R. Such spaces are called inner product spaces. Orthogonality in a general inner product space V is defined analogously to the orthogonality of vectors in R, The notion of orthogonality and orthogonal projection makes inner product spaces extremely convenient for the study of partial differential equations and finite element methods. After discussing the...

## Linear Spaces

In the first section let us refresh the knowledge of linear algebra and show its application to finite-dimensional spaces of functions. A linear space V usually is defined over a general commutative body (field) B. However, the real line R and the complex plane C are the only commutative bodies of practical importance for most applications related to partial differential equations and numerical methods. Definition A.l (Linear space) Let B R or B C. A nonempty abstract set V endowed with two...

## A U

Figure 5.2 Carle David Tolme Runge (1856-1927). Figure 5.2 Carle David Tolme Runge (1856-1927). C.D.T. Runge contributed significantly to the fields of differential geometry, interpolation, and numerical solution of algebraic and ordinary differential equations. He also was active in experimental physics, where he investigated the wavelengths of the spectral lines of elements. Nowadays his explicit second-order method (5.36) is widely used in the slightly more general form Yk+1 Yk + Atfc (l -...

## A uf

Cannot be determined since its coefficient matrix is negative definite. In such cases it is customary to multiply the equation by ( 1) so that Definition 1.1 can be applied. Moreover, notice that Definition 1.1 only applies to second-order PDEs. Later in this text we will discuss two important cases outside of this classification hyperbolic first-order systems in Section 1.5 and elliptic fourth-order problems in Chapter 6. Remark 1.2 Sometimes, linear second-order PDEs are found in a slightly...

## Ai M da Jn oxf xix x

+ ( (Mi 1 1 + Ml2u2) + (M12vi + M22u2) dS Jdn--v-' ux 1 v-v-< dx2 The transformation relations (6.77) yield Inserting this relation into (6.91), we obtain - Mn-4 + 2M12 + A 22 d (6.92) + f M + Ml s dS. - QvipdS ftp dx. Jon dv ds J rM. Jn Last, in order to prepare the weak formulation for the incorporation of the boundary data V* Qt + dM*s ds (part of the traction and soft support boundary conditions), we need to include the quantity dM s ds into the boundary integrals explicitly. This...

## Appendix A Basics Of Functional Analysis

This chapter presents elementary linear functional analysis which is needed for a first course in PDEs and modern numerical methods. Linear spaces are presented in increasing order of complexity, as shown in Figure A. 1. This text is not a traditional course in functional analysis. It assumes less at the beginning and does not address all abstract concepts of a standard functional-analytic course. On the other hand, topics needed for the study of PDEs and numerical methods, such as the Lp and...

## Appendix B Software And Examples

This chapter is devoted to the discussion of selected topics related to finite element software. In Section B. 1 we present an efficient way of connecting the packages PETSc, Trilinos and UMFPACK to a finite element solver. Section B.2 gives a brief description of the highperformance modular finite element system HERMES. The full manual posted on our web page offers additional technical details. At the end of Section B.2 we present numerical results obtained with HERMES, where the efficiency of...

## Example Diffraction problem

The last example taken from 76 is concerned with an electromagnetic diffraction problem in the domain ft ( 10,10)2 (0,10) x ( 10,0) with reentrant corner. The Maxwell's module of HERMES (see Paragraph B.2.3) is employed to discretize the time-harmonic Maxwell's equations by means of hierarchic hp edge elements. The edge elements use the same ip-FEM kernel as the elliptic module that was described in Paragraph B.2.2. The technology of the hierarchic edge elements is slightly different from the...

## Example Insulator problem

This time it is our goal to calculate the distribution of the electric field induced by an insulated conductor in the vicinity of a point where the conductor leaves the wall. The computational domain Q C R2 corresponding to this axisymmetric problem is depicted in Figure B.8. Figure B.8 Computational domain (all measures are in millimeters). Figure B.8 Computational domain (all measures are in millimeters). The wall itself, where we are not interested in the solution, is not included in the...

## Example Spherecone problem

The next problem also deals with electrostatics. A metallic sphere of the radius 200 mm carries an electric potential 100 kV. The distance of the sphere to the ground is 1000 mm. There is a metallic cone 100 mm above the sphere with zero electric potential. The cone is 500 mm high and its bottom has the radius 100 mm. The axisymmetric computational domain fi is depicted in Figure B. 14 (notice that the figure describes the boundary conditions used). We solve equation (7.25) in cylindrical...

## Interfacing with PETSc

The PETSc solver package (Portable, Extensible Toolkit for Scientific Computation) was developed at Argonne National Labs. It can be downloaded from the web page http www-unix.mcs . anl .gov petsc petsc-2 index.html where also installation instructions, documentation and an extensive amount of additional information can be found. Rather than trying to give another description of the package here, let us present the source code sMatrix_PETSc . cpp. This is a PETSc version of the sMatrix utility,...

## Modular structure of HERMES

HERMES is a modular object-oriented FEM system designed to facilitate the portability of the hp-FEM technology to various PDE models in engineering and science. The system consists of two main modules FEM ftp-FEM Module containing the finite element discretization technology, such as the mesh processing algorithms, interior mode elimination algorithms, assembling algorithms, a-posteriori error estimation algorithms, etc. Algebraic Module with a variety of solvers for systems of linear and...

## Sparse Matrix Solvers

Efficient solvers for sparse systems of linear algebraic equations are key ingredients of finite element programs. Nowadays an engineer or researcher hardly can afford developing matrix solvers on his her own, and thus public domain software packages play an increasingly important role. Moreover, as the finite element software becomes more complex, the question of efficient simultaneous interfacing to multiple matrix solver packages matters. Every matrix solver comes with its own unique...

## The elliptic module

The system of nonlinear PDEs is entered via the definition of nonzero components of the vectors and matrices in the matrix equation A (pl(x.u, VU)g )) + A (ft(x.u. VU) (B.2, d d + -7T- (P5(x. u, Vu)u(x)) + - (P(i(x, u. Vu)u(x)) Here (cc) (ui(x). U2(x) UNtq(x))1 is the unknown solution with the components tij 6 Hl(Q), i 1, Nrq. The matrix parameters P .P2 Pi of the type Ncq x N,,q may be arbitrary, some of them can be zero, but each equation in the system must stay a scalar second-order elliptic...

## The Highperformance Modular Finite Element System Hermes

Nodal elements (such as the Lagrange, Whitney, or Nedelec elements) are naturally suited for meshes where all elements have the same polynomial degree. This is why they are most suitable for problems with nice solutions. However, many problems in computational engineering and science exhibit significant local behavior in the form of steep gradients, singularities, boundary and or internal layers, etc. These phenomena can be most efficiently resolved by means of hierarchic finite element methods...

## The Maxwells module

The time-harmonic Maxwell's equation (7.62) is considered in the two-dimensional form V ( ijT'V-E) K2erE F in ft, (B.5) where E is the complex phasor of the electric field (i.e., the underlined quantity in (7.61)). In the two-dimensional setting, the relative permeability ir nT(x) is a scalar in 2D, whereas the relative permittivity er e, (x) is a 2 x 2 tensor. By to and k tu poeo we denote the frequency and wave number, respectively. Here, as in Chapter 7, the symbol V (d dx2, d dxi)T stands...

## Beam And Plate Bending Problems

In Chapters 2-4 we considered second-order PDEs whose weak formulations took place in the Sobolev space H1(Q). This space required the piecewise-polynomial finite element approximations to be globally continuous (i.e., continuous across element interfaces). Now we are going to study fourth-order problems with the weak formulations in H2(Q). Finite element approximations conforming to H2(Q) are required to be once continuously differentiable. Since the fourth-order PDEs are encountered in...

## Bending Of Elastic Beams

There are two basic one-dimensional models for the bending of elastic beams The Euler-Bernoulli model consisting of one fourth-order PDE, and the Timoshenko model based on a pair of coupled second-order equations. The Timoshenko model is simpler to solve in the sense that standard if1-conforming elements can be used for its discretization, and it is known to better capture the purely three-dimensional behavior of the structure (such as large deformations). On the other hand, the higher-order...

## BT ejA es T Rs

In this case generally one cannot access the vectors Zi. Substituting (5.95) into (5.93), one obtains yfc+l yk + Atk bi< f> (Yk + 9l, tk + ClAtk). Thus the vectors can be used instead, but at the price of s additional evaluations of the right-hand side < & . The classical Newton's iteration Let us now turn our attention to the solution of the nonlinear algebraic system (5.96) for the vectors , i 1,2, , ,s. For the sake of clarity, let us define the vector and write the system (5.96) in...

## C W Wk

Here u,n is a vertex of Kf such that x (v, ) xand the shape functions were defined in Paragraph 6.6.3. Analytical formulae for the unknown real coefficients are obtained from the conditions The symbols e,f, and g stand for the edges of K such that e Xa'(bi), f Xk c2), and g x c (63), and the points x, xf, and xg are the midpoints of the edges e, , and g, respectively. Now we combine the technique developed for the Hermite elements with the trick introduced in the previous step. For every grid...

## Classification of Sobolev spaces

We already know that all Lp-spaces are Banach spaces and, moreover, that the space L2 is a Hilbert space. Let us see about the Sobolev spaces Theorem A.22 (Wfc'p(fi) is a Banach space) Let Q c Rd be an open set, k > 1 aninteger number and p 6 1, oo . The Sobolev space is a Banach space. Proof We need to show that every Cauchy sequence i C Wk'p(Q.) has a limit Wk-p(Q.). It follows from the Cauchy property of fn Li in the Wk'p-norm that for every multi-index a < k the sequence Dufn 1 C LP(Q)...

## Continuous Elements For D Problems

After learning about the general concept of nodal finite elements in Chapter 3, the reader should know how to design general nodal finite elements of the form (K, P, E), and be able to perform the following operations check the unisolvency of the element (K, P, E), construct the unique set of nodal shape functions 8 , 02, , Qnp satisfying the delta property (3.4), use the set of degrees of freedom E and the nodal shape functions 0 , 02, , to construct the local interpolant la, construct the...

## Dc

U(x,T) is the solution corresponding to the time t T. Notice that the coefficients cne ' converge to zero very fast as the time grows, and therefore after a sufficiently long time T the solution will be very close to zero in il. Hence, the heat transfer problem evidently is a well-posed in the sense of Hadamard. Now let us reverse the time by defining a new temporal variable s T t. The backward heat transfer equation has the form We consider an initial condition uo(x) corresponding to s 0,...

## Divuvdx uVvdx u vvdS

Exercise A.61 Show that the sequence of domains i2n in Example A.51 is uniformly bounded, i.e., that there exists a bounded domain fl such that Qn C fl for all n. Exercise A.63 Show in detail that (A.89) defines an inner product in Hk (f2). Exercise A.64 Consider a function f C( 1,1 ) defined as f(x) x3 for 1 < x < 0 and f(x) x2 for 0 < x < 1. Find the largest k for which f Hk( 1,1). Exercise A.65 Consider a domain fi B(0, R) C R2, 0 < R < 1 e....

## Du o du tfOu

W + (1 dx + a'2> dx2 XlX2JhT f1'1) 1)' and decide if (and where in R2) it is elliptic, parabolic, or hyperbolic. Exercise 1.10 In R3 consider the equation OX i.i'2 OX 2.X3 and decide if (and where in R2) it is elliptic, parabolic, or hyperbolic. Exercise 1.11 In R3 consider the equation

## E E I w aoW f ij

Write the relations of the coefficients a,ij, Cj, a.Q and alJ,bl,Ci, oq. Exercise 1.4 Use Definition 1.1 to show that equation (1.8) from Example 1.1 is elliptic. Exercise 1.5 Use Definition 1.1 to show that equation (1.9) from Example 1.2 is parabolic. Exercise 1.6 Use Definition 1.1 to show that equation (1.10) from Example 1.3 is hyperbolic. Exercise 1.7 Verify that the function u(x.t ) defined in (0, tt) by the relation (1.17) is the solution of the heat-transfer equation (1.15) with the...

## E M S Vi

-4 j i> b bub2j ,bs)T, c (ci,c2, ,cs)r. The relation (5.87) represents the implicit RK method (b,c,A) defined in (5.44). The consistency condition for explicit RK methods (5.41), extends naturally to implicit RK methods via (5.88), 3 1 3 1 Moreover, we have the following lemma Lemma 5.9 The coefficients of an implicit RK method (6, c, A) defined by collocation satisfy the conditions

## Edge Elements

In the early era of finite element methods for the Maxwell's equations it was generally assumed that 7i1(n l) rf was the correct space for the discretization of the electric field E. However, the globally continuous discretizations exhibited spurious waves and other unwanted phenomena, the origin of which was not known (see, e.g., 61, 85 and 115 ). Later it was realized that the space ii(curl. ft ,) was larger than H1 (ft ,) '' The space H(curl, ft ,) admits discontinuous functions and...

## Ehpx

Where z (z , z2, , zn)'1 e CN is the vector of unknown complex coefficients. The usual finite element discretization yields a system of complex-valued linear algebraic equations of the form Here A is an N x N complex matrix and a complex vector. Only the complex version of UMFPACK can handle the complex arithmetics. However, the code is written in such a way that also real solvers can be employed. This is done by representing the complex system as a 2N x 27V real system. All matrix solvers...

## Ei o xKlx xeKi

The reader does not have to worry about the inverse reference maps in (2.64), since in the element-by-element procedure they are never evaluated explicitly. This will be explained in detail in Paragraph 2.4.9. Notice that forp 1 relation (2.64) yields the hat functions (2.25). An example of a quadratic Lagrange nodal vertex basis function is shown in Figure 2.22. The bubble functions are local to element interiors. There are p 1 bubble basis functions per every element Km G T .p, defined as (92...

## Equations For The Field Vectors

In the following we turn our attention to the original Maxwell's equations expressed in terms of the field vectors E and H. For our purposes it is not practical to cover the material properties of the involved media in their most general form (such as anisotropy, dependence on state variables, etc.). We assume that the permittivity e, permeability n and the electric conductivity 7 are scalar functions depending on the position in space only. 7.3.1 Equation for the electric field Consider...

## Equivalence Of Nodal Elements

Let us return to the one-dimensional Lagrange nodal elements for a moment again. Assume the reference domain Ka ( 1,1) and p + 1 disjoint points l yi< y2< < yp+1 1. The polynomial space has the form P Pp(Ka) and the degrees of freedom are defined as Li(g) g(yi), i 1, 2, ,p + 1 for all g e P. Consider another interval K C R connected with Ka through the affine reference map (2.37), Xfc Ka > K. We construct the Lagrange element (K,P, E) by defining new nodal points yi, y2, , Vp+i G K,yt x...

## Ev cav iM CQM gf

Let e() inf (u) v W and let J Lj be a minimizing sequence, i.e., ca Ih'n - vm 2 < a(vn vm, vn - vm) 2a(vn,vn) + 2a vm,vm) - a(vn + vm,vn + vm) where + vrn) W thanks to the convexity of W. Now E(vn), E(vm) eg implies Ih'n vm > 0 as n, m oo. Thus is a Cauchy sequence in V and there exists a limit u e V, vn > u. Since W is closed, we also have u e W. The continuity of E implies Let us show that the solution u G W is unique. Suppose that both m and u2 are solutions. Clearly the sequence U2,...

## Exercises

Exercise 5.1 Use Definition 5.5 to prove that a solution of (5.11) that is asymptotically stable in the forward direction, is necessarily unstable in the backward direction. Exercise 5.2 Consider the ODE (5.11) with the right-hand side (5.10). Extend the procedure of construction of the linear system (5.48) to the general (b, c, A) RK method. Exercise 5.3 Use Theorem 5.3 to prove that the stability function R(z) 1 + zbT(I z_4)-1l of every explicit higher-order RK method (b,A) is polynomial....

## F

Fi l(vi) vl x)dx -hi + -hi+i -(hi + hi+l). (2.29) Hence the system of linear algebraic equations (2.27) has a tridiagonal form illustrated in Figure 2.3. Figure 2.3 Tridiagonal stiffness matrix Sn for piecewise-affine approximations in ID. It follows from Corollary 2.1 that the stiffness matrix Sn is invertible, and thus there exists a unique vector Yn containing the coefficients of the approximate solution un Vn. This model situation is particularly interesting, since the linear algebraic...

## Firstorder Hyperbolic Problems

This section is devoted to first-order hyperbolic problems of the form These equations differ from the previously studied second-order PDEs significantly and methods other than FEM are usually used for their numerical solution. PDEs of the form (l .103) are referred to as conservation laws, and they play an important role in the continuum mechanics and fluid dynamics. The (generally nonlinear) flux function ( 1 2, , fd)T, where d is the spatial dimension, consists of d directional fluxes fi Rm...

## Ftfc nffc wnvk wvi

And thus u ,f e Since , V _i, it is j i,f 0, and we can define For example, the Legendre polynomials can be constructed via the Gram-Schmidt procedure EXAMPLE A.44 (Legendre polynomials) Let us use the monomial basis of the space V L2( 1,1 ),Bmnn u> i, (> 2, '3, '4, 1, x, x2, x3 and the Gram-Schmidt process to create an orthonormal basis of the space V. In the first step normalize wj, Ln(x) ('1 u'i u'i l 2. and define Vi span i> i . Next define the element 6 Vi by 0. Thus u< 2 w2 u'2 x....

## G t t a ai ao ai a

Similarly the condition g (0,1)T const, on the edge e3 yields 62 0, and the last condition g ( 1, l)7 const, on the edge e2 means that a2 61. Hence any polynomial g P has the form Leu0(g) Le2fi(g) Le3io g) 0 g 0 for all g P. (7.96) is a basis in P. Now any g P can be expressed uniquely as and the left-hand side of the implication (7.96) can be written as aiLei gi) + ot2Leu 0(32) + a3ieu o(g3) 0, aile2,o(ffi) + a2Le2t0(g2) + a3Le2fi(g3) 0, ai e3,o(5i) + a2Le3fi(g2) + a3Le3fi...

## Hermite Elements In D

In contrast to the one-dimensional case, Hermite elements do not conform to the Sobolev space H2 in 2D and 3D. Nevertheless, they find important application in nonconforming approximations to fourth-order problems, computational geometry, surface reconstruction, and elsewhere. We focus on triangular elements, since the construction of Hermite quadrilaterals can be done easily using the product geometry of the reference domain Kq. The cubic Hermite element (Kt, Ps(Kt), E) on the triangular...

## Higherorder Elements

In Section 2.2 we constructed a Galerkin sequence Vi C V2 in the Sobolev space V by subdividing selected mesh elements into subelements of the same polynomial degree ( i-refinement). Sometimes, much faster convergence can be achieved by increasing the polynomial degree of the elements instead (p-refinement). Such approach usually is more efficient in elements where the solution is very smooth, without singularities, oscillations, or boundary internal layers. An illustrative example is given in...

## Higherorder Hermite Elements In D

This section is devoted to the design and properties of higher-order Hermite elements. The cubic Hermite elements are extended to arbitrarily high polynomial degrees in both the nodal and hierarchic fashion in Paragraphs 6.3.1 and 6.3.2. The conditioning properties of the nodal and hierarchic higher-order shape functions are compared in Paragraph 6.3.3. The basis of the space VhtP is constructed in Paragraph 6.3.4 and the integrals from the weak formulation (6.21) are transformed to the...

## Higherorder Irk Methods

It was discovered in early 1970s that general (implicit) RK methods could be generated by combining classical collocation methods with higher-order numerical quadrature rules. Several RK methods derived via the traditional Taylor expansion techniques turned out to actually be collocation methods. A classical book on higher-order IRK methods is 60 . Let us return to the (nonautonomous) ODE system (5.11), (5.12) resulting from the MOL, The collocation constructs the approximate solution X(t) Y(t)...

## Higherorder Nodal Elements

In this section we extend the lowest-order Q1- and P1 -elements to higher-order Lagrange elements. The quadrilateral case, based on the product Gauss-Lobatto points, is described in Paragraphs 4.3.1 and 4.3.2. The quality of the Lagrange interpolation is discussed in Paragraph 4.3.3. Higher-order triangular elements are constructed using the Fekete points in Paragraphs 4.3.4 and 4.3.5. The basis of the space VftiP for regular hybrid quadrilateral triangular meshes is presented in Paragraph...

## Hpxj hPxj xj j M

Where gh,p Km 6 Pl Km) for all Krn 6 7)liP, as illustrated in Figure 2.30. interpolation on piecewise-affine elements. interpolation on piecewise-affine elements. Higher-order case On a general higher-order finite element mesh 'Th r> , as the reader may guess, the interpolation problem is decoupled by subtracting the piecewise-affine vertex interpolant gvh p from the interpolated function g. The function g g vanishes at all grid points, and can be projected locally onto the polynomial spaces...

## Dgs

The formula (A.34) is still used in modern digital calculators due to its high efficiency. Its origin sometimes is attributed to Heron of Alexandria who described it in his famous Metric, but as a matter of fact the formula was already known to the old Babylonians. The concept of contractive operators and the Banach fixed point theorem play an important role in the analysis and numerical solution of nonlinear problems. Let us mention the basic results and show their applications Definition A.39...

## J dT J

E q dx LbJ 9K(E)) for all j 1,2, , k- l)k. 7.5.4 Transformation of weak forms to the reference domain In this paragraph we transform the integrals involved in the weak formulation (7.77) of the model problem to the reference domain, as required by the element-by-element assembling procedure. We focus on triangular elements, but we will point out where the quadrilateral case differs. Recall from Paragraph 7.4.3 the weak formulation Find E V such that a(E, E) 1(F) for all F e V, (7.129) where the...

## IKIlQ up N

On the other hand, gives the maximum difference of u and We use the Hl-norm in what follows. The problem was solved twice, using the piecewise-affine FEM and the hp-FEM. In both cases it was our goal to attain the best possible accuracy using as few degrees of freedom as possible. The approximate solution, its gradient, finite element meshes, and a-posteriori error estimate e ,.p uref ui,.p based on a very accurate reference solution ,.,. , are shown in Figures B.3-B.7. The efficiency of the...

## Imi cuy

Since X is linear, it is continuous if and only if it is bounded (see Lemma A.24). The following definition generalizes the Lipschitz continuity of functions and introduces the Holder spaces Definition A.59 (Holder continuity, Holder space Ck'i3(fl)) We say that a function f C(fl) is Holder-continuous with the exponent 3 > 0 if there exists a constant Cf such that ( Ci) - f(x 2) < Cf xi -x2f for all x ,x2 e ft. The space Ck'0 Q) consists of functions whose ath partial...

## Index

Abelian group. 172 algorithm adaptive quadrature. 63 RK5(4) method, 183 assembling P1 dements jn 2D, 136 PP Qr-elements in 2D. 163 element-by-element in ID, 56 Hermite elements in ID, 231 vertex-by-vertex in ID, 56 enumeration of bubble DOF in 2D, 162 DOF for Hermite elements in 1D, 230 DOF for Lagrange elements in ID, 78 edge DOF in 2D. 161 edges in 2D, 159 vertex DOF in 2D, 134 area moment of inertia, 212 array Butcher's, 180 connectivity ID, 78 2D, 134 CSR, 84 autonomization of RK methods,...

## Info

1 + + z2 + + , exp(c) R z) + 0(z2 We see that the function R(z) is a consistent approximation of exp(s). The stability domain of the function is SR C ze C R(z) < 1 C z e C 1 - z < 1 , which is the complement of the closed complex circle with the center at 1 + Oi and radius 1. In particular, and therefore the implicit Euler method is stable for all values of A 6 C and all time steps At (we say that it is absolutely stable). 5.3.4 Stability functions for general RK methods The stability...

## What Are Lobatto Shape Functions

Figure 2.23 Performance of an iterative matrix solver on two differently conditioned matrices of the same size and sparsity structure. The horizontal axis represents the number of iterations and the vertical one shows the norm of the residuum of the approximate solution. Definition 2.2 (Condition number) Let M be a nonsingular n x n matrix. The product where . j is some matrix norm, is called condition number of the matrix M ( with respect to the norm . j. One may use, for example, the standard...

## Interpolation On Finite Elements

Assume some restricted set of functions C (such as, for example, polynomial, piecewise-polynomial or trigonometric polynomial functions) in a linear space V and a function g 6 V that does not belong to C. The prototype approximation problem is to find a suitable function gc C (approximation of g) such that gc is in some sense close to g. The measure of the quality of the approximation (abstract distance of gc from 3), can be defined as an error estimate, the norm g gc v if the space V is...

## Interpolation On Nodal Elements

The interpolation on finite elements is a procedure that takes a function g V(Oh) and produces its suitable piecewise-polynomial representant in the finite element space gh,p 6 Vh,p Qh)- Here by Qh we mean a suitable open polygonal domain that approximates the domain fl for the purposes of the finite element discretization (more about this will be said in Chapter 4). Various interpolation techniques with different quality and computational cost can be used, ranging from the fastest fully...

## J A

Says that the total dielectric flux P out of any (simply-connected) volume V with a sufficiently regular boundary S is equal to the total electric charge Q contained in the volume V. The total electric charge Q is defined by where q is the electric charge density. The symbol u x) stands for the outer normal vector to the surface S at a point x S. Finally, Gauss' law for magnetism, postulates that the magnetic flux out of any volume V with a boundary S is zero, or, in other words, that the...

## J nJ n

And the linear form I e V' is defined by 2.2.2 Finite-dimensional subspace Vn C V The Galerkin procedure assumes a sequence of finite-dimensional subspaces of the infinite-dimensional Hilbert space V, satisfying (2.18), Let n > 1 be a natural number. Consider a partition of the interval fl (a, b), and define the finite element mesh are called finite elements, and the value h(n) max (x_.n' ) 1 < i< M 1 l U Remark 2.3 For historical reasons the subscript h h(n) is often used instead of the...

## J o ds

Suppose, moreover, that the given twisting moment M*s(s) has jumps at the corner points s sj, where 0 < Sj < b, j 1, 2, , p (p < Nc, b < , s (0, b) on rir). Then fb fj * p - I - d,s _ (Sj + 0) - AC (Sj - 0msj ) + M*sip b Writing M*s(sj + 0) M*s(sj 0) hj, we see that the corner discontinuities in Mus and A *s. produce the following terms, (A g(si + 0) - Mva(Si - 0))i> (st) - hM*i)> which are not present when the boundary dfl is smooth. Taking this fact into account, from (6.92) we...

## JcJ oJo

The difference of the electric potentials at points A and D is called voltage and denoted by uab- If the loop C is closed, the following holds (fields with this property are called conservative). Point sets with the same potential (curves in 2D and surfaces in 3D) are called equipoten-tials. In 2D an equipotential curve Ccl2 starting from a point A g R2 can be constructed easily (numerically) using the relation 9,.(C(s)) const V ,.(C(.s)) C' s) 0 o E C .s)) C'(.s) 0 (7.23) i.e., C is...

## Jdn os Jm dsJTtr ds

However, in reality the boundary dQ often has corners. Suppose that the boundary dQ is defined parameterically by means of the parameter s (0,I), and that there are Nc corner points 0 < Si < I, i 1, 2, , Nc. In general the function Mvs(s) has jumps at these points since the normal vector v varies discontinuously. The integration by parts with a piecewise continuous function Mus(s) gives (we assume Af s(0) Mvs(l)y.

## L

P I Q E)-oK,jtjdZ I E-ij , j 1,2,3, (7.109) where oatJtj is the unit tangential vector to the edge a,j of K oriented compatibly with tj through the map xk- The tangential vectors t3 are transformed by xk according to the relation Let the matrix T of the type 2x2 represent the transformation . Then, after substituting for OKjtj from (7.110), relation (7.109) becomes f J E-ijdt, J 1,2,3, (7.111) To satisfy equation (7.111), the matrix T has to have the form Thus finally the correct transformation...

## L fdxl qmds l mdsghmai

The boundary integral over r.,s only contains A *, because ij> 0 on r.SiS, and no boundary integral over Tr is present since both ip and dip dv are zero on Tci. The Dirichlet lift G(x) is implemented into the left-hand side in the usual way, by decomposing the moments A ,_, into M,j(w) MtJ(W + G) Mjj(W) + M,j(G). ij 1.2, and leading all integrals containing MLj(G) over to the right-hand side. The weak formulation of equation (6.90) reads Find a function W e V(S1) such that a(W. ip) l t) for...

## List Of Figures

1.1 Jacques Salomon Hadamard (1865-1963). 6 1.2 Isolines of the solution u(x, t) of Burger's equation. 7 1.3 Johann Peter Gustav Lejeune Dirichlet (1805-1859). 14 1.4 Maximum principle for the Poisson equation in 2D. 27 1.5 Georg Friedrich Bernhard Riemann (1826-1866). 41 1.6 Propagation of discontinuity in the solution of the Riemann problem. 42 1.7 Formation of shock in the solution u(x, t) of Burger's equation. 44 2.1 Boris Grigorievich Galerkin (1871-1945). 46 2.2 Example of a basis...

## Lowestorder Elements

Let Q C Rd, where d is the spatial dimension, be an open bounded set. If the Hilbert space V consists of functions defined in 0 and the Galerkin subspaces Vn C V comprise piecewise-polynomial functions, the Galerkin method is called the Finite element method (FEM). Let us begin with the exposition of the simplest case of piecewise-affine elements in one spatial dimension. V (aj Vw) + agu (ail )' + aou , (2.20) where e L2(ii), in a bounded interval fi (a, 6) C R, equipped with the homogeneous...

## Lowestorder Hermite Elements In D

For the first exposition of the Hermite elements let us consider problem (6.7) with the clamped boundary conditions (6.11 ) from Paragraph 6.1.4. The weak formulation for this case was derived in (6.14). Consider a subdivision a xq < x < < xm b of the domain Q. For i 1, 2, , M denote Ki (xj_x, Xi). Recall from Paragraph A.4.2 that if1-functions are continuous in one spatial dimension, Applying (6.18) to the derivative of w v' 6 one obtains Any approximate solution to problem (6.14) has to...

## Lv hP

And inserting this linear combination into (6.51), one obtains a system of pt 3 linear algebraic equations, A P fflP E A fc)(x) dx 0, k 4, 5, ,Pi, (6.52) for the unknown coefficients a . Transformed to the reference domain Ka, with the orthogonal hierarchic shape functions (6.36) this linear system simplifies substantially to Aujm(OAujk(Od A(s-sj .p)( )AWfc(Od ) (6.53) 4 A(5- ,p)(0Awfc(0d , (6.54) where k 4, 5, ,pi. Hence no system of linear algebraic equations needs to be solved. Here g(xK (Q)...

## M

Some remarks are in order before we introduce concrete data structures and algorithms. Generally, data structures differ from implementation to implementation. A safe way to avoid criticism for the complicatedness or inoptimality of one's data structures and algorithms is not to expose them. On the other hand, the presentation of the data structures and algorithms may be of considerable help to beginners. Therefore let us try to be concrete, without claiming that our data structures or...

## Method Of Lines

In this section the reader will need some basic facts about second-order parabolic problems from Chapter l. Recall the general second-order parabolic equation (1.76), where Lisa second-order elliptic operator of the form (1.1), fi c r2 is a bounded domain with Lipschitz-continuous boundary, T > 0, Qt x (0, T) is the corresponding spacetime cylinder, and 6 C(Qr)- The classical regularity assumptions are weakened after the problem is stated in the weak sense. For example, in the special case...

## Mg f x dx

We look for coefficients Pi, 02,03 such that Applying A to the basis B, by (A.8) and the delta property fi(vj) 6 j, we have Calculating the integrals of the basis functions 2 and v3, we obtain Pi 1 3, 02 A 3,03 1 3, and thus Thus by (A.9) we can integrate all quadratic polynomials using their function values at -1,0 and 1, A(g) h a) + h g) + h g) + (o) + this is the Simpson's rule for the interval ( 1,1) . Exercise A.l Consider the set _Mnxn of all real n x n matrices. 1. Show that _Mnxn is a...

## N

Next let us calculate the angle of arbitrary geometrical sequences u and i> n f 1 in the space I2 equipped with the 2-product (A.67). Consider geometrical sequences given by the parameters u uor vn vqs1, 0 < r. ,s < 1, 0 < uQ, v0 R. We have 4. Last consider an open set Q (-1,1) c R and the Hilbert space L2(Q). The angle of the functions f(x) 1 and g(x) x is (u.v) ( I ,.r d.r a arccos r arccos -

## O V u

Where the conductivity 7 is a function of the spatial variable. Boundary conditions The boundary Oil can be split into two open (not necessarily connected) disjoint parts Fp and T j. We consider the perfect conductor boundary condition (7.57), Et 0 onFp, and the impedance boundary condition (7.71), Here t t(x) is the positively-oriented unit tangential vector t ( 1*2, vi)1, where v (vi,V2)T is the unit outer normal vector to the boundary < 9fi. The impedance A A(a ) > 0 was defined in...

## Oc

P(v) Wi)vw, for all u& V (A.77) defines a unique orthogonal projection operator P e C(V, V), P2 P, (v - Pv, w) 0 for all v V and w e P(V). Moreover, P 1 and V P(V) (I - P)(V) is an orthogonal direct sum. The assumption of closedness of the subspace Vj is essential and we will discuss it in more detail in Remark A.5 and Example A.48. On the contrary, the assumption of the existence of an orthonormal basis By, is not necessary, since one always can have such basis by Theorem A. 12....

## Og

Where M is the matrix norm (A. 15), Let u, v V. Substituting D l(L + U) for M and u v for w in (A.39), we obtain IIFw-F IU IILrHL + t Xii-tOlloc < IID-'iL + UMu- vlU. Hence the operator F is a contraction if D 1(L + U) < 1. Looking at the structure of the matrix D_1(L + ( ), we have 1 n 1 n D l L + U) max V ly- + uVJ max , - V a0- . l< < n dn 1 < i< n au . ' Thus a sufficient condition for the operator F to be contraction is a > ciij for all 1 < < n. (A.40) Every matrix A with...

## On

Corollary 1.2 (Comparison principle) Let fi C Rd be an open bounded set and L an elliptic operator of the form (1.71). Suppose that functions u,v G C2 (fi) O C(Q) solve the equations Lu fu and Lv fv, respectively, and fu < fv in n, u < v on < 9fi. Proof Apply the minimum principle to w v u. Corollary 1.3 (Continuous dependence on boundary data) Let Q. c Ud be an open bounded set and L an elliptic operator of the form (1.71). Suppose that u and U2 solve the...

## Oos G

Figure B.21 Approximate solution of the micromotor problem. Top electric potential ph.p (zoom 1 and 6). Bottom left detailed view of the singularity of -E ,.(> I V jat a corner of the electrode (zoom 1000). Bottom right Error estimate based on a reference solution (zoom 1000). Figure B.22 The hp-mesh (zoom 1,6, 50, 1000). Large sixth-order elements are used far from the electrodes and small quadratic elements are placed at the reentrant corners. The piecewise-affine mesh had geometry...

## Ox Ox

Some of the above-defined quantities are depicted in Figure 6.23. Figure 6.23 The transversal force, shear resultant, and bending and twisting moments. Further details on the derivation of the thick plate model can be found, e.g., in 124 , In the next paragraph, equations (6.68), (6.69) will be used to deduce the Kirchhoff thin plate model. In addition to the hypotheses (P1)-(P4) of the Reissner-Mindlin plate model, the Kirchhoff model imposes the normal (Love's, Kirchhoff's) hypothesis, which...

## P

Using the value co in (A.51 ), we obtain ry. Using the value co in (A.51 ), we obtain Let us remark that in the space lp(R), equipped with the discrete p-norm (A. 12) or the discrete maximum norm (A. 10), the Holder inequality (A.50) attains the following form Let u,v Rn, and 1 < p, q < oc such that l p + l q 1. Then This inequality sometimes is called discrete Holder inequality. The Minkowski inequality carries the name of a German mathematician Hermann Minkowski. Although he was...

## Partial Differential Equations and the Finite Element Method

The University of Texas at El Paso Academy of Sciences of the Czech Republic A JOHN WILEY & SONS, INC., PUBLICATION Copyright 2006 by John Wiley & Sons, Inc. All rights reserved. Published by John Wiley & Sons, Inc., Hoboken, New Jersey. Published simultaneously in Canada. No part of this publication may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, recording, scanning, or otherwise, except as permitted...

## Preface

Rien ne sert de courir, il faut partir point. Jean de la Fontaine Many physical processes in nature, whose correct understanding, prediction, and control are important to people, are described by equations that involve physical quantities together with their spatial and temporal rates of change (partial derivatives). Among such processes are the weather, flow of liquids, deformation of solid bodies, heat transfer, chemical reactions, electromagnetics, and many others. Equations involving...

## Qjt

V x ( T1 Vx ) -Vx . (7.43) If the partial derivatives of H are continuous, they can be interchanged and we obtain V x ( i1 V x E) - (V x H). (7.44) Substituting for V x if from Ampere's law (7.3) and using the constitutive relations (7.7) and (7.9), we can write dE d 2E 0Ja V x ( r1 V x ) +7 - . (7.46) In practice the third term on the left-hand side sometimes is neglected when dealing with lower frequencies (typically less than 1 MHz). 7.3.2 Equation for the magnetic field Taking the curl of...

## Secondorder Parabolic Problems

Next let us turn our attention to linear parabolic problems (the notion of parabolicity was introduced in Definition l.l). Let Q. C Rd be an open set with Lipschitz-continuous boundary. We will study a class of linear parabolic equations where t is the time, u u(x, t), f f(x, t) and L is an elliptic operator of the form (l.l) with time-independent coefficients. The equation (l .76) is considered in a space-time cylinder QT ft x (0, T), where T > 0. 1.3.1 Initial and boundary conditions...

## Selected General Properties

Second-order PDEs (or PDE systems) encountered in physics usually are either elliptic, parabolic, or hyperbolic. Elliptic equations describe a special state of a physical system, which is characterized by the minimum of certain quantity (often energy). Parabolic problems in most cases describe the evolutionary process that leads to a steady state described by an elliptic equation. Hyperbolic equations describe the transport of some physical quantities or information, such as waves. Other types...

## Selected Time Integration Schemes

There exist many excellent papers and books on the numerical solution of ODEs, and numerous sophisticated ODE packages can be downloaded from the Internet. However, one should not think that all important problems in the theory and numerics of ODEs have been solved. On the contrary Significant progress has been made recently in the development of new methods and in understanding of the existing ones, and the numerical solution of ODEs continues being a very active research area. The...

## The source code of sMatrix Trilinos cpp

implementation of Trilinos solvers under the sMatrix interface include sMatrix.h include < Epetra_ConfigDefs.h> include < Epetra_SerialComm.h> include < Epetra_CrsMatrix.h> include < Epetra_Map.h> include < Epetra_Vector. h> include < Epetra_LinearProblem.h> include < AztecOO.h> include < sstream> include < fstream> include < iomanip> include < cstdio> include < cstdlib> fprintf(stderr, s n, msg) fflush(stderr) exit(0) sMatrix sMatrix(int iRank,...

## Ti

The proof is finished by realizing that u - ur N(f). EXAMPLE A. 10 (Linear operators) 1. Consider the spaces V R3 and W R3, along with a map V > W defined by For every u,v e V it is f(u + v) 2(u + v) 2u + 2v f(u) + f(v) and moreover f(au) 2au a2u af(u) for every u G V and a R. Therefore is a linear operator. It is easy to see that N(f) ((), 0, 0)T . Moreover, since for every vector iv W we can define a vector v w 2 G V such that f(v) w, it is R(f) W- Therefore is a bijection. 2. Next consider...

## Tl T i

It is left to the reader as an easy exercise to verify that Accordingly the set of degrees of freedom contains three linear forms, E ei, o, e2, o, e3,o . where P R is defined as the integral of the tangential component of the field E on the edge er Lemma 7.7 (Unisolvency) The finite element (K,. P, E) is unisolvent. Proof According to Definition 3.2 we have to show that the following implication holds Let us find some basis in the space P first. A general polynomial g 6 P1(Af) 2 has the 1,6)...

## Transient Problems And Ode Solvers

The nature changes at every time instant, and the numerical simulation of evolutionary processes plays an important role in applied sciences and engineering. Transient problems can be very complicated and the spectrum of numerical methods for their solution is accordingly wide. At the introductory level it is natural to begin with the Method of lines (MOL), which has a prominent position due to its ability to add temporal evolution to all numerical methods for stationary PDEs without altering...

## Ux xXlX ox

U3(x) w(xi,x2), and (6.66) yields the momentum components > (1 _ 7) A f d2w + dJL ( + d2yj) Substituting the shear resultants Qi,Q2 from (6.73) into (6.68), we finally obtain the Kirchhoff thin plate model in its well-known form, In the following we introduce several types of boundary conditions for equation (6.74), derive its weak formulation, and prove the existence and uniqueness of the weak solution. The boundary conditions for plates are typically prescribed in local coordinates, where v...

## J

In what follows, by the symbol I we denote the n x n identity matrix, I diag(l, 1, , 1) (the dimension n will always be clear from the context). Definition A.14 (Nonsingular, singular, and inverse matrix) Let M be a real or complex n x n matrix. Then M is said to be nonsingular if its n columns are linearly independent vectors. Otherwise AI is singular. The matrix M l such that MM l M 1M I is said to be the inverse of M. In Definition A. 14 one can use linearly independent rows instead of...

## Vector formulation With the notation

The problem (B.2) can be rewritten to div (_4.Vit) + (P5u) + (P6u) + P7u F in Q, Ox oy ( .4Vu v) + (P uvi + P6w2)i 9n.i on IV*. Spatial discretization The polynomial degrees 1 < p < 10 on the finite elements can be defined either via a data file or by means of a function in the code that assigns the polynomial degrees to elements based on the coordinates of their vertices. The polynomial degrees may differ from element to element. The solver constructs the corresponding hierarchic basis of...

## Veru ip

The (p l)(r 1) element-interior nodes (bubble nodes) can be sorted in any unique way, for example as With this point set in hand, the Lagrange Q '-element is constructed as follows 4.3.2 Lagrange-Gauss-Lobatto Qp'r-elements It is natural to construct the master Q -element on the reference domain Kq first, and then to extend it to an arbitrary convex quadrilateral domain K. qp< r-element on the reference domain Kq In the sense of Definition 3.1 the master element is a triad (Kq, Qp'r(Kq), g),...

## W V

Also in the space L2(ft) the Cauchy-Schwarz inequality, u(x)v(x)dx < w(x)v(a ) dx Jn Jn < jju(x) 2dxj JJv(x) 2dxJ H 2 H 2, is a consequence of the Holder inequality (A.50) with p q 2 and the triangular inequality. We know from Lemma A.32 that every inner product induces a norm. Conversely, there are norms which induce an inner product Lemma A.33 (Parallelogram rule) Let V be a real normed space. If the norm satisfies the parallelogram rule + t-ll2 + m - r 2 2...

## WkpQ Wkptt A Traces of Wfcpfunctions

Since the Sobolev space Wk'p(il) always is a subset of the corresponding space Lp(ft), the W 'fc p-functions are only defined almost everywhere in ft. Since the boundary dQ is a zero-measure subset of ft, it might seem that the boundary values (traces) of W -functions never can be well defined. However, the notion of trace is associated with the whole class of VKfc'p-equivalent functions, and it is defined using a representant that is continuous up to the...

## Wx ux

This problem belongs to the class of first-order hyperbolic conservation laws that will be studied in Section 1.5. There the reader will learn how to derive the exact solution to (1.98), (1.99) in the form uo(x ct) + uo(x + ct)--UAx ct) --U (x + ct) where U x) is a primitive function to U (x) Exercise 1.19 Can equation (1.89), when equipped with a Neumann boundary condition on the whole boundary dfl, have a unique solution How would this change in the stationary case Lu f Exercise 1.20...

## On rd

The solution w(x) is sought in the form where the unknown function W e V(ft) satisfies the homogeneous boundary conditions It is advantageous to develop the weak formulation from the equation d2Mll + 2 d2Ml2 + d2M22 (69Q) which is equivalent to (6.74) through (6.68) and (6.76). We multiply (6.90) with a test function i> (x) V(fl), whose minimum regularity will be determined later from the final weak forms, and integrate over ft, Jn dxi V dxi dx2 J Ox2 dxi dx2 J Jn

## [fvi j lvdi [ftp J lpdx

J Dlfpdx - f(0 )< p(0 ) + (0+M0+) Thus (A.83) holds and the above-defined function D uf is the weak derivative of . EXAMPLE A.54 (Weak differentiability in one dimension II) Discontinuous functions in ID are not weakly differentiable Consider the function By the same token as in Example A.53, the only candidate for the weak derivative DlJ is the zero function (with an arbitrary value at x 0). If zero is the weak derivative of , for all ip e Cjf3 ( 1,1) we have Dlfip dx - j ftp' dx - j...

## Distributions and weak derivatives

The following compact notation is practical for operations with partial derivatives Definition A.53 (Multi-index) Let d be the spatial dimension. Multi-index is a vector (ai, 2, ,Qd) consisting of d nonnegative integers. By ct Q> we denote the length of the multi-index a. Let f be an m-times continuously differentiable function. We define the ath partial derivative of f by Note that D f for a (0, 0 0). To give at least two other examples, we have for a (1,0, 0 0), and for a (1,1 1) one...

## Interfacing with Umfpack

UMFPACK is a set of routines for solving nonsymmetric sparse linear systems by means of the Unsymmetric MultiFrontal method. The software was developed at the University of Florida at Gainesville. Being a direct solver, UMFPACK is substantially different from the iterative solvers Trilinos and PETSc. We use it successfully for indefinite problems arising in the discretization of the time-harmonic Maxwell's equations, where the iterative solvers do not perform well. Documentation and source...

## Example Lshape domain problem

The first numerical example deals with a problem whose exact solution u is known, and thus the error function Ch,p u utip can be calculated exactly. We consider an L-shape domain fi C R2 with a reentrant corner, shown in Figure B.2. Figure B.2 Geometry of the L-shape domain. Considered is the equation Au 0 in il with the Dirichlet boundary conditions u x) R(x)2 i sin(20(sc) 3 + tt 3) for all x dil. Here R(x) and 0(x) are the standard spherical coordinates in the plane. The exact solution has...