Parameter estimation for Euler equations with uncertain inputs Sergiy Zhuk Tigran Tchrakian Abstract— The paper presents a new state estimation algorithm for 2D incompressible Euler equations with periodic boundary conditions and uncertain but bounded inputs and initial conditions. The algorithm converges (in L2 -sense) to a least squares estimator given incomplete and noisy observations. It is also shown experimentally that the proposed algorithm applies to several conservation laws developing shock discontinuities. The results are illustrated by numerical examples. I. I NTRODUCTION Data assimilation algorithms represent a backbone of modern cyber-physical systems. These algorithms allow one to optimally combine a priori knowledge encoded in equations of mathematical physics with a posteriori information in the form of sensor readings. Weather forecasting is one of many examples where data assimilation is applied to improve predictions generated by hydrodynamical models. The divergence-free Euler equations provide the most basic model of incompressible flows of homogeneous inviscid fluids [7], yet this model may provide insights in studies of turbulence. We refer the reader to [2] for a detailed overview of various mathematical questions related to Euler equations. of the solution of 2D Euler equations globally (in time) [7, Cor.3.3, p.116]. Assuming periodic boundary conditions, we apply Fourier-Galerkin (FG) approximation: we project the vorticity equation onto a subspace generated by {einx }|n|≤N and obtain an ODE for the projection coefficients, a FG model. Note that Fourier-Galerkin approximation possesses a spectral convergence rate provided the solution of the Euler equation is smooth [1]. Although the main results of this work may be derived for the case of bounded domains with non-penetration boundary conditions or unbounded domains (R2 ), we focus on smooth periodic flows to simplify the convergence analysis. We refer the reader to Section II, where other types of boundary conditions and weak solutions of Euler equations are discussed. In this paper we design a new state estimator for 2D Euler equations subject to uncertain but bounded input and unknown initial conditions. We rely on a vorticity-stream formulation of the Euler equations. The resulting vorticity equation describes the rotation of the vorticity of the fluid velocity field (see Section II). We stress that the homogenous vorticity equation possesses a nice property: the L2 -norm of the vorticity is conserved over time. The latter, in fact, allows one to prove existence and uniqueness Given an ODE representing projection coefficients, we design a discrete-in-time FG model such that the L2 -norm of the discrete vorticity is conserved over time (see Section III-B). This property allows us to prove the convergence of the approximations provided by the discrete FG model to a solution of the continuous FG model. We then derive the state estimator (in the form of a minimax filter) for the discrete FG model and note that the combination of the mentioned convergence proof with results of [6] may be used to derive a continuous formulation of the minimax filter (see Section IIIC). To the best of our knowledge, this result is new. A similar approach has been used to design data assimilation algorithms for scalar macroscopic traffic flow models [9] and flood models (StVenant equations) [11]. We stress that the proposed method shows very good performance not only for smooth periodic flows but also for conservation laws with shock discontinuities and non-periodic boundary conditions (see Section IV). IBM Research, Dublin, {sergiy.zhuk,tigran}@ie.ibm.com Data assimilation for systems of hyperbolic conservation laws based on the calculus of variations Ireland, was proposed in [3], where the authors adopt the strategy: “optimize”, then “discretize” so that the estimate of the initial density does not depend on a discretization method. In contrast, the present paper solves the filtering problem; that is, the state estimate at time instant, t, depends on the observation at the same time instant and the previous estimate only. A comparison of classical filtering algorithms (extended Kalman filter and ensemble Kalman filter) for scalar conservation laws with quadratic non-linearity may be found in [4]. Adaptive parameter estimators for hyperbolic equations were considered for instance in [5]. where ū ∈ R2 is a given vector, ~u0 ∈ C 2 (Ω)2 is a 1 smooth 1-periodic R vector-function and f ∈ C (ΩT ) has zero mean, Ω f (x, t)dx = 0. It is not hard to prove (see [7, Prop 2.4,p.50]) that ~u = (u1 , u2 )> ∈ C 2 (ΩT )2 is a smooth 1-periodic solution of the incompressible Euler equation on ΩT : This paper is organized as follows. The next section presents a brief overview of Euler equations. Section II-A presents discrete in time FG model for Euler equations; the state estimation algorithm is shown in section III-C. Numerical experiments are given in section IV. Section V contains concluding remarks. We note that Euler equation (2) has the unique smooth solution ~u for Ω = R2 provided f~ = 0 and ~u0 is a smooth function such that div(~u0 ) = 0 and curl(~u0 ) ∈ L1 (R2 ) (see [7, L.3.2, p.93 and Cor.3.3, p.116]). For the case of bounded domains with smooth boundary one may also derive existence of smooth solutions for (2) provided ~u0 and f~ are smooth and ~u0 · ~n = 0 on ∂Ω (see [10, p.356]). Weak solutions (or Sobolev space solutions) for (2) were constructed in [12] provided | curl(~u0 )| is bounded, and f~ ∈ C(0, T, Lp (Ω)) is so that | curl(f~)| is bounded. Solutions of L2 (0, T, L2 (R2 )2 )-class corresponding to so called vortex sheets were discussed in [7, p.303]. Notation. Let Ω denote a subset of R2 with boundary ∂Ω, and ~n is a normal vector for ∂Ω poining outside, ΩT := Ω × [0, T ]. C s (Ω) denotes a space of continuously differentiable functions on Ω (up to order s), L2 (Ω) is the space of square-integrable functions on Ω, H 1 (Ω) is a Sobolev space of functions with weak derivatives of L2 (Ω)-class. div(~u) = ∂x1 u1 +∂x2 u2 , curl(~u) = ∂x1 u2 −∂x2 u1 , ∇u = (∂x1 u, ∂x2 u)> , ∇⊥ u = (−∂x2 u, ∂x1 u)> ~u · ~v denotes the canonical inner product of R2 , H 2 := H × H denotes the cartesian product of H withR itself. L2 (t0 , t1 , H) := {f : f (t) ∈ T H and 0 kf (t)k2H dt < +∞}. We write u = v a.e. on Ω if u(x) = v(x) for almost all x ∈ Ω. d~u + (~u · ∇)~u + ∇p = f~ , dt div(~u) = 0 , ~u(0) = ~u0 + ū , (2) where f~ ∈ C 2 (ΩT ) and curl(f~) = f , and the pressure p is a function of ~u and f~. A. Fourier-Galerkin approximation Define Z (~u · ∇w)vdx . b(~u, w, v) := Ω Assume that div(~u) = 0 and ~u, w, v are smooth 1-periodic functions on Ω. We find integrating by parts: II. E ULER EQUATIONS b(~u, w, v) = −b(~u, v, w) . Assume that ω verifies the vorticity-stream formulation of the Euler equation: Now, we apply Fourier spectral method to construct a finite dimensional approximation of (1). Let s ∈ iπr −1 s·x Z2 and define φs (x) := e 2r . It is known that {φs }s∈Z2 is a total orthonormal system in L2 (Ω) and (φs , φk )L2 (Ω) = δs1 ,k1 δs2 ,k2 where si , ki are components of s, k. Now, we multiply (1) by φs ∂t ω + ~u · ∇ω = f , −∆ψ = ω , ⊥ ~u = ū + ∇ ψ , ω(0) = curl(~u0 ) , 2 (x, t) ∈ ΩT := [−r, r] × [0, T ] , (1) (3) and integrate over Ω. We obtain a weak vorticity formulation: d (ω, φs )L2 (Ω) + b(~u, ω, φs ) = (f, φs )L2 (Ω) . dt P Let us now set ω N (x, t) := |ki |≤N ωk (t)φk (x) with ωk (t) := (ω(t), φk )L2 (Ω) and plug ω N into the above formulation. We get the following ODE for the coefficients: X ω̇s + b(~u, φk , φs )ωk (t) = (f, φs )L2 (Ω) . |ki |≤N Define ω ~ (t) := {ωs (t)}|s1,2 |≤N , f~(t) {(f (t), φs )L2 (Ω) }|s1,2 |≤N and set := B(~ ω ) := {b(~u, φk , φs )}|s1,2 |,|k1,2 |≤N . Now we use this notation to rewrite the finite dimensional vorticity formulation in the vector form: d~ ω + B(~ ω )~ ω = f~ with initial condition ω ~ (0) = dt ω ~ 0 where ω ~ 0 := {ωs (0)}|s1,2 |≤N . We stress that B(~ ω ) = −B > (~ ω ) as b(~u, φk , φs ) = −b(~u, φs , φk ) by (3) so that B(w) ~ is skew-symmetric. A. Problem statement Let ω ~ solve ω ~ (0) = ω ~0 , kω N (·, t) − ω(·, t)kL2 (Ω) ≤ eg(t) N −2l kω(·, 0)k2H l (Ω) + eg(t) N 2−l max kω(·, τ )kH l (Ω) , g(t) ≥ 0 . 0≤τ ≤t (7) This so called spectral convergence rate justifies our choice of the state equation. Namely, (4) is a standard finite dimensional Fourier-Galerkin model with uncertain input f~. Now, following the idea of [14], [15] we may incorporate the effect of the unresolved modes (the projection error) by simply adding another model error term e and introduce an additional algebraic equation to filter out inadmissible e (see [8], [13]) so that the exact projection coefficients of ω will be among the solutions of (4). On the other hand, the above convergence rate estimate suggests that the effect of e is negligible for reasonably large N and l (so for smooth ω) and it may be therefore absorbed by f~ (by increasing the size of the ellipsoid (6)). B. Discrete Fourier-Galerkin model III. M AIN RESULTS d~ ω + B(~ ω )~ ω = f~ , dt Let us stress that if ω ~ is the solution of (4), ω N 1 corresponds to ω ~ and ω solves (1) then according to [1, T.5.1] we have: (4) and assume that a vector-function ~y is observed in the following form (|k1,2 | ≤ M ): X ~yk (t) = (Hk , φs )L2 (Ω) ωs (t) + ηk (t) , (5) |s1,2 |≤N where Hk ∈ L2 (Ω) and ~η = {ηk }|k1,2 |≤M is a measurable vector-function modelling noise. Our aim is to construct a state estimate ω̂(T ) for ω ~ (T ) given data ~y and assuming that Z T ω ~ 0> S~ ω0 + f~> Qf~ + ~η > R~η dt ≤ 1 (6) 0 provided S, Q, R are symmetric positive definite matrices of appropriate dimensions. Assume f = 0. If we multiply (1) by ω and integrate over Ω we get: 21 ∂t kωk2L2 (Ω) = (f, ω)L2 (Ω) as b(~u, ω, ω) = 0 by (3). Hence, L2 (Ω)-norm of ω is conserved. We stress that (4) has the same property, namely k~ ω k2R2N +1 is conserved as B(~ ω ) is skewsymmetric. In what follows we propose a method which approximates ω ~ (jh) by w̃j := w̃(jh), j = T 0, m, h := m , m ∈ N and the norm of w̃j is conserved. By Newton-Leibniz formula we get: Z (j+1)h ω ~ ((j + 1)h) = ω ~ (jh) − B(~ ω (τ ))~ ω (τ )dτ jh Define Bj := B(~ ω (jh)). Approximating the integral by the trapezoidal rule one gets: h h ω ((j+1)h) = (I− Bj )~ ω (jh)+O(h3 ) . (I+ Bj+1 )~ 2 2 1 ω N (x, t) P := |ki |≤N ωk (t)φk (x) where the coefficients ωk are components of ω ~ d~ ω By noting that Bj+1 = Bj + hB( (jh)) + O(h2 ) dt we can simplify the above equation compromising the order of the approximation: specifically, approxd~ ω imating Bj+1 by Bj + hB( (jh)) we reduce the dt order down to O(h2 ) (locally); if we simply take Bj we get O(h)-approximation. In what follows we stick to the latter and define w̃j as a solution of the linear system: (I + h h Bj )w̃j+1 = (I − Bj )w̃j , w̃0 = ω ~ 0 . (8) 2 2 Note that (I + symmetric and h 2 Bj ) Kj := (I + is invertible as Bj is skewh h Bj )−1 (I − Bj ) 2 2 is the Caley transform of the skew-symmetric matrix Bj . Hence, Kj is an orthogonal matrix and so kw̃j k2R2N +1 = kw̃0 k2R2N +1 . For j = 1, . . . , m − 1 and jh ≤ t < (j + 1)h we define the following functions: w̃j+1 +w̃j U (m) (t) := , V (m) (t) = w̃j and 2 (m) W (t) = w̃j − (t − jh)B(w̃j )U (m) (t). Since kw̃j k2R2N +1 = kw̃0 k2R2N +1 it follows that kV (m) kL2 (0,T ) , kU (m) kL2 (0,T ) ≤ C1 < +∞ and so kW (m) kL2 (0,T ) ≤ C2 < +∞. Therefore, the sequences of piecewise constant functions, {V (m) }m∈N , {U (m) }m∈N contain weakly convergent subsequences in L2 (0, T ). The same holds true for the sequence {W (m) }m∈N . We denote the convergent subsequences by {V (mi ) }, {U (mi ) } and {W (mi ) } respectively and let V ∗ , U ∗ and W ∗ be their limiting functions. We claim that all the three mentioned sequences converge strongly and V ∗ = U ∗ = W ∗ . dW (mi ) Indeed, = −B(V (mi ) (t))U (mi ) (t) and dt dW (mi ) so is bounded in L2 (0, T ). Hence, dt {W (mi ) } weakly converges in H 1 (0, T ) that implies strong convergence in L2 (0, T ). Now, (mi ) (mi ) 3 2 1 2 1 2 k ≤ h m kw̃0 k2R2N +1 C2 P2N +1 2 j=1 kB(ej )k , ej is j-th kV −U where C := canonical basis L2 (0,T ) vector. Hence, by taking weak limit (m → ∞) and using weak lowersemicontinuity of L2 -norm we get: V ∗ = U ∗ . On the other hand kW (mi ) − V (mi ) k2L2 (0,T ) ≤ Pm h3 (mi ) (t)k2R2N +1 ≤ j=1 kB(w̃j )U 3 3 h 4 3 Ckw̃0 kR2N +1 . By the same argument we get: W ∗ = V ∗ . Taking weak limits in dW (mi ) = −B(V (mi ) (t))U (mi ) (t) we find dt dW ∗ that: = −B(W ∗ )W ∗ (t). dt Now, recalling the spectral convergence rate estimate (7) and noting that there exists unique smooth solution of (1) we deduce that any convergent subsequence of U (m) , V (m) and W (m) has the same limit. The latter proves that the entire sequences U (m) , V (m) and W (m) are weakly convergent and share the same limiting function W ∗ which is the unique solution of (4). Let us summarize the above results. Proposition 1: If f = 0 then (i) (4) has a unique solution ω ~ such that ω (t)kR2N +1 , k~ ω (0)kR2N +1 = k~ (ii) a sequence of piecewise constant functions V (m) (t) = w̃(jh), jh ≤ t < (j + 1)h converges to ω ~ in L2 (0, T ) provided w̃(jh) T solves (8) and h := m , m ∈ N. We note that it is not hard to generalize point (ii) of the latter proposition to the case f 6= 0. We omit this generalization here for the sake of space. C. State estimator Following [6] we introduce the following system of linear Hamiltonian equations: > I− h 2 Bj h −1 2Q > h 2 H RH I+ h 2 Bj Uj+1 Vj+1 = > I+ h 2 Bj h −1 2Q > h 2 H RH I− h 2 Bj (9) where Pj := Uj−1 Vj for j > 0 and P0 = S −1 , and Bj := B(ω̂j ), and ω̂j solves the following system I Pj , of linear equations: (I + h h Bj + Pj+1 H > RH)ω̂j+1 2 2 h h (10) = (I − Bj − Pj H > RH)ω̂j 2 2 h Pj+1 H > R~yj+1 + Pj H > R~yj + , 2 2 (a) Estimate at t = 0.05 (b) Truth at t = 0.05 (c) Estimate at t = 1.00 (d) Truth at t = 1.00 (e) Estimate at t = 3.00 (f) Truth at t = 3.00 where H := {(Hk , φs )L2 (Ω) }|k1,2 |,|s1,2 |≤N . By combining the idea of the proof of Proposition 1 with a well-known stabilization effect brought by Pj , the solution of approximated Riccati equation (see [6] for the details), one may prove the following proposition. Proposition 2: Let Uj , Vj , Pj and ω̂j be defined T , m ∈ N. A sequence as above and h := m of piecewise constant functions V (m) (t) = Vj , U (m) (t) = Uj , ω̂ (m) (t) = ω̂j for jh ≤ t < (j +1)h converges to U , V and ω̂ in L2 (0, T ). Moreover, U , V and ω̂ represent the unique solution of the following system of equations: dω̂ = −B(ω̂)ω̂ + P H > RH(~y − H ω̂) , ω̂(0) = 0 , dt dU = B > (ω̂)U + H > RHV , U (0) = I , dt dV = Q−1 U − B(ω̂)V , V (0) = S −1 , P = V U −1 . dt (11) In fact, the latter proposition presents a continuous version of the minimax filter (equation (11)) and equations (9)-(10) define an approximation for the filter. In the following section we illustrate the convergence properties of this approximation on numerical examples. IV. N UMERICAL EXPERIMENTS A. Euler equations For the divergence-free Euler system, we look for solutions in the space, span{φk (x) := φk1 (x)φk2 (y)}k1 =0...Nx ,k2 =0...Ny , where ψk1 (z) = sin(k1 πz/L), L = 2r is the length of each side of the spatial domain, and Nx and Ny are the numbers of basis functions in each of the spatial Fig. 1: Estimates and truth for different times dimensions. The choice of basis satisfies the boundary conditions, ω = 0 on ∂Ω. The corresponding Fourier-Galerkin model (4) is subject to the perturbation vector f~ which has only two non-zero entries, meaning that in effect, we perturb only two modes: one low, and one high frequency mode. For the filter, we first generate observations ~y by running (4) the Euler system forward in timewith initial condition represented by a translation and scaling of Gaussian distributions. The filter is initialized to zero, meaning we assume no knowledge of the initial condition. Observations are taken on an evenly spaced 15 × 15 grid away from the boundary, where the state is fixed at zero. For both the generation of observations and the filtering, we choose Nx = Ny = 60. Figure 1 shows the estimate and truth at different times. The sparsity of the observations is apparent from Figure 1a, n=1, m=1 n=2, m=2 1.4 1.4 Truth Estimate time = 1.00 Truth Estimate 1.2 1.2 1.2 Estimate Truth 1 Observations Spectral coefficient Spectral coefficient 1 0.8 0.6 0.8 Perturbed observations 1 0.6 0.4 0.8 0.4 0.2 0.2 0 0 0.5 1 1.5 2 2.5 Time 3 3.5 4 4.5 −0.2 5 u 0.6 0 0 0.5 1 1.5 2 2.5 Time 3 3.5 4 4.5 5 0.4 (a) Projection coefficients for (b) Projection coefficients for k1 = k2 = 1 k1 = k2 = 2 n=3, m=3 0.2 1.4 0.6 Truth Estimate 0 1.2 Relative L2−norm of estimation error 0.4 Spectral coefficient 0.2 0 −0.2 −0.4 0 1 2 3 4 5 6 x 0.8 0.6 Fig. 3: True traffic state versus estimate at t=1 0.4 0.2 −0.6 −0.8 1 0 0 0.5 1 1.5 2 2.5 Time 3 3.5 4 4.5 5 0 0.5 1 1.5 Time 2 2.5 3 time = 3.00 (c) Projection coefficients for(d) Relative L2 -error of the k1 = k2 = 3 estimate 1.2 Estimate Truth Observations 1 Perturbed observations 0.8 B. Lighthill-Whitham-Richards model (LWR) In this section we apply the above method to scalar conservation laws. Recall from [9] that the standard equilibrium traffic-flow model, LWR model consists of a scalar conservation law, ∂t u(x, t) + ∂x f (u(x, t)) = 0, (12) with initial data u0 (x) = u(x, 0) where u : R × R+ → R is the traffic density, x ∈ R and t ∈ R+ are the independent variables, space and time respectively, and f : R → R is the flux function. A typical flux function is f (u) = uVm 1 − uum where the constants, Vm and um , are the maximum speed and the maximum density respectively. We impose periodic boundary conditions on the interval (0, 2π). Unlike the Euler equation presented above, 0.6 u which is from an early point in the estimation. In Figure 1c, we see that by t = 1.00, the estimate does still not fully mimic the truth, but that the flow is being captured. In Figure 1e, the estimate appears to capture the truth quite well. Figures 2a - 2c show the estimated and true projection coefficients for some of the low wave numbers. 0.4 0.2 0 0 1 2 3 4 5 6 x Fig. 4: True traffic state versus estimate at t=3 this model develops shock discontinuities even subject to periodic boundary conditions and smooth initial condition. The filter has been applied to the Fourier-Galerkin model2 which has been used to generate observations. The filter has been initialized to 0. The estimation results are presented on the figures 3-5. We refer the reader to [9] for the further details. C. StVenant equations This final example shows that the proposed state estimator may handle systems of conservation laws 2 The model had artificial viscosity term activated on higher order modes to allow for correct shock tracking time = 6.00 1.2 Estimate Truth Observations 1 Perturbed observations 0.8 u 0.6 0.4 0.2 0 0 1 2 3 4 5 Fig. 7: h, hu and the estimate after 20 time steps 6 x Fig. 5: True traffic state versus estimate at t=6 Fig. 8: h, hu and the estimate after 120 time steps Fig. 6: Initial conditions for h, hu and the estimate V. C ONCLUSION with non-periodic boundary conditions. The standard equilibrium flood model consists of a system of scalar conservation laws: ∂t h + ∂x (hu) = 0 , ∂t (hu) + ∂x (hu2 + gh2 )=0 2 (13) with boundary conditions u(0, t) = ul (t) and h(0, t) = hl (t) on (0, 1), where h is the fluid depth, u is the averaged velocity and g is the gravitational constant. The Discountinuos Galerkin method has been applied to the above equation to generate observations. There was no perturbation, so f~ = 0. The initial conditions for the model and filter are presented at Figure 6. The estimation results are presented on the figures 7-8. We refer the reader to [11] for the further details. The paper presents the discrete and continuous versions of the minimax state estimator for Euler equations. The estimator is derived for a FourierGalerkin model which approximates smooth solutions of the Euler equation with spectral convergence rate. A very curious topic for the future research is to develop the idea of [9] for weak solutions of Euler equations, that is to combine the presented state estimation approach and vanishing viscosity method [1] to design a state estimator for L∞ solutions of Euler equations in vorticitystream formulation. We stress that the presented method and convergence results apply to NavierStokes equations in dimension 2 without major modifications. R EFERENCES [1] C. Bardos and E. Tadmor. Stability and spectral convergence of fourier method for nonlinear problems: on the shortcomings of the 2/3 de-aliasing method. Numerische Mathematik, 129(4), 2015. [2] C. Bardos and E. Titi. Euler equations of incompressible ideal fluids. Uspekhi Matematicheskikh Nauk, 62(3):5–46, 2007. [3] Claude Bardos and Olivier Pironneau. Data assimilation for conservation laws. Methods and Applications of Analysis, 12(2):103–134, 06 2005. [4] S. Blandin, A. Couque, A. Bayen, and D. Work. On sequential data assimilation for scalar macroscopic traffic flow models. Physica D, 2012. [5] M. Demetriou and H. Banks. Adaptive parameter estimation of hyperbolic distributed parameter systems: Nonsymmetric damping and slowly time varying systems. ESAIM: Control, Optimisation and Calculus of Variations, 3:133–162, 1998. [6] J. Frank and S. Zhuk. Symplectic möbius integrators for lq optimal control problems. In Proc. IEEE Conference on Decision and Control, 2014. [7] A. Majda and A. Bertozzi. Vorticity and incompressible flow. Cambridge Univ. Press, 2002. [8] V. Mallet and S. Zhuk. Reduced minimax filtering by means of differential-algebraic equations. In 5th Int. Conf. on Physics and Control, 2011. lib.physcon.ru. [9] T. Tchrakian and S. Zhuk. A macroscopic traffic data assimilation framework based on Fourier-Galerkin method and minimax estimation. IEEE Transactions on Intelligent Transportation Systems, (99):1–13, 2014. special issue. [10] R. Temam. Navier-Stokes equations: Theory and Numerical Analysis. AMS Chelsea Publishing, 2001. [11] S. Tirupathi, S. Zhuk, T. Tchrakian, and S. McKenna. Data assimilation for 1d shallow water equations: a minimax approach. in preparation, 2015. [12] V. Yudovich. Non-stationary flow of an incompressible liquid. Zh. Vychisl. Mat. Mat. Fiz., pages 1032–1066, 1963. [13] S. Zhuk. Kalman duality principle for a class of illposed minimax control problems with linear differentialalgebraic constraints. Applied Mathematics and Optimisation, 68(2):289–309, 2013. [14] S. Zhuk. Minimax projection method for linear evolution equations. In Proc. IEEE Conference on Decision and Control, ieeexplorer.ieee.org, 2013. [15] S. Zhuk, J. Frank, I. Herlin, and R. Shorten. Data assimilation for linear parabolic equations: minimax projection method. SIAM J. Sci. Comp., 2015. to appear.