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
(1) | ||||
(2) | ||||
(3) |
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 , :
(4) |
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
(5) |
for the grid points in a distance to the boundary. The symmetric difference approximation leads to the symmetry condition and yields the MOL-equation
(6) |
In a similar way, the Robin boundary condition is approximated by the equation
(7) |
which may be solved for by
(8) |
Thus, may be eliminated from the equation (5) with yielding
(9) |
with
(10) |
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:(11) | ||||
(12) |
where is the -th unit vector. The initial conditions are simple evaluations of the function on the grid,
(13) |
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
(14) |
where are the first non-negative solutions of the equation
(15) |
with the convention that , for . The corresponding normalized eigenvectors have the components
(16) |
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).
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
(17) |
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
(18) |
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
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
(20) |
with . Optima may be computed by using some multiplier function for the ODE restriction (11) and considering the Lagrangian
(21) |
The partial derivatives of the Lagrangian with respect to and recover (11), (13) and the other ones are
(22) | ||||
(23) | ||||
(24) |
Hence, the Karush-Kuhn-Tucker conditions, in (22)–(24), show that the control may be eliminated by
(25) |
and a necessary condition for the optimal solution is that it solves the following boundary value problem:
(26) | ||||
(27) |
The homogeneous differential equation (27) for has the simple solution
(28) |
with the matrix exponential
(29) |
With given by (25), this solution may be used in (19) to yield the solution for (26) with the coefficient functions
(30) |
Considering the function satisfying for and , the integral may be written as
(31) |
The result (30) may be used in different ways. The first application is computing the solution for a given target profile .
Lemma 2.3
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,
(34) |
with coefficients belonging to some reasonable form of the control . Then, by the boundary condition in (27), the corresponding target profile has the form
(35) |
where, by (30), the coefficients of are given by
(36) |
We will use this construction in our numerical example.
![]() |
![]() |
![]() |
![]() |
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].
![]() |
![]() |
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.
References
-
[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.