A Stiff MOL Boundary Control Problem for the 1D Heat Equation with Exact Discrete Solution

Method-of-lines discretizations are demanding test problems for stiff integration methods. However, for PDE problems with known analytic solution the presence of space discretization errors or the need to use codes to compute reference solutions may limit the validity of numerical test results. To overcome these drawbacks we present in this short note a simple test problem with boundary control, a situation where one-step methods may suffer from order reduction. We derive exact formulas for the solution of an optimal boundary control problem governed by a one dimensional discrete heat equation and an objective function that measures the distance of the final state from the target and the control costs. This analytical setting is used to compare the numerically observed convergence orders for selected implicit Runge-Kutta and Peer two-step methods of classical order four which are suitable for optimal control problems.


Introduction
One main area of application for stiff integration methods is semi-discretizations in space of time-dependent partial differential equations in the method-of-lines approach.In order to test new methods in this area one may rely on PDE test problems with known analytic solution or reference solutions computed by other numerical methods.However, both approaches have its drawbacks.For PDE problems the accuracy is limited by the level of space discretization errors and in computing reference solutions one has to trust the reliability of the used code.This background was our motivation to develop the current test problems with exact discrete solutions for a finite difference semi-discretization in space with arbitrarily fine grids.
It is known that, in contrast to multi-step type methods, one-step methods may suffer from order reduction if applied to MOL systems especially with timedependent boundary conditions, see [6,7].Motivated by our recent work [4,5] on Peer two-step methods in optimal control the present example is formulated as a problem with boundary control.
The paper is organized as follows.In Section 2, we apply a finite difference discretization with a shifted equi-spaced grid for the 1D heat equation with general Robin boundary conditions and derive exact formulas for the solutions of the discrete heat equation and an optimal boundary control problem.These analytical solutions are used in a sparse setting to study the numerically observed convergence orders in Section 3 for several one-step and two-step integration methods which are suitable for optimal control.

A discrete heat equation with boundary control 2.1 Finite difference discretization of the 1D heat equation
We consider the initial-boundary-value problem for a function Y (x, t) governed by the heat equation where Ψ(x) and u(t) are given functions.The homogeneous Neumann condition at x = 0 may be considered as a short-cut for space-symmetric solutions Y (−x, t) ≡ Y (x, t).The coefficients of the general Robin boundary condition are nonnegative, β 0 , β 1 ≥ 0 and nontrivial (β 0 , β 1 ) = (0, 0).The equation ( 1) is approximated by finite differences with a shifted equi-spaced grid with stepsize ξ = 1/m, m ∈ N: For the approximation of the boundary conditions, also the outside points x 0 = −ξ/2 and x m+1 = 1 + ξ/2 will be considered temporarily.In the method-of-lines approach with central differences, approximations y j (t), j = 1, . . ., m are defined by the differential equations for the grid points in a distance to the boundary.The symmetric difference approximation 0 != ξY x (0, t) ∼ = (y 1 − y 0 ) leads to the symmetry condition y 0 ≡ y 1 and yields the MOL-equation In a similar way, the Robin boundary condition is approximated by the equation which may be solved for y m+1 by Thus, y m+1 may be eliminated from the equation ( 5) with j = m yielding with Hence, we have θ = 3 for the Dirichlet condition and θ = 1 for the pure Neumann condition.Collecting all equations (5), ( 6) and ( 9), the following MOL system for the vector y(t) = y j (t) j=1,...,m is obtained: where e m is the m-th unit vector.The initial conditions are simple evaluations of the function Ψ on the grid, The basis of our construction is that the eigenvalues and eigenvectors of the symmetric matrix M are known, which is well known for special values of θ, at least.Lemma 2.1 For m ≥ 2, the eigenvalues of the matrix M ∈ R m×m from (12) are given by where ω k , k = 1, . . ., m, are the m first non-negative solutions of the equation Proof: In the main equation ( 5), the ansatz v = e iωx j m j=1 gives In the first equation, we have with the same factor λ := −(4/ξ 2 ) sin 2 ωξ/2 .In order to satisfy the eigenvalue condition in the last component, we consider the equation 0 since cos(ωξ) = 1 − 2 sin 2 (ωξ/2).The last grid point is x m = 1 − ξ/2 and with the trigonometric formulas for cos(ω − ωξ/2), sin(ω − ωξ/2) and the identity sin(ωξ) = 2 sin(ωξ/2) cos(ωξ/2), we may proceed with Hence, the different versions of θ in (10) verify the condition (15) for m = 1/ξ.Rearranging (15) as tan(ω) = β 0 /(2mβ 1 ) cot(ω/(2m)), for β 0 > 0 it is seen that exactly m solutions exist in (0, mπ) since the function ω → cot(ω/(2m)) is monotonically decreasing and positive.Finally, the vector norms are computed for ω = 0.
Abbreviating ω/m =: Ω and using cos which leads to the value of the normalizing factor ν in (16).
The functions f k are monotonically decreasing in ω and an iteration with initial value ω = max{1, (k − 1)π} converges to the desired solution ω k at least for 2m >

Exact solution of the discrete heat equation
Knowing the eigenvectors and eigenvalues of the linear problem (11), the computation of its solution is straight-forward.The representation y(t) = k=1,...,m η k (t)v [k]  leads to Since the matrix M is symmetric, the inner product with v [j] yields the decoupled equations η j (t) = λ j η j (t) + γv m u(t), which can be solved easily leading to the following result.
Lemma 2.2 With the data from Lemma 2.1, the solution of the initial value problem (11), ( 13) is given by Remark 2.2 The presence of the terms v m indicates that simple sparse solutions with only a few terms in (19) may not exist due to the inhomogeneous boundary condition (2).

Exact solution of an optimal control problem
The inhomogeneity u(t) inherited from the boundary condition (2) may be considered as a control to approach a given target profile ŷ ∈ R m at some given time T > 0. In an optimal control context, controls are searched for minimizing an objective function like with α > 0. Optima may be computed by using some multiplier function p(t) for the ODE restriction (11) and considering the Lagrangian The partial derivatives of the Lagrangian L with respect to p(t) and p(T ) recover (11), (13) and the other ones are Hence, the Karush-Kuhn-Tucker conditions, ∂ (•) L = 0 in ( 22)-( 24), show that the control u(t) may be eliminated by and a necessary condition for the optimal solution is that it solves the following boundary value problem: The homogeneous differential equation ( 27) for p has the simple solution with the matrix exponential With u given by ( 25), this solution may be used in (19) to yield the solution for (26) with the coefficient functions Considering the function ϕ 1 (z) = 1 0 e zt dt satisfying ϕ 1 (z) = (e z − 1)/z for z = 0 and ϕ 1 (0) = 1, the integral may be written as The result (30) may be used in different ways.The first application is computing the solution for a given target profile ŷ.
Lemma 2.3 Let ŷ ∈ R m be given.Then, the coefficient vector η(T ) of the solution y(t) = V η(t) of the boundary value problem (26), ( 27) is given by the unique solution of the linear system with the positive semi-definite matrix Q = (q k ) m k, =1 having the elements Proof: At the end point T , the formula (30) simplifies to This equation may be reordered to the form given in (32) with the matrix elements (33).Finally, we consider the quadratic form of the matrix Q with some vector w = (w j ), obtaining This means that Q is semi-definite and I + Q definite and the system (32) always has a unique solution.
In general, solutions computed with (32) will not be sparse, i.e., they will have m nontrivial basis coefficients in the state y and the Lagrange multiplier p. Due the special inhomogeneity in (26), sparse solutions for the state y probably do not exist.However, by adjusting the target profile ŷ, one may simply start with a sparse multiplier p(t) with, for instance, two terms only, p(t) = δ 1 e λ 1 (T −t) v [1] + δ 2 e λ 2 (T −t) v [2] , with coefficients δ 1 , δ 2 belonging to some reasonable form of the control u.Then, by the boundary condition in (27), the corresponding target profile has the form where, by (30), the coefficients of y(T ) are given by We will use this construction in our numerical example.

Test case: Dirichlet boundary control problem
To illustrate an application of the derived expressions for the exact discrete solutions of the linear heat equation equipped with different boundary conditions, we T ∈ R m , state vector y(t) ∈ R m , and M as defined in (12) with θ = 3.We set δ 1 = δ 2 = 1/75 in (34) and compute the target profile ŷ ∈ R m from (35) with coefficients for y(T ) defined in (36).We will compare numerical results for four time integrators of classical order four: the symmetric 2-stage Gauss method [3,Table II.1.1],the symmetric 3-stage partitioned Runge-Kutta pair Lobatto IIIA-IIIB [3,Table II.2.2] and our recently developed two-step Peer methods AP4o43bdf and AP4o43dif [5].The two one-step methods are symplectic and therefore well suited for optimal control [3,8].Two  The initial value for the multiplier is set to p(T ) = y h (T ) − ŷ with y h (T ) being the approximation of y(T ) with time step h.In this case, the Karush-Kuhn-Tucker system decouples and only two systems of linear ODEs have to be solved.In the second scenario, the optimal control problem is solved for all unknowns (y, p, u) by a gradient method with line search as implemented in the Matlab routine fmincon, see e.g.[1,9] for more details, and the errors for the control are discussed.The control variable u(t) is approximated at the nodes t ni = t n + c i h, i = 1, . . ., s, chosen by the time integrators on a time grid {t 0 , . . ., t N } with step size h and s stages [2,5].In both test cases, we use m = 250 and m = 500 to also study the influence of the system size.The number of time steps are N = 2 k with k = 4, . . ., 11.In Fig. 1, results for the first test scenario are shown.Not surprisingly, the serious order reduction for the symplectic one-step Runge-Kutta methods are clearly seen.This phenomenon is well understood and occurs particularly drastically for time-dependent Dirichlet boundary conditions [7].This drawback is shared by all one-step methods due to their insufficient stage order.Note that the number of affected time steps increases when the system size is doubled.In contrast, the newly designed two-step Peer methods for optimal control problems work quite close to their theoretical order four for the state y and the adjoint p.The order reduction for the one-step methods is also visible for the more challenging fully coupled problem.The results plotted in Fig. 2 show a reduction to first order for the approximation of the control, whereas the two-step methods perform with order two for this problem.We refer to [2] for a discussion of the convergence order for general ODE constrained optimal control problems.Once again, the range of the affected time steps depends on the problem size.It increases for finer spatial discretizations.
Figure 2: Dirichlet heat problem with m = 250, 500 spatial points, solved by gradient descent for (y, p, u): max n,i |u(t ni )−u h (t ni )| test scenarios are considered.First, the accuracy of the numerical approximations for y(T ) and p(0) are studied, where the exact control u(t) = −γp m (t)/α is used.