Assignment 1: Ill-posedness and Regularization

Computational and Variational Methods for Inverse Problems

CSE 397/GEO 391/ME 397/ORI 397

Assignment 1: Ill-posedness and Regularization due September 20, 2021 In this assignment we will explore the nature of ill-posedness for inverse problems involving the three canonical equations of mathematical physics: the heat equation (parabolic), Poisson’s equation (elliptic), and the wave equation (hyperbolic). We will also study regularization methods for mitigating ill-posedness. For the heat equation inverse problem we discussed in class, you will derive expressions for the eigenvalues and eigenvectors of the discretized parameter-to-observable map and compare them to the eigenvalues and eigenfunctions for the continuous operator we saw in class. Next you will repeat the analysis in the Jupyter notebooks we covered in class, but now for Poisson’s equation (with an elasticity interpretation). Finally, there is an optional question (for extra credit) on analyzing the ill-posedness of an initial condition inverse problem for the wave equation. Problem 1 Recall the inverse problem for the 1D heat equation discussed in class and in the accompanying Jupyter notebook 1.Ill posedness.ipynb1 in which we seek to reconstruct the initial temperature field from the observed temperature field at a later time T. As we have seen, this problem is severely ill-posed (it is called exponentially ill-posed, due to the exponential decay of the eigenvalues). In class we stated expressions for the eigenvalues/eigenfunctions of the continuous operator F governing this problem. In addition to confirming these are correct, in this assignment you will derive the eigenvalues/eigenvectors for the discrete operator F, which helps us understand the role of discretization in ill-posedness of the inverse problem. Let u(x, t) denote the temperature field and u(x, 0) := m(x) the initial temperature profile. Given the length of the rod L, the thermal diffusivity k > 0, the final time T, and homogeneous Dirichlet boundary conditions (u(0, t) = u(L, t) = 0), the parameter-to-observable map F(m) can be written as F(m) := u(x, T), where for a given m(x), u(x, T) denotes the solution at the observation time T of the heat equation    ∂u ∂t − k ∂ 2u ∂x2 = 0 0 < x < L, 0 < t ≤ T, u(x, 0) = m(x) 0 < x < L, u(0, t) = u(L, t) = 0 0 < t ≤ T. We discretize the problem using 2nd order centered finite differences in space and implicit Euler in time. We denote with h = L nx the mesh size and with ∆t = T nt the time step used in the discretization. We would like to study the nature of the ill-posedness of the discretized inverse problem Fm = d, where m is an nx − 1 vector representing the discretized parameter field we wish to invert for (i.e., the initial condition), and d is the corresponding data vector containing the temperature observations at the grid points. The (nx −1)×(nx −1) matrix F, arising from the discretization of the parameter-to-observable map F, takes the form F = (I + ∆tK) −nt , 1The Jupyter notebook is available for downloading on Canvas under Class 2 at 1.Ill posedness.ipynb where the discretized diffusion operator K is an (nx − 1) × (nx − 1) matrix given by K = k h 2         2 −1 −1 2 −1 −1 2 −1 . . . . . . . . . −1 2 −1 −1 2         , a) Confirm that the eigenvalues λi and corresponding eigenfunctions vi(x) of the continuous operator F are given by λi = e −kT( πi L ) 2 , vi(x) = r 2 L sin  πi x L  , i = 1, 2, . . . In other words, show that F(vi) = λivi . b) Derive expressions for the eigenvalues and eigenvectors of the discrete parameter-to-observable operator F using the fact that the eigenvalues µi of K are given by µi = k 4 h 2 sin2  πi 2nx  , i = 1, 2, . . . , nx − 1, and the jth component (j = 1, 2, . . . , nx − 1) of the corresponding eigenvector ui is given by [ui ]j = r 2 L sin  πi jh L  , i = 1, 2, . . . , nx − 1. c) Set T = 1, L = 1, nx = 100, and nt = 100, and plot the decay of the eigenvalues of the continuous (F) parameter-to-observable operator as a function of i. Do this for the following values of k: 0.0001, 0.001, 0.01, 0.1 (all on the same plot). Explain the behavior you see in the plot. [When plotting here and in Problem 1d, set the lower limit of the vertical scale to 10−16.] d) Set L = 1, k = 0.01 and T = 0.1, and plot the decay of the discrete eigenvalues as a function of i for different space and time discretizations, and compare to the exact eigenvalue spectrum. Use (nx, nt) = (20, 20),(40, 40),(80, 80),(160, 160). What do you observe as you increase the resolution? Comment on the regularizing nature of discretization. Problem 2 Here we wish to study ill-posedness and regularization for an inverse problem governed by Poisson’s equation, specifically the axial deformation of an elastic rod of length L. The inverse problem is to infer the distributed body force (per unit length) m(x), from observations of the axial displacement u(x) everywhere along the length of the rod (even this problem is ill-posed!). Given m(x), we can find u(x) by solving the forward elasticity problem  −k ∂ 2u ∂x2 = m(x) 0 < x < L, u(0) = u(L) = 0, where k is the elastic modulus of the rod, and homogeneous Dirichlet boundary conditions (i.e., the rod is fixed) are imposed at both ends. The parameter-to-observable map can thus be written as F(m) := u(x) for a given m(x). a) Derive expressions for the eigenvalues λi and eigenfunctions vi(x) of the continuous operator F. Hint: the eigenfunctions are the same as those for the heat inverse problem above (but the eigenvalues are not!). b) Discretize the scaled Laplacian −kuxx using the same central finite difference method as for the heat inverse problem above, so that K is defined exactly as above. The discrete parameter-toobservable operator can thus be written as F := K−1 . Compare the eigenvalues of F (making use of your knowledge of the eigenvalues of K given in Problem 1b) with those of F for this problem by plotting the decay of the eigenvalues for both operators as a function of i, for the case k = L = 1. For the discrete operator, take nx = 100. c) We shall attempt to solve the inverse problem without regularization directly from m = F −1d. To generate synthetic data, let’s take the true body force to be2 mtrue = max(0, 1 − |1 − 4x|) + 100 x 10(1 − x) 2 . To generate the data, add i.i.d. normally distributed noise3 η with mean zero and standard deviation σ = 10−4 . The resulting noisy observation of the displacement field is d = Fm + η. Experiment with several different choices of noise, mesh size nx, and elastic modulus k (similar to the experiments we did with the heat equation inverse problem in TikhonovRegularization.ipynb) to study the stability of the inversion. d) Using Tikhonov regularization with α = 10−7 , 10−6 , 10−5 , 10−4 , 10−3 , 10−2 , and for L = 1, k = 1, and nx = 200, compute the regularized least squares solution mα for the data generated in part (c) above. Which value of α gives the best inversion in the “eyeball” norm? e) Determine the (approximate) optimal value of the regularization parameter α in the Tikhonov regularization according to the L-curve criterion. f) Determine the (approximate) optimal value of the regularization parameter α in the Tikhonov regularization according to Morozov’s discrepancy criterion, i.e., find the largest value of α such that kFmα − dk ≤ δ where δ = kηk and mα is the solution of the Tikhonov-regularized inverse problem with regularization parameter α. g) Plot the L2 norm error in the reconstruction, kmtrue −mαk, as a function of α, where mα is the Tikhonov regularized solution. Which value of α (approximately) minimizes this error? Compare the “optimal” values of α obtained in parts (d), (e), and (f) of this problem and comment on any differences. Problem 3 (Optional; for extra credit) Here we will study the stability of the inverse problem for the wave equation through a model 1D problem. The forward problem is to find the transverse displacement u(x, t) of a cable of length L, tension k, 2This is the same function we used in the Jupyter notebook for the heat inverse problem TikhonovRegularization.ipynb. Feel free to use modify this notebook to answer parts (c)–(g) of this problem. 3Python provides the function randn through the numpy.random library. and mass density ρ, fixed at both ends (u(0, t) = u(L, t) = 0), that is initially at rest (u˙(x, 0) = 0) and is plucked with an initial displacement of u(x, 0) = m(x). In other words, given m(x), solve    ∂u2 ∂t2 − c 2 ∂ 2u ∂x2 = 0 0 < x < L, 0 < t ≤ T, u(x, 0) = m(x) 0 < x < L, ∂u ∂t (x, 0) = 0 0 < x < L, u(0, t) = u(L, t) = 0 0 < t ≤ T for the displacement field u(x, t), where c := p k/ρ is the wave propagation speed. The goal of the inverse problem is to infer the initial displacement m(x) from observation of the cable’s position u(x, T) at a later time t = T. Thus the parameter-to-observable map is given by F(m) := u(x, T). a) Derive expressions for the eigenvalues λi and eigenfunctions vi(x) of the continuous operator F as a function of L, c, T, i. What do these expressions tell you about the ill-posedness (i.e., stability) of the inverse problem? How would you choose the observation time T to provide for the perfect inference of m(x)? b) What would happen to your ability to reconstruct m(x) if you were limited to observing at T = L 2c ? Can you think of some other quantity to observe that would address this problem? c) In the real world (i.e., for a real cable), what might limit your ability to reconstruct m(x)?