Efficient methods for nonlinear parabolic PDEs

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content EFFICIENT METHODS FOR NONLINEAR PARABOLIC PDES Copyright 2003 by Dong-jin Kim A Dissertation Presented to the FACULTY OF TH E GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements of the Degree DOCTOR OF PHILOSOPHY (MATHEMATICS) May 2003 Dong-jin Kim Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UMI Number: 3103916 Copyright 2003 by Kim, Dong-jin All rights reserved. ® UMI UMI Microform 3103916 Copyright 2003 by ProQuest Information and Learning Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. ProQuest Information and Learning Company 300 North Zeeb Road P.O. Box 1346 Ann Arbor, Ml 48106-1346 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UNIVERSITY OF SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIVERSITY PARK LOS ANGELES, CALIFORNIA 90089-1695 This dissertation, written by q - f t L . w i __________ * under the direction o f h isS dissertation committee, and approved by all its members, has been presented to and accepted by the Director of Graduate and Professional Programs, in partial fulfillment of the requirements for the degree o f DOCTOR OF PHILOSOPHY Director Date August 1 2 , 2003 Dissertation^Committee Chair \ y ^ A f / \ i i) \ Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Dedication To My Family, Especially my incredible soul-mate, Susan Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Acknowledgements My sincere thanks are due to the members of my committee: Dr. Wlodek Proskurowski (Chairperson), Dr. Firdaus E. Udwadia (CoChairperson), Dr. Ed ward Blum and Dr. Igor Kukavica for their assistance and encouragement in the completion of this Ph.D. dissertation. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Contents D edication ii A cknow ledgem ents iii A bstract ix I An efficient approach for solving a class of nonlinear 2D parabolic PDEs 1 1 Introduction 2 2 M odel problem: porous m edium equation 5 3 Linearization approach 8 3.1 Standard linearization m e th o d ............... 9 3.2 Operator splitting m e th o d ............................... 13 4 Porous m edium equation w ith strong absorption 16 5 Num erical experim ents 18 5.1 Efficiency of the splitting m e th o d .............. 19 5.2 Finite extinction phen o m en o n..................... 23 6 Concluding remarks 35 7 Appendix: Linear stab ility analysis 36 II Approxim ate Inertial Manifold with Perturbation Tech nique 39 8 Introduction 40 iv Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9 Prelim inaries 43 10 A pproxim ate inertial m anifold w ith perturbation technique 48 10.1 Perturbation tech n iq u e................................ 48 10.2 Construction of approximate inertial manifolds ....................... 55 11 N um erical experim ents 65 11.1 Euler-Galerkin method vs. our perturbation m ethod....................... 74 12 C oncluding remarks 92 III Inertial m anifold w ith perturbation technique 94 13 Inertial manifold w ith perturbation technique 95 13.1 Inertial manifold .......................................................... 95 13.2 Reaction diffusion equations with a polynomial nonlinearity . . . 100 14 C onvergence results 103 B ibliography 113 v Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. List o f Figures 5.1 Comparison of the efficiency............................................................... 29 5.2 Comparison of the a c c u ra c y ........................................ ... 30 5.3 Comparison of the accuracy vs. efficiency..................... ..... 31 5.4 Comparison between A and P ............................................................ 32 5.5 Time evolution of the maximum height .................... 33 5.6 Finite extinction p h en o m en o n ................................. 34 11.1 Convergence result of Example 1 for M = 4 0 ....................... 78 11.2 Convergence result of Example 1 for a = 1 .................................... 79 11.3 Accuracy of Example 1 for M — 4 0 ....................... 80 11.4 Accuracy of Example 1 for M — 8 0 .................................................. 81 11.5 Nonlinear term and its first and second derivatives of Example 1 82 11.6 Time evolution of the numerical solution of Example 1 83 11.7 Convergence result of Example 2 for M = 4 0 ................................. 84 11.8 Convergence result of Example 2 for b = 1 .................................... 85 11.9 Accuracy of Example 2 for M — 4 0 .................................................. 86 vi Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11.10 Accuracy of Example 2 for M — 8 0 ...................................... 87 11.11 Nonlinear term and its first and second derivatives of Example 2 88 11.12 Time evolution of the numerical solution of Example 2 . . . . . 89 11.13 The importance of P-mode for the better accuracy of our P2 . . 90 11.14 Comparison of the accuracy vs. efficiency............................. 91 vii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. List of Tables 5.1 Comparison of the efficiency for Example 1 ................................. 22 5.2 Comparison of the a c c u ra c y ........................................................... 23 5.3 Comparison of computational costs (per temporal s te p ) ................ 25 5.4 T* for fixed h = .0 4 ......................... 26 5.5 T* for fixed k — l e - 5 .................................. 26 5.6 Modified algorithm, T* for fixed k = le-5 after T ° .......................... 27 11.1 Comparison of the accuracy for Example 1 (non-polynomial) . . . 76 11.2 Comparison of the accuracy for Example 2 (polynomial) ............ 76 11.3 Comparison of the efficiency for Example 1 (non-polynomial) . . . 77 11.4 Comparison of the efficiency for Example 2 (p o ly n o m ial)............ 77 viii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A bstract In Part 1 we consider a class of nonlinear 2D parabolic equations that allow for an efficient application of an operator splitting and a suitable linearization of the discretized problem. We apply our scheme to study the finite extinction phe nomenon for the porous medium equation with strong absorption. A comparison of accuracy and computational efficiency of the resulting algorithms for several test problems are presented. In P art 2 we present an explicit form for the con struction of an Approximate Inertial Manifold for a class of nonlinear dissipative partial differential equations by using a perturbation technique. We investigate two numerical examples of the reaction diffusion equation with polynomial non- linearity and non-polynomial nonlinearity. We show a comparison of accuracy and computational efficiency for the Euler-Galerkin method and our perturbation method. The Inertial Manifold and its convergence results are presented in Part 3. We introduce an Inertial Manifold in the form of a power series expansion by using a perturbation technique. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Part I An efficient approach for solving a class of nonlinear 2D parabolic P D E s 1 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 1 Introduction We study an efficient use of the linearization approach for a class of nonlinear 2D parabolic PDEs where the nonlinearity has the form of powers of the solution function. We apply the operator splitting technique in our approach to improve the efficiency of the scheme. This algorithm is then used to numerically compute the extinction time to high accuracy for the porous medium equation with strong absorption. We consider a finite difference method for nonlinear 2D parabolic PDEs, and we use an implicit scheme of the Crank-Nicolson type in view of its superior stability properties. The discretization of nonlinear PDEs leads to systems of nonlinear algebraic equations. These can be solved using Newton’s type methods, which is costly. Instead, one can use a linearization approach that eliminates the need for iterations at each temporal step. Solving the arising large linear system 2 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. of algebraic equations straightforwardly using direct methods would be extremely inefficient. Thus, the linearization by itself does not solve the main difficulty, the high computational cost, in dealing efficiently with two-dimensional problems. To remedy this, to the linearized equations we apply the operator splitting technique that allows for computationally more efficient solution processes. The nonlinearity in the PDEs considered has the form of powers of the solution function. This is quite common situation in many applications, especially, in the porous medium type equations, see [1], [2], for which various numerical schemes have been introduced, see for example, [2], [6], [7]. We apply our scheme in the porous medium equation, and we use it to study the finite extinction phenomenon for the porous medium equation with strong absorption. Throughout the paper we use the following finite difference operator nota tion. A continuous function u = u(x, y, t) is discretized on an equidistant spatio- temporal grid (Xi,yj,tn) = (ih ,jh ,n k ), 1 < i,j < M , 0 < n < N, where h = A x = Ay, k = A t are step sizes in spatial and temporal directions and n = 0 is for an initial function. The discrete function at the point (Xi,yj,tn) is denoted by uithn. The forward difference operator in each direction is defined as D+t ■ D .±.tU n — ^n+l 'U /n k D+x • D+xUi — D + y : D + yU j Uj+1 Uj h 3 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. To simplify the notation, indices over variables that do not change under the particular operator (are temporarily frozen) are suppressed. The backward difference operator is defined analogously as n p, _ — U i-\ ^ D— x * T'— rg U i — , etc. h These one-sided operators are first order accurate approximations of the deriva tives, i.e., (D+x — D)ui = Oih), etc. Additionally we use the second order opera tor: 7 v p . n n ui+i 2Ui + iq_i - b ' + s - b ' —x'U'i J-'— x -C '+ x ^ i — 2 ; which is second order accurate, i.e., (D+XD -X — D 2)ui — 0 {h 2). This paper is arranged as follows: In Chapter 2, we present our model problems in some detail. Chapter 3 is devoted to the analysis of computational efficiency of the different algorithms for porous medium equation. Chapter 4 contains the analysis for the porous medium equation with strong absorption. Numerical ex periments for several problems with various conditions are reported in Chapter 5, while Chapter 6 contains the concluding remarks. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 2 M odel problem: porous m edium equation Let us consider the following initial-boundary value problem as our model problem: ut — A (um) — cuP in x I, (2.1) «|anxi = 9, (2.2) u(x,y, 0) = uq(x , y), (2.3) where A denotes the Laplace operator in j R2, fl C R2 is a bounded spatial region, and / = (0, T), T < oo, is a time interval. We will use this model problem to describe our method. If the absorption term cuP is absent, i.e., c = 0, and if m > 2, the equation (2.1) is well-known as the porous medium equation. From the application point of 5 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. view, the object of the porous medium equation can represent the scaled density of the gas [1]. Consider an ideal gas flowing isentropically in a homogeneous porous medium. The flow is governed by the following three laws: 1. Equation of state: p = PoPa, (2.4) where P = P (x,y ,t) is the pressure, p = p(x,y,t) is the density, and a E [1, oo) and P0 E R + are constants. 2. Conservation of mass: ypt + div(pv) = 0, (2.5) where v = v(x, y, t ) is the velocity vector and rj E R + is the porosity of the medium. 3. Darcy’s law: i/v = — pV F , (2.6) where v E R + is the viscosity of the gas and p, E R + is the permeability of the medium. If we substitute P and v from (2.4), (2.6) into (2.5), we obtain ypt + d i v ( - ^ p V p a ) = 0, i.e., = ^ ^ d iv (p V p Q). U T ] 6 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. then we obtain the porous medium equation ut — A(um), m > 2. This porous medium equation commonly occurs in many other applications [2]. The nonlinear partial differential equation (2.1) with ra >2, 0<p<l,c>0 is called the porous medium equation with strong absorption. It appears in various applications, most frequently to describe a phenomenon of thermal propagation in an absorptive medium, where the solution u stands for temperature. In other applications u is a concentration and the process is described as diffusion with absorption. The solution u is required to be nonnegative, and thus we consider nonnegative and bounded initial condition u0(x,y). Furthermore, in the case of zero boundary conditions, i.e., g(x) = 0, it is known (see [4]) that the solution u vanishes identically as time t goes on and extinction in finite time arises, i.e., there exists a time T* > 0 such that u(-, t) is not identically zero for 0 < t < T* but u{-, T*) = 0. Here, T* is called an extinction time of a solution u and the behavior is known as finite extinction phenomenon. 7 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 3 Linearization approach Consider the model porous medium equation (without absorption): ut — A (um), m> 2. Applying the standard Crank-Nicolson scheme (which is 0 ( k 2 + h2) accurate) one obtains D+tun = - ^(D+XD -X + D +yD -y)u™ + (D+XD -X + D+yD - y)u™ +- ^ . Rearranging the un and un+i terms gives k k u n + 1 — —(D+XD ^X + D+yD_y)u™ + 1 = un + - ( D + XD^-X + D +yD_y)u™. (3.1) This yields a system of nonlinear algebraic equations that can be solved using Newton’s iterative method at each temporal step. At (n + l)th temporal step, this system of nonlinear equations is F{u) = u - ^ (D +XD„X + D+yD^y)um — 6 = 0, 8 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where u = (ui.i, w2 li,- ■ b = un + §( D+XD-X + D +yD^y)u™, and un is the computed solution at the previous time level, a known value. Then Newton’s iterative method is given by J ( u (,))tu('+1) = ~F (u{l)), where w/^+1) = u^l+ 1^ — u®, the Jacobian m atrix J(u®), M 2 by M 2 is ,/(„«>) = l _ ^ ( D+XD . X + D+yD-sJtuW)™-1, and I is the M 2 by M 2 identity matrix. The attraction of the Crank-Nicolson scheme is its unconditional stability, see Appendix, for all values of A (in practice though one uses moderate values of A , to avoid slowly decaying oscillations). The main problem one needs to deal with in solving (3.1) is its computational cost. The Jacobian matrices, M 2 by M 2, easily become very large, and thus computational cost is extremely high. One can somewhat improve the efficiency of the scheme by exploiting the spaxsity of the Jacobians, see [10]. 3.1 Standard linearization m ethod Here, we follow the idea of Richtmyer [8] who applied it to one dimensional prob lems ut — (um)xx (some call it the Richtmyer’s linearization method, see [9]). Recently, this idea also was applied to Korteweg-de Vries equation, see [5]. 9 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. For m > 2, u™ + 1 = u™ + mu™ ^rtn+i - un) + 0 ( k 2), or <+1 = < ' + " » < 1“n+1, where un+i « wn+i - n„. Substituting (3.2) into (3.1) and rearranging we obtain: [J - ~ y {D+xD ^x + D + y D -y ^^ju jn + i = k(D+xD_x + D+yD_y)u™, (3.3) which induces a linear system in the unknown u n+\. Eqn.(3.3) needs to be solved at each time level n = 1,2, ■ • •, N. In matrix form it can be written as (3.2) A l u . n -fl (3.4) where A is M 2 by M 2 block tridiagonal matrix / T\ D2 A D i T 2 D$ D m - 2 T m - i D ,M V / (3.5) Dm - i Tm each blocks Tj s and D,f's are M by M tridiagonal and diagonal matrices, respec tively. They are not constant diagonal (Toeplitz) but depend on the solution u at the previous time level n: T, = ! + ?£■%, 10 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. T- = ± 3 4i/m— 1 — 7 /m_1 i a l ,j,n U2,j,n -U m — 1 1 ,j,n 4:<- m —1 — U‘ m —1 3 ,j,n -u m —1 M —2,j,n mX u m —1 1 ,j,n U . m — 1 2 ,j,n \ for 1 < j < M and A = k /h 2. The M 2 vector b in (3.4) has the form: -U m —1 M,j,n M ,j,n J \ U m — 1 M —l,j,n U m —1 M ,j,n j (3.6) A A where u™\dn is M 2 vector with the contribution from the boundary, 11 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. | __ I . ..m l uk |an - uk |anx + uk |an„ u 1,0, k V, 2,0 ,k UM - 1,0,f e u M,0,k n.m Ul,M + h k U, 2,M+l,fc U y u M ,M + l,k j where indices 0 and M + 1 denote the boundary values, B is M 2 by M 2 block tridiagonal matrix, I is the M by M identity matrix, and T is M by M tridiagonal + “ 0,1, k 0 m"1 u 0,2,fc U U, 0 ,M,k 0 y w M + l,M ,k j matrix 12 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. T I -4 1 I T I -4 1 B = A , T = I T I -4 1 V / T / V / It should be noted that the linearized equations (3.3) are identical to those obtained by using only one iteration of the Newton’s method with the initial guess = un, the computed solution at the previous time level. This can be verified simply by inspection. The cited papers failed to make this observation. Again, the straightforward attem pt to solve system (3.4) is unacceptably costly. The half bandwidth of A is M, and thus applying a band matrix solver would cost 0 { M 2) flops per unknown, far from optimal. 3.2 O perator sp littin g m eth od In order to remedy the high computational expenses in solving system (3.4) we propose the operator splitting, each in only one coordinate direction. We add a benign (because it does not affect the 0 (k 2 + h2) accuracy of the Crank-Nicolson scheme) term with a truncation error 0 ( k 2): T T l^ k2 - j - i 13 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. to the left hand side of eqn.(3.3). As a result we can factorize the left hand side of eqn.(3.3) as follows: ( / - ^ D +,D _ « » r 1 ) ( / - = b, (3.7) where b — k(D+xD^x + D+yD -y)u™. By splitting eqn.(3.7) into two simpler systems we have: { i - ~ D +xD ^ - ' - y n+i = b, (3.8) = “V ), (3-9) where ujn+i is an intermediate value and u n + 1 = un + u > n+i- When the bounded spatial region 0 is a square, [a, 0\ x [a, (3\, the intermediate boundary values can be given from (3.9) by ( T T tJ z \ I ~ ~ ^~ D + y D -ygn J idn+l ~ 9n)j where dQy = {(x,y)\x = a,p, a < y < /?}, u (x ,y ,t) = g(x,y,t), for (x,y,t) € SO x J. The structure of the matrices in (3.8) and (3.9) is identical. They differ in the ordering of variables, horizontal (along x-axis) in (3.8) and vertical (along y-axis) 14 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. in (3.9). The matrices in equations (3.8-3.9) are block diagonal with tridiagonal blocks: for eqn.(3.8) Tx A Tm- i T- \ mA - / + 2 T „ J-M t y j,n U2,j,n \ m - i l,j,n “ 2,j,n —?/m _1 3 j,n —?/m _ 1 UM -2 ,j,n 9?/m — 1 ,im ~ 1 a M,j,n V —n m ~ l M —l,j,n ZU M,j,n y for 1 < j < M. For eqn.(3.9) the indices in 7} are transposed: u instead of uT,flni f°r 1 < 'i-,J < M. Since the entries depend on the solution u the matrices in (3.8) and (3.9) are identical only if the solution is symmetric in x and y. Each tridiagonal block can be solved independently of each other at the cost proportional to its size, i.e., O(M). Since there are 2 M blocks to be solved the total cost is 0 (M 2), or G( 1) flops per unknown, which is optimal. 15 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 4 Porous m edium equation w ith strong absorption In this chapter we study the finite extinction phenomenon for two-dimensional porous medium equation with strong absorption. ut = A(um) - cup, (4.1) where m > 2 , 0 < p < 1 and c > 0 . As in chapter 2 , eqn.(4.1) can be discretized using the Crank-Nicolson scheme: 0 ( (D+xD— x D-\.yD_y)un -b (^D+xD— x "b D+yD— y)uinjr- ^ ) - | « + < + 1 ), or k Un+ 1 _ ~^(.D+XD — X + D +yD -y)U n+l = un + ~(D+XD -X + D+yD^y)un — + w ^ +1). 16 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Since m > 2 , applying the Richtmyer’s linearization method to u™ + 1 and rear ranging gives ( l - + D+,D .„ ) u r ‘) ^ + i = *(£>+,£>., + D + .B - .K - y « + < + i), where a > r a + 1 = w n + 1 — w n. Applying the operator splitting we obtain a system of two equations, similar to those in (3.8-3.9): (/ - = b, (4.2) ( TV) n \ I - — D+yD ^ ,v ^ - 1Jujn+1 = ujn+i, (4.3) where u> n+i is an intermediate value and K f* b = + o + .o -.K ' - y « + <+i)- Here, however, in the presence of the additional nonlinear zero order term cup, the first system (4.2) is nonlinear (because of the term u^+1). To solve it we use Newton’s method, again with the Jacobian matrix that is block diagonal and thus each tridiagonal block can be solved independently. As a consequence, the computational cost can significantly be reduced as in the case of the porous medium equation (without absorption). 17 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 5 N um erical experim ents In this chapter we consider two different numerical experiments. First, we inves tigate the efficiency of our linearization with splitting method for solving two- dimensional problems. We compare three different methods: two iterations of Newton’s method (denoted by N in Tables and Figures), and the linearization method, first in its standard form (denoted by L), then combined with the op erator splitting (denoted by S), in terms of the computational efficiency (cost). We also comment on the accuracy of the schemes. Second, we study the finite extinction phenomenon for the porous medium equation with strong absorption employing our scheme. We present the experimental results of the numerical extinction time values for various spatial and temporal step sizes. All our pro grams, written in Matlab, are implemented on a PC with Pentium III processor at 933MHz. 18 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5.1 E fficiency o f the sp littin g m eth od We consider problems defined in a unit square spatial domain divided equidistantly into M grid points in the x- and y- directions. Thus the spatial step size of the uniform grid becomes h — To study the computational efficiency (cost) of the schemes (per temporal step) we compare the CPU time and the number of Mflops for M=20, 40 and 80 with a fixed value of A — p- = 1 , where k is the temporal step size. We report on the cost reduction factor as well as the power, p, of the cost’s growth model: cost = (D(MP) (computed as the slope of the least- squares linear polynomial of logarithmically scaled points). The total cost is then increasing linearly with the number of temporal steps, i.e., proportional to the value of It should be noted that the Crank-Nicolson scheme allows for the use of higher values of A than of the order of 1 , see Table 5.2. The normalized L2-norm of the error is defined as: lerrorll = ~- 1 1 1 M M (Ui,j,n~u(ih ,jh ,T ))2, ij =1 where u ^ n is the numerical solution and u(ih, jh,T ) is the analytical solution at the time t = nk = T, and h and k are the spatial and temporal step sizes, respectively. We investigate two numerical examples which are different initial-boundary value problems of the same porous medium equation ut — A (u5). 19 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. • Example 1. We choose initial and boundary conditions corresponding to the exact solution of ut = A (u5) which is defined as [ 1 0 ] u(x,y,t) = c/^(2f + x + y), (5.1) for 0 < x, y < 1 , 0 < t < 1 . • Exam ple 2. Similarly, we consider the initial-boundary value problem corresponding to the exact solution 1 J x 2 + y2 u{x, y, t) = x_ t , (5.2) for 0 < x,y < 1, 0 < t < 0.5. The exact solution of the Example 2 (given analytically by Barenblatt and thus called the Barenblatt solution) is a radially symmetric self solution u > 0, defined on {(x, i)|x 6 RN, 0 < t < T} ([1], [3]). It is of the general form / \ l/(ro — 1) u(x, t) = ( T - (A T K + £ T ||x ||7 ( T - t f - K \ where K — C = ^ = n 5 an< ^ A > 0. For the Example 2, we choose N = 2, T = 1, m = 5 and A — 0, thus obtaining (5.2). The computational cost per temporal step of solving the linear system repre sented by system (3.4) is the critical component of the total computational ex penses. A brute force band Gaussian Elimination of the M 2 by M 2 system with the half-band width M would cost 0 (M A ) flops, a prohibitively high expense. 20 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The pentadiagonal (or block-tridiagonal) matrices (3.5) have the same struc ture as the discrete Laplacian, although they are not constant diagonal (Toeplitz) and thus do not allow for the use of FFT-based Fast Solvers. One possibility of reducing the cost is to apply nearly optimal reordering before the Elimination. Such a tool is provided in Matlab (as a default option) using the backslash (\) operation in the sparse mode; it employs the Minimum Degree ordering algorithm. Because of this the flop count for solving (3.4) decreases to about 0 (M 36), see Fig.5.1. At the same time, for the splitting method, see Fig.5.1 the computational cost for solving (3.8-3.9) is of about 0 ( M 2) flops, or 0(1) per unknown, which is optimal. Thus, our analysis from section 3.2 is confirmed by the experiments. One should point out that the cpu growth factors are significantly closer, of about 0 ( M 3) and 0 (M 2'5), respectively. This measure of efficiency is much more computer and language dependent. In this implementation the difference in cpu time between the two methods is less pronounced than the flop count. Neverthe less, for M = 80 the actual cost (in cpu) reduction ratio from Newton’s method to our operator splitting method is about 8 times (see Table 5.1). The computational cost for both Examples 1 and 2 is almost identical, therefore we provide the data only for one of them. We also provide the comparison of the accuracy of the methods considered. One should remark that the accuracy is heavily problem dependent as the dis cretization error depends not only on the numerical scheme and the parameters of 21 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 5.1: Comparison of the efficiency for Example 1 M h Number of Mflops CPU time (msec) N L S N L S 2 0 .0476 .53 .26 .03 48 24 1 2 40 .0244 5.13 2.56 . 1 0 281 138 65 80 .0123 78.13 39.06 .41 3,362 1,669 442 discretization but also on the smoothness of the solution. Because of this we pro vide the two Examples. The computed solution in Example 1, Fig.5.2a, is much more accurate, in the range of 1 0 ~ 7 to 1 0 “ 8 (except for the standard linearization method), than in Example 2, Fig.5.2b (10~ 5 to 10-6) for the given range of grid sizes, M. Here, p is the power of the discretization error model: error = 0 ( M P) with the theoretical value for smooth enough functions of p = — 2. The exper iments indicate that the accuracy of our method is almost the same as the one after two iterations of Newton’s method, also see Table 5.2. We also verified the influence of larger values of A on accuracy, see Table 5.2. Only for the coarse grids the slowly decaying oscillations affect the error. To bring yet another perspective we plot the comparison of the accuracy vs. efficiency of the considered algorithms, see Fig.5.3. It clearly shows the superiority of the splitting method. As a consequence, the results presented in the next Section were obtained only with the latter. 22 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 5.2: Comparison of the accuracy Example 1 Example 2 M A = 1 A = 10 A = 1 N L S N L S N L S 2 0 2.3e-8 2.1e-7 2 . 1 e- 8 2.0e-7 5.9e-4 1.5e-7 5.0e-6 5.1e-6 4.9e-6 40 5.9e-9 1.9e-8 5.8e-9 7.4e-9 1.3e-6 6.3e-9 1 .8 e- 6 1 .8 e- 6 1 .8 e- 6 80 1.5e-9 2.3e-9 1.5e-9 1.6e-9 8 .8 e - 8 7.0e-10 6.5e-7 6.5e-7 6.5e-7 5.2 F inite extin ction p h enom en on To study the finite extinction phenomenon numerically we consider the initial- boundary value problem (2.1 - 2.3) with m = 2, p = 0.5, c = 5 and zero boundary condition (g = 0). The presented results were obtained with the splitting method. • Exam ple 3. We choose the initial function u0(x,y) given as (see [17]): , 1 if (x, y) — (0 , 0 ), U Q (x,y) = \ (5.3) where O = [ - 1 .2 , 1 .2 ] x [ - 1 .2 , 1 .2 ]. We denote M to be the number of spatial steps in the x- and y- directions, i.e., (M — 1) is the number of grid points in each direction. Thus the step size of the uniform grid is h = The second equation of our system, (4.3), is linear. We solve it just as described in Section 5.1. On the other hand, the first equation (4.2) is nonlinear. To solve 23 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. it we use Newton’s iterative method: at the (n + l)th temporal step, the system of (M — 1) decoupled nonlinear equations in the x-direetion is F(u) — u — kD+xD -xUnU + 2.5 k^/u — 6 = 0 where u = • • •, uitM-i)T, i — 1 , 2 , ...,M — 1 and 6 — 'U ’n ~ t" 2.5k-\Jun with the known value of un. The Jacobian matrix J(u) is J(u) = I — kD+xD-xun + 1.25 k/y/u, where I is the (M — 1) by (M — 1 ) identity matrix. The Newton iterations are performed until the stopping criterion < t, s = 0 , 1 , • • • is satisfied. Here, the initial guess is u ^ h — un for n — 0,1, ■ • •, and r is the given tolerance. To prevent the singularity of the Jacobian matrix, we use regularization, i.e., Jnew — •/ + using the machine precision e ~ 2.2e— 16. Moreover, we force the solution of (4.2 - 4.3) to remain non-negative at each tem poral step. On the average, the algorithm uses only about two Newton’s iterations to solve the nonlinear equation (4.2). As a consequence, the computational cost per temporal step is only about twice that of the cost discussed in Section 5.1 (without the absorption term ), see Fig.5.4 and Table 5.3. Here, we denote by P the porous medium equation, Example 2 , and by A the porous medium equation with absorption, Example 3. Note that in order to make the comparison we 24 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. modify Example 2 and run it with the same the number of the spatial steps, M — 20,40,80 and equal exponents, m = 2. Table 5.3: Comparison of computational costs (per temporal step) No. of Kflops CPU time (msec) M P A P A 2 0 23.5 68.3 12.9 25.4 40 96.6 284.5 58.3 128.0 80 391.5 1,147.8 362.1 735.1 In Fig.5.5 we plot the time evolution of the maximum height, H of the com puted solution, u of (5.3) in standard and log-log scales. Note th at the log of the solution is decreasing rapidly in the vicinity of the extinction time, T*. Addition ally, in Fig.5.6 we plot traces of the computed solution, u on the xy-plane at 4 different points in its time evolution. The goal here is to accurately compute the extinction time, T*, the time at which the solution becomes identically zero, u(x,y, T*) = 0. Thus, one must accurately compute both the solution, u, and the time T*. The latter task imposes additional restriction: the temporal step k must be smaller or equal the required precision in determining the value of T*. This point is well illustrated by the following Table 5.4: 25 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 5.4: T* for fixed h=.04 k .01 .001 .0001 .00001 A 6.25 .625 .0625 .00625 r j - i* .26 .263 .2634 .26343 It is clear th at the number of significant digits in the computed value of T* is limited by the — log1 0 k. It doesn’t mean though that it is sufficient to take small values of k (the temporal step size) as of course the spatial discretization error also plays an important role, see Table 5.5. Table 5.5: T* for fixed k = le-5 h .08 .04 A .0015625 .00625 T* .26360 .26343 In these experiments with k = le-5 we limited the refining of the spatial grid because of large computational costs. To remedy this we designed a modified algorithm. From the point of view of accuracy one could conduct the experiments with much larger values of A (see Table 5.2) than those in Table 5.5 if it was not for the resolution requirement that k be sufficiently small at the extinction time T*. The time evolution (in the log-log scale) of the maximum height, H of the computed solution, u, see Fig.5.5, shows an abrupt downward turn at some point, 26 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. T° in the vicinity of the extinction time, T*. The modified algorithm employs a much larger time step, k° for t < T°, and a small k, corresponding to the required resolution, afterwards. In the experiments reported below, see Table 5.6, we have chosen an ad hoc value k° = lO O fc, and the program was adaptively changing the temporal step size a t a point when the slope, p of the maximum height of u became p < p°. Again, ad hoc we have chosen p° = — 20. The reported A in Table 5.6 is that of A 0 = p , i.e., th at for t < T°. Table 5.6: Modified algorithm, T* for fixed k = le-5 after T° M 30 60 1 2 0 240 h .08 .04 . 0 2 . 0 1 A° .15625 .625 2.5 1 0 T* .26360 .26343 .26337 .26336 For a fixed spatial step size, h the total computational cost is proportional to the number of temporal steps, -p for the original algorithm and ^T° for the modified one. For example, for h = .04 the computed value of T° turned out to be T° = .241. In this case the reduction of the number of temporal steps was from 26,343 to 2,484, or more than tenfold. We have not made an attem pt to optimize the speed up. 27 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Although we do not know the exact value of the extinction time we can obtain a very rough estimate of the power, p of the asymptotic error model using the Richardson extrapolation method on the values in the last table. We obtain: 2P = II3 4 3 - 2 II3 7' = ^6- We can thus conclude th at the accuracy of this computation is of about 0 ( h 15), which is consistent with the experiments in Section 5.1. 28 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Computational cost N, p=3.60 O 1 0 ' L, p=3.60 3 1 0 " 1 S, p=1.96 ,-2 2 10' N, p=3.06 = > 1 0 ' Number of grid points, M Figure 5.1: Comparison of the efficiency Comparison of the efficiency (measured in Mflops and CPU time (msec)) for the Newton (N) method (two iterations), standard linearization (L) method and our splitting (S) method as a function of the number of grid points, M for Example 1. Here, p is the power of the cost’s growth model: cost — 0 (M P). 29 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Accuracy i_ 10 o t T 3 , o - 1 0" L, p=-3.27 N, p=-1.98 S, p=-1,92 Example 1 1 0 1 0 10" o 1nr6 3 L, p=-1.48 N, p=-1.47 S, p=-1.46 Example 2 10" 1 0 Number of grid points, M 1 0 Figure 5.2: Comparison of the accuracy Comparison of the accuracy for the Newton (N) method (two iterations), standard linearization (L) method and our splitting (S) method as a function of the number of grid points, M for A = 1. Here, p is the power of the error’s growth model: error = 0 ( M P). 30 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Accuracy vs. efficiency Example 1 .-7 L, p=-1.05 ^ * ' * ^ N’ P="0-64 Example 2 S, p=-0.57 L, p=-0.48 io - 7 CPU time (msec) Figure 5.3: Comparison of the accuracy vs. efficiency Comparison of the accuracy (L 2 -error) vs. efficiency (CPU time) for the Newton (N) method (two iterations), standard linearization (L) method and our splitting (S) method. 31 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Efficiency; Example 3 vs. modified Example 2 & ° 3 g 1 0 1 1 0 s 1 0 A, p=2.04 P, p=2.03 10 ' 1 0 1 0 ° 8 m E, | 1 0 2 •H D 0 . o 10 ’ 10 ' A, p=2.43 P, p-2.40 Number of spatial steps, M 1 0 Figure 5.4: Comparison between A and P Comparison between the porous medium equation with strong absorption (A) and the porous medium equation (P) in terms of the number of Kflops (Top) and CPU time (Bottom), as a function of the number of spatial steps, M for A = 1. Here, p is the power of the cost’s growth model: error = 0 ( M P). 32 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Extinction phenomenon 0.8 ® 0.6 0.4 0.2 0.05 0.15 0.25 0.2 0.3 0.35 Timet Iog(t) Figure 5.5: Time evolution of the maximum height H of the numerical solution, u with the spatial step size h — .04 and temporal step size k = .0001. Bottom plot is the log-log plot of the top one. 33 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.6: Finite extinction phenomenon Finite extinction phenomenon of the numerical solution for Example 3 with the spatial step size h = .04 and temporal step size k = .0001. Here, H is the maximum height of the numerical solution at the time t. 34 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 6 Concluding remarks For a class of nonlinear parabolic PDEs in 2D rectangular domains we have con structed an operator splitting algorithm of optimal efficiency. We have verified experimentally for the porous medium equation that the computational cost of our scheme is 0 (1 ) flops per unknown per temporal step while the accuracy remains the same as for two Newton’s iterations. In the presence of an additional zero order term (strong absorption term) the asymptotic efficiency remains unchanged, 0 (1 ), with the leading constant only twice larger. We have modified the algorithm to properly change the temporal step size. This allows computing the extinction times extremely accurately and with significant computational savings. 35 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 7 Appendix: Linear stability analysis Strictly speaking the linear Fourier analysis applied to nonlinear equations to show the stability of the scheme cannot be rigorously justified. Nevertheless it has been found to be effective in practice. For an example of such analysis applied to Korteweg-de Vries equation see [5]. First of all, let us assume that the solution function u is bounded in the given spatio-temporal region and so let v = max |um-1|. 36 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Substituting it in eqn.(3.3) the corresponding linear equation to which we apply the von Neumann analysis becomes yyifc \ I 2 v(D-\-xD— x 3 ~ D+yD— y) = kv[D+xD— x 3 ~ D-^yD— y)unj i.e., mA , ^ p ,q ,n + 1 — ^ V \W p —l ,q ,n + l “ I - 'W p+ x,q,n+ l 3“ ^ p , q —l ,n + l 3“ ^ p ,q + l ,n + l 4 :tV p ,q ,n + l) Av(Up— x,q,n+l 3“ U p^.x,q,n+1 3“ Uptq— X,n+ 1 3” M p ,g -f-l,ra+ l ^'Up,q,n+l) • (7-1) Let (7.2) where (3 — Kxh, 7 = nyh and (3 , 7 e [— 7r, 7r], Substituting (7.2) into (7.1) and then dividing by £nei/ 3pe * 79, we obtain (£ - 1) - — ^~v(e~if3 + eip + e - *7 + ei 7 - 4)(£ - 1) -*P JP - 1 - o— * 7 _ L e* 7 Av(e~lP 3- el(i 3- e ^ 7 + e1 7 - 4), which can be written as 1 — T ^ - v { 2 cos j3 + 2 cos7 — 4) )(£ — 1) = A w (2cos/? + 2 cos 7 — 4). Hence, — 4Awfsin2( |) + sin2(7)N ] £ - 1 = - —s 1 4 - 2mAw^sin2( |) + sin2( |) 37 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Therefore, the amplification factor £ is 1 + 2(m - 2)Au^sin2 (f) + sin2(^)^ 1 + 2raAu^sin2 ( f ) + sin2(^)^ Since m > 2, it is clear that 0 < £ < 1 for all positive A and all (3, 7 . Thus the Crank-Nicolson method for the porous medium equations is unconditionally linearly stable. 38 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Part II A pproxim ate Inertial M anifold w ith Perturbation Technique Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 8 Introduction The aim of this paper is to present an explicit form for the construction of ap proximate inertial manifolds (AIM) for a class of nonlinear dissipative partial differential equations by using a perturbation technique. Partial differential equations can be treated as infinite-dimensional dynami cal systems on a suitable Hilbert space H [15]. For certain classes of dissipa tive dynamical systems (including reaction-diffusion equations, pattern formation equations, and fluid dynamic equations such as the Navier-Stokes equations) it is possible to ensure that the asymptotic behavior of the system is confined to a finite-dimensional subspace. This approach is usually referred to as inertial manifold theory [11],[12]. An inertial manifold is a finite-dimensional, exponentially attracting, positively invariant Lipschitz manifold. It allows for the reduction of the infinite dimensional 40 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. dynamics to a finite dimensional system of ordinary differential equations referred to as inertial form which shares all the long term dynamics of the original problem. Existence results of inertial manifolds can be found in Foias, Sell, and Temam [12] and Mallet-Paret and Sell [22] (see the references therein). However, there are still many dissipative partial differential equations for which the existence of inertial manifolds is not known. In the construction of an AIM, one seeks to relax the restrictive requirements of invariance, and exponential attraction, so that they hold only up to a neigh borhood of the manifold [23]. From a practical point of view then, AIMs can be applied regardless of the existence of an actual inertial manifold. Even for equations possessing inertial manifolds, AIM can be useful, since they provide approximations of the attractor by finite dimensional smooth manifolds which are explicitly defined. We consider here a class of nonlinear evolutionary equations of the type dn — j£ + Au F(u) — 0, (8-1) where A is a suitable linear, unbounded, self-adjoint and positive operator on a Hilbert space H, while F is a nonlinear operator satisfying some conditions. Here we assume that an inertial manifold exists for the class of equations (8.1) (see [12]). Instead, we approximate an inertial manifold as the graph of a smooth function < 3 > : PH — ► QH, where P is an orthogonal projection with 41 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. finite dimensional range in H and Q = I — P. The function $ is approximated as an explicit form modified by the approximate solution of Qu where u is the solution of equation (8.1). This approximate solution of Qu is obtained by using perturbation technique. This paper is organized as follows: In Chapter 9, we present some background material about dissipative dynamical systems. In Chapter 10, we present a finite series for an approximate solution of Qu by using a perturbation technique. And we construct Approximate Inertial Manifolds. Numerical results are presented in Chapter 11, while Chapter 12 contains the concluding remarks. 42 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 9 Prelim inaries Let ft denote a regular bounded set of RF (1 < n < 4). Assume we are given an infinite dimensional Hilbert space H = L2(ft) with an inner product (•, •) and the corresponding norm | ■ |. The nonlinear evolutionary equation which we will study has the form H it ~ + Au + F(u) = 0, (9.1) where A, with the appropriate boundary condition (Dirichlet, Neumann, periodic), is a suitable linear, unbounded, self-adjoint and positive operator on H with dense domain D(A) C H, while the nonlinear operator F is at least a continuously differentiable function from R into R . Furthermore, if F satisfies the growth conditions Ci|s|^ — C 2 5: F(s)s < c3\s\k — c2, Vs e R, with k> 2, (9.2) 43 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. F'(s) > — c4, Vs E i?, (9.3) where the q ’s are positive constants, then eq.(9.1) has a unique solution for any positive time given an initial condition u0 E H, see [21]. Suppose further that A rx is compact. As a result, H has an orthonormal basis {e-j}JL1 consisting of eigenfunctions of the operator A, with corresponding eigenvalues Xj, i.e., Acj A j Cj, j 1,2, •" ■ , where the eigenvalues satisfy 0 < Ai < A 2 < • ■ •, Xj — ► oo as j — * co. Furthermore, in terms of this basis the operator A can be represented by OO Au = ^2X j (u,ej)ej . (9.4) i=i Since we have a representation of A by (9.4), we can use this to define its fractional powers by the natural expression OO Aau = Y l X < j ( u’ej)ei- (9-5) l=i The domain of Aa in H, D(Aa), is OO OO D(Aa) = {u :u = Y dcjej ,Y^ \cj\2\ f < oo}. l= i l= i We can make D(Aa ) a Hilbert space by using the inner product (u,v)D(A*) = {Aau,Aav), 44 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. which gives rise to a corresponding norm IMId(A“) = \\Aau\\. In the sequel, we denote by || • || the norm of D(A1^2) and | • \p the norm of LP(Q), 1 < p < oo (| • |2 = | ■ I), while the norm of any other Banach space B is denoted by || • ||s. W ith the assumption of the existence of ^-dimensional inertial manifold A4, we let Pjv denote the orthogonal projection onto the span of the first N eigenfunc tions {ei, e2, ■ • •, eN} of the operator A and let QN = I — Pn- We write simply P = Pn and Q — Qn - If u = u(t) is a solution of (9.1) we define p = p(t) and q = q(t) by p — Pu and q = Qu. By projection of (9.1) on P H and QH, we find that p and q are solutions of the coupled system of equations ^ + Ap + PF(p + q) = 0, (9.6) -Jj; + Aq + QF(p + q) = 0, (9.7) where Ap = APu = PAu and Aq = AQu = QAu. In the following, we will use the following useful properties (see [17],[18]) : for all a > 0, e — 1 / Xn+i, i,j € N, \Aaq\2 < e\Aa+1/2q\2, Vg G QD(Aa+1/2), (9.8) 45 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. \q\ < cq0e, \ft > to, dt — C qi€, Vi > ti, ^ ^'uj• V i ^ t j , (9.9) (9.10) (9.11) where each tim e i0,ij, tj depends on the data (f2, A, F ) and on R when an initial value u(0) = uq satisfies |uo| < R. If q G QD(A), (9.8) gives \q\2 < e\A1/2q\2 \A1/2q\2 < e\Aq\2 where each inequality is for a = 0 and a = 1/2, respectively. Therefore, we obtain M < AM \, Vg G QD(A). (9.12) We assume that, for every u0 G D(A), the initial value problem (9.1), with initial value tt(0) = uq, has a unique solution S(t)uo defined for all i > 0 and S(t)uo G D(A) for all i > 0. The solution operator S(t) : H — > H has the following properties: 5(0) = /, S(t)S{s) = S(s)S(t) = S(s + i), S(t)u0 is continuous in u0 and t. The solution S(t)u0 is uniformly bounded in time in H, and if u0 G D(A) then S(t)uo is uniformly bounded in D(A). See [20]. 46 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. As in Foias et al. [12], we define an inertial manifold A4 for (9.1) to be a subset of H satisfying the following properties: (i) Ad is a finite dimensional Lipschitz manifold, (ii) M is invariant, i.e., S{t)M. C M , for all t > 0, (iii) M. attracts exponentially all solutions of (9.1), i.e., dist(S(t)u0,M ) — > 0 as t — > oo (9.13) for every uQ £ H and the rate of decay in (9.13) is exponential, uniformly for U q in bounded sets in H. Property (iii) implies that an inertial manifold must contain the global attractor. 47 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 10 A pproxim ate inertial manifold w ith perturbation technique In this chapter we will be looking for an approximate inertial manifold Afapp which is constructed as the graph of a smooth function. This function < f > app is defined by introducing a suitable approximation of the result obtained by using a perturbation technique. 10.1 P ertu rb ation technique Consider the split system of equations p + Ap + PF(p + q) = 0, (10.1) q Aq -f QF(p + q) = 0, (10.2) 48 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where the over dot denotes the derivative with respect to time. Here, the objective is to solve the equation for q (10.2). This can be accom plished by the perturbation technique. Let us define a small positive number e = 1/Ajv+1 and the operator eA on QD(A ) as the linear operator L. Then the linear operator L is constant order with respect to e, th at is, Lq — eAq ~ 0(1). This is because LQ = L (q,ei)ei i=N+1 oo = eA (< ? > ei)ei i=N+l 1 0 0 = V H \(q , ei)ei A iV +i oo ^ i(Qi & i)& i> A j > 1 * = A T + 1 and OO OO Lq\ = X) A i(q,ei)ei > i= JV + l i=N+1 kl- Multiplying eq.(10.2) by e we obtain a perturbed equation eg + Lq + eQF(p + q) = 0, (10.3) where operators L and Q are linear and the overdot denotes a derivative with respect to the independent variable t. Here, we apply the perturbation technique for the perturbed equation (10.3) for q. This perturbed equation (10.3) is a regular perturbation problem because 49 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. its solution , say q(x,t,e), converges as e — ► 0, (i.e., iV — > oo uniformly with respect to the independent variables x and t) to the zero function which is the solution of the limiting problem Lq — 0. So the idea is to expand the solution q(x, t, e) in a power series in e, q(x, t, e) = $o(x, t ) + e0i(z, t) + e202(^, t) + • ■ ■ (10.4) and use a partial sum of the infinite series (10.4) as an approximation of the exact q(x, t, e), as k ?(£>e) = (10.5) i— 0 where the spatial index is suppressed to simplify the notation. (Hereafter we suppress the spatial index, i.e., q(t,e) = q(x,t,e) and = 4>i(x,t).) And the unknown functions <pi(t) are solved recursively, i.e., in the order ( j> 0, - , < p k. By substituting the finite series expansion (10.5) into our perturbed equation (10.3) and doing a derivative with respect to t we obtain + ecj) i + • • • + + L{(j)Q + ecf>i + ■ ■ • + e f e < ^ f c ) + eQF(p + 0o + e0i + • • • + 4 > k) — 0, (10.6) where (f> i s are functions of x and t. Collecting coefficients of equal powers of e in eq.(10.6) gives L(j)o + e[0o 4“ L(f> i + QF(p -f- 0o)] ~f ■ * * = (1QA) 50 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where only terms up to O(e) have been retained. Here, we assume that QF(p+(j>o) is constant order with respect to e, i.e., Q F (p + 0 o) = O (l). Equating the coefficient of e° to zero in (10.7) yields 0o = 0, 00 = 0 . Therefore, equating the coefficient of e1 to zero in (10.7) we obtain L(f> i + QF(p) = 0, (10.8) i.e., 0! = —L QF(j>) Moreover, our assumption (10.8) can be rewritten as (10.9) QF(p) = 0(1). (10.10) To obtain the higher-order terms up to 0(ek) in the power series solution (10.5), the nonlinear term F satisfies the condition for Taylor’s theorem to use Taylor’s formula for QF (p + £i= i e*0i) about p. If F : R —> R be of class Ck, we can consider Taylor’s formula for F (p + Ya=i e*0i): k £ i=1 r k F(p + J2e% ) - F{p) + F'{p) £ e V i + ’ • • + - ^ { p ) b=i J J■ r k Li=l ' k i=1 (10.11) 51 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where the prime denotes the derivative with respect to p and the value of p lies between p and p + X)i=i e*0*- We now expand (10.11) in the form of a power series in e. Since we are only interested in terms up to 0(e k), we only keep terms to 0(ek) inside the brackets in (10.11), that is, E e* 0* t=i 'k-j+l E *=1 + 0(ek+1), for j = 1,2,” • ,k. (10.12) Thus by using the relationship (10.12) we can rewrite Taylor’s formula (10.11) as follows: F ( p + E £i« i=1 ■i=1 1 £01 + 0(ek+1) k-j+l E e V i F0 + eFi + • • • + H h ekF^ + 0(ek+1), (10.13) F0 = F(p), for jf = 1,2, • ■ •, fc - 1, = F ( p ) ^ + i r ( p ) ^ . 0*1 0 * 2 * l + * 2 = j 0 ) 0 * 1 0* 2 0 ! * l + * 2 + * 3 = V *3 + E 0*i 0*2''' 0v _ * 1 + * 2 H ---h * 7 = j (10.14) for 0 < ii, * 2 , • • ■ , ij < j < k, and = F(p)0f c + -F"(p) + ^ < 3 ) w E l 0*1 0*2 0*3 * l + * 2 + * 3 = f c 52 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. For instance, first three terms Fi, F-i and F 3 are f t = F\p)4> 1 , F2 ■ F'(p)<k + l^F"(p) 53 ^ * 1 & 11+ 12= 2 * 2 f t = ft(p )* i + %F"(p) ii+i2=3 + ^ r (3)(p) &h4’te4> h ii+h+is— S F'{p)'h + | f t » [0i02 + + | ^ (3)(p) [0i0i0i] ft'(p)03 + ft"(p)0102 + ^yft’t3)tP)0 Since the operator Q is linear, Taylor’s formula for F defined by (10.13) can be used directly to induce a Taylor expansion for QF (p + e< pi + ■ ■ • + ef c < /> * ,) about p ( p + b QF Ip + e ’&J = Oft + f<3ft + 0Q ft + ••■ + e Oft + 0(£*+1). (10.15) Using (10.15), we can rewrite eq.(10.6) as e24 > i + es4 >2 H ------- 1 - ek+l(f)k + eL(j) 1 + e2L (j> 2 + • ■ • + ^ L < f> k + eQF0 + e2QF1 + --- + efe+ 1 <3Ff c = 0 . (10.16) Only if each QFi, i = 0,1,2, • • ■ , k, is constant order with respect to e, i.e., QFi — 0(1) i = 0,1,2, • • ■ , fc , (10.17) 53 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. collecting coefficients of equal powers of e in eq.(10.16) gives L(f> i + Q F q + e 0i -F L (f> 2 + QF\ 0fc-1 + L (j> k + 0, where terms up to 0(ek) have been retained. Equating the coefficient of each power of e to zero, we obtain the explicit for mula for each 0 ,- as follows: 0! = - L ~ xQF q = —L~1QF(p), = -L-'iQ F, + fc) - -L "1 (<?*»* + + ), 0fc = — L 1(QFk-i + 0fc-i) = - i - 1(O F '(p )4 _ , + iQ ii" ( !,) X ) [* .* » ] + ■ ii+i2=fc— 1 4 ' 3 f |3) (p ) E I* . + • • • + 7 r 3 W « F<‘ “ 1 ) (p ) ^ _1 + F - i ^ 11+42+13=^— 1 1 h Thus we can get approximate solutions q(x,t,ey s, say <& (£), as finite sums of 0i(t)’s> i.e., k Qk(t) = ' £ e i< j) i(t), k = 1 ,2 ,---. (10.18) 4=1 Our approximate inertial manifolds Affe are closely related to these approxima tions qFs. 54 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10.2 C on stru ction of approxim ate inertial m an ifolds The aim of this section is to present the construction of the sequence (Mk)keN of approximate inertial manifolds in detail. That is, the approximate inertial manifolds M k are constructed as the graphs of the functions : PD{A) — > ■ QD(A) which will be suitable approximations of the previous perturbation results under the assumption th at the nonlinear operator F has the Taylor expansion given by (10.11). Considering k = 1 in the equation (10.18) gives the first manifold M i when F is a differentiable function satisfying the condition (10.17). In practice, we assume that F is a Cl function which satisfies F{p) = 0 ( 1) with |F '| < K n (10.19) where a positive constant K n is constant order with respect to e. The approximate solution q(x,t,e) of (10.2), say qi(t), is qi(t) = e< l> i{t) — -eL~lQF(p{t)) — —A~1QF(p(t)). (10.20) Since q\ is a function of t through p(t), the related function dq, the function of p, can be defined by qi itself, i.e., $x(p) = - A ~ lQF(p). (10.21) 55 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The graph of the function $ 1 is the approximate inertial manifold defined in C. Foias, 0 . Manley and R. Temam [16] in the context of the Navier-Stokes equation and subsequently used in many other articles, both theoretical and numerical, (see [17], [18], [24]) According to [16],[17],[18], all solutions of eq.(9.1) are attracted by a thin neighborhood of M i as shown in the following theorem. T h eo rem 1 For t sufficiently large, t > t*, any orbit of (9.1) remains at a dis tance in H of M i bounded by ce2; c is an appropriate constant which depends on the data and t* depends on the data and on R when the initial function u0 satisfies | « o | < R - Proof. Let u = p + q be an orbit of the given nonlinear equation (9.1). For every t > 0, because of the relation (10.21), p(t) + (p(t)) lies in M i which is the graph of $ i. Therefore, we have dist(«(t),M i) = iniv€Ml\u(t) - v\ < \q(t) - #i(p(i))|. By the definitions of q and 'fq in (10.2) and (10.21) respectively, we can evaluate the norm as follows since \Q\ < 1: A{q{t)-$i{p{t))) = ~[QF(p + q ) - QF(p) + q] A{q{t) - $ i( p ( i) ) ) 1 < \Q(F(p + q) - F(p))\ + |^| < \F(p + q)-F (p)\ + \q\ < K u\q\ + 1 ? | 56 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where the last inequality comes from the mean value theorem using (10.19). Using (9.9,9.10) for j — 0,1, we obtain from (10.22) |A(q(t) - (t)) | < {Ku cqo + cqi)e L = Cje, for t > m ax(t0,U)- (10.23) Since q — G QD(A), using a known property (9.12), we can conclude that |< z (i)-* i(t)| < e\A(q(t)-Mt))\ < c1e2. (10.24) Q.E.D. Next we construct another manifold M 2 if the nonlinear operator F is a C2 function which satisfies the sufficient conditions of (10.17): F(p) = 0{ 1), w ith |F '| < K2i and |F"| < K 22, (10.25) where positive constants K2i and K 22 are constant order with respect to e. The orbits of (9.1) become closer to M 2 than to the manifold M \ by a factor e. This approximate inertial manifold M 2 is the graph of the function < 3 > 2(p) which is obtained by considering an improved approximation of q. In (10.18), q2(t) is rewritten as follows: Q 2(t) = e0i(f) + e2(f> 2(t) = - e l r 1 (QF(p(t)) + QF'{p(t))e<t>i(t) + e^. (t)^ - - A ~ x (QF(p(t)) + QF’(p(t))qi(t) + qxit^j (10.26) 57 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where the overdot and the prime denote the derivative with respect to time and p, respectively. Here, at the second step k — 2, we will update q\ in (10.26) using $ !, the best approximation of q that we have obtained. That is, the nonlinear term F(p) in (10.21) is replaced by F(p + $i(p)). And we can rewrite qi(t) as follows: Qi(t) = ~ ^ A - lQF{p + ^ 1{p)) = - A - 1g F '(p + $ 1(p))(p + $ 1(p)) (10.27) where the second equality (10.27) comes from the chain rule and the overdot denotes the derivative with respect to time. In (10.27) p + $i(p) can be approximated by p and p given by (10.1) is better approximated by p = -A p - PF(p + (p)). Hence, we can derive a better value for q\ in (10.26) and denote it by <l>i(p) = A ^Q F'lp + §i(p))^Ap + PF(p + §i(p))^. (10.28) Thus the function $ 2 is defined by $ 2(p) = - A - 1( q F{p) + Q F’(p)$1(p) + $ !(p )), (10.29) where is given by (10.28). 58 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The manifold M 2, the graph of < J > 2 (p), provides a higher order approximation to the orbits than M i. This is stated in the following theorem. T h eo rem 2 For t sufficiently large, t > t*, any orbit of (9.1) where (9.2) and (9.3) hold remains at a distance in H of M 2 bounded by ce3; c is an appropriate constant which depends on the data and t* depends on the data and on R when the initial function uq satisfies |w o| < R. Proof. Let u = p + q be an orbit of (9.1) with a boundary condition. As in the proof of Theorem 1, in order to evaluate the distance in H to jM2, it suffices to estimate the norm of < 3 > 2(p(i)) ~ ?(£)■ Using the definitions of q and $ 2 in (10.2) and (10.29) respectively, we have A ($ 2 — q) = QF(p + q)~ QF(p) - QF'{p)§ 1 + q - $ x. So \A ( $ 2 -q ) \ < IQF(p + q)~ QF(p) - QF'(p)$i \ + \q - . (10.30) Since \Q\ < 1 we get \QF(p + q ) - Q F ( p ) - Q F ’(p)$1\ = \Q (F(p + q) - F(p) - F'(p)^)j < |F(p + q) - F(p) - F'(p)$i| . (10.31) Consider the Taylor approximation of F(p + q) in p as follows: F(p + q)= F(p) + F\p)q + \F"{p)q\ (10.32) 59 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where p is between p and p + q. Substituting the Taylor approximation (10.32) into (10.31) gives the following inequality: |F(p + q)~ F(p) - F 'W * ,! < |F'(p) (? - $ ,)| + i | f "(p)?2|. Using (9.9), (10.24), and (10.25), we have ! * » (< f - * i)| + l \ F " ( p ) q 2\ < K n I? - $ i| + ^ p |? 2| < C\K2l£1 + 2C ,0-^22e2 Let ~ 2 = c2ie . Thus the first member on the right in (10.30) satisfies that |QF(p + q) - QF(p) - QF'(p)$i\ < c2ie2. (10.33) For the second member in (10.30), we will estimate A(4>i — q). The time derivative of eq.(10.2) and the definition of <fq in (10.28) give A ( $ i - q ) = q + QF’(p + q)(p + q) + QF'(p + $ i) [Ap + PF(p + = q + QF'(p + q)q + QF'(p + $ x) (p + Ap + PF(p + < f> i) +Q (V (p + q)~ F \p + < F x))p. (10.34) Now we need to estimate the different four terms in the right-hand side of (10.34) for large time t. For the first two terms we can use the property regarding time derivatives (9.10) and the conditions of F (10.25). For large t \q\ < cq2e (10.35) 60 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. | QF'{p + q)q\ < | F'{p + q)q\ < K2i\q\ < K2lcqle. (10.36) For the third term in the right-hand side of (10.34), we use (10.1) and substitute p -I- Ap as —PF(p + q) to give, QF\p + 4>i) (p + Ap + PF(p + 4b QF'{p + $ j) PF(p + $ 0 - PF(p + q < K 2i\F(p+ $ i) - F(p + q)\ < K 2 21\ ^ - q \ < h K l e 2, (10.37) (10.38) (10.39) where the first two inequalities (10.37,10.38) come from the mean value theorem, and equation (10.24) from Theorem 1 is used for the last inequality, (10.39). The fourth term in (10.34) can be approximated as follows: \Q (F>(p + q )-F '(p + < £> ,))p\ < \(F/(p + q ) - F l(p + ^ 1))p\ < K22\(q - $i)p| < K22\q - UIpU where the last inequality is given by Holder inequality. Since Q is a bounded smooth domain in R4, we can use the Sobolev embedding theorems, i.e., IO(P(P + 9)--P '(p + 4'1)}p| < K22\ q ~ ^ \ i IpU 61 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where C is a positive constant. Moreover, the H 2 norm is equivalent to the norm \A • | (see [21]). Therefore, we can obtain the following inequality \Q(F'(p + q)-F'(p + $ 1))p\ < C K 2 2\\q — $ i||h 2IIpI|h2 < C"X22| ^ % - $ i ) P p | (10.40) where C' is a positive constant. Here, \Ap\ — \APu\ = \PAu\ < \Au\ < cui (10.41) where the last inequality comes from the bounded property of u (9.11). Hence, using the result of Theorem 1 (10.23) and (10.41), we infer from (10.40) \Q{F'(p + q)-F'(p + $ 1))p\ < C ' K ^ A i q - ^ W A p l < C ' c ^ K ^ e . (10.42) Thus, from (10.35), (10.36), (10.39), and (10.42), we obtain that for large time t |H ($i — q) | < (cq 2 + K 21cql + diK^e + C'c\cuiK 2 2)t Let = c22e. Since 4>i — q € QD(A), using (9.12), the second member in (10.30) is approx imated as |$ i — q\ < e|H($! — q)| < c22e2. (10.43) 62 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Finally, we conclude from (10.33) and (10.43) that | A (4> 2 - q) | < ( 6 2 1 + c22)e2 =* c2 e2 (10.44) Hence, because of (9.12), we have I $2 ~ q \ < c2e3 (10.45) Q.E.D. Now we aim to generalize the idea to construct A4k, (k > 3), recursively under the assumption th at the nonlinear operator F in the evolution equation (9.1) is a Ck function which satisfies the corresponding sufficient conditions of (10.17). We start by rewriting in (10.18) < ? f c eL 1 ( QFq + eQF\ + - • • + ef c iQF^— x + tl4 > i + • • ■ ef c l< f> k - 1 k-il + V,QF"(p) E 'V . L i= 1 Jfc-2 + QF^k- 1 \p) e0i fc-i + Q F '(p)q^ + + 1 + f t- i )+ 0 (el+1). (10.46) We need to obtain better approximation for <jk-i by using the known best approximation 4h--i of q up to this stage. That is, we define dW i by §k-i(p) = A QF'{p + 4 > f c -i(p)) Ap + PF(p + (10.47) 63 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Thus we can construct the manifold M k with the equation q = $ k(p), : PD(A) — > QD(A). Here, &k(p) can be defined recursively as (10.48) where $k-i(p) is defined in (10.47). In the following chapter we will consider first two approximate inertial man ifolds, $ 1 and $ 2 ) only. That is, we will show that our approximate inertial manifold with two terms, < 3 > 2, gives more accurate results than other nonlinear Galerkin methods through numerical experiments. 64 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 11 Num erical experim ents In this chapter we present some numerical experiments to support the results of the previous chapter. We will show that the higher order approximate inertial manifold $ 2(p) gives more accurate results than < E q (p) and the Euler-Galerkin method. We will compare accuracy of different methods. We investigate two numerical examples of the reaction diffusion equation subject to homogeneous Dirichlet boundary conditions with different nonlinear terms: 1. non-polynomial nonlinearity and 2. polynomial nonlinearity. For the purpose of computational demonstrations, we consider the reaction diffusion equation with a general nonlinear operator ut — vuxx + F{u) = 0, x e (0,7r), t > 0 subject to homogeneous Dirichlet boundary conditions, u(0, t) = u(7 T , t) — 0, t > 0 (11.2) 65 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where the diffusivity v is a, positive constant. We note th at the 1-D Laplacian operator — dxx with the Dirichlet boundary condition (11.2) is a linear, unbounded, self-adjoint and positive operator possessing a compact inverse. We apply semi-finite difference method to obtain a system of ordinary differen tial equations. A continuous function u = u(x, t ) is discretized on an equidistant spatial grid Xj = jh, 1 < j < M, where h = is the step size in spatial direction. Discretizing the boundary valued problem (11.1), (11.2) in the spatial variable it can be reduced to a system of ordinary differential equations where M is the number of unknowns and indices 0 and M + l denote the boundary Let us consider a column vector u : u(t) = (ui, • ■ • ,UM)T whose components are real valued functions of t. In matrix form eq. (11.3) can be w ritten as u0 — 0 — Um+i duj v . n . — (uj+ 1 — 2Uj + Uj- 1 ) + F(v,j) — 0, for j = 1, • • ■ , M (11.3) values. — + Au + F(u) = 0, (11.4) where A is a tridiagonal matrix of order M : A = tridiag /l2’2 h2’ h2 66 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. and F(u) is a column vector of the length M with the contribution from the homogeneous boundary condition (11.2) : / \ F M - F(u2) F(u) V / \ F(u i) F(um)- (11.5) / F(uM- i) F{um) - 1 This ODE system (11.4) can be solved using spectral methods by projecting the set of equations in (11.4) onto the eigenbasis of A (see [15]). Let E be the matrices whose columns are the eigenvectors of A, ordered in terms of increasing eigenvalues 0 < Ai < A 2 < • • • < Am - Let U = E~xu, i.e., u = EU. Substituting EU into u and taking E~l in eq. (11.4), we obtain dU dt + AU + G(U) — 0, (11.6) where A = E -'A E diag(Ai, A2, • • •, Xm ) Let G(U) - E~1F{EU) (G\, G2, • • •, Gm )t ■ 67 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Eq. (11.6) is the starting point for the computational analysis of approximate inertial manifolds since it is suitable for decomposition into dominant and enslaved modes. We assume th at the dominant p-mode is Np. Then we can decompose eq. (11.6) into dominant and enslaved modes, i.e., < ^ + A pp + Gp(p, q) = 0 (11.7) ^ + A qq + Gq(p,q) = 0 (11.8) where A p = diag(Ai,A2,---,A Np), A q — diag(Ajvp+i, A jv p+2, • • ■ , Am), Gp = (Gi, G2, • • •, Gnp)t , Gq — ( G jV p + i, G n p+ 2, ■ • ■, G m ) t , V = (U1,U2,---,UNp)Th= {pi,p2,---,PNp)T, q — (f/jVp+i, U n p+ 2, • • •, U m )T = = (#1 , Q 2 , • • ■ , qNq)T, (p,g) = (U1,U2,--- ,UM)T = (Pi,-" ,PNp,Qir-- ,QNq)T- The construction of the approximate inertial manifold through the result of the previous chapter can be given by $1 (p) = - ^ G q i p ) (11.9) $ 2 (p) = - A g 1(Gq(p) + G'q(p)^1(p) + ^i(p)j (11.10) $ 1 (p) = Ag1^ ( p ,$ i( p ) ) ( A p p + G p(p,$i(p)))^ (11.11) 68 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where the overdot and the prime denote the derivative with respect to time and p , respectively. The approximate inertial manifold (11.9), q = $i(p), is one of the simplest ver sions of the nonlinear Galerkin methods which was introduced by Foias, Manley, and Temam [16]. Another simple nonlinear Galerkin method is the Euler-Galerkin AIM introduced by Foias, Sell, and Titi [13] and subsequently used in many stud ies. (see [14],[15]) The Euler-Galerkin AIM (say, $>(p)) is given by q = $(p) = — r(I 4- TAq) ~ 1G q(p) (11.12) where I is the identity matrix and the time scale r is equal to l/A^rp+i for the sake of simplicity, (see [14]) The standard Galerkin approximation of (11.6) using Np modes can be ex pressed by setting q = 0 in (11.7), i.e., + App + Gp(p) = 0. (11.13) We will compare the accuracy of these different methods using the solution obtained by the semi-finite difference method with sufficiently fine discretization as our exact solution. We will consider four different methods : 1. the standard Galerkin approximation (denoted by SG in Tables and Fig ures), 2. the Euler-Galerkin approximation of equation (11.12) (denoted by EG), 69 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3. the perturbation method with the first term only (11.9) which was intro duced by Foias, Manley, and Temam [16] (denoted by PI), 4. the modified perturbation method with two terms (11.10) (denoted by P2). Three versions of the Galerkin’s methods for solving (11.1) with some initial and boundary conditions were previously investigated : SG, EG, P I. We compare them with our AIM obtained by the modified perturbation method using two terms ffq and $ 2 in (11.10). All our algorithms are implemented in M atlab and the numerical simulations were performed with the ODESUITE stiff stable integrator ODE15s. We investigate two numerical examples of the reaction diffusion equation sub ject to homogeneous Dirichlet boundary conditions for positive constants a, b: • Exam ple 1. ut - \ u xx - e«/d+““) u(x, 0) lt(0 , t) — u(7T, t) • E xam p le 2. ut — 0.16 uxx + us — bu u(x, 0) u(0, t) = u(lt, t) T O = 0 = 0, x G (0,7 r) = 0, t > 0. = 0 = sinx, xG (0,7t) = 0, t > 0. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. First five figures (Fig.11.1 through Fig.11.5) are related to the numerical results of Example 1, the non-polynomial nonlinearity example. Fig.11.1 and Fig.11.2 show convergence diagrams with errors measured at time T — 1. Figure 11.1 is for different values of the parameter a but for fixed number of unknowns M = 40. Figure 11.2 is for different M but fixed a = 1. In these figures, we display the L2 errors as a function of the p-mode, Np. The value Np is represented on the horizontal axis and the corresponding L2 error is on the vertical axis. The value of Np grows by 4 in the range [4,16] which corresponds to 10% through 40% of the total number of modes M. In each of Figs.11.3,11.4, we show four plots in which the L2 errors are compared for different methods as a function of a for different NpS when M — 40 and 80, respectively. In both Fig.11.3 and Fig.11.4, p-mode Nv has same four different values 4, 8 , 12, and 16 which correspond 10% through 40% of M = 40 in Figure 11.3 and 5% through 20% of M — 80 in Figure 11.4. In each plot, the value a is represented on the horizontal axis and the corresponding L 2 error is on the vertical axis. Here, a grows as powers of 2 in the range [ 2 -4 , 24]. In Figure 11.5, we show four plots of different nonlinear terms and their first and second derivatives. Each plot has different values of a — 1/16, 1/2, 2, 16. We note that the function F(u) = eu/(1+au) is increasing on [0, oo) with F(Q) = 1 , .F(oo) = e1/® . If a > 1/2, the function is concave for u > 0. If a < 1/2, the function is convex for u < ( 1 — 2a)/a2 and then concave for u > ( 1 — 2 a)/a2. 71 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The time evolution of the numerical solution of Example 1 is in Figure 11.6 for different values of a = 1/16, 1 / 2 , 2 , 16. Next six figures (Fig.11.7 through Fig.11.13) are related to the numerical re sults for Example 2, the polynomial nonlinearity example. Specifically, Figure 11.7 through Figure 11.11 correspond to Figure 11.1 through Figure 11.5, respec tively. Here, L2 error is measured at time T — 40. The final time T = 40 was chosen well into the steady state by demanding th at |n(T) —u{T—1)| < 10~n . In each Fig.11.9 and Fig.11.10, b has four values of 1, 2, 3, and 4. The time evolution of the numerical solution of Example 2 is in Figure 11.12 for different values of b = 1, 2, 3, 4. Figure 11.13 shows the importance of choosing sufficiently large value of the p-mode to achieve better accuracy of our perturbation method with two terms (P2). Here, the L2 -errors for the Euler-Galerkin method (EG) and the perturbation method with two terms (P2) as a function of b for Example 2 when M = 40 are presented. In Fig.11.1 and Fig.11.7, for both non-polynomial and polynomial nonlinearity examples, our perturbation method with two terms P2 is superior to the others for all considered values of Np (10% through 40% of M). We note that changing parameter a in Example 1 and b in Example 2 contributes to significant changes in the derivative terms F'{u) or F"{u). This allows us to study their influence on the performance of all algorithms. For all a and b we considered, P2 is superior to the others when Np — 12 and 16. But for smaller values of nodes Np = 4 and 72 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 8 , P2 could not improve the accuracy when a — 1/16 and 16 for Example 1 and b = 4 for Example 2. This is because P2 needs the conditions (10.25) related to the nonlinear term. When Np = 4, i.e., e = I/A 5 , the first or the second derivative of the nonlinear term is increasing too much to satisfy the conditions (10.25) as u is increasing. (See Figure 5 and 10.) This means that Np = 4 is too small to satisfy the conditions (10.25). We conclude that the conditions (10.25) for P2 must be satisfied in order to gain some improvement in accuracy. Furthermore, we can see from Fig. 11.2 and Fig. 11.8 that the accuracy depends on e i.e., on the value of the p-mode Np not on the percentage of p-mode (i.e., ratio Np/M). For example, Np = 8 is 40% of M — 20, 20% of M — 40 and just 10% of M = 80, i.e, same p-mode but different percentage. All these cases show very close L2-errors in Fig.1 1 . 2 and Fig.11.8. This is much clearer in figures, Fig.11.3, Fig.11.4, Fig.11.9, and Fig.11.10. We have the same pattern of accuracy when we preserve the same e for different total value of the mode M = 40 and 80 in the both Example 1 and Example 2. Summarizing the experiments, we see that P2 improves the accuracy of P I and that P 2 is much more accurate than SG and EG for both polynomial and non polynomial nonlinearity whenever the value of Np satisfies the condition (10.25). If Np is too small to satisfy the condition (10.25), we cannot gain any improvement in accuracy by using our perturbation method P2. (See Figure 11.13) 73 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11.1 Euler-G alerkin m ethod vs. our perturba tion m eth od In this section we concentrate on comparing our perturbation method P2 with the Euler-Galerkin method which was introduced by Foias, Sell, and Titi [13] and subsequently used in many studies, (see [14],[15]) We investigate the efficiency as well as the accuracy of these methods in detail. We use the same examples to the previous subsection but with parameters a and b fixed as 1. The total number of modes M is 40. The final time T is 1 for both examples. We add the result of accuracy for Example 2 at T = 40. W ith similar computation to the previous experiment, we compare L2-errors for accuracy. To study the computational efficiency (cost) of the schemes we compare the CPU time and the number of Mflops for Np = 4, 8 , 12 and 16. For the non-polynomial nonlinearity example (Example 1), the error of P2 is at least about 500 times smaller than the error of EG, see Table 11.1. This Table 1 1 . 1 corresponds to Figure 11.3. Furthermore, the ratio of error (EG /P2) is increasing up to 10000 times as Np is increasing from 4 to 16. However, for the polynomial nonlinearity example (Example 2), P2 could not improve the accuracy as much as for the Example 1. The ratio of error (EG/P2) is less than 81, see Table 11.2. When T = 40, the solution arrived at the steady state, the improvement of the accuracy is better as Np = 16. Thus, we can conclude about the accuracy that 74 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. P2 gives a better accuracy than EG for all reasonable values of Np — 4 through 16 when total mode M = 40. Moreover, its improvement in accuracy is increased as Np is increased. On the other hand, P2 is computationally costlier than EG, see Table 11.3 and Table 11.4. The number of Mflops of P 2 is about 2.5 times more than the one of EG and P 2 uses about 2 times more cpu time than EG for both Example 1 and Example 2. This is because P2 has one additional term, the first derivative term in (11.10), in comparison with EG. Therefore, the cost of P 2 doubles that of EG. Here, we note that the computational cost for Np = 16 is less than the cost when Np — 12, see Table 11.4 and Figure 11.14. This is because the higher accuracy when Np = 16 allows coarser step sizes of time integration in ODE solver, ODE15s, even if the number of unknowns is increased. To bring another perspective, we plot the comparison of the accuracy vs. efficiency of the considered algorithms, see Figure 11.14. It clearly shows that P 2 has better accuracy but is more costly than EG. 75 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 11.1: Comparison of the accuracy for Example 1 (non-polynomial) Np e EG P2 EG /P2 4 . 1 2 2 9.00e-2 1.826-4 495 8 .039 1.48e-2 7.64e-6 1937 1 2 .019 5.136-3 8.85e-7 5797 16 . 0 1 2 2.41e-3 2.08e-7 11587 Table 11.2: Comparison of the accuracy for Example 2 (polynomial) Np e T=1 T=40 EG P2 EG /P2 EG P2 EG/P2 4 .253 2.85e-2 6.31e-3 5 5.59e-2 9.33e-3 6 8 .080 7.71e-4 2.33e-5 33 1.76e-3 2.74e-5 64 1 2 .040 1.83e-5 2.26e-7 81 5.54e-5 2.21e-7 251 16 .025 4.35e-7 2.54e-8 17 1.79e-6 2.84e-9 630 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 11.3: Comparison of the efficiency for Example 1 (non-polynomial) Number of Mflops CPU time (sec) Np EG P2 P2/EG EG P2 P2/EG 4 14.0 35.2 2.5 .359 .679 1.9 8 24.6 58.5 2.4 .638 1.112 1.7 12 29.9 73.2 2.5 .735 1.382 1.9 16 35.7 80.1 2.2 .865 1.508 1.7 Table 11.4: Comparison of the efficiency for Example 2 (polynomial) Number of Mflops CPU time (sec) Np EG P2 P2/EG EG P2 P2/EG 4 13.0 29.7 2.3 .342 .586 1.7 8 17.3 43.6 2.5 .468 .839 1.8 12 20.5 52.0 2.5 .512 .981 1.9 16 19.1 45.7 2.4 .482 .866 1.8 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. a=1/16 Fixed M ,-5 «> 10' SG - B - P1 .-10 a=1/2 a=2 .-1 0 Np Figure 11.1: Convergence result of Example 1 for M = 40 The L2 -errors committed by the four different methods with different P-mode Np and the parameter a for the reaction diffusion equation with the non-polynomial nonlinearity. 78 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. M~20 Fixed a M=40 .-5 SG ^ EG - e - P1 ~e~ P2 .-10 .-5 .-10 M=60 M=80 .-5 .-■10 24 Np CM .-10 Np Figure 11.2: Convergence result of Example 1 for a = 1 The L2 -errors committed by the four different methods with different P-mode Np and the number of unknowns M. Each Np is corresponding to 10% through 40% of four different M ’s. 79 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. M -4 0 Np=4 (10%) .-2 S G -<t E G -B - P1 - e - P2 ,-2 1 0 ° 10" 1 0 “ 10" 1 0 " Np=8 (20%) epsilon=0.0385 1 0 " 10“ 1 0 .-2 b . O t a Np=16 (40%) 1 0 - 2 1 0 icf4 i-8 1 0 epsilon=0.0120 ,-8 10 ,0 ,-2 ,2 10 1 0 1 0 ' Figure 11.3: Accuracy of Example 1 for M — 40 Accuracy of Example 1 as a function of a for different Np’s when M — 40 : a grows as powers of 2 in the range [2~4,24] for each different value of Np = 4,8,12,16 which is corresponding to 10% through 40% of M = 40. Here, epsilon is e = 1/A j v p +i- 80 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. M=80 Np=4 SG E G - e - P2 & ? 3 , - G epsilon=0.1204 « Np=8 (1 0 % ) .-4 k -6 .-8 .-2 Np=16 (2 0 % ) ,-2 .-6 epsilon=0.0108 .-2 a Np=12 (15%) epsilon=0.0l81 .-8 1 0 a Figure 11.4: Accuracy of Example 1 for M — 80 Accuracy of Example 1 as a function of a for different Np’s when M = 80 : a grows as powers of 2 in the range [2~4,24] for each different value of Np = 4,8,12,16 which is corresponding to 5% through 20% of M = 80 . Here, epsilon is e = iAjVp+i- 81 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. a=1/16 a=1/2 2 0 1 0 0 1 0 DF DDF -20 1 2 4 0 3 n a-16 snsn 0.1 0.2 0.3 0.4 Figure 11.5: Nonlinear term and its first and second derivatives of Example 1 F(u) is the nonlinear term of Example 1. F, DF and DDF stand for F(u), F'{u) and F"(u), respectively. 82 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. u (x ,t) u(x,t) a=1/16 a=1/2 Figure 11.6: Time evolution of the numerical solution of Example 1 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. b=1 Fixed M b=2 — ^— SG -<h EG -B - P1 -©- P2 .-5 f j 10" & , 1 .2 1 0 ' 1 0 1 0 “ b = 3 .0 10 ' .-5 10 ,0 .1 .2 10 10 ' 1 0 ‘ Np Figure 11.7: Convergence result of Example 2 for M = 40 The L2 -errors committed by the four different methods with different P-mode Np and the param eter b for the reaction diffusion equation with the polynomial nonlinearity. 84 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. M=20 Fixed b u 0 l m i. 0 ) 1 CM - J SG - M r " E G - a - P1 - e - P2 .-10 2 4 6 8 M=40 k -10 M=60 M=80 -10 24 Np -10 24 Np Figure 11.8: Convergence result of Example 2 for b = 1 The L2-errors committed by the four different methods with different P-mode Np and the number of unknowns M. Each Np is corresponding to 10% through 40% of four different M ’s. 85 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Np=4 (10%) M=40 Np=8 (20%) SG EG - a - P1 - © _ p2 3 ,-10 epsilon=0.0803 ,-10 Np=12 (30%) Np=16 (40%) epsilon =0.0250 ,-10 b epsilon=0.0402 ,-1 0 b Figure 11.9: Accuracy of Example 2 for M = 40 Accuracy of Example 2 as a function of b for different iVp’s when M = 40 : b has four values of 1, 2, 3, and 4 for each different value of Np = 4,8,12,16 which is corresponding to 10% through 40% of M — 40. Here, epsilon is e = l/Ajvp+i- 86 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Np=4 (5%) M=80 ,-5 CM .-10 Np=8 (10%) .-5 .-10 Np=12 (15%) Np=16 (20%) .-5 .-1 0 b .-5 .-10 b Figure 11.10: Accuracy of Example 2 for M = 80 Accuracy of Example 2 as a function of b for different Np’s when M = 80 : b has four values of 1, 2, 3, and 4 for each different value of Np — 4,8,12,16 which is corresponding to 5% through 20% of M = 80. Here, epsilon is e = 1 / A j v + i. 87 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. b=1 15 DF 10 5 0 -5 1 1.5 0.5 2 b=2 15 10 5 0 -5 0 0.5 1 1.5 2 b=3 b=4 15 10 5 0 ■ 5 0.5 1 1.5 2 0 u . -5 0.5 1 1.5 2 0 u u Figure 11.11: Nonlinear term and its first and second derivatives of Example 2 F(u) is the nonlinear term of Example 2. F, DF and DDF stand for F(u), F'(u) and F"(u), respectively. 88 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 11.12: Time evolution of the numerical solution of Example 2 89 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Np=4 (10%) 1 0 . Np=8 (20%) 0.08 o 0.06 j 0.04 0.02 b Figure 11.13: The importance of P-mode for the better accuracy of our P2 The L2 -errors committed by the Euler-Galerkin method (EG) and the perturba tion method with two terms (P2) as a function of b for Example 2 when M = 40 and T — 40. 90 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Accuracy vs. efficiency eg - e - P 2 .-2 ® 1 0 Example 1 .-6 -8 1 0 1 0 "1 .0 10 o Example 2 -10 .0 1 0 ' CPU tim e (sec) Figure 11.14: Comparison of the accuracy vs. efficiency Comparison of the accuracy (L2 -error) vs. efficiency (CPU time) for the Euler- Galerkin method (EG) and the perturbation method with two terms (P2) for each a, b = 1 and M — 40. 91 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 12 Concluding remarks In this paper we have presented an approximate inertial manifold (AIM) by us ing a perturbation technique. Our AIM is constructed by the modification of the perturbation result for numerical purposes. The perturbation result is ex pressed explicitly as a finite sum in e which is an approximate solution of induced perturbed equation from the split q equation. Compared with the other versions of the nonlinear Galerkin methods such as Euler-Galerkin AIM (EG) and Foias, Manley, and Temam AIM (elsewhere P I), our AIM induced by the modified perturbation method with two terms (denoted by P2) is more accurate for both polynomial and non-polynomial nonlinearity whenever the number of p-modes, Np, satisfies the conditions (10.25). The nu merical result shows the superior accuracy of our AIM obtained by considering higher values of Np. At the same time computational cost almost doubles. This 92 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. is because our AIM includes the derivative of the nonlinear term as well as the nonlinear reaction term itself using the perturbation technique. 93 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Part III Inertial manifold w ith perturbation technique Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 13 Inertial m anifold w ith perturbation technique Here we assume that an inertial manifold exist for the class of equations (9.1). Instead, we search for an inertial manifold M , the graph of a smooth function < 3 > mapping P H into QH, in the form of power series expansion by using a pertur bation technique. 13.1 Inertial m anifold The nonlinear evolution equation is again du . , . + Au -T F(u) — 0, (13.1) 95 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where F together with its derivatives are uniformly bounded (say, by K ) , i.e., \F ^ ( u ) \< K , i = 0,1,2,..-. (13.2) We assume that q = < & (p) is a smooth function of p. Then we have l - T - T I - n '-' m where the prime and the overdot denote derivatives with respect to p and t, respectively. By Using (13.3), we can rewrite eq.(10.2) as §'(p)p + A$(p) + Q F(p + $(p)) = 0. (13.4) From (10.1) p can be substituted as —Ap — PF(p + <h(p)) and thus eq.(13.4) is rewritten as $'(p) [-A p - PF (p + $(p))] + A $(p) + QF(p + $(p)) = 0. (13.5) In this case, the dependent variable is q = <h(p) and the independent variable is p. Introducing the same operator L as the previous section and multiplying eq.(13.5) by e we obtain a perturbed equation e$'(p) [— Ap - P F (p + <&(p))] + L${p) + eQF(p + $(p)) = 0, (13.6) where operators L, P, and Q are linear and the prime denotes a derivative with respect to the independent variable p. 96 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. We can use the same idea as the previous section to solve the perturbed equa tion (13.6), i.e., we expand the solution d> (p) in a power series in e, $(p) = 0o(p) + e0i(p) + e202(p) + ---. (13.7) and need to solve for the unknown functions 0»(p) recursively. We substitute the formal power series expansion (13.7) into our perturbed equation (13.6) and do a formal derivative with respect to p we obtain e |^0o + e0i + e202 ' ■ ‘j j ^ —Ap — P F (p-\- < f > o + £0i + £202 + ■ • •) -\-L + e202 + •■•) + eQF (p -f- 0o + e0i + e20 2 + • • = 0 (13.8) where 0*’s are now functions of p. With the similar method to the previous section, collecting coefficients of equal powers of e in eq.(13.8) gives j L0o + e 0 o [—A p — PF (p + 0o)] + L < j) i + QF(p + 0 o ) + •■• = 0, (13.9) where only terms up to O(e) have been retained. Equating the coefficient of e° to zero in (13.9) yields 0o = 0, which is the asymptotic limit of the regular perturbation problem (13.6). Fur thermore, 0o = O . 97 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Therefore, equating the coefficient of e1 to zero in (13.9) we obtain L < j > i + QF(p) = 0, i.e., fa = ~ L ~ 1QF(p). (13.10) To obtain the higher-order terms, we use Taylor series for P F (p + c'fa) and Q F (p + X 4 S 1 e% fa -\-----) about p. These formula can be induced using the same notation as the previous section: i.e., P F (lP + f l = P F o + + £ 2 p F 2 + • • ■ + enPFn + 0 (en+1), (13.11) QF [P + J2 eV fl = QFq + eQF1 + e2QF2 + ••■ + enQFn + 0(en+1). (13.12) i— l Using Taylor approximation (13.11) and (13.12), we rewrite eq.(13.8) such as b '0 + efa + e2 fa + • • -j ]^~Ap — P F (p + f a e f a + e2 fa + • • •) + L ^fa + efa + e2fa + • • • j + eQF (p + fa + efa + e2 fa + • • = 0 e2fa + e3fa + ■ ■ • 11- A p - P F 0 - ePF1 - e2PF2 + ■ • eLfa + e2L(j>2 T ^ L<fa -f- + g QFq + e2QF\ + e3QF2 + + 0. (13.13) Collecting coefficients of equal powers of e in eq.(13.13) gives Lfa + QF0 + e Lfa + QFi — fa[Ap + P F q] 98 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. L fa + QF2 - fa[A p + PF 0] - fa P F i L fa + QF-i — fa [Ap + PF$\ — fa P F i — fa P F 2 + e 0 L f a + QFn_i - (t>fai[Ap + P F 0] - 4>fa2P F i ----------fa P F n- 2 where terms up to 0 (en) have been retained. Equating the coefficient of each power of e to zero, we obtain the explicit for mula for each fa as follows: = - L ^ Q F q = - L - ' Q F f a ) , = - L - 1(Q F l - f a [ A p + P F 0\ fa fa L - 1[Q F ,{ p )< j> 1 - fa [ A p + PF{pl = - X T 1 ( q f ? - & { A p + W 0] - # P F ^ - i '1 (QF'ip)^ + y QF"(p)tf - < j,'2]Ap + PF(p)} - ^PF'lp)^ ), -XT1 (q F 3 - fo[Ap + PF0] - <&PF1 - # PF2 -L-1(Q F '(p)<h + faF''(p)-El‘ t > i < t > i] + faF''(l>) E [MsM V i+j=3 i+j+k=3 M A P + PP(p)l - &PF'(p)h - < f[ \PF'(p)<h + f a p " ( p ) E [Ms] i+j=2 -L -1 ( QF'{p)fa + Q F"{p)fafa + ^Q F "(p )fa 99 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -fo[Ap + PF(p)\ - & PF'(p) * - ft [iW (p)fc + = — L [QFn., - [Ap + PF0] - 4>'„-iPFi - & ^ P F 2 - = ~ L - 1(QF'(p)4,r , . l + ^ Q F "(p ) Y , [ « ( ] + ••• + i+ j= n — 1 - t i - M p + PF(p)] - K-,PF'(p)<h < - 3 [p p 'W * + ^ P P " (p )^ S ]--------- # [p f '(p ) ^ 2 + i + j = n - 2 Z J. J / 13.2 R eaction diffusion equations w ith a p oly nom ial nonlinearity In this section we consider a general reaction diffusion equation with a polynomial nonlinearity (see [11], [18], [24]). This is a good example to apply our inertial manifold with perturbation method. Let Q be an open bounded set of R n with boundary F and g be a polynomial of odd degree with a positive leading coefficient 2p— 1 g(s) = 5Z fys-7, b2 p -i > 0. 3= 0 We consider the following problem involving a real valued function u = u(x,t)\ du dt dA u + g(u) = 0, in fi x R +. (13.14) 100 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The equation is supplemented with the initial and the Dirichlet boundary conditions u(x,0) — u0(x), i G f l, (13.15) u — 0, x £ F. (13.16) For the mathematical setting of this problem let us introduce the operator Au = —d A u on H = L2{V L ) and nonlinear operator F(u) = g(u) and define the semigroup S(t) [11] (see the references therein) S(t) :u 0 e H -» ■ u(t) £ H. Then it is known (See [11]) that the semigroup S(t) associated with the boundary-value problem (13.14-13.16) possesses a maximal attractor A which is bounded in H q(Q), compact and connected in L 2(Q). Moreover, we can get the same property above for the boundary-value problem (13.14,13.15) with one of the following boundary conditions : On Neumann — = 0 on T (13.17) ov Periodic il = (0, L)n and u is f2-periodic (13.18) where v is the unit outward normal on T. Let uQ £ B(0, R) C H be given where B(0, R) C H is a ball centered at 0 of radius R. Then there exists a time to \\u\\Lo o < K 0, Vt > t0- 101 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Thus we obtain that 3 K s.t. |F*(s)| < K , Vi — 0,1,2, • • •, |s| < K 0 which satisfies our condition (13.2) for nonlinear term. For the long-term dynamics of the reaction diffusion equation with one of three boundary conditions (13.16-13.18) we can apply our inertial manifold with perturbation method. 102 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 14 Convergence results In the previous chapter we have found the exact form of a smooth function $ : P D (A ) — ► Q D (A) whose graph is an inertial manifold. This function is obtained by a formal power series expansion in e through the perturbation technique: 00 $ (p) = (14.1) i= 1 Under the assumption th at the nonlinear operator F is uniformly bounded and infinitely differentiable every function 4>i(p) in the infinite series is uniformly bounded smooth function, too. This is because each 4>i(p) depends on only - Thus if the infinite series converges pointwise and Y ^ i con verges uniformly then term-by-term differentiation in the perturbation technique (13.8) can be defined. That is, the infinite series < I> (p ) — el0*(p) can be the exact solution of the q-mode equation (10.2). 103 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Let < L n be the partial sum = n = l ,2 , i— 1 (14.2) We will show that the sequence {<& „} is convergent. Since our solution space H is a Hilbert space, it is enough to show that the sequence { < 3 > n} is a Cauchy sequence by the Cauchy criterion. We will introduce three lemmas which are useful to prove the convergence. Lem m a 14.1 Y ,e iFi = F (p + $ n) + 0 (e n+1), i=0 for n = 1,2,3, Proof. It is easy to see th a t for k = 2,3, • • •, n, because «£ = K - h -i + O (e"+1) , = ft +••• + £ ” 1 + e" k+2< p n~k+2 + + €n k+2<j>n-k+2 + '' ' + 4 -tt+1 + O O (£ "-‘+ 2) ® !U +1 + O (e1- 1 ) O (e”- ‘+2) + O (e”+1) . (14.3) + 4 > n 104 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Now let us consider Taylor formula for F(p + Tn) about p of order n : F(p + $ n) 1 + 0(en+1) (14.4) for some p in the neighborhood of p in the first equality. Substitute the formula (14.3) into Taylor expansion (14.4) then F(p + < f > n) can be rewritten as F(j> + *„) = F (p) + F '(p )$ „ + • • • + + • • • + K > + ole"*1) n\ (14.5) Thanks to the definition of Fj’s in eq.(10.13) we can see that eq.(14.5) is equal to eq.(10.13): (14.5) F(p) + F'(p) + • • • + „ F W (P) 1 601 k\ + 0 (en+1) n— k+l H e<& by(10.13) i= 1 F0 + eFi -i + ehFk + h e”Fn + C?(en+1) (10.13). Thus F(p + $„) = F0 + eF1 + • • • + ekFk + ■ ■ • + enFn + £?(en+1). 105 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Lemma 14.2 = -e L ~ 1(Q F(p + ^ 1) - &n-i[A p + P F (p + $„_,)]) + 0 (en+l), for n = 2 ,3, • • •. Proof. Using each explicit formula of 0*’ s we have obtained in the previous section: 01 = -L -'Q F o , 02 = - L - ^ Q F r - ^ A p + P F ^ , 03 = - L - 1 (Q F2 - (ff2[Ap + PF0] - f t P F ^ , 0n = - L - 1 ( Q F n_ ! - + PFo] - < l > , n_ 2P F 1 - 0 U ^ 2 ------------ 0 / l P P n - 2) the partial sum 4 > n is written as follows: = £01 + £202 + ' ' ' + £n0n - - e L - ' Q F o - e L - 1 (e Q F i - e0'[^4p + P F 0j ) - e l T 1 ( e 2Q F 2 - (?<ff2[Ap + P F 0] - e ^ e P F ^ - e L - l [en~l QFn^ - [Ap + PF0) ----------e ^ e ^ P F ^ ( n — 1 n — 2 n — 3 £ e’Q P - £ e ^ A p + PPFi\ - £ e2( j ) '2[Ap + elPF^ i=0 i=0 t=0 106 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. By Lemma 1 the first term in the equation (14.6) can be approximated as follows: n — 1 E ^ Q Fi = QF(P + *n-i) + 0 ( e n). (14.7) i= 0 Moreover, for k = 0,1,2, • • •, n — 2, n — l k J 2 ^ P F i = Y / eiPFi + ek+1P F k+1 + --- + en~1PFn j= 0 i= 0 k = E e iF F i + 0 (efc+ 1 ) i= 0 Therefore k n — l E = E + C0(efc+1) i= 0 i= 0 L e^al p i?(p + $ n_1) + (9 (en) + C ,(efc + 1) = P F (p + $ n_1) + 0 (e fc+ 1). (14.8) By using(14.8) we can simplify the second term in (14.6) as follows: n — l n — j - l E E + 1=1 i= 0 n — l n — j — 1 - E^>'[^+ E ^ 3 jf= l 1=0 by(H.8) g + + + ] 1=1 n — l - E + P F (P + $ n-i)] + 0 ( e n). (14.9) i=i 107 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Substituting approximations (14.7) and (14.9) into the equation (14.6) the partial sum becomes n —l $ n = —eL~l [QF{p + - J 2 + P F (P + + 0(en+1) = - e l T 1 ( q F(p + < f > n _ x) - ^ [ A p + PF(p + $ . _ ! ) ] ) + 0 (en+l) which completes the proof. L em m a 14.3 The operator norm || L 1 || is at most 1, i.e. || L 1 ||< 1. Proof. L ^ u OO i= 1 OO since L = eAQ L 1 J 2 (u » ei)ei i= N + 1 o o \ ^JV +1 A . E i= N + l ' H (u, ei)ei ~ ( \ N + l Y u \ ,2 = L - r - IKei)l i= N + 1 V A* / OO < j r \(u,ei)I 2 i= JV + l oo < £ |( u ,e i )|2 i = 1 = U Now we prove the main theorem regarding the convergence of the infinite series. 108 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. T heorem 3 {$„} is a Cauchy sequence. P ro o f. We will show that |$ n — $ m| — > 0 as n, m — ■ > 0 0 . Let n > m. Since F(u) is uniformly bounded smooth function of u in D (A), the first-order Taylor formula for F(p + $ n-i) at p + gives F(p + $ n_i) — F(p + $m_i) + F /(ft)($„_i — $ m_i) (14.10) = F(p + $ m_i) + F'(u)(em(fim + em+1( j)m+i H h en V n - 1 ) (14.11) for some u in the neighborhood of p+^m-i in the first equality (14.10). Moreover, since P and Q are linear projections the approximation (14.11) becomes P F (p + $„_!) = P F (p + ^m -i) + O (em) , (14.12) Q F(p + $ n_x) - Q F { p P ^ m_l) p O { e m). (14.13) By Lemma 2 and Lemma 3 we obtain that eL -1(Q F(p + $„_,) - K - M p + PF(p + $„_!)]) + 0 (e n+1) Lemma2 eL -1 QF(p + [Ap + PF(p + ^ m^ )} + 0 (e m+1) < e I I L - 1 +e I I L - 1 Q F(p + 4 > r a _ 1 ) — QF(p + $ m_i) ^ [ A p + P F (p + $„_,)] - C - i [Ap + P F (p + $ m_x)] 109 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Lemma3 < QF(p + ^ 1) - Q F ( p + ^ m_1) (14.14) +6 K - l i AP + P F (P + $ «-l)] - &m-l[AP + P F (P + $m-l)] (14.15) + 0 (em+1) . In the last inequality the first term (14.14) becomes QF(p + ^ n_1) - Q F ( p + ^ 1) O (em+1) (14.16) which is because of the equation (14.13). Similarly, thanks to eq.(14.12) we can rewrite the second term (14.15) as fol lows: by(14.12) < K - l l AP + P F iP + - $ m -l\A P + P F (P + $m-l)] ^ - 1 [ A p + P F ( p + $„,_!) + O (£"•)] - V ^ l A p + P F ( p + «•„„!)] - < _ , ] [A p + P F ( p + # „ _ ,)] + O (£-) K - i - »m -i A p + P F ( p + * m_!) + O (e”* 1 ) . (14.17) In this inequality (14.17), since F is uniformly bounded by K , we can show that Ap + P F (p + $ m-i) is bounded for all p in PH: Ap + PF(p + < b T O — 1 ; < \Ap\ + \PF(p + $m_i)| < Aat • \p\ + |.F| 110 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. < XN ■ sup|p| + K Let M. (14.18) Furthermore, since every function fa is uniformly bounded smooth function, we obtain n — l (14.19) Finally, the equations (14.14) through (14.19) can be summarized as follows: by(14.14,14.15) < e l^n ^m| Q F(p + $n_x) - QF(p + $m_i) + e |&n_x[Ap + PF (p + $„_!)] - Q'm-AAp + PF(p + $ m_i)] +0 (em +1) by(14.16— 14.19) , x . N < O (em+1) + eMO (em) + O (em+1) O (em+1) — » 0 as m — > oo. This completes the proof. Q.E.D. Now let $ be the limit of the power series OO 4 > = E £ ‘&' i—1 We need to show th at i Pfa is uniformly convergent in order to define term- by-term derivative in the perturbation technique (13.8). With the similar idea 111 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. to Theorem 1 we can use the Cauchy criterion to show th at the sequence { < E > (J converges uniformly. That is, it is enough to show that the sequence is a Cauchy sequence because our solution space H is complete. Here is the Theorem for the uniform convergence of the sequence {$',}. Theorem 4 {<f^} is a Cauchy sequence. However, we are not going to give the proof in detail because the equation (14.19) itself can prove this theorem exactly. 112 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Bibliography [1] Aronson, D.G.; The porous medium equation, in: Fasano, A. and Primice- rio, M.(Eds.), Nonlinear diffusion problems, Lecture Notes in Mathematics, v.1224, Springer-Verlag, (1985) [2] Betelu, S.I., Aronson, D.G. and Angenent, S.B.; Renormalization study of two-dimensional convergent solutions of the porous medium equation, Physica D, 138, 344-359, (2000) [3] Benilan, P., Chipot, M., Evans, L.C. and Pierre, M.(Eds.); Recent advances in nonlinear elliptic and parabolic problems, Longman Scientific & Technical, (1988) [4] Bandle, C., Nanbu, T. and Stakgold, L; Porous medium equation with ab sorption, SIAM J. Math. Anal., v.29, 1268-1278, (1998) [5] Feng, B.F. and Mitsui, T.; A difference method for the Korteweg-de Vries and the Kadomtsev-Petviashvili equations, J. Comput. Appl. Math. 90, 95- 116,(1998) [6] Jager W. and Kacur J.; Solution of porous medium type systems by linear approximation schemes, Num. Math. 60, 407-427, (1991) [7] Mikula, K.; Numerical solution of nonlinear diffusion with finite extinction phenomenon, Acta Math. Univ. Comenianae, v.LXIV, 2(1995), 173-184 [8] Richtmyer, R.D. and Morton, K.W.; Difference methods for initial-value prob lems, Interscience Publishers, New York, (1967) [9] Smith, G.D.; Numerical solution of PDEs: finite difference methods, Claren don Press, Oxford, 3rd edition, (1985) [10] Sommeijer, B.P. and Van der Houwen, P.J.; Algorithm 621, Software with low storage requirements for two-dimensional nonlinear parabolic PDEs, ACM Trans. Math. Software, v.10, 378-396, (1984) 113 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [11] R. Temam; Infinite-Dimensional Dynamical Systems in Mechanics and Physics, 2nd ed. Springer Verlag, Berlin, 1988. [12] C. Foias, G.R. Sell, R. Temam; Inertial Manifolds for Nonlinear Evolutionary Equations, J. Differential Equations 73, 309-353 (1988) [13] C. Foias, G.R. Sell, E.S. Titi; Exponential tracking and approximation of inertial manifolds for dissipative nonlinear equations , J. Dyn. Differential Equations 1, 199-243 (1989) [14] C. Foias, M.S. Jolly, LG. Kevrekids, G.R. Sell, E.S. Titi; On the computation of inertial manifolds, Phys. Lett. A, no. 131, pp. 433-436, 1988. [15] L. Russo, A. Adrover, G. Continillo; Dynamic behavior of a reaction/diffusion system: wavelet-like collocations and approximate inertial manifolds, COG 2000, ST.Petersburg, Russia, 2000 IEEE. [16] C. Foias, O. Manley, R. Temam; Modelling of the interaction of small and large eddies in two dimensional turbulent flows, Math. Modelling Numer. Anal. 22, 93-118 (1988) [17] M. Marion; Approximate inertial manifolds for reaction-diffusion equations in high space dimension, J. Dynamics Differential Equations 1, 245-267 (1989) [18] A. Debussche, M. Marion; On the Construction of Families of Approximate Inertial Manifolds, J. Differential Equations 100, 173-201 (1992) [19] M.S. Jolly; Explicit Construction of an Inertial Manifold for a Reaction Dif fusion Equation, J. Differential Equations 78, 220-261 (1989) [20] R. Temam; Navier-Stokes Equations and Nonlinear Functional Analysis, CBMS Regional Conferences, SIAM, 1983. [21] J. Robinson; Infinite-Dimensional Dynamical Systems, Cambridge University Press, 2001. [22] J. Mallet-Paret, G. R. Sell; Inertial manifolds for reaction diffusion equations in higher space dimensions, J. Amer. Math. Soc. 1 (1988), 805-866. [23] C. Foias, M.S. Jolly; On the numerical algebraic approximation of global at tractors, Nonlinearity 8, 295-319 (1995) [24] J. Novo, E.S. Titi, S. Wynne; Efficient methods using high accuracy approx imate inertial manifolds, Numer. M ath. 87, 523-554 (2001) 114 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.

Asset Metadata

Creator Kim, Dong-jin (author)

Core Title Efficient methods for nonlinear parabolic PDEs

Contributor Digitized by ProQuest (provenance)

School Graduate School

Degree Doctor of Philosophy

Degree Program Mathematics

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag mathematics,OAI-PMH Harvest

Language English

Advisor Proskurowski, Wlodek (committee chair), Udwadia, Firdaus (committee chair), Blum, Edward (committee member), Kukavica, Igor (committee member)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c16-368512

Unique identifier UC11339375

Identifier 3103916.pdf (filename),usctheses-c16-368512 (legacy record id)

Legacy Identifier 3103916.pdf

Dmrecord 368512

Document Type Dissertation

Rights Kim, Dong-jin

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

Linked assets

University of Southern California Dissertations and Theses