## 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...

## Acknowledgments

I acknowledge with gratitude the assistance and help of many friends, colleagues and students in the preparation of the manuscript.1 Tom s Vejchodsky (Academy of Sciences of the Czech Republic) read a significant part of the text and provided me with many corrections and hints that improved its overall quality. Martin Zitka (Charles University, Prague, and UTEP) checked Chapter 2 and made numerous useful observations to various other parts of the text. Invaluable was the expert review of the...

## 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...

## AoxCminQ infi

The weak formulation of the problem (1.26), (1.48) is derived as follows Assume that u C (0) n C1 (O). Multiply (1.26) with a test function v C n) n C1 ), integrate over SI, and use the Green's theorem to reduce the maximum order of the partial derivatives. The boundary integrals do not vanish as they did in the homogeneous Dirichlet case, and we get an extra boundary term, (cii Vw Vu + a0uv) dx - f aigvdS f fvdx. J n Jan ov Jn Here v is the unit outer normal vector to dQ and dujdv Vu-v....

## 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...

## Ba

The potentials A and ip are not unique Many different pairs of A and generate the same fields B and E. It is easy to verify that equations (7.29) and (7.30) are invariant under the transformations < P + fr + C> < 7-31) where C is an arbitrary real constant. Transformations (7.31) and (7.32) are called gauge transformations. Coulomb and Lorentz gauges While the constant C may be made unique by requesting < p(x) > 0 as cc > oo, various uniqueness conditions may be imposed on A and ip....

## 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)...

## Lvv

The following theorem presents an important abstract existence and uniqueness result for operator equations Theorem 1.4 (Existence and uniqueness of solution for strongly monotone operators) Let V be a Hilbert space, f 6 V and L 6 (V,V') a strongly monotone linear operator. Then for every f e V' the operator equation Lu f has a unique solution u V. Proof According to Lemma 1.1 the operator L satisfies the stability estimate (1.21). Moreover, if v 6 fl(L)1, then (Lv, v) 0 and CIMI2 < (Lv. v)...

## 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...

## D

Figure A.31 Parallelogram ABCD in R2. Exercise A.45 Consider a normed space V where the norm satisfies the parallelogram rule (A.62), u + I'll2 + 1 - v 2 2 2 + 2 2 for all u. v V. Show that the identity + v + II2 - II + v - 2 II + w 2 - u - 'II2 + II + tt-H2 - II - wf Exercise A.46 Prove Lemma A.35 ffS is a nonempty subset of an inner product space V, then S1- is a subspace ofV. Exercise A.47 Consider the Hilbert space V R3, endowed with the standard Euclidean inner product, and its subspace W...

## 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,...

## Ddg

The diagonal blocks Bu are defined as Bu I- Atkau (Yk + g,,tk + cAtk). and the nondiagonal blocks By, i j, have the form BtJ -A tkaij Yk+gj,tk + c lAtk). The standard convergence theorem for the Newton's method implies quadratic convergence for sufficiently small Atk. But one has to realize that this convergence rate comes at a high price At each step of the loop (5.100) (5.101) the complete Jacobi matrix (5.102) of the size Ns x Ns has to be reconstructed....

## Dg

The simplified Newton's method assumes the form (J-AtkA Jk)AGn -*(G), (5.103) This iterative procedure is economical, since a single LfZ-decomposition of the matrix I - AtkA Jt only is needed at each time step, i.e., an order of 2(Ns)s 3 operations. Termination criterion for the simplified Newton's method Suppose that the time step Atk is sufficiently small so that the iteration (5.103)-(5.104) converges linearly, i.e., that there exists some contraction factor 0 < ui < 1 so that Let G be...

## 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 3.2 Prove Proposition 3.2. Exercise 3.3 Prove Proposition 3.3. Exercise 3.4 Prove Proposition 3.4. Exercise 3.5 Write the explicit form of the affine inverse map xj from Proposition 3.3. Exercise 3.6 Construct the explicit formula for the inverse of the biaffine map x a Kq > K, corresponding to a convex quadrilateral with straight edges, whose nonpolynomial character was discussed in Corollary 3.1. Exercise 3.7 Consid ra regularfinite element meshT pOver...

## 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....

## Fuxtt gEiiiiuxii

Since x'(t) u(x(t),t) is the slope of the characteristics and u(x(t), t) is constant, also in this case the characteristics are straight lines. A characteristic curve passing through (x,q, 0) has the slope (xo,0) ito(xo)- When two different characteristic curves, carrying two different values of the solution on them, intersect, a discontinuity (shock) is born. This is illustrated in Figure 1.7. Figure 1.7 Formation of shock in the solution u(x, t) of Burgers' equation. Nonlinear hyperbolic...

## 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...

## Amv

Figure B.9 Approximate solution piof the insulator problem. Figure B.10 Details of the singularity of Ei,.,, V ,.,, at the reentrant corner, and the discontinuity along the material interface (zoom 1000). Figure B.ll The Ap-mesh, global view. Large fifth-order elements are used far from the singularity and material interface, small quadratic elements are placed close to the reentrant corner and the material interface The .p-mesh, details of the reentrant corner (zoom The .p-mesh, details of the...

## 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...

## Skp

Figure A.21 Graph of the function (A.42). To define the Riemann integral, cover the interval (a. b) with a partition a a-y < .( i < < x b such that x, -.r,_i < (b-a) n for all 1 < i < n. The Riemann integral is defined as (R)f g(x) dx lim V - x,-,). (A.43) where is an arbitrary point in the subinterval t. .r,). Clearly the limit (A.43) does not exist, because always can either be minus one or one (and thus the result of the Riemann integration can be anything between (b a) and b a)....

## 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

Figure A.7 Linear operator in R2 (rotation of vectors). ( (w2))cu- ((-sin a, cos o)t)bu- (-sina,cosa)T. Therefore is represented by the matrix Choosing now an arbitrary two-dimensional vector v (a,b)T e V, for f(v) we have a cos q b sm a a sin a + b cos a 2. Let us return to the derivative operator D 6 C V, W), D(u) u', from Example A. 10. In the spaces V P4 ( 1.1) and W P3 ( 1,1) we consider the monomial bases By 1, x, x2, x3, x4 3, 4, 5 and Bw l,x,z2,x'J wi, W21 VJ3, W4 , respectively. It is...

## 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.

## Jk

Next assume an edge a of the element K that lies on the impedance boundary F . Let e* be the edge of Kt such that a xk(et). It follows from (7.128) that As before, the symbol tx stands for the unit tangential vector to the edge e,, oriented as shown in Figure 7.7. It follows from (1 .Ml) and (7.128) that the tangential components are transformed via the relation e xK(C)) ta(xK t)) ti(i). (7.132)

## 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...

## Lm

Exercise 3.11 Consider a nodal finite element (Kq, P, E) on the reference square domain Kg with P span l, xi, X2, xf, X1X2, x , xix , x x2 andT, Li,L2, , Lg , where Li(g) g -1,-1), L2(g) g( 1,-1), L3(g) ff(l,l), L4(g) g(-1,1), Ls(g) J g(-l,x2)dx2, L6(g) J g(l,x2)dx2, L7 g) J g x , 1)dx , L (g) J g(x1,l)dx1. 1. Check the unisolvency of this finite element. 2. Construct the corresponding nodal basis (if relevant). 3. Write the formula for the local element interpolantXKq (g), and apply it to the...

## 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

The variational formulation (5.5) is approximated using the sequence of approximations listed in Paragraph 4.1.2 The domain Q is replaced with a simpler domain fi , suitable for meshing, boundary conditions are moved from dil to the new boundary coefficients and data are extended to if fif, < . Q, the space V is replaced with a piecewise-polynomial space Vh,p built on the finite element mesh, exact integration is replaced with the Gaussian quadrature, etc. After inserting the construction...

## 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...

## S

This holds for all right-hand sides if the coefficients satisfy cs 1, bs 0, a-sj bj for all j 1,2, , .s. The remaining coefficients are determined routinely (see, e.g., 25 ). The total number of evaluations of the right-hand side in n steps of the algorithm is only n (s 1) + 1 instead of n s. Therefore we speak about an effectively (s l)-stage method, of the type RKp(p-l). Among numerous known embedded RKp(p 1) methods, the most mature were constructed by J. R. Dormand and P. J. Prince (see 42...

## 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...