Sparse Dynamics for Partial Differential Equations
Hayden Schaeffer, Stanley Osher, Russel Caflisch, Cory Hauck
aa r X i v : . [ m a t h . NA ] D ec Sparse Dynamics for Partial Differential Equations
Hayden Schaeffer, Stanley Osher, and Russel CaflischUniversity of California, Los Angeles, Department of MathematicsCory HauckOak Ridge National Laboratory
Abstract
We investigate the approximate dynamics of several dif-ferential equations when the solutions are restricted to asparse subset of a given basis. The restriction is enforcedat every time step by simply applying soft thresholdingto the coefficients of the basis approximation. By reduc-ing or compressing the information needed to representthe solution at every step, only the essential dynamics arerepresented. In many cases, there are natural bases de-rived from the differential equations which promote spar-sity. We find that our method successfully reduces thedynamics of convection equations, diffusion equations,weak shocks, and vorticity equations with high frequencysource terms.
In this work, we investigate the approximate dynamics ofvarious PDEs whose solutions exhibit behaviors on multi-ple spatial scales. These scales may interact with one an-other in a non-linear manner as they evolve. Many physi-cal equations contain multiscale (as well as multiphysics)phenomena, such as the homogenization problems frommaterial science and chemistry and multiscale systems inbiology, computational electrodynamics, fluid dynamics,and atmospheric and oceanic sciences. In some cases, thephysical laws used in the model can range from molecu-lar dynamics on the fine scale to classical mechanics onthe large scale. In other cases, the equations themselvescontain high wave number oscillations that separate intodiscrete scales, on top of the smooth underlying behaviorof the system.The main source of difficulty in multiscale computationis that accurate simulation of the system requires all phe- nomena to be fully resolved. The smaller spatial scalesinfluence the global solutions, thus they cannot be ignoredin the numerical computation. In some cases it is possi-ble to derive an analytic equation for the effect of smallscales on the solution (see [1, 13]). In practice, how-ever, it may not be possible to derive a simple expres-sion that represents the fine scale behavior. Many problemdependent methods have been proposed in the literature,while a few provide a general methodology for model-ing the macroscopic and microscopic processes that yieldmultiscale models. For example, some general meth-ods include the heterogenous multiscale method [7], theequation-free method [10], multiscale methods for ellip-tic problems [12], and the sparse transform method [4].For an overview of general multiscale approaches see [6].A key difference between our method and the methodsin [7,10,12] is that we are directly resolving all of the sig-nificant scales in the solution. By contrast, the methodsof [7, 10, 12] directly resolve only the coarse scales of thesolution, and they separately “reconstruct" the fine scalesolution (as well as its effect on the coarse scales).From the perspective of mathematics, multiscale meth-ods began with representation of a function using a globalbasis, such as Taylor series or Fourier series. More so-phisticated bases have appeared; for example, any oneof the many wavelet bases used in imaging and compu-tational physics. The key to the basis approximation isthat each basis element represents behavior on a specificscale, therefore the coefficients of the basis provide com-plete information about the underlying function. This isalso the principle behind multiresolution and decomposi-tion methods.As the methods of multiscale and multiphysics mod-eling developed over the past few decades, so did cor-responding methods in imaging and information science.One of the fundamental ideas in imaging is that of spar-1ity. Sparse data representation is used throughout imag-ing from compression to reconstruction. Early advancesin sparse techniques were, e.g., in [3, 5], which presenteda convex minimization approach to the computationallychallenging sparse basis pursuit problem. Many mod-els have been proposed which use sparsity to produceboth more efficient numerical methods and better qual-ity solutions. Some applications of sparsity to imaginginclude: compressive sensing, reconstruction of imagesfrom sparse data [15, 16] , and recovery of images usingsparse regularization [9, 14]. The underlying principle ofsparsity is that images can be approximately representedby a small number of terms with respect to some basis.Inducing sparsity, creating effective bases, and developingefficient computational algorithms have been intensely ac-tive fields in information science.For imaging and information science, one of the rea-sons for the success of sparse methods is their abilityto resolve drastically different phenomena with a smallamount of information. This is also a principal goal ofmultiscale modeling. In this work, we transfer sparsitymethodology, which was developed for information sci-ence, to multiscale nonlinear differential equations andshow that it can be an effective tool for accurately com-puting solutions using less information.In particular, we propose solving PDEs with the con-straint that the approximate solution resides on a sparsesubspace of a basis. In this way, the complexity of themethod will depend on the number of basis terms retainedand ( nearly ) independent of the grid size. In the followingsections we will discuss the general problem and the op-timization method used to induce sparsity in the solution.Also, the general numerical method will be explained aswell as results of numerical experiments. The method istested on an advection equation with oscillatory velocity,a parabolic equation with oscillatory coefficients, a con-servation law with oscillatory diffusion, and the vorticityequations with high frequency source terms. We concludewith a discussion on the proposed work and implications.
In general, the problem can be stated as follows. Assum-ing x ∈ R n and t > , let u ( x, t ) : Ω → R be theapproximate solution of ∂u∂t = F ( u ) (1)subject to the constraint: u ( x, t ) = P j b u j ( t ) φ j ( x ) ,where the number of non-zero c j at a given time step is sparse. The operator F ( · ) can be non-linear, non-local,and dependent on the derivative of u . The basis terms φ j are assumed to exist on separate scales, which is true ofmost bases (for example, Legendre polynomials, Fourier,Wavelets, etc. ). In this way, the basis terms represent dif-ferent global behaviors. The method involves two steps:evolve the PDE forward in time and project the updatedsolution onto a sparse subset. We first address sparsityinduction through soft thresholding. At a given time step, the problem of projecting the up-dated solution onto a sparse subset is equivalent to fittinga solution u n with corresponding coefficients { b u nj } at the n th time step to a solution u whose corresponding coef-ficients { b u j } are sparse. This can be written as a con-strained least squares fit as follows. min u || u − u n || L s.t. (2) u = X j b u j φ j & { b u j } is sparseExpanding with respect to the basis and assuming that thebasis is orthonormal, this constrained optimization prob-lem is related to the following unconstrained problem:(L0) b u = argmin b u λ || b u || + 12 || b u − b u n || L (3)where b u is the vector of coefficients. The “norm” || · || is the number of non-zero coefficients in equation 2. Thismakes equation 3 both non-convex and difficult to solve.By replacing the L norm by the L norm, we get thefollowing convex relaxation of equation 3:(L1) b u = argmin b u λ || b u || + 12 || b u − b u n || L (4)Note that since b u ∈ C , the L norm is || b u || L := P j | b u j | ,where | b u j | = qb u j b u j . The solution of equation 4 is givenby the following equation: b u j = S λ (cid:0)b u nj (cid:1) = max (cid:0) | b u nj | − λ, (cid:1) b u nj | b u nj | (5)In general, this can be computed for a non-orthonormalbasis, which is equivalent to a basis pursuit problem with2he L norm as a sparse regularizer. In that case, the so-lution must be found by an iterative method rather thanthe simple shrinkage provided here as an example. Theresulting minimizer b u is a proximal solution which lies ona sparse subset of the original coefficient domain [2]. Thiscan be used to show that the solutions form a contractionmap in the L norm. Alternatively, we can simply applythe soft thresholding on the coefficients directly in orderto induce sparsity in this way. Assuming u ( x, t ) is periodic in the domain Ω ⊂ R n , onenatural basis is the Fourier basis, whose coefficients arethe Fourier transform of u ( x, t ) . This is appropriate forthe examples shown here. For the rest of this work, wewill use the Fourier basis; however, the overall method-ology presented here is independent of the correspondingbasis.Taking the Fourier transform of the PDE from equa-tion 1 and discretizing the resulting differential equationin time yields a multistep scheme. Since our method doesnot depend on the choice of numerical updating, we canassume that the scheme takes the following form: b v = Q (cid:0)b u n − q , ..., b u n (cid:1) (6)The updated solution b v may be sparse depending on boththe PDE and the update operator Q , but in general willhave non-trivial values everywhere depending on the ap-proximation and implementation. The auxiliary variable v is projected onto a sparse subspace by the shrinkage op-erator: b u n +1 = S λ ( b v ) (7)Altogether, the update in the spatial domain is simply: u n +1 = X j S λ (cid:0) Q (cid:0)b u n − q , ..., b u n (cid:1)(cid:1) φ j (8)Unlike traditional projections, this is non-linear and adap-tive. Rather than sorting the coefficients and retaining afixed number of large amplitude terms or keeping termswhose wavenumbers are below some cutoff, the shrink-age allows the number and choice of non-zero coefficientsto evolve over time. Also, this is not the same as hardthresholding the solution at every step ( i.e. keeping onlythe terms larger than a fixed value) since the coefficientsthat remain have decreased their magnitude by λ . Mostimportantly, the projection does not favor any particularpart of the spectrum; instead the amplitude of the coef-ficient determines if it remains. In terms of the Fourier basis, the importance is placed on the amplitude not thewavenumber.For general convergence, as long as λ = Cdt p for p larger than the accuracy of the scheme used to update thevariable in time, then the shrinkage operation does notchange the spatial accuracy of the original method and themethod will still converge as dt → . For example, dis-cretize using the forward Euler method, and then expandthe shrinkage operator to get b u n +1 = b u n + dt \ F ( u n ) + O ( λ ) . (9)Therefore, to have convergence as dt → , the shrinkageparameter must be λ = Cdt α . In general, the shrink-age operator is non-expansive in each coefficient, hencenon-expansive in coefficient norms. This may help withobtaining a general convergence result. In this section, we discuss the application of the proposedsparse method to several equations with different numeri-cal schemes.
The convection equation we consider is the following: ∂ t u = a ( x ) ∂ x u (10)where the coefficient a ( x ) is highly oscillatory.Let k be the wavenumber and use spectral Leap Frog asthe updating for equation 10 to obtain b v n +1 = b u n − + 2 dt b a ∗ ( i k b u n ) (11)in which ∗ is the convolution operator over frequency.The time step is O ( dx ) to preserve the stability condi-tion in equation 10. In Figure 1, the coefficient is chosenas follows. a ( x ) = 14 exp (cid:18) . . x )1 + 0 . x ) (cid:19) This choice of a ( x ) exhibits both fast and slow modes, butthe particular structure is not directly needed.Figure 1 illustrates the performance of the sparse solu-tion method on this example by comparison of the sparsesolution, the true solution produced using a standard fullyresolved method, and a “low frequency solution" pro-duced by solving equation 11 for wavenumbers k in theinterval | k | ≤ K in which K is the number of modes in3 a) True (black) and Sparse (blue) Solution in x space (b) True (black) and Sparse (blue ’x’) Solution in x space,zoomed in(c) Sparse (blue) Solution and low frequency (red) Solu-tion in x space (d) True (black) Solution and Sparse (blue) Solution in b u domain Figure 1: Convection with highly oscillatory coefficients. The solution is shown in x space, zoomed in version in x space, and in the b u space. The solutions are shown on a 512 grid, with dt = 2 e − , dx = 1 . e − , and λ = 5 e − .sparse solution. In Figure 1 (a-b) the sparse solution pro-duced by our method and the true solution at t = 1 areplotted in the spatial domain at a given time. In Figure 1(c) the sparse and low frequency solutions are displayed.The low frequency solution contains parasitic modes com-mon to the Leap frog scheme. In Figure 1 (d) the sparseand true spectra are plotted. The sparse spectrum capturesthe largest amplitude coefficients throughout the domain.In fact, out of the 512 coefficients used in the true solu-tion, only 27 are retained in the sparse one (about . ). The parabolic equation we consider is the following: ∂ t u = ∂ x ( a ( x ) ∂ x u ) (12)where the diffusion coefficient a ( x ) is highly oscillatory.The coefficient is assumed to be bounded, i.e. A M ≥ a ( x ) ≥ A m > . This is also related to the ellipticcase ∂ x ( a ( x ) ∂ x u ) = f , since an elliptic equation canbe solved by taking a parabolic scheme to steady state.Alternatively, the corresponding parabolic scheme can beiterated forward for a small number of time steps in or-der to find a partial solution to the elliptic problem. Thenby using the partial solution, the locations of the non-zerocoefficients can be extracted and the elliptic problem canbe solved by a Galerkin method on these coefficients [4].The updated scheme we use for equation 12 is forwardEuler. b v n +1 = b u n + dt i k b a ∗ ( i k b u n ) (13)The time step is O ( dx ) to preserve the stability con-dition as well as the highly oscillatory nature of the coef-ficient a ( x ) in equation 12. In Figure 2, the coefficient is4 a) True (black) and Sparse (blue) Solution in x space (b) True (black) and Sparse (blue ’x’) Solution in x space,zoomed in(c) True (black) Solution and Sparse (blue) Solution in b u domain (d) Diffusion coefficient in x space Figure 2: Parabolic diffusion with highly oscillatory coefficients. The solution is shown in x space, zoomed in versionin x space, and in the b u space. The solutions are shown on a 2048 grid, with dt = 1 . e − , dx = 3 . e − , and λ = 2 . e − .chosen as follows. a ( x ) = 110 exp (cid:18) . . x )1 + 0 . x ) (cid:19) In Figure 2 (d) the highly oscillatory diffusion coeffi-cient a ( x ) is plotted in space. In Figure 2 (a-b) the sparsesolution produced by our method and the true solution at t = 1 are plotted in the spatial domain at a given timeand are nearly indistinguishable. The high frequency in-formation is near the scale of the grid size, which can beseen in the zoomed in plot. In Figure 2 (c) the true andsparse spectra are displayed. The sparse spectrum cap-tures the largest coefficients throughout the domain andnot just the low wavenumbers. In fact, out of the 2048coefficients used in the true solution, only 53 are retainedin the sparse one (about . ). In time, this number ofnon-zero coefficients as well as the identities of the non-zero coefficients will change in order to capture variousbehaviors. To investigate the sparse dynamics of conservative lawswith diffusion, we use the viscous Burgers type equation. ∂ t u + 12 ∂ x (cid:0) u (cid:1) = ∂ x ( a ( x ) ∂ x u ) (14)The LHS of equation 14 is the standard Burgers advectionterm and the RHS is diffusion related to equation 12. Theequation exhibits three separate phenomenon: 1. Smoothlarge scale behavior from the diffusion, 2. Small scaleoscillations induced from the high frequencies in the co-efficient a ( x ) , and 3. Non-linear interactions between fre-quencies from the advection term. The update scheme intime is the standard total variational diminishing Runge-Kutta 2. b u = b u n + dt i k (cid:16) b a ∗ ( i k b u n ) − \ F ( u n ) (cid:17) v n +1 = 12 ( b u n + b u ) + dt i k (cid:16) b a ∗ ( i k b u ) − \ F ( u ) (cid:17) F ( u ) = u . As before, we have the stability con-dition dt is O ( dx ) .For Figure 3, the diffusion coefficient is chosen as: a ( x ) = 0 .
075 exp (cid:18) .
65 + 0 . x )1 + 0 . x ) (cid:19) The convolutions in the diffusion and non-linear terms aredone in the spectral domain, rather than by other meth-ods such as the psuedo-spectral method. The various dy-namics can be seen in the spatial and in the spectral plots(see Figure 3). The true, sparse, and low frequency solu-tions are plotted in the space in Figure 3(a-b). The lowfrequency projection is done by thresholding any coeffi-cients outside of a particular range. Specifically, the num-ber of low wavenumbers retained is the same as the sparsesolution, although their identities are dramatically differ-ent. The sparse solution captures the local and global be-haviors of the solution more accurately than the low fre-quency projected solution. In Figure 3 (c-d), the spectrumof the true solution is compared to the sparse and low fre-quency spectra, respectively. The local peaks in the spec-tra are related to the wavenumbers in the diffusion coef-ficient a ( x ) and the harmonics induced by the non-linearadvection term. Notice that in this case, each of the dis-tributions in the spectral domain do not decay as rapidlyas in the parabolic case. The sparse solution contains 130coefficients out of a total possible 1024, about . . The vorticity equation we consider is derived from the (2-d) incompressible Navier-Stokes equation [11]: ∂ t u + ( ∇ ⊥ ∆ − u ) · ∇ u = γ ∆ u + f (15)where u is the vorticity (not the velocity). Similarlyto equation 14, equation 15 exhibits three separate phe-nomenon: 1. Smoothness from the diffusion term on theRHS, 2. Small scale oscillations induced from the highfrequencies in the source term f , and 3. Non-linear inter-actions between frequencies from the advection term onthe LHS. However, since the operator ∇ ⊥ ∆ − is smooth-ing (in some sense), the advection term can be viewed asless non-linear than the one found in equation 14. In termsof the numerical method, the operator ∇ ⊥ ∆ − dampensthe coefficients by a factor which acts as | k | − .For the numerical implementation, the diffusion termis discretized using Crank-Nicolson while the advectionterm is lagged. Since the operators are diagonalized in thecoefficient basis, the steps can be invertible and lead to asimple updating scheme. b v n +1 = 2 dt γdt | k | (cid:16) ik ⊥ (cid:0) | k | − u n (cid:1) ∗ ik b u n + b f (cid:17) + 2 − γdt | k | γdt | k | b u n For Figure 4, the forcing term is chosen to be: f ( x, y ) = 0 .
025 sin(32 x ) + sin(32 y )1 + 0 .
25 (cos(64 x ) + cos(64 y )) The standard stability condition is used for choosing thetime steps in order to insure capture of all small scale be-haviors. In Figure 4 (a-b), the true and sparse solutionsare plotted in the spatial domain. Notice the oscillationsintroduced by the source term interact with the two vor-tex patches and thus contribute to the global behavior ofthe solution. The spectra of the true and sparse solutionis plotted in Figure 4 (c-d), where the low wavenumbersare located in the middle of the domain. The sparse solu-tion retains about . of the coefficients. In the sparsespectrum, the coefficients are located throughout the do-main, including the highest frequency itself (seen on theboundary of the spectral domain). In Figure 4 (e), the L and L ∞ error are shown to decrease as the resolution in-creases. This sparse solutions, as well as the other exam-ples presented here, converge as the spatial discretizationgoes to zero.It was observed that as the dimension increases, thesparsity of the solution also increases (since it is propor-tional to the product of the sparsities in each dimension).Thus the method scales well with dimension.It is worth noting that in a related work, wavelet hardthresholding was used to separate coherent and incoherentstructures in turbulent flows [8]. In the examples from the previous section, the PDEs con-tained a mixture of multiscale properties with a diffusionterm. The combination of non-linear and oscillatory termscreated a large range of wavenumbers in the solution,while the dynamics produced a range of amplitudes. Thisgives the necessary structure for sparsity with respect tothe Fourier basis. In general, the highest order derivativewill determine the appropriate basis in which the solutionscould be sparse.If the spectrum is more localized, i.e. non-zero regionsin the low frequency regime, then the proposed model canbetter condition the numerical method. Empirically, the6 a) True (black), Sparse (blue), and Low Frequency(red) Solution in x space (b) True (black), Sparse (blue), and Low Frequency(red) Solution in x space, zoomed in(c) True (black) and Sparse (blue) solutions in b u do-main (d) True (black) and Low Frequency (red) solutions in b u domain Figure 3: Viscous Burgers equation with highly oscillatory diffusion. The solution is shown in x space, zoomed inversion in x space, and in the b u space. The solutions are shown on a 1024 grid, with dt = 7 . e − , dx = 6 . e − ,and λ = 6 . e − .shrinkage operator acts as a non-linear filter on the coef-ficients. It was observed that for a fixed C and p , where λ = Cdt p , the numerical updates presented here with dt larger than theoretically and numerically possible in theoriginal scheme will converge. In the case of the vorticityexample, dt can be taken much larger when soft thresh-olding is applied than in the standard scheme. Also, thenon-linear filter seems to reduce parasitic modes and spu-rious oscillations found in spectral approximations for lin-ear and for non-linear slightly viscous hyperbolic equa-tions (see Figures 1 (b) and (c)) and 3 (g).One key point is that our method works by fully re-solving the solution. Its efficiency is gained by omittingmodes that are insignificant. This requires that λ is smallenough that the filter does not annihilate essential fea-tures. For example, if the initial data is smaller than λ for a particular unstable mode, then our approximate so-lution will not match the true dynamics. As the grid isrefined, the mode will be captured (since λ decreases as ∆ x decreases).In terms of complexity, each iteration is dominated by the convolution step. The convolution in the coef-ficient domain (spectral domain) can be done explicitlyover the n s ( t ) -sparse vectors rather than transformingback onto the spatial grid, which is O (cid:0) n s ( t ) (cid:1) at eachstep. When n s ( t ) << N log N , convolving in the spec-tral domain rather than transforming back and forth be-tween domains decreases the computational cost of thealgorithm. Knowing a priori the maximum sparsity, i.e. n s,max = max t n s ( t ) , faster routines and transformscould be optimized for specific problems and applications.For example, one can optimize the transform between thespatial and coefficient domains knowing the given sparsityat the current step and the non-trivial coefficients’ identi-ties. In the linear cases, as in equation 13, the operation b a ∗ can be stored as a large but sparse matrix, reducingthe updates to a sparse matrix - sparse vector operationat every iteration. Our goal in this work is to formulatea PDE solver that promotes sparsity. In future work wewill present a study of the computational complexity andspeed.When the dynamics are dominated by a linear term, for7 a) True Solution in x space (b) Sparse Solution in x space (c) Low Frequency Solution in x space(d) True Solution in b u domain (e) Sparse Solution in b u space (f) L (solid) and L ∞ (dashed) Errors verse gridstep size Figure 4: Vorticity equations with high frequency source term. The solution is shown in x space, zoomed in version in x space, and in the b u space. The solutions are shown on a 256 by 256 grid, with dt = 0 . , dx = 0 . , γ = 0 . ,and λ = 0 . .example high viscosity, the identities of the non-trivialcoefficients settle over time. This was also observed inthe non-linear cases, but over a longer time period. Thiswould enable creation of a sparse basis for elliptic equa-tions (e.g., with oscillatory coefficients).In many of the cases here, it is possible to get hyper-sparse solutions (those with 1% or fewer coefficients) atthe cost of accuracy. In this work, we have proposed a method to fully resolvethe solutions of multiscale PDEs while only retaining im-portant modes. The reduced dynamics created by thesparse projection properly captures the true phenomenaexhibited by the solution. The sparse projection amountsto a shrinkage of the coefficients of the updated solution atevery time step. There exist many possibilities for usingthe sparsity to optimize individual algorithms and create faster, more efficient computational routines.
The research of H. Schaeffer was supported by the De-partment of Defense (DoD) through the National DefenseScience and Engineering Graduate Fellowship (NDSEG)Program. The research of S. Osher was supported byONR: N00014-11-1-719. The research of R. Caflisch wassupported by DOE: DE-FG02-05ER25710. The researchof C. Hauck is sponsored by the Office of Advanced Sci-entific Computing Research; U.S. Department of Energy.The work was performed at the Oak Ridge National Labo-ratory, which is managed by UT-Battelle, LLC under Con-tract No. De-AC05-00OR22725. Accordingly, the U.S.Government retains a non-exclusive, royalty-free licenseto publish or reproduce the published form of this con-tribution, or allow others to do so, for U.S. Governmentpurposes.8 eferences [1] V. I. Arnol’d.
Mathematical Methods of ClassicalMechanics . Springer, May 1989.[2] J.F. Cai, S. Osher, and Z. Shen. Split Bregman meth-ods and frame based image restoration.
Multiscalemodeling & simulation , 8(2):337–369, 2009.[3] E.J. Cand„es, J. Romberg, and T. Tao. Robust uncer-tainty principles: Exact signal reconstruction fromhighly incomplete frequency information.
Informa-tion Theory, IEEE Transactions on , 52(2):489–509,2006.[4] I. Daubechies, O. Runborg, and J. Zou. Asparse spectral method for homogenization multi-scale problems.
Multiscale Modeling & Simulation ,6(3):711–740, January 2007.[5] D.L. Donoho. Compressed sensing.
InformationTheory, IEEE Transactions on , 52(4):1289–1306,2006.[6] W. E and B. Engquist. Multiscale modelingand computation.
Notices Amer. Math. Soc ,50(50):1062–1070, 2003.[7] W. E, B. Engquist, and Z. Huang. Heterogeneousmultiscale method: a general methodology for mul-tiscale modeling.
Physical Review B , 67(9):092101,2003.[8] M. Farge and K. Schneider. Coherent vortex simu-lation (cvs), a semi-deterministic turbulence modelusing wavelets.
Flow, Turbulence and Combustion ,66(4):393–426, 2001.[9] T. Goldstein and S. Osher. The split Bregmanmethod for L1-regularized problems.
SIAM Journalon Imaging Sciences , 2(2):323–343, 2009.[10] I.G. Kevrekidis, C.W. Gear, J.M. Hyman, P.G.Kevrekidid, O. Runborg, and C. Theodoropoulos.Equation-free, coarse-grained multiscale computa-tion: Enabling mocroscopic simulators to performsystem-level analysis.
Communications in Mathe-matical Sciences , 1(4):715–762, 2003.[11] A.J. Majda and A.L. Bertozzi.
Vorticity and incom-pressible flow , volume 27. Cambridge UniversityPress, 2001.[12] J. Nolen, G. Papanicolaou, and O. Pironneau. Aframework for adaptive multiscale methods for el-liptic problems.
Multiscale Modeling & Simulation ,7(1):171–196, 2008. [13] G. Papanicolau, A. Bensoussan, and J.-L. Lions.
Asymptotic Analysis for Periodic Structures . Else-vier, January 1978.[14] L.I. Rudin, S. Osher, and E. Fatemi. Nonlinear totalvariation based noise removal algorithms.
PhysicaD: Nonlinear Phenomena , 60(1):259–268, 1992.[15] H. Schaeffer and S. Osher. A low patch-rank inter-pretation of texture.