Log In Sign Up

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.


page 1

page 2

page 3

page 4


A linear implicit Euler method for the finite element discretization of a controlled stochastic heat equation

We consider a numerical approximation of a linear quadratic control prob...

Numerical analysis of a Neumann boundary control problem with a stochastic parabolic equation

This paper analyzes the discretization of a Neumann boundary control pro...

Error Inhibiting Schemes for Initial Boundary Value Heat Equation

In this paper, we elaborate the analysis of some of the schemes which we...

Discrete Adjoint Implicit Peer Methods in Optimal Control

It is well known that in the first-discretize-then-optimize approach in ...

Sparse Representations of Solutions to a class of Random Boundary Value Problems

In the present study, we consider sparse representations of solutions to...

Adjoint-based optimization of two-dimensional Stefan problems

A range of optimization cases of two-dimensional Stefan problems, solved...

1 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 time-dependent 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.

2 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 governed by the heat equation


where and are given functions. The homogeneous Neumann condition at may be considered as a short-cut for space-symmetric solutions . The coefficients of the general Robin boundary condition are nonnegative, and nontrivial .

The equation (1) is approximated by finite differences with a shifted equi-spaced grid with stepsize , :


For the approximation of the boundary conditions, also the outside points and will be considered temporarily. In the method-of-lines approach with central differences, approximations are defined by the differential equations


for the grid points in a distance to the boundary. The symmetric difference approximation leads to the symmetry condition and yields the MOL-equation


In a similar way, the Robin boundary condition is approximated by the equation


which may be solved for by


Thus, may be eliminated from the equation (5) with yielding




Hence, we have for the Dirichlet condition and for the pure Neumann condition. Collecting all equations (5), (6) and (9

), the following MOL system for the vector

is obtained:


where is the -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

are known, which is well known for special values of , at least.

Lemma 2.1

For , the eigenvalues of the matrix from (12) are given by


where are the first non-negative solutions of the equation


with the convention that , for . The corresponding normalized eigenvectors have the components


with constants .

Proof: In the main equation (5), the ansatz gives

In the first equation, we have

with the same factor . In order to satisfy the eigenvalue condition in the last component, we consider the equation , i.e.,

since . The last grid point is and with the trigonometric formulas for and the identity , we may proceed with

Hence, the different versions of in (10) verify the condition (15) for . Rearranging (15) as , for it is seen that exactly solutions exist in since the function is monotonically decreasing and positive. Finally, the vector norms are computed for . Abbreviating and using , we get

which leads to the value of the normalizing factor in (16).   

For later use, we introduce the diagonal matrix

and the unitary matrix

satisfying .

Remark 2.1

The well-known frequencies for Dirichlet boundary conditions are and , for Neumann conditions. For general values the equation (15) may be rewritten in fixed point form


The functions are monotonically decreasing in and an iteration with initial value converges to the desired solution at least for , since .

2.2 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 leads to


Since the matrix is symmetric, the inner product with yields the decoupled equations , 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 indicates that simple sparse solutions with only a few terms in (19) may not exist due to the inhomogeneous boundary condition (2).

2.3 Exact solution of an optimal control problem

The inhomogeneity inherited from the boundary condition (2) may be considered as a control to approach a given target profile at some given time . In an optimal control context, controls are searched for minimizing an objective function like


with . Optima may be computed by using some multiplier function for the ODE restriction (11) and considering the Lagrangian


The partial derivatives of the Lagrangian with respect to and recover (11), (13) and the other ones are


Hence, the Karush-Kuhn-Tucker conditions, in (22)–(24), show that the control 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 has the simple solution


with the matrix exponential


With given by (25), this solution may be used in (19) to yield the solution for (26) with the coefficient functions


Considering the function satisfying for and , 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 be given. Then, the coefficient vector of the solution of the boundary value problem (26), (27) is given by the unique solution of the linear system


with the positive semi-definite matrix having the elements


Proof: At the end point , 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 with some vector , obtaining

This means that is semi-definite and 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 nontrivial basis coefficients in the state and the Lagrange multiplier . Due the special inhomogeneity in (26), sparse solutions for the state probably do not exist. However, by adjusting the target profile , one may simply start with a sparse multiplier with, for instance, two terms only,


with coefficients belonging to some reasonable form of the control . Then, by the boundary condition in (27), the corresponding target profile has the form


where, by (30), the coefficients of are given by


We will use this construction in our numerical example.

Figure 1: Dirichlet heat problem with spatial points and given exact control : (left) and (right).

3 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 consider the following ODE-constrained optimal control problem with an incorporated boundary control of Dirichlet type:

with , , , , state vector , and as defined in (12) with . We set in (34) and compute the target profile from (35) with coefficients for 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].

Figure 2: Dirichlet heat problem with spatial points, solved by gradient descent for :

Two test scenarios are considered. First, the accuracy of the numerical approximations for and are studied, where the exact control is used. The initial value for the multiplier is set to with being the approximation of with time step . 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 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 is approximated at the nodes , chosen by the time integrators on a time grid with step size and stages [2, 5]. In both test cases, we use and to also study the influence of the system size. The number of time steps are with .

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 and the adjoint . 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.


  • [1] J.T. Betts.

    Practical Methods for Optimal Control and Estimation Using Nonlinear Programming

    SIAM, 2010.
  • [2] W.W. Hager. Runge-Kutta methods in optimal control and the transformed adjoint system. Numer. Math., 87:247–282, 2000.
  • [3] E. Hairer, Ch. Lubich, and G. Wanner.

    Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations

    , volume 31 of Springer Series in Computational Mathematics.
    Springer, Heidelberg, Berlin, 2006.
  • [4] J. Lang and B.A. Schmitt. Discrete adjoint implicit peer methods in optimal control. J. Comput. Appl. Math., 416:114596, 2022.
  • [5] J. Lang and B.A. Schmitt. Implicit A-stable peer triplets for ODE constrained optimal control problems. Algorithms, 15:310, 2022.
  • [6] Ch. Lubich and A. Ostermann. Runge-Kutta approximation of quasi-linear parabolic equations. Math. Comp., 64:601–627, 1995.
  • [7] A. Ostermann and M. Roche. Runge-Kutta methods for partial differential equations and fractional orders of convergence. Math. Comp., 59:403–420, 1992.
  • [8] J.M. Sanz-Serna. Symplectic Runge–Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more. SIAM Review, 58:3–33, 2016.
  • [9] J.L. Troutman. Variational Calculus and Optimal Control. Springer, New York, 1996.