Discretization Issues

In order to have a multigrid solver for the full optimization problem with high efficiency, one has to guarantee that certain operators on the discrete level are h-elliptic. In particular, the Hessian which is a symmetric positive definite operator should have a discretization which is h-elliptic. The Hessian may be a pseudo differential operator and its symbol may no longer be a trigonometric polynomial. So we have to change slightly the definition of h-ellipticity to include that case. Thus,

$$\mbox{\tt h-ellipticity:} \qquad \hat {\cal H}^h (\theta) \geq C \sum _j \vert sin (\theta _j /2) \vert ^m$$ (15)

$$\mbox{\tt quasi-ellipticity:} \qquad \hat {\cal H}^h (\theta) \geq C \sum _j \vert sin (\theta _j) \vert ^m$$ (16)

To demonstrate the difficulty that may arise with discretization we consider the following example.

Example I: A Control Problem Let $\phi$ be the state variable, $\alpha$ the design variable and consider the minimization problem,

$$\min _{\alpha} \frac{1}{2} \int _{\partial\Omega} (\frac{\partial \phi}{\partial n}-d)^2 ds$$ (17)

where $\phi$ satisfies

$$\begin{array}{lr} \Delta \phi = 0 & \Omega \\ \phi = \alpha & \partial\Omega \end{array}$$ (18)

and where $\Omega = \{ (x,y) \vert y >0 \}$. A simple calculation shows that the necessary (optimality) conditions are given by

$$\begin{array}{lr} \Delta \phi = 0 & \Omega \\ \Delta \lambda = 0... ...Omega \\ \frac{\partial \lambda}{\partial n}= 0 & \partial \Omega. \end{array}$$ (19)

Following a calculation as in a previous lecture we obtain that the symbol of the Hessian is given by

$$\hat {\cal H} ({\bf k} ) = \vert {\bf k} \vert ^ 2,$$ (20)

which is an elliptic symbol. Consider now the following discretizations.

Discretization I: h-elliptic Hessian. A uniform grid with vertex discretization is used with the usual 5-point Laplacian

$$L^h = \frac{1}{h^2} \left[ \begin{array}{rrr} & 1 & \\ 1 & -4 & 1 \\ & 1 & \end{array}\right].$$ (21)

The discrete approximation to $\phi$ are parameterized as $\phi _{l,k}$ where $(l,k)$ is a grid point. For discretization of the normal derivative at the boundary we use a ghost point with a central discretization

$$\frac{\partial \phi}{\partial n}\approx \frac{1}{2h} ( \phi _{l,1} - \phi _{l,-1} )$$ (22)

where $k=0$ stands for the boundary, and $k=-1$ represent ghost points. The above formula is used in discretization of the cost functional. The actual implementation is to add one equation for the ghost point, by approximating the interior equation also on the boundary. Thus, relating the ghost values to other values in the interior.

The design variables $\alpha$ are taken at the grid vertices. A Fourier change in $\alpha$ by $\exp (i \theta l)$ introduce the following changes

$$\phi _{l,j} = \exp (i \theta l) d(\theta)^j$$ (23)

$$\lambda _{l,j} = A \exp (i \theta l) d(\theta)^j$$ (24)

where $d(\theta) \leq 1$ is the solution of

$$d(\theta) + 1/d(\theta) - 4 + 2 cos(\theta) = 0.$$ (25)

Note that if $d (\theta )$ is a solution, then so is $1/ d (\theta )$. The first correspond to a bounded solution in half space while the other is unbounded. Define

$$\mu (\theta ) = \frac{1}{2h} [ 1/d( \theta) - d ( \theta ) ]$$ (26)

and observe that $\mu (\theta )$ is the symbol for normal derivative at the boundary. A simple calculation, as was done previously, shows that the symbol for the discrete Hessian is

$$\hat {\cal H}^h (\theta) = \mu ^2 (\theta ).$$ (27)

To see that this is indeed h-elliptic notice that $d( \pi ) = (6 - \sqrt{32})/2$, and thus it is bounded away from zero.

Discretization II: Quasi-elliptic Hessian. Here we consider a cell-centered scheme, i.e., the variables approximating $\phi$ are located at the center of the cells and the control variables $\alpha _l$ are located at the boundary grid points. The discretization for the interior points is the same as before. The boundary condition also here uses ghost points with a central discretization for Neumann boundary condition. If we denote the grid points by $(l,k)$ then the $\phi$ variable are $\phi _{l+\frac{1}{2}, k+\frac{1}{2}}$ The contribution of $\alpha$ to the boundary condition at a general point in the bottom boundary is

$$\frac{1}{2} ( \phi _{l+\frac{1}{2} ,\frac{1}{2} } + \phi _{l+\frac{1}{2} ,-\frac{1}{2} } ) = \frac{1}{2} ( \alpha _{l} + \alpha _{l+1} )$$ (28)

where $k=0$ are the boundary points and $\phi _{l, -1/2}$ are the ghost variables. The Neumann derivative is approximated by

$$(\frac{\partial \phi}{\partial n})_{l+1/2,0} \approx \frac{1}{h} ( \phi _{l+1/2,1/2} - \phi _{l+1/2,-1/2} )$$ (29)

and the corresponding symbol is

$$\mu _c (\theta ) = \frac{1}{h} ( \frac{1}{\sqrt{d(\theta )}} - \sqrt{d(\theta )} ).$$ (30)

Defining

$$C(\theta ) = \frac{1}{\sqrt{ d(\theta ) }} +\sqrt{ d(\theta )}$$ (31)

the solution is given by

$$\phi _{l,j} = \frac{\cos( \theta /2)}{C( \theta ) } \exp (i \theta l) d(\theta)^j$$ (32)

$$\lambda _{i,j} = \mu _c( \theta ) \frac{\cos( \theta /2)}{C( \theta ) } \exp (i \theta l) d(\theta)^j.$$ (33)

Again a simple calculation shows that symbol of the Hessian is

$$\hat {\cal H} (\theta ) = \mu _c^2 ( \theta ) \frac{\cos( \theta /2)}{C( \theta ) }.$$ (34)

Since $cos ( \theta /2)$ vanishes at $\pi$ the symbol is not h-elliptic, but only quasi-elliptic. This means that no local relaxation of the design variables will result in smoothing. That is, high frequency errors in the design variable will not be damped fast. A standard multigrid method will not work for this discretization. In this example one could guess the quasi-ellipticity of the Hessian from the fact that perturbation in the design variable of the form $(-1)^l$ do not cause any changes in $\phi$.

Example II: Consider the optimization problem

$$\min _{\alpha} \frac{1}{2} \int _{\partial\Omega} (\phi -d)^2 ds$$ (35)

$$\vspace{2mm}$$ (36)

subject to the equation,

$$\begin{array}{lr} \Delta \phi = 0 & \Omega \\ \frac{\partial \phi}{\partial n}= \alpha _x & \partial \Omega \end{array}$$ (37)

where $\Omega = \{ (x,y) \vert y >0 \}$. A simple calculation shown that the necessary conditions are given by

$$\begin{array}{lr} \Delta \phi = 0 & \Omega \\ \Delta \lambda = 0... ...n}= \phi - d & \partial\Omega \\ - \lambda _x = 0 & \partial\Omega \end{array}$$ (38)

and the symbol of the Hessian is

$$\hat {\cal H} ({\bf k} ) = 1.$$ (39)

Discretization I: h-elliptic Hessian. We use here cell-centered discretization for the state and costate and the design variables are at the boundary nodes. The approximation for the tangential derivative is

$$( \frac{\partial \alpha}{\partial x} ) _{l+1/2} \approx \frac{1}{h} (\alpha _{l+1} - \alpha _l )$$ (40)

whose symbol is given by

$$s(\theta ) = \frac{2i}{h} \sin (\theta /2)$$ (41)

The normal derivative is the same as in discretization II of example I. A change in the design variable by a Fourier component $\exp (i \theta l)$ results in the following changes in state and costate

$$\begin{array}{ll} \phi _{l,k} = \frac{s(\theta ) } { \mu _c (\the... ...\frac{s(\theta )}{ \mu _c^2 (\theta )}\exp (i \theta l) d(\theta)^j \end{array}$$ (42)

and the symbol of the Hessian is

$$\hat {\cal H}^h (\theta ) = \frac{4 \sin ^2 {(\theta /2)}}{\mu _c ^2 (\theta )}.$$ (43)

We gave this example to show that although the symbol of the Hessian for the differential level was just the constant 1, on the discrete level the Hessian is more complicated. Taylor expansion of this symbol at $\theta = 0$ shows that it indeed approximate the PDE. Note however, that the Hessian has a condition number which is independent of h. The ratio of largest to smallest value that the symbol attains is independent of $h$. This means that this problem presents no special difficulties even for very large number of design variables.

Discretization II: Quasi-elliptic Hessian. Here we discretized the problem using cell-vertex variables for the state and adjoint variables. The design variables are given at the boundary nodes. A simple calculation shows that

$$\hat {\cal H}^h (\theta ) = \frac{\sin ^2 (\theta )}{\mu ^2 (\theta )}$$ (44)

and this is a quasi-elliptic symbol, due to the $sin ^2 (\theta)$ term, which vanishes at $\theta = \pi$.