1 / 36

Lecture 3

Lecture 3. A Review of ADI and Operator-Splitting Methods for the Solution of Initial Value Problems. 1. Preliminary Remarks Operator-splitting methods have a controversial reputation. Some practitioners see them as particular time- discretization methods , while others view them as

hamal
Download Presentation

Lecture 3

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Lecture 3 A Review of ADI and Operator-Splitting Methods for the Solution of Initial Value Problems

  2. 1. Preliminary Remarks Operator-splitting methods have a controversial reputation. Some practitioners see them as particular time- discretization methods, while others view them as iterative methods. They are clearly both, even if in this presentation we focus on the first point of view. Incidentally, let us mention that various numerical methods and algorithms are disguised operator-splitting methods as we will see on a particular example. Their main feature is to take advantage of decompositionproperties in the structure of the problem to be solved, some of these properties being obvious, while some others are not.

  3. 2. Time-Discretization of Initial Value Problems by ADI and OS Schemes Let us consider the following IVP (flow for the Dyn. Syst. community): (IVP) dφ/dt + A(φ) = 0 in (0, T), φ(0) = φ0. Operator A: V →V and can be multivalued. Suppose that (DP) A =ΣJj = 1 Aj holds. A natural question is: Can we take advantage of (DP) to solve (IVP) ? As we all know here: the answer is yes!

  4. 2.1. ADI type schemesWe suppose for the moment that J = 2. A first candidate to take advantage of (DP) is provided by the Peaceman-Rachford scheme, namely (with τ > 0 a time-discretization step): (1) φ0 = φ0; forn ≥ 0, φnbeing known, solve (PRS) (2) (φn + ½ – φn)/(½τ) + A1(φn + ½) + A2(φn) = 0, (3) (φn + 1 – φn + ½ )/(½τ) + A1(φn + ½) + A2(φn +1) = 0. (PRS) is of the backward-Euler type for A1 and of the forward-Eulertype for A2 on the time interval [tn, tn+½], the situation being reversed on [tn+½, tn+1].

  5. A second candidate is the Douglas-Rachford scheme (of the predictor-corrector type); if J = 2, the DR-scheme reads as follows (1) φ0 = φ0; forn ≥ 0, φnbeing known, solve (DRS) (2) (φn + ½ – φn)/τ + A1(φn + ½) + A2(φn) = 0, (3) (φn + 1 – φn)/τ + A1(φn + ½) + A2(φn +1) = 0. These schemes, which have been around for about 50 years, have motivated ahuge literature, either as methods to approximate the solution of time dependent problems or as iterative methods to solve linear and nonlinear steady state problems in finite or infinite dimension. In the second case a strategy of periodically variableτ is advocated (E. Wachpress and others). Further remarks are in order; among them:

  6. ● If V is an Hilbert space and A1 and A2 are monotone operators (possibly multivalued) both schemes are unconditionally stable. ● The most general convergence results in the context of monotone operators in Hilbert spaces are those of P.L. Lions-B. Mercier (1979). ● The D-R scheme can be generalized to decompositions with more than two operators (J > 2). ● The particular case A = – 2, A = A1 + A2 with A1= – ∂2/∂x12 and A2 = – ∂2/∂x22 explains the terminology ADI, but the applicability of (PRS) and (DRS) goes much beyond the above type of decomposition since Aj can be for example the subgradient of a convex functional. ● The above schemes are O(τ), at best, in general; however if A1 andA2 are both linear and commute, (PRS) is O(τ2); T. Arbogast, J. Douglas et al have found a way to make (DRS) second order accurate via a slight modification (if A1, A2 are smooth enough)*. ● (PRS) and (DRS) are not stiff A-stable; indeed, (PRS) is very close to the Crank-Nicolson scheme. ● Several French scientists have contributed to ADI; among them, J. Lieutaud, D. Gabay, P.L. Lions, E.Godlewsky, M.Schatzmann.

  7. On the Russian side important contributions by Dyakonov (who passed away recently). What about Spain? Important contributions byE. Fernandez-Cara et al, Bermudez, Moreno, Pares,….. ● In order to improve the A – stability properties of (PRS) we introduced (in the mid-1980s) the so-called (by us) θ-scheme(the German scientists call it the Fractional Stepθ-scheme). The idea is quite simple: (i) Consider θ (0, ½ ) (0 < θ < 1/3, in practice).(ii) Denote (n + )τby tn+ .(iii) Decompose [tn, tn+1] as [tn, tn+1] = [tn, tn+θ][tn+θ, tn+1– θ][tn+1– θ, tn+1]. Use the above decomposition of the time interval [tn, tn+1] to time-discretize the initial value problem (IVP) by the following variant of the Peaceman-Rachford scheme (where θ* = 1 – 2θ):

  8. Description of the θ - Scheme (1) φ0 = φ0; forn ≥ 0, φnbeing known, solve (2) (φn+θ– φn)/(θτ) + A1(φn+θ) + A2(φn) = 0, (θ-scheme) (3) (φn+1–θ– φn+θ)/(θ*τ) + A1(φn+θ) + A2(φn+1–θ) = 0, (4) (φn+1 – φn+1–θ)/(θτ) + A1(φn+1) + A2(φn+1–θ) = 0.

  9. Some Properties of the θ-Scheme ● Properties are more problem dependent than for P-R and D-R; basically, for 1/4 < θ < 1/3the scheme is stiff-A stable. It is only O(τ), but there are situations where θ= 1 – 1/√2makes it O(τ2) (actually, nearly O(τ3); it would be O(τ3) if √2= 1.5). ● Quite good to capture steady state solutions (compared to (PRS) and (DRS)). ● Very popular in Germany (R. Rannacher, S.Turek, E. Bänsch, V. Heuveline and others) for the simulation of incompressible viscous flow. Also used in various places (Australia among them) for simulating non-Newtonian fluid flow.

  10. 2.2. LIE’S, STRANG’S and MARCHUK-YANENKO Schemes Co-existing with ADI, there have been another family of Operator- Splitting schemes going back (cf. Chorin, Hughes, Marsden et al, CPAM 1978) to S. Lie and widely used by CFD, Physics, Chemistry, Russian, etc,,, people. It is related to the Trotter formula in Semi- Group Theory. The basic idea is very simple: suppose that A is linear in (IVP), namely that (IVP) dφ/dt + Aφ = 0 in (0, T), φ(0) = φ0. We have then φ(t) = e – Atφ0, and therefore:

  11. (FSGR) φ(t + τ) = e – Aτφ(t),  t ≥ 0. We also have (FOSR1)e – A2τe – A1τ– e – (A1 + A2)τ= ½ (A2A1– A1A2)τ2 + O(τ3), and (FOSR2)e – ½A1τe – A2τe – ½A1τ– e – (A1+A2)τ= O(τ3)( = 0 if [A1,A2] = 0). Relations (FSGR) and (FOSR1) leads to the following OS scheme, known by some as the Lie’s scheme:

  12. Descriptions of the Lie’s Scheme (I) Condensed form: (1) φ0 = φ0; n ≥ 0, φn→ φn+½ → φn+1 as follows: (LIE’S SCHEME) (2) φn+½ = e – A1τφn, (3) φn+1 = e – A2τφn+½ . The Lie’s scheme is exact if A1 and A2commute.

  13. Descriptions of the Lie’s Scheme (II) Practical and Generalized Form: ( with 0 ≤ ,  ≤ 1,  +  = 1) (1) φ0 = φ0; n ≥ 0, φn→ φn+½ → φn+1 as follows: dφ/dt + A1[φ, tn + (t – tn)] = 0 on (tn , tn+1), (2) φn+½ = φ(tn+1), φ(tn) =φn, dφ/dt + A2[φ, tn+ + (t – tn)] = 0 on (tn , tn+1), (3) φn+1 = φ(tn+1). φ(tn) =φn+½,

  14. Some Properties of the Lie’s Scheme ● Easily generalizable to J > 2 (when simulating particulate flow we may have J ≈ 10). ● Unconditionally stable if the Aj are monotone operators (possibly multivalued). Indeed this scheme is quite robust. ● First order accurate at most, in general. ● Different time and space discretizations can be used to discretize the subproblems (2) and (3) (including closed form solutions). ● The LIE’S scheme is asymptotically inconsistent; when applied as an iterative solver to compute a steady state solution, in general φn and φn+½ converge to different limits whose distance to the exact solution is O(τ) at best [it may happen that ½(φn + φn+½) has better convergence properties]. This makes some practitioners worry (I was one of them) if steady state solutions are required.

  15. Description of the Strang’s Scheme (I) Relations (FSGR) and (FOSR2) leads to the following OS scheme, known as the Strang’s scheme: (1) φ0 = φ0; n ≥ 0, φn→ φn+½ → φ#n+½ → φn+1 as follows: dφ/dt + A1(φ, t) = 0 on (tn , tn+½), (2) φn+½ = φ(tn+½), φ(tn) =φn, dφ/dt + A2(φ, tn+ ½) = 0 on (0, τ), (3) φ#n+½ = φ(τ), φ(0) =φn+½,

  16. Description of the Strang’s Scheme (II) dφ/dt + A1(φ, t) = 0 on (tn+½, tn+1), (4) φn+1 = φ(tn+1). φ(tn+½) =φ#n+½, Properties of the Strang’s Scheme: ● Unconditionally stable if the Aj are monotone operators. ● O( 2) if the Aj are smooth enough. ● Generalizable to J > 2. ● Asymptotically inconsistent but provide steady state solutions whose distance at the exact solution is O( 2). ● The sub-problems (2)-(4) have to be solved by schemes which are themselves second accurate (at least).

  17. The following implicit Runge-Kutta scheme of ordertwo (due to J. Cash) is well-suited to such a task: When applied to the solution of dX/dt + f(X, t) = 0 on (t0, tf), (IVP) X(t0) = X0 the scheme reads as follows: (1) X0 = X0; then for m ≥ 0

  18. Xm+θ+ θΔt f(Xm+θ, tm+θ) = Xm (2) Xm+1 – θ = (1 – θ)/θ Xm+θ + (2θ – 1)/θ Xm(3) Xm+1+ θΔt f(Xm+1, tm+1) = Xm+1 – θ(4) The above scheme is stiff A-stable and 2nd order accurate if θ = 1 – 1/√2 (actually, almost 3rd order accurate).

  19. ● There exist (e.g., M. Schatzman) variants of the Strang’s scheme which are O( 4) but they are not unconditionally stable. Description of the Marchuk-Yanenko scheme Suppose that to implement the Lie’s scheme we discretize the sub- initial value problems using just one step of the backward Euler scheme. We obtain then the following scheme: (1) φ0 = φ0; n ≥ 0, φn→ φn+1/J ……→ φn+(J – 1)/J→ φn+1as follows: (2)j (φn+j/J – φn+(j – 1)/J)/τ+ Aj (φn+j/J, tn+1) = 0, for j = 1, 2, ……., J.

  20. The above scheme is known as the Marchuk-Yanenko scheme. If the Aj are are monotone the scheme is unconditionally stable. It is O(τ) at best, robust and relatively easy to implement. 3.ADI and Augmented Lagrangian algorithms (i) Consider the following functional J : V → R; assume that J is differentiable and that (DP)J = J1+ J2, J1and J2 being both differentiable. (ii) Consider now the minimization problem (MP) Minv  V J (v).

  21. (iii) Suppose that (MP) has a solution denoted by u. We have then: (OC) J’(u) = J1’(u) + J2’(u) = 0. (iv) Problem (MP) is clearly equivalent to (EMP) Min {v, q}  Wj (v, q), with W = {{v, q}| {v, q}  V × V, v – q = 0} and j (v, q) = J1(v) + J2(q). (v) Observe that {u,u} is solution of problem (EMP).

  22. (vi) Assuming that V is a real Hilbert space, we associate to (EMP) the following saddle-point problem Find {{u, p}, λ}  (V × V) × V such that (SDP) Lr(u, p; μ) ≤ Lr(u, p; λ) ≤ Lr(v, q; λ),  {{v, q}, μ}  (V × V) × V, where r > 0 and (LA)Lr(v, q; μ) = j(v, q) + ½ r ||v – q||2 + (μ, v – q). (vii) If {{u, p}, λ}is solution of (SDP) thenp = uwhere uis solution of (MP).

  23. (viii) To solve (SDP) we use the following Relaxation/Uzawa algorithm ((ALG2) in various related references): (1) {u– 1, λ0} is given inV × V; n ≥ 0, {un – 1, λn} → pn →un → λn + 1 as follows (2) J2’(pn) + r (pn –un – 1) – λn = 0, (3) J1’(un) + r (un –pn) +λn = 0, (4) λn + 1=λn + r (un –pn). Comparing (3) and (4) shows that λn + 1 = – J1’(un) which implies in turn that:

  24. AL  ADI (a) r (pn –un – 1) + J2’(pn) + J1’(un – 1) = 0, (b) r (un –un – 1) + J2’(pn) + J1’(un ) = 0. We have recovered thus the Douglas-Rachford scheme with  = 1/r. One can recover the Peaceman-Rachford scheme by updating λn a 1st time between (2) and (3). 4. A striking application. Let Abe a d × d matrix, symmetric and positive definite. Wedenote by λ1the smallest eigenvalue of A; we have then:

  25. (EVM)λ1 = Minv SAv.v with S = {v| v Rd, ||v|| = 1}. (EVM) is equivalent to (EVM-P) Minv Rd[½ Av.v + Is(v)] with 0 if v  S, Is(v) = i.e., Is(.) is the indicator functional of S. + ∞ if v Rd \S If usolves (EVM) we have the following optimality condition:

  26. (NOC) Au + ∂IS(u) = 0, with ∂IS(u) a generalized differential of IS at u. To (NOC), we associate the following flow that we time-discretize by the Marchuk- Yanenko Scheme; we obtain thus: u(0) = u0, (NOC-F) du/dt + Au + ∂IS(u) = 0 and then (1) u0= u0; (M-Y) n ≥ 0; un→ un+ ½ → un+1 as follows

  27. (2) (un+½ – un)/ + Aun+½ = 0, (3) (un+1– un+½)/ + ∂IS(un+1) = 0. Eq. (2) implies that: (2)’ un+½ = (I +  A)–1un Observing that IS = IS, Eq. (3) implies that: un+1 = Arg max v  Sun+½.v , namely (3)’ un+1= un+½/||un+½||. We have reinvented thus the inverse power method with shift.

  28. 5. Another application: Bingham Flow in a Pipe Participating in the Fall of 2005 to a conference on visco-plasticity in Banff (BC), we had the pleasant surprise to discover that the ‘visco- plastic’ community was quite in favor of the Augmented Lagrangian (AL) approach for the numerical simulation of visco-plastic flowwith stress yield, like Bingham’s. We are going to discuss thus the AL solution of a simple Bingham flow problem as an illustration of the above methodology: Let Ω be a bounded domain of R2; we denote by Γ the boundary of Ω. We consider then the following problem from the Calculus of Variations u  H10(Ω), (BFP) J(u)  J(v),  v  H10(Ω),

  29. with J(v) = ½μ∫Ω|v|2dx + g ∫Ω|v|dx – C∫Ωvdx. The idea here is to uncouple nonlinearity and derivatives; to do that we will treat v as an additional unknown q and force the relation v – q = 0 by penalty and Lagrange multiplier. To implement the above idea, we proceed as follows: (i) Introduce Q = (L2(Ω))2, W = {{v,q}| v  H10(Ω), q Q, q = v }, j(v,q) =½μ∫Ω|v|2dx + g∫Ω |q|dx – C∫Ωvdx. (ii) Observe that (BFP) is equivalent to: (BFP-E){u, p}  W, j(u,p) j(v,q),  {v,q}  W.

  30. (iii) Introduce the following augmented LagrangianLr(with r > 0)from (H10(Ω) × Q) × Q into R: Lr ({v,q},m) =j(v,q) +½r ∫Ω|v – q|dx + ∫Ωm.(v – q)dx and observe that if {{u,p},l}is a saddle-point of Lr over the space (H10(Ω) × Q) × Q(i.e., {{u,p},l}  (H10(Ω) × Q) × Q, (SDPP) Lr ({u,p},m)  Lr ({u,p},l)  Lr ({v,q},l),  {{v,q},m}  (H10(Ω) × Q) × Q, then {u,p}solves (BFP-E) which implies that usolves (BFP) and p = u.

  31. (iii) In order to solve (BFP-E) we advocate the following algorithm(a disguised Douglas-Rachford ADI algorithm): (1)u – 1andl0are given in H10(Ω)and Q. For n ≥ 0, un – 1andlnbeing known, solve • pn  Q,Lr ({un – 1, pn }, ln)  Lr ({un – 1, q}, ln), q  Q, then (3) un  H10(Ω),Lr ({un, pn }, ln)  Lr ({v, pn }, ln), v H10(Ω), and finally (4) ln+1 = ln + r ( un– pn ). The convergence follows from, e.g., RG-Le Tallec, 1989.

  32. The sub-problems (2) and (3) are simpler that what they look like since (a)1/r (1 – g / |Xn(x)|) if |Xn(x)| > g, pn( x ) = 0if |Xn(x)| g, withXn= r un – 1 + ln. (b) Problem (3) is equivalent to un  H10(Ω), (LVP) (μ + r)∫Ωun.vdx= C∫Ωvdx +∫Ω(r pn – ln).vdx,v H10(Ω), a ‘simple’ linear problem indeed.

  33. 6. More applications

  34. Above, we have visualized the behavior of the mixture of solid sphericalparticles with a Newtonian incompressible viscous fluid in a rotating cylinder, the number of particles being 160, here. When the angularvelocity is sufficiently large, the particles cluster in 3 sub-populations, of equal size, approximately (this reminds the formation of Bénard rolls in heated flow). The simulation relies on the Lie’s scheme, with J of the order of 10. We could have illustrated our presentation with many more results from numerical experiments related to many areas of applications. Indeed, new applications of ADI/OS are discovered almost every day (e.g., Monge-Ampère type equations) and we are not sure that the Scientific Community is fully aware of the capabilities of the ADI/OS methods. 7. Some References [1] GLOWINSKI, R. and P. LE TALLEC, Augmented Lagrangians andOperator –Splitting Methods in Nonlinear Mechanics, SIAM, Philadelphia, PA, 1990.

  35. [2] GLOWINSKI, R., Finite element methods for incompressible viscous flow. In Handbook of Numerical Analysis, Vol. IX, P.G. Ciarlet and J.L. Lions, eds., North-Holland, Amsterdam, 2003, pp. 3-1176. [3]MARCHUK, G.I., Splitting and alternating direction methods. In Handbook of Numerical Analysis, Vol. I, P.G. Ciarlet and J.L. Lions, eds., North-Holland, Amsterdam, 1990, pp. 197-462. [4] DEAN, E.J. and R.GLOWINSKI, An augmented Lagrangian approach to the numerical solution of the Dirichlet problem for the elliptic Monge-Ampère equation in two dimension, Electronic Transactions in Numerical Analysis, 22, (2006), 71-96.

  36. Additional references on ADI and OS can be found in refs.[1]-[3]. In particular, the Chapters 2 and 6 of ref. [2] are almost self-contained introductions to ADI and OS methods. A recent reference concerning inverse problems (4th order) of elliptic nature is: [5] F. DELBOS, J.CH. GILBERT, R. GLOWINSKI & D. SINOQUET, Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach, Geophys. J. International, 164, (2006), 670-684.

More Related