Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
A finite element approach on singularly perturbed elliptic equations
(USC Thesis Other)
A finite element approach on singularly perturbed elliptic equations
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
A FIN ITE ELEMENT APPROACH ON SINGULARLY PERTURBED ELLIPTIC EQUATIONS by Andreas Heinrich Fugmann A Thesis Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree MASTER OF SCIENCE (Mathematics) November 1994 UNIVERSITY O F SO U TH ER N CALIFORNIA TH E G RA D U A TE SC H O O L U N IV ER SITY PARK LOS A N G ELES. C A L IFO R N IA 9 0 0 0 7 This thesis, •written by under the direction o f hAJSt Thesis Committee, and approved by all its members, has been pre sented to and accepted by the Dean of The Graduate School, in partial fulfillm ent of the requirements fo r the degree of jA o ^ 'c v Dtan D a te..Jf9^ 9^ 99.39.-1 .. 199.9... THESIS COMMITTEE Acknowledgements I would like to thank Professor Zhiqiang Cai for fruitful discussions and helpful ideas. I also want to thank the other members of my Thesis Com m ittee Professor Edward K. Blum and Professor Chunming Wang. For proof reading and moral support as the work proceeded special thanks to Terri R. Kang. And last but not least thanks to m y housemates Stefan Beifuss, M artin Eberhardt and Christian Molt. They were always there when I needed a break. Contents A cknow ledgem ents List O f Figures ii 1 Introduction 1 2 Transform ation 5 2.1 Different approach w ith Galerkin M e th o d .................................................... 5 2.2 Numerical treatm ent of the p ro b lem ............................................................. 10 2.3 Choice of Basis fu n c tio n s.................................................................................. 16 3 A nalysis of the Galerkin m ethod 21 3.1 The reduced p ro b le m ......................................................................................... 22 3.2 The error in rank one update ....................................................................... 27 4 Suggestions to avoid the error in scaling 33 4.1 The proper b a s i s ............................................................................................... 33 4.2 O ther possibilities............................................................................................... 35 5 C onclusions 38 A p pendix A Program s testing the new a p p ro a c h .............................................................................41 A .l The new m e t h o d ............................................................................................... 41 A.2 Working on the reduced p r o b le m ................................................................ 43 A ppendix B Program s to avoid the error in scaling ................................................................. 45 B .l Program using the new b a s is ........................................................................... 45 B.2 Using an a priori estim ate for the integral ................................................. 46 ii List Of Figures 1.1 Sample p ro b le m .................................................................................................. 2 2.1 Linear basis functions........................................................................................ 18 2.2 Averaged basis f u n c t io n .................................................................................. 20 3.1 Linear basis on reduced problem ......................................... 23 3.2 Averaged basis on reduced p ro b le m .............................................................. 24 3.3 Stepbasis so lu tio n ................................................................................. 27 4.1 Solution with new basis .................................................. 35 4.2 Solution with approximation of in teg ra l....................................................... 37 iii Chapter 1 Introduction We w ant to consider th e problem —eu" + au‘ — / , on fi = (0,1), ( i.i) it(a;) = 0, on c?fi, where e > 0 is a given small perturbation param eter and th e coefficient a is a sufficiently sm ooth arbitrary function. This problem can be treated by the standard Galerkin m ethod, where hat functions are used on a uniform m esh as test and trial space. The error estim ate for this m ethod, shows convergence of order 0 ( j ) . T h at m eans the approxim ation is difficult for problems with small perturbation e. We would have to choose mesh size of 0 (e ), to get a good approxim ation. Take for exam ple the problem: -eu"-\-u> = 1, on fi = (0,1), (1.2) u(a;) = 0, on with e = 0.01, and an uniform mesh size h = 0.1. The solution obtained by the stan dard Galerkin m ethod fails to capture the boundary layer at the outflow boundary x ~ 1. The approxim ate solution is not close to the exact solution, as it oscillates in the whole domain (F ig .:l.l). Figure 1.1: Sample problem 1.25 0.75 0.25 A large variety of different approximation schemes have been suggested to over come these difficulties in the boundary layer. We want to have a stable m ethod that 2 converges uniformly in e. That would allows us to choose the meshsize indepen dent of e. One of the latest advances on this problem was made by E. Stynes, M. O ’Riordan [2]. They developed a Galerkin-Method for * —eu" + afaiju'fa;) 4- 6(a:)u(a;) = /(x ), on = (0,1), < (1.3) u(x ) = 0, ond fl, . with / , 6 6 ( ^ [0 ,1], a € C 2[0,1] and a(®) > a > 0. It approximates / , b and a by piecewise constant functions, and solve this prob lem using a Galerkin-Method with L-splines. This method is O(hi) accurate in dependent of e. The restriction a(a;) > o: > 0 limits this approach to problems without turning points. As equations of this type frequently arise from stationary convection-diffusion problems by linearization, this restriction rules out many prac tical applications. To solve those flow problems, we want to develope a method that allows a more general choice of a. However, in this work we only try to establish a m ethod working for any constant a. This method later would have to be generalized to arbitrary or at least piecewise constant a. In chapter 2 we transform problem eqroriginal to a form, where we have V- ellipticity of the associated bilinear form for an arbitrary constant a. Then we discuss the numerical treatm ent of this problem. Finally we test this method with linear and averaged basis functions. Chapter 3 analyses the disappointing results and points at the sources of the error. This is picked up in chapter 4, where two 3 possibilities to fix the problem are presented. Although the m ethod works with the basis given in chapter 4, the possibility to generalize the algorithm for an arbitrary function a, is doubtful at best. The success of the approximation relies on the fact, th a t this basis allows an interpolation of the exact solution with an accuracy of order e independent of the meshsize h. Our new basis functions are tailored to fit the two boundary layers of the average of the solution for the problem (1.2) and its adjoint. They would probably fail to approxim ate interior boundary layers. 4 Chapter 2 Transformation 2.1 Different approach with Galerkin Method T he m ain reason for the restriction on the coefficient a in M. Stynes, E. O ’Riordans approach [1] is the necessity to prove coercivity of the bilinear form associated with (1.3). Since either a > 0 or —a > 0 either the problem (1.1) almost satisfies the restrictive condition a > 0, so we try to exploit this taking an average of the solution of these two problems. Let u+ be the solution of { -e u " + au'+ = / , on fi = (0,1), ^ u+(a:) = 0, Va; G dfl. 5 And let u_ be the solution of its adjoint: —€.v!!_+au'_ = f, on ft = (0,1), u_(s) = 0, Vs G c?ft, (2.2) and define u = |(u + + u_) and u = ^(u+ — «_). Then the sum of the equations (2.1) and (2.2) leads to -eu" -f aid = / , o n fl = (0,1), u(s) = u(x) = 0, Vs € dt1, and the difference yields (2.3) —eu', + aul = 0, on ft = (0,1), u(s) = u(x) = 0, Vs G dft. We use the latter, solve for u" and for s G ft we integrate over [0, s] to get (2.4) u" — u € f uH (t)dt = f —u\t)dt Jo Jo € u'(x) — ^(O) = -(u (s ) — u(0)) u'{x) : - u (s ) + u'(0). 6 To find the initial value u'fO), we integrate again f u'(x)d; Jo u (l) — u(0) 0 «'(0) x = f — u (s) + u'fOWx Jo £ — f u(x)dx -f- u'(0) £ Jo — f u(x)dx + u'(0) £ Jo f u(x)dx. £ Jo Now we are able to express u‘ in term s of w, we get v!(x) = — u ( x ) f u(x)dx e £ Jo . This result is used to replace v! in (2.3), and we obtain our m ain integral-differential equation in u only -£u" -f- ~ ( u — fo udx) = / , on = (0,1), u(aj) = 0, Va; G dCl. (2.5) We m ultiply this equation w ith a test function v € V = {v G i f 1(fi)|u(0) — u (l) = 0} and integrate over ft. Using integration by parts, we get the associated bilinear form B(u,v) :~ /*1 rl I —£U,r(x)v(x) H u(a;)u(a;) — f u(t)dtv(x)dx Jo £ Jo r l a 2 r l a 2 r l r l I u' (x)v* (x)dx H f (u(x)v(x))dx I u(x)dx I v(x)d\ Jo £ Jo £ Jo Jo a 2 r 1 — e(u', v1 )-\ ((« ,« ) — I u(x)dx I v(x)dx), € JO Jo w here (.,.) is the usual £ 2 innerproduct. So, the appropriate variational form ulation is Find u € V such th at B(u, v) = (f,v ) \/v ,e V . Obviously, the bilinear operator B is sym m etric. B is bounded since |2 ? (tt,i> )| = |e(u /,v1 ) + — (u,u) — — f udx f vdx\ e c Jo Jo < e\(u',v')\ + — \(u,v)\ + — \ [ udx\\ f vdx\ t t Jo Jo B is also V-elliptic, < z 2 f 1 \B(u,u)\ = |e(u/, u7 ) H ( ( « ,« ) — I udx I udx)\ € Jo Jo ~ \e(u\u') + ~ ( f («(a;) ~ f u{t)dtdx) £ Jo Jo (2.6) Using the boundedness of ( /,.) , the Lax-M ilgram-Theorem proves the existence and the uniqueness of the solution for (2.6). 8 Once we have the solution u for problem (2.6) is solved, th e solution it can be com puted as u(x) = f udt — x f udt. Jo Jo A lternatively u can be com puted in th e sam e way as uy because of th e following tedious com putation: ~eu" + ait' = f u" = 1 r - u -----/ e c I I o 1 N . 3 ~(it(x ) - 'u(O)) - - [ fd t e £ Jo = a „ . . 1 f x P , —u ( x ) -----I fd t £ £ Jo U^x) = u'(0) -f- — u(x) — — f fd t 6 £ JO f u((x)dx = Jo f u'{Q)dx H — f udt — — f / fd td x Jo £ Jo £ Jo Jo «(1) — u(0) ~ u '(0 )(l — 0) + — f udx — f f fd td x £ Jo Jo Jo 0 = u'(0) H — f udx — f f fdtdx £ Jo Jo Jo u ' ( 0) = ----1 udx + j 1 fdtdx £ Jo Jo Jo u '(x ) = — — f udt + — f f fdsdt + - f — u(x) - £ Jo £ Jo Jo £ (2.7) 9 Again we use this result to replace u' in (2.4), which gives us the following integral- differential equation for u, - e u " + *£(u - /o1 udx) = f(/o1 fd x - fdtdx), on ft = (0,1), u(:r) = 0, Va: G £?fi. k Except for the right hand side, this is th e sam e equation as (2.5). Once we found a m ethod to solve (2.5) for arb itrary functions / , we could use the sam e m ethod here to get u. Finally, u+ is easily recovered as u + u. 2.2 Numerical treatment of the problem In th e following, let us take uniform m esh on ft = [0,1], of w idth h = Let {a;t -}£_0 be th e nodes, w ith x, = ~, for i = and let be the basis of test and trial space Vh in our G alerkin m ethod. Furtherm ore, we assum e th a t our basis functions have th e following properties, Vi = 1, ...,rc: • 4 > ,‘(ajj) = 5{j, • $,• has local support, i.e. = 0 Vx 0 [si-i,aJi+i], • $i(x) + $ i+i(x ) = 1 Vrr € [a;,-,a;i+i]. 10 The last property makes sure th a t the bases functions interpolate at least constant functions exactly in [x i,x n], for any meshsize h. In the common case, where the bases functions are ju st translations of each other, $ ,(x + a?,--!) = $ i(x ) Vx £ f?, i G {1, ...,n — 1}, we have Vx € [x,-, x,-+i]: <&,(x) = <&i(xf_i + A x), w ith A x = x — x ,„i, = 3>:+i(£(i+i)-i - h + A x) = 1 — $ ,(x , — h + A x) = 1 — $,-(x — h). (2.9) T h at allows us to com pute rl rxi fxi+i f <&i(x)dx — I 3?j(x)dx + I $j(x)rfx JO Jxi— j Jxi JfXi fxi+l ' $ ,(x )d x + I 1 — $ ,(x — h)d; ar,-_i Jxi t> Xi rxi — I $i(x)dx + I 1 — $ ,(x — h)dx Jxi—i J x _ i fXi — I 1 dx = h. J v i-l 11 Following th e standard procedure of finite elem ent m ethods, we discretize the variational problem as ! ■ Find Uh € Vh such th a t B(uh,v) - (f,v ) V w € Vi. Since Uh(x) = f°r some v = ( t y , u n_i) € R n_1, B ( ., .) is bilinear and it is sufficient to check the basis functions of V/,, we have to solve the following problem: { Find v € f?n-1such th at EIL-11 ®j) = ( /, * i) Vi € { 1 , n - 1}. This linear system of equations is usually w ritten as Av = 6, with A = ( B ( $ ;,$ j) ) n_ ixn_i 6 = ((f,* j))n- lxi. In our case the stiffness m atrix A is sym m etric and positive definite, because B ( ., .) is sym m etric and V-elliptic, see [3]. T hat proves the regularity of A. Hence a unique solution of Av — b exists. Before we solve the linear system of equations, we have to com pute the stiffness m atrix. Let A = (ay) € with 2 2 i i ais = B( 4i , = £(«;,*;) + dx [ ijd x . t € JO JO 12 W hile the first two sum m ands on the right hand side are zero for m ost of the indices, £ (< & ;., * ,) = 0 Vjz - j\ > 1. (2.10) 2 j 2 The last term — fQ <&idxf0 <&jdx does usually not vanish. So A is a full m atrix. T hat would m ake the com putation of a solution expensive. B ut the special structure of A allows a different treatm ent. The stiffness m atrix A can be w ritten as A ~ A + D, w ith A = D = (dij) € where «y = €($;,;•) + a2 t 1 f 1 d{j = I $idx I $jdx. e Jo Jo Obviously A is the stiffness m atrix for a reduced variational problem , Find u £ V such th at < ( 2 . 11) 2 ? ( u ,v ) « ( /,v ) Vu G F , k with a2 B(u,v) := e(u',v') H -----(«,«)• 13 B is also a symmetric, V-elliptic and bounded bilinear form, as B is. Therefore, A is again positive definite and symmetric. In addition, (2.10) shows th at A is tridiagonal. Using Gauss algorithm, we need only O(n) flops to solve Av = b. D, on the other hand, can be written as D = — ^ := (Jo $idx) € R n~l . D is only a m atrix of rank 1. Therefore, the solution of Av = b can be computed as a rank one update of the solution for the reduced system Av = b. By the Sherman Morisson formula [5], the solution v of Av = (A + ( - ^ ) $ $ T)v = b (2.12) is v — v — f3w, where v = (3 = w = A~lb, i " 1# . (2.13) So the system Av = b can be solved efficiently for any b. We only have to choose a set of basis functions and pick a quadrature formula to compute b = ((/, $ j))n-ix i = (fo f$ jd x )n- ixi- The basis functions will be our basic concern in Chapter 3. Here, let us put these problem aside for a moment, and compute an approximate solution for u using the approximation U /4 for u, i.e. uh(x) = 53?=/ Uj$,(x). For the 14 convenience of notation, define $o = 0 and < 3 > n = 0. Then an approximation Uh for u can be computed as follows, for x E £*]: f x r l Uh{x) = — {J Uhdt — x j Uhdt) fit rx rl n~i — “ ( / Y V i $ i d t - X / Y , v&idt) C Jo g Jo fri = ~ ( J 2 Vi f $idt+ f vk-i$ k -id t+ f vk$hdt - x Y 'v ; f $idt) e Jo Jxk— 2 J*t-1 Jo This simplifies further, for basis functions with fydx = h , for x E [a;*-!,®*]: a f x fx uh(x) = —( ^ hvi 4- I vk~ i$ k -id t+ vk$kdt - x22(vih)) e , ‘ =1 J x k - 2 ^ x k - 1 ,'=i a fx f x = ~ ( h £ Vi + uf c _i / $ k-idt + vk ^kdt — x h ^ v i ) . £ f = l J x k~ 2 •'Xk - 1 j Consequently the approximate solution u+ is for x G [a:jfc_i, a jjfc ]: '» * f it , f^k— 2 = uf$ ,(s) + - ( 5 3 v{ / $idt+ vk- i $ k-idt t= rl ^ i=l fx y 1 + I vk$kdt - x J Z v i / <M<). •^*-1 ,-_ x Jo In the special case, /q $;dx = h, for x E a:*]: w+(ar) = uh(x) + u A (®) 15 2.3 Choice of Basis functions The easiest choice for basis functions would be to pick hat functions, i.e. *,■(*) = V x€ V x€ [x;,xi+1], 0, otherwise. These functions satisfy all properties required, including jjf $;dx = h. As all these basis functions are translations of each other, i.e. $,(x) = $ i(x — x,-_j), we only need to consider $ 1 and < J > 2 to get the stiffness m atrix A which is tridiagonal. The diagonal elements are f 1 M '.d x + — f 1 dx Jo e Jo 1 1 , a2 fh x x . = 2 el — — da; + 2— / — — dx Jo fin e Jo h n a2 x3 2k + 2 'e 3h2 |A _ o C , 2a2fl l* = ° - z h + ~3T- The off-diagonal elements are 2h —1 1 , a2 f2h 2h — x x — h rin I or t = Cl ~h~h + T A _ h ^ a2 rh h — x x £ £ Jo h h h h dx dx h a2 rh = ------1 ---- I e e Jo h2 -dx £ a 2, x 2 x 3 h a2 f h hx — x 2 - i h t K2h 3 h2)lx=0 £ a2 ,h h . = - h + T < 2 “ 3> _ £ a2 h h £ 6 Since /q1 $idx = h , we have $ = (h,...,h)T. So for (3 in th e Sherm an M orrison Form ula we get w ith w = v4-1 3 > . Now we try this m ethod on the sam ple problem (1.2), where / = 1 and a = 1. T hen 6 = (/q l$>jdx)n- i xi — (h,...,h)T = 3> . T h at simplifies our rank one u p d ate form ula further, because 6 = < & im plies w = v. Hence a _ SElLi1 P ~ £ i v - ^ n - l - t + and finally, v = v — j3w = (1 - P)v 17 = ( 1 “ — i JU c_ h So the solution for Av = b can be obtained by scaling the solution of the reduced problem w ith th e factor 1 (2-14) 1 - a U Z i ® i Rem ark: T he only property of the basis functions we used in above com putations is /q1 <& idx = h . Therefore, this form ula for th e scaling factor can be used for any choice of basis functions th a t has this property. Figure 2.1: Linear basis functions l 0.75 0.5 0.25 0 - 0 . 2 5 - 0 . 5 18 It is not surprising th a t our m ethod does not work properly w ith this basis, since it already failed for the direct approach of (1.2) using Galerkins m ethod, as long as (h > e). Inspired by the results of [eol] a different and m ore promising choice for of basis- functions is average of th e functions $ ,■ and $,■, where is solution of — e$"(x) + a<&f(a;) < t> ,(x ) 0, Vx € (a?i_i, Xi) U(a?.', ® i+i), 1, 0, Vx G (x,_i, x,-+i) j (2.15) and is its analogue for the adjoint problem, -e tfj'(x ) - flWj(x) Vi(x) 0, Vx G (*,-_!, Xi) U(a:* ) ®i+i), 1, 0, Va; G (s-t— l j •Pi+i) ■ (2.16) So for our exam ple (1.2) on the uniform mesh, $ ,(x ) = i Vx G [xf,x i+1], 0, otherwise, (2.17) 19 and 1 - e\ - eheH, Va; € [xi-i, a:;], 1 “ " Z f i * > Va; € ® i+i], 0, otherwise. (2.18) Using Vh = {|($> + ^ = 0.1 and e = 0.01. As figure 2.2 shows, the approx- Figure 2.2: Averaged basis function - 0.2 -0 . 4 im ate solution is oscillating at the boundaries and it is off by a shift. Experim ents show this holds, unless h < e. 20 Chapter 3 Analysis of the Galerkin method A fter trying two different sets of basis functions m ore or less unsuccessfully, let us use th e results to analyze the errors at different stages of th e com putation. As we have organized it, the com putation consists of three basic steps: 1. Solution of th e reduced problem 2. R ank one u p d ate of the solution 3. C om puting u by integration As the last p a rt depends basically only on th e accuracy of th e previous steps and contains no num erical difficulty in itself, because th e antiderivative-s of th e basis functions are available, only th e first two steps are discussed bellow. T he observa tions are based on th e exam ple problem (1.2) treated above. 21 3.1 The reduced problem We first consider the approximation of the reduced problem (2.11) by linear basis functions (3.1) and averaged basis functions (3.2). Both approximations seem to be reasonable accurate as long as the scaling is considered . T hat is exactly what can be expected for a purely elliptic problem like this. Unfortunately, we get oscillations, which are due to the two boundary layers. Recall, we computed u as the average of u+, the solution of (2.1), and u_, the solution of (2.2). Since u+ has a boundary layer at the outflow boundary x = l, has a boundary layer at its outflow boundary - the inflow boundary of u+, because it is the solution to the adjoint problem. We shall see later, that the solution of the reduced problem u differs from u only in a scaling factor. Therefore, the shape of «is the same as the shape of u. The double boundary layer we observe is a result of our construction, to get a V-elliptic bilinear form. Therefore, we had to expect the trouble with the linear basis functions since those already showed oscillations in (Fig.: 1.1. The hope is for averaged basis functions to do better. Because they are developed through the same process as u, from basis functions that approximate and u_ properly. Our hopes are disappointed by the experiment (Fig.:3.2). The cause maybe th at the averaging does not reflect the fact, th at the adjoint solution dominates as we approach 0, where the solution dominates closer to the boundary 1. The problem of both basis functions is the inability to capture the boundary layers. Since this is exactly the same problem that is faced in [1]. The basis functions 22 Figure 3.1: Linear basis on reduced problem 0 . 014 0.012 0.01 0.008 0.006 0.004 0 . 002 used there to approxim ate u become candidates also. They turn out to work, as well as their analogue {’ P,•}” _!, for the adjoint problem. They approximate the solution w ithout oscillations. But each of them fails at one boundary layer, where the shape of the basis function prohibits appropriate fitting. The same thing can be said about similar basis functions, which are m ore closely related to our reduced problem. Let be the solution of = 0, V® € < M Xi) = t 3 ,1 ) = 0, V® < = [®,-_i,®i]c, 23 Figure 3.2: Averaged basis on reduced problem 0.014 0.012 0.01 0.008 0.006 0.004 0.002 And define the basis functions by $ ,(x ) = 4>i(x), Vx e [xi_i,x;], 1 - - h), Vx e [®i, Zi+i], 0, otherwise, (3.2) An a posteriori estim ate, which proves this, can be obtained as follows. Let uj be the interpolant of th e solution u 71— 1 ui(x) u(x,)$,-(x). »=i 24 Let L be the differential operator J 2 _ Lu = —cu" H u. e Looking at th e associated bilinear operator B , we get B{u — Uh-,u — Uh) — B(u — — ui) + B(u — Uh,ui — Uh) = B(u — Uh, u — u j ) + B {u,uj — Uh) — B(iih, tt/ — u/i) = B(u - u h, u - u /) + ( / , u j - fiA ) - ( / , « / - « fc ) = B{u,u — ui) — B{ufau — U{) = ( / , M “ « / ) “ { L U f a U — U j ) < | | / | | | | u - t i / | | + ||Lfifc||||ti - « / | | . (3.3) Rem ark: A t this point we would like to have L uh = or b e tte r L< f?j = 0, which defines a possible choice of basis functions. U nfortunately these functions do not in terpolate constant functions as required. W ith the basis functions (3.2) we, have t Jo {=l = L h ^ u x i ^ d x ) ^ e «=1 J0 a2 n_1 i = t (X > ? a)* e ,=i a 2 < — maaifu,-) (3-4) € Since max(vi) and f are bounded, we get an upper bound using (3.3) and (3.4), a2 B{u — u h ,u ~ u h ) < Ci\\u — ui\\ + C2— l|w — w/||. (3.5) e Using the definition of B , we get a lower bound, a2 — II^-w aH L ^ B { u - u h, u - u h) (3.6) 2 We com pare upper and lower bound and divide by > 0, u - uh\\l3 < - ui 1 1 + C2\\u - 'S/ll CL Then as ||w — u/|[ < C^/i, by a m axim um principle argum ent: a x Ifi-fifcllLa < Ch.*. (3.7) 26 If we look at the formula for the error propagation computing the scaling factor (3.11) we would need to have C3 = 0(e) or C2 = O(e) to succeed in getting a good approximation for u in (2.6). Figure 3.3: Stepbasis solution 0.014 0.012 0.01 0.008 0.006 0.004 0.002 0 3.2 The error in rank one update In our special case with a and / constant, the rank one update turns out to be only a scaling factor. This approach can be justified even for the not discretized version of the variational problem (2.6): 27 Let u be the solution of the reduce variational problem (2.11). And assume the solution u can be found by scaling u, i.e. 3c: u = cu then Vu € Vh ' B{u,v) = (f,v) {cu', u') H -----(cu,v) — j cudx j vdx = ( /, v) G ? fl fl e{u',v')-\-----c(u,v) — c l udx I vdx ~ (/?*>) c Jo Jo c? r1 f 1 c(/, v ) -------c I udx I vdx — ( /, u) £ JO J o c f (1 ju) — — c f udx(l,v) = f(l,v) £ ■ JO c { f I udx) = f £ i/Q C = 1 „ a2 fo Ux 1 c S (3.8) Comparing this result to the discretized version (2.14), the factors differ only in fa udx, which is replaced by its approximation h Zl"=1 The correct factor would be 50, where its approximation is only 9.9909, if we use the basisfunctions (3.2), h = 0.1 and e. The explanation for this behavior is obvious, when we look at the propagation of the error. We try to compute the scaling factor c from the value of the integral I = fo udx as 28 For our exam ple (1.2), w ith e = 0.01, the scaling factor is c (/) = 50. A nd the error in I propagates as A c « c '(/)A 7 = — (1 - — )-2 t P e f } ~ c ( I ) 2A I = 100 502A / (3.9) So the absolute error gets magnified by a factor of 250000. T he situation is a little better, if we look at th e relative error. Since / w e: A c 7 A / 7 7 cC / ) _ r C(I) t f I » c ( I ) ~ = 5 0 ^ . (3.10) So the relative error for th e approxim ation of u gets m ultiplied by the scaling factor. This is, at least for small e, c(f) = (1 — — f udx)-1 £ *J 0 = (1 - \ f e(l - £ J o 1 + e e + (i - (i - e< + 1 J_ 2e (3.11) 29 So the approxim ation of the integral has to have relative error th a t is considerably less than — , which cannot be achieved in general w ith a m esh th a t has e < h. A basis th a t works with our special problem is given in C hapter 4. A nother different look at the problem is obtained by a change in functions for the variational problem. Define v = u ~ r Jo udx. T hat leads to a couple of equation for v : v = u — f udx, Jo v = u", v { 0 ) « (!) J ' vdx o = u (0) — I udx = — f udx, Jo Jo (1) — f udx = — f udx = 4 ” u(0) = u (l), Jo Jo f (u — f udt)dx f udx — f udt = 0. Jo Jo Jo Jo (3.12) Using these equations, we can transform 2.5 transform s into the differential equation —cv + —v = u(0) fo vdx f in D = (0,1) «(!)» 0. 30 The associated variational formulation, with V = f f 1(S2) is Find v G V such that B(v, u) v(0) fo vdx (/,« ), V uGV 0, which is similar to the reduced problem (2.11). Only this tim e we have free peri odic boundary conditions and an integral condition to fix the free boundary value f{ J vdx = 0. The value of the boundary values are related to the solution u by (3.12). After solving this problem, u could be recovered as u (s) = v(x) — i>(0) or u(:r) = u(a;) — u (l). The discretization of the problem is pretty similar to the discretization for the reduced problem, since only the boundary conditions changed. The m ain part of the stiffness m atrix is still tridiagonal, but we need to add some basis functions at the boundaries 1 - $i(ar), V s € [0,h], 0, otherwise, $ o(a0 = < and $ „ (a ;)= 1 — ^ n(cc), Va; € [1 — h, 1] 0, otherwise 31 The integral equation then translates to n r1 Y . Vi I ${(x)dx = 0 J° and the periodic boundary condition remains v0 — Vi = 0. Although the method yields good results for given boundary values, we get ex actly the same approximation for u as with the method discussed before. The integral condition determines the boundary conditions in a way, such that the integral for the approximate solution is zero. In other words the integral over the positive part is the same as the integral over the negative part. If we look at the exact solution, it increases rapidly in a boundary layer around 0 and decreases rapidly in a boundary layer around 1. These two small regions, in fact their width is 0(e ), represent the negative part of the integral. These boundary layers are hard to captured by basis functions, if the meshsize is bigger than e. That means the negative part of v is bound to be approximated badly, and so are the boundary values u(0) and u(l). Exactly the values we depend on to get u from v by a correct shift. 32 Chapter 4 Suggestions to avoid the error in scaling 4.1 The proper basis O ne possibility to get a b e tter approxim ation is the use of a basis th a t approxim ates b o th surges perfectly. Since th e functions < £ , ■ (3.2) fit perfectly for th e boundary layer at 1. T he functions — x ) (4.1) fit th e boundary layer a t 0, because u is sym m etric to and < £ , • fits u a t 1. W e choose * . ' = < A , V i € { l f J + l r c - 1 } , /i, f o r i = | s j . 33 This basis approxim ates u well in the boundary layers, while fx still has to be speci fied. Looking at the conditions on our basis functions, th e interpolation of constant functions should be exact, with k = |_f J , (i{x) = 1 - tj ) k(x), Var € [aub-i,arjfc], 1 Va; € [a:*, a;fc+i], 0, otherwise , The sym m etry relation of and ipi (4.1) implies, th a t there is no change in elem ents of th e stiffness m atrix by changing from < f> i to rfri, for * = 1, • ••> [§ J — 1 V ’. ) = cr = e f ■ tft> ?(x)dx H f ij}f(x)dx JO 6 JO rl = e I <f/?(2xi — x)dx H ------- / <f>f(2xi — x)dx Jo e Jo J rZxi — l f 1 ~4^{x)dx -) j —tf>j(x)dx 2X i € . Jo = e —<j>'?(x)dx + — J^ — <$(x)dx (4.2) T he analogous result holds for £?(<&, < ^ ,- +i). The only change occurs, for the middle elem ent of the stiffness m atrix: &kk — —2 + ^ . 2h 2 h e £ e ( l - e « ) 2 34 The solution u com puted by th at stiffness m atrix is alm ost exact as seen in figure (4.1). The nodal values of the solution are indicated by ”o”, the So there is a set of Figure 4.1: Solution w ith new basis basis-functions, th a t work for this particular problem. 4.2 Other possibilities In both m ethods discussed above, the error was due to an im proper approxim ation of the integral udx. If we could somehow get a good approxim ation for th a t 2 1 integral, we could treat the reduced problem replacing / by / + f0 udx. This 35 new problem then would be easy to solve. We would get a proper approxim ation in L < z sense according to (3.7). One possible estim ate would be the integral of the solution of (1.1) with e = 0 and boundary condition u(x) = 0 only for the inflow boundary. This change has only an im pact in the boundary layer of w idth 0 (e) [5]. T he integral estim ate has an error of 0 (e), which is not bad, because we only consider sm all e. The problem is, th a t this error does not depend on the mesh size, i.e. the approxim ation can not be im proved by refining the mesh. Nevertheless for small e an h > e the results are b etter than the results obtained by th e first approach. Taking the results of this first step to com pute an improved approxim ation for / 0 X udx and repeating the sam e algorithm we get an even b etter approxim ation. T he nodal values obtained using the m ethod once are indicated, by ”0 ” in figure (4.2), for the second step, using th e first one to approxim ate the integral, they are indicated by ”x” . Further iterations lead to now im provem ent, the solutions converge to th e unique approxim ate solution 36 Figure 4.2: Solution with approximation of integral 37 Chapter 5 Conclusions Since we only discussed the simple case that a and / are constant, we are not able to draw any final conclusions, about the new approach of the original problem (1.1). The new problem (2.5), replacing the original has two boundary layers, instead of one. We have to choose a special set of basis functions, that capture both boundary layers. W ith this basis, we get an accurate approximation of the exact solution for h > t. That suggests, it maybe critical for the success of a more general method, to know all boundary layers, including interior ones in advance. The m ain problem of our method is the error propagation in the rank one update. We showed, th at the method is very sensitive to errors made in the approximation of integral of the solution. The relative error is magnified by — in our special case of the rank one update. The new method can be used in cases, where we are able to obtain a reasonable a priori estim ate for the integral of the solution to (2.5). Then we can reduce the problem to purely elliptic one, which can be solved easily. 38 Over all, the new approach already yields a significant amount of problems for the easy case discussed here. Therefore it is questionable, if we are able to generalize this method to cases of arbitrary or piecewise constant functions a and / . 39 Reference List [1] M. Stynes and E. O’Riordan, Finite element method for elliptic convection dif fusion problems, BAIL V, Proc. Fifth Internat. Conf. on Interior and Boundary Layers, (J.J.H . Miller,ed.), Boole Press, Dublin, 1988, p.65-76. [2] M. Stynes and E. O ’Riordan, An analysis of a singularly perturbed two-point boundary value problem using only finite element techniques, M ath. Com puta tion 56, (1991), 663-675. [3] C. Johnson, Numerical solution of partial differential equations by the finite element method, Cambridge University Press, 1990. [4] H. Goering, A. Felgenbauer, G. Lube, H.-G. Roos and L. Tobiska, Singularly perturbed differential equations, Akademie Verlag Berlin 1983. [5] G. Dahlquist, A. Bjork, Numerical methods, 2nd Edition, Prentice Hall 199x (Preprint). 40 A p p en d ix A Programs testing the new approach For the testing we used M athematica Version 2.2 A .l The new method (* approximation to the original problem using the new method *) (* sind linear basis functions *) epsilon=0.01; length=l; n=10; h=l/n; a = l ; (* upward flank of the hat function *) Phi[x_]:=x/h; (* alternatively, the averaged basis functions are used *) (* phi[x_] =(1-Exp[l/epsilon*x])/(1-Exp[l/epsilon*h]); *) (* psi[x_] = (1-Exp[-l/epsilon*x])/(l-Exp[-l/epsilon*h]); *) (* Phi[x_]=l/2*(phi[x]+psi[x]); *) (* downward flank of the hat function *) Phim[x_] :=1-Phi[x] ; (* set of basis functions *) bphi[x_,i_]:=Which[x < (i-l)*h,0, x < i*h ,Phi[Mod[x,h]], x < (i+l)*h,Phim[Mod[x,h]] , True ,0 ] (* bilinear form for the reduced problem *) B[u_,v_]:=epsilon* 41 Integrate [D[u[x] ,x]*D[v[x] ,x] ,-[x,0,h}] + a"2/epsilon+Integrate[u[x]*v[x],{x,0,h}] (* linear form (.,f) for f=l *) b[u_] :=Integrate[u[x] ,{x,0,h}] ; (* fill in part produced by bilinear form *) fillln=b [Phi]+b[Phim]; (* diagonal for tridiagonal stiffness matrix *) Diag=B[Phi,Phi]+B [Phim,Phim]; (* off diagonal for the tridiagonal stiffness matrix *) offDiag=B[Phi,Phim]; (* stiffness matrix for the reduced problem *) A=Table[Switch [Abs[i-j],l,offDiag,0,Diag,_,0], {i,n-l>,{j,n-l>]; (* old version, computed stiffness matrix for full problem *) (* same results as the more efficient reduced+scaling method *) (* A=A“a“2/epsilon*Table[fillIn~2,-Ci,n-l}-,-Cj ,n-l}-] ; *) (* the load vector b for f=l *) bf=Table[fillIn,-Ci,n-l>] ; (* solution of linear system of equations *) sol=LinesLrSolve[A,bf] ; (* scaling that represents rank one update *) sol= 1/(l-h/epsilon*Sum [sol[[i]],{i,1,n-l}])*sol; (* boundary values by boundary conditions *) extsol=Append [Prepend[Append[sol,0],0],0]; (* antiderivative of basisfunction *) inter[x_]=Integrate[Phi[x],x]; (* approximate solution for the averaged problem *) appSol[x_]=Sum[sol[[i]]*bphi[x,i] ,-(i, 1 .n-l}] ; (* approximate solution of the original problem *) corr=Sum[sol[[i]] ,-[i,l,n-l>] ; corSol[x_]=appSol[x] + a/epsilon*( h*(Sum[extsol[[i]],{i,1,Floor[x/h]+1}]-x*corr) +(inter[Mod [x,h]]-inter[0] )* ext sol [ [Ceiling [x/h]+l]] -(inter[h-Mod[x,h]]-inter[0])* extsol[[Ceiling[x/h]]]); (* exact solution of the sample problem *) Clear[eps]; exactSol=DSolve[-(-eps*u ’1 [x] +u1 [x] ==1, u [0] ==0,u[1] ==0},u[x] , eps=epsilon; exaSol=u[x] /. exactSol[[1]] ; (* plot output for comparison *) Plot [-[corSol [x] ,exaSol)-,-Cx,0,1},PlotRange->-{-{0,1} ,{-0 .7,l}}] A.2 Working on the reduced problem (* working on the reduced problem (using a=l,f=l) *) (* upward flank of the hat function *) Phi[x_]=x/h; (* alternatively: *) (* averaged basis functions *) (* phi[x_] = (l-Exp[l/epsilon*x])/(l-Exp[l/epsilon*h]) ; *) (* psi[x_] = (l-Exp[-i/epsilon*x])/(l-Exp[-l/epsilon*h]); *) (* Phi[x_]-l/2*(phi[x]+psi[x]) ; *) (* or new basis functions, that produce no oscillations *) (* Phi[x_]:=(Exp[(h-x)/epsilon]-Exp[(h+x)/epsilon])/ *) (* (1-Exp[2*h/epsilon]); *) (* downwards flank of the test function *) Phim[x_]=l-Phi [x] ; (* definition of basis functions *) bphi[x_,i_]=Which[x < (i-l)*h,0, x < i*h ,Phi[Mod[x,h]] , x < (i+l)*h,Phim[Mod[x,h]], True ,0 ] ; (* bilinear form for the reduced problem *) B[u„,v_]= epsilon* Integrate[D[u[x] ,x]*D[v[x] ,x] ,{x,0,h}]+ l/epsilon*Integrate[u[x]*v[x] ,{x,0,h}]; (* linearform (.,f) for f=l *) b[u_]=Integrate[u[x] ,{x,0,h}] ; load=b [Phi]+b [Phim] ; (* numerical values for the elements of the stiffness matrix *) epsilon=l/100; n=10; h=l/n; Diag=N[B [Phi,Phi]+B[Phim,Phim],40] offDiag=N [B [Phi,Phim],40] A=Table[Switch[Abs[i-j] ,1 ,offDiag,0,Diag,_,0] ,{i,n-l},{j ,n-l>] ; (* numerical vale for the load vector *) bf=Table[load,{i,n-i>]; C* solution of the linear system of equations *) sol=LinearSolve[A,bf]; (* compute the solution as linear combination of basis functions *) corSol[x_]:=Sum[sol[[i]] *bphi[x,i],{i,1,n-l}]; (* compute the exact solution *) Clear[eps]; exactSol=DSolve[{-eps*D [u[x],{x,2}] + l/eps*u[x]==1, u[0]==0,u[l]==0},u[x] ,x] ; eps=epsilon; exaSol=u[x] /. exact Sol [[1]] ; (* plot of exact result and approximation *) Plot[{corSol[x],exaSol>,{x,0,1}, PlotStyle->{{GrayLevel [0] }■, {Dashing [{0.02,0 .02)-] }}-, PlotRange->{{0,1} ,{0,1. 5*epsilon»] ; 44 A p p en d ix B Programs to avoid the error in scaling We used MATLAB Version 4.1 for the programs in this section B .l Program using the new basis t constants ep=0.01; % perturbation parameter h= 0.1; % meshsize exphep=exp(h/ep); * / , i index of test-function i=[ 1 1 2:n 2:n 2:n (n+l)*ones(size(l:(n+1))) ]; 7 , j index of basis-function j = [ 1 (n+1) 1: (n-l) 2:n 3: (n+1) l:(n+l) ]; 7 , index of middle element middle=ceil(n/2); delta=-ep/h+exphep/(exphep-l) % diagonal element of the stiffness matrix Diag=h/ep* (1+exphep-2) / (1-exphep) “ 2; % off diagonal element of the stiffness matrix offDiag=h/ep*(-exphep)/(1-exphep)“2; % building the stiffness matrix as a sparse matrix s=[ 1 -1 offDiag*ones(size(l:(n-l))) Diag*(ones(size(2:n))) offDiag*(ones(size(3:(n+1)))) 1-delta ones(size(l:(n-l))) 1-delta ]; A=sparse(i,j,s); A(middle,middle)=-2+2*h/ep*exphep~2/(1-exphep)"2; A(n+1,middle)=2*delta; 45 ' / , load vector for f=l b=[0 h*ones(size(l:(n-l))) 0]; b(middle)=2*h*delta; % solve linear system of equations v=A\b1; * 1 execute the shift y=v-v(l); * / , plot the exact solution and nodal values x=0:h:1; xs=0:1/1000:1; ys=sol(xs); plot(x,y,1o1,xs,ys); B.2 Using an a priori estimate for the integral 5i c o n s t a n t s a=l; h = l / n ; exphep=exp(h/ep); X using standard Galerkins method with hat functions % on a the original problem and its adjoint with epsilon=0 ' / , and reduced boundary conditions. ' / , i index of test-function Bi=[ 2:n l:n-l ] ; * / , j index of basis-function Bj = [ 1 :n-l 2:n ] ; Bs=[ -0.5*ones(size(2:n)) 0.5*(ones(size(l:n-l))) ]; B=sparse(Bi,Bj,Bs); Bplus=B; Bplus(n,n)=l/2; Bminus=-B; Bminus(l,l)=l/2; * / . treat the reduced problem * / . i index of test-f unction Ai=[ 2:n-l l:n-l l:n-2 ] ; 46 i j index of basis-function Aj=[ l:n-2 i:n-l 2:n-l ]; i , diagonal element of the stiffness matrix Diag=a/ep*(4*exp(a*h/ep)*ep-a*h+a*exp(2*a*h/ep)*h)/ (exp(2*a*h/ep)-l); i t off diagonal element of the stiffness matrix offDiag=2*a*exp(a*h/ep)/(1-exp(2*a*h/ep)); 7. construction of the stiffness matrix As=[ offDiag*ones(size(2:n-l)) Diag*(ones(size(l:n-l))) offDiag*(ones(size(l:n-2))) ]; A=sparse(Ai,Aj,As); b=[h*ones(size(l:n-l))]; v=A\b’; y=[0 v' 0]; % solve the original (eps=0) bp=[b h/2] ; vp=Bplus\bpJ; yp-[0 vp’]; 7. solve the adjoint of original (eps=0) bm=[h/2 bj ; vm=Bminus\bm*; ym=[vm' 0] ; i compute the correction factor with the integral of the i , solution estimated as the average of the integal of the 7, two (eps=0) solutions above corr=l+(sum(yp)+sum(ym))*h/(2*ep); it solve the reduced problem v=A\b*; y= [o V1 0] ; it scale the solution of the reduced problem vc=v*corr; % take the improved estimate for the integral % and solve the reduce problem with new left side bc=b+h'‘2*sum(vc)/ep; vcl=A\bc>; yc=[0 vc' 0]; 47 y c l = [ 0 v c l J 0] ; 5 1 p l o t t h e s o l u t i o n s t o g e t h e r w i t h t h e e x a c t s o l u t i o n x = 0 : h : i ; x s = 0 : 1 / 1 0 0 0 : 1 ; y s = s o l ( x s ) ; p l o t ( x , y c , ' o ' , x , y c l , ' x ' , x s , y s ) ; 48 INFORM ATION TO U SE R S This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter face, while others may be from any type of computer printer. The quality of this reproduction is dependent upon the quality o f the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproduction. In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand comer and continuing from left to right in equal sections with small overlaps. Each original is also photographed in one exposure and is included in reduced form at the back of the book. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6" x 9" black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. A Bell & Howell information Company 300 North Zeeb Road. Ann Arbor. M l 48106-1346 USA 313/761-4700 800/521-0600 UM I Number: 1376453 UMI Microform 1376453 Copyright 1995, by UMI Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. UMI 360 North Zeeb Road Ann Arbor, MI 48103
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Implementation aspects of Bézier surfaces to model complex geometric forms
PDF
Independent process approximation for the coupon collector's problem
PDF
Trifluoromethanesulfonates (triflates) for organic syntheses
PDF
The analysis of circular data
PDF
Complementarity problems over matrix cones in systems and control theory
PDF
Structural equation modeling in educational psychology
PDF
Diminishing worlds?: the impact of HIV/AIDS on the geography of daily life
PDF
A General Approach To One-Step Iterative Methods With Application To Eigenvalue Problems
PDF
Pulse oximetry failure rates
PDF
Altered interaction of human endothelial cells to the glycosylated laminin
PDF
The study of temporal variation of coda Q⁻¹ and scaling law of seismic spectrum associated with the 1992 Landers Earthquake sequence
PDF
Structural geology of the Chiwaukum schist, Mount Stuart region, central Cascades, Washington
PDF
Hand function in older adults: the relationship between performance on the Jebsen Test and ADL status
PDF
Performances of orthodox and liberal Jewish children in third grade on praxis subtests of the SIPT: a cross cultural study
PDF
How casual attribution influences behavioral actions taken by self-aware individuals
PDF
The Japanese demonstratives ko, so and a
PDF
"It was in the air": Adolph Gottlieb, the pictographs, New York, and the Zeitgeist of the 1940s
PDF
The effects of dependence among sites in phylogeny reconstruction
PDF
A Kalman filter approach for ionospheric data analysis
PDF
Air terrorism and the international aviation regime: the International Civil Aviation Organization
Asset Metadata
Creator
Fugmann, Andreas Heinrich
(author)
Core Title
A finite element approach on singularly perturbed elliptic equations
School
Graduate School
Degree
Master of Science
Degree Program
Mathematics
Degree Conferral Date
1994-11
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
mathematics,OAI-PMH Harvest
Format
masters theses
(aat)
Language
English
Contributor
Digitized by ProQuest
(provenance)
Advisor
Cai, Zhiqiang (
committee chair
), Blum, Edward K. (
committee member
), Wang, Chunming (
committee member
)
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c18-1272
Unique identifier
UC11357956
Identifier
1376453.pdf (filename),usctheses-c18-1272 (legacy record id)
Legacy Identifier
1376453-0.pdf
Dmrecord
1272
Document Type
Thesis
Format
masters theses (aat)
Rights
Fugmann, Andreas Heinrich
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA
Tags
mathematics