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

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 $Y (x, t)$ governed by the heat equation

$\partial_{t} Y (x, t) =$	$\partial_{x x} Y (x, t), (x, t) \in [0, 1] \times [0, T],$	(1)
$\partial_{x} Y (0, t) =$	$0, β_{0} Y (1, t) + β_{1} \partial_{x} Y (1, t) = u (t),$	(2)
$Y (x, 0) =$	$Ψ (x),$	(3)

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) \equiv Y (x, t)$ . The coefficients of the general Robin boundary condition are nonnegative, $β_{0}, β_{1} \geq 0$ and nontrivial $(β_{0}, β_{1}) \neq (0, 0)$ .

The equation (1) is approximated by finite differences with a shifted equi-spaced grid with stepsize $ξ = 1 / m$ , $m \in N$ :

x_{j} = (j - \frac{1}{2}) ξ, j = 1, \dots, m .

(4)

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, \dots, m$ are defined by the differential equations

y_{j}^{'} = \frac{1}{ξ^{2}} (y_{j - 1} - 2 y_{j} + y_{j + 1}), j = 2, \dots, m - 1,

(5)

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} \equiv y_{1}$ and yields the MOL-equation

y_{1}^{'} = \frac{- y_{1} + y_{2}}{ξ^{2}} .

(6)

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

β_{0} \frac{y_{m} + y_{m + 1}}{2} + β_{1} \frac{y_{m + 1} - y_{m}}{ξ} = u (t),

(7)

which may be solved for $y_{m + 1}$ by

y_{m + 1} = \frac{2 β_{1} - β_{0} ξ}{2 β_{1} + β_{0} ξ} y_{m} + \frac{2 ξ}{2 β_{1} + β_{0} ξ} u (t) .

(8)

Thus, $y_{m + 1}$ may be eliminated from the equation (5) with $j = m$ yielding

y_{m}^{'} =

\frac{1}{ξ^{2}} (y_{m - 1} - θ y_{m}) + γ u (t)

(9)

with

θ = \frac{2 β_{1} + 3 β_{0} ξ}{2 β_{1} + β_{0} ξ} = 3 - \frac{4 β_{1}}{2 β_{1} + β_{0} ξ}, γ = \frac{2}{(2 β_{1} + β_{0} ξ) ξ} .

(10)

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, \dots, m}

is obtained:

	$y^{'} =$	$M y + γ e_{m} u (t),$		(11)
	$M =$	$\frac{1}{ξ^{2}} ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} - 1 & 1 1 & - 2 & 1 ⋱ & ⋱ & ⋱ 1 & - 2 & 1 1 & - θ \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠,$		(12)

where $e_{m}$ is the $m$ -th unit vector. The initial conditions are simple evaluations of the function $Ψ$ on the grid,

y (0) = ψ, ψ = {(Ψ (x_{j}))}_{j}^{m} .

(13)

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 \geq 2$ , the eigenvalues of the matrix $M \in R^{m \times m}$ from (12) are given by

λ_{k} = - 4 m^{2} {sin}^{2} (\frac{ω_{k}}{2 m}), k = 1, \dots, m,

(14)

where $ω_{k}, k = 1, \dots, m,$ are the $m$ first non-negative solutions of the equation

tan (ω) tan (\frac{ω}{2 m}) = \frac{β_{0}}{2 m β_{1}},

(15)

with the convention that $ω_{k} = (k - \frac{1}{2}) π$ , $k = 1, \dots, m,$ for $β_{1} = 0$ . The corresponding normalized eigenvectors $v^{[k]}$ have the components

v_{j}^{[k]} =

ν_{k} cos (ω_{k} \frac{2 j - 1}{2 m}), j = 1, \dots, m,

(16)

with constants $ν_{k} = 2 / \sqrt{2 m + sin (2 ω_{k}) / sin (ω_{k} / m)}$ .

Proof: In the main equation (5), the ansatz $v = (R e^{i ω x_{j}})_{j = 1}^{m}$ gives

\frac{1}{ξ^{2}} (v_{j - 1} - 2 v_{j} + v_{j + 1}) =

\frac{1}{ξ^{2}} R e^{i ω x_{j}} (e^{- i ω ξ)} - 2 + e^{i ω ξ}) = - \frac{4}{ξ^{2}} {sin}^{2} (\frac{ω ξ}{2}) v_{j} .

In the first equation, we have

	$\frac{1}{ξ^{2}} (- v_{1} + v_{2}) =$	$\frac{1}{ξ^{2}} (- cos \frac{ω ξ}{2} + cos (3 \frac{ω ξ}{2})) = \frac{4}{ξ^{2}} ({cos}^{3} (\frac{ω ξ}{2}) - cos \frac{ω ξ}{2})$
		$= - \frac{4}{ξ^{2}} {sin}^{2} (\frac{ω ξ}{2}) v_{1},$

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 = e_{m}^{T} (M v - λ v)$ , i.e.,

	$0! =$	$v_{m - 1} - (θ + λ ξ^{2}) v_{m} = cos (ω (x_{m} - ξ)) - (θ + λ ξ^{2}) cos (ω x_{m})$
	$=$	$(cos (ω ξ) + 4 {sin}^{2} (\frac{ω ξ}{2}) - θ) cos (ω x_{m}) + sin (ω ξ) sin (ω x_{m})$
	$=$	$(1 + 2 {sin}^{2} (\frac{ω ξ}{2}) - θ) cos (ω x_{m}) + sin (ω ξ) sin (ω x_{m}),$

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

	$0 =$	$(1 + 2 {sin}^{2} (\frac{ω ξ}{2}) - θ) (cos (ω) cos (\frac{ω ξ}{2}) + sin (ω) sin (\frac{ω ξ}{2}))$
		$+ 2 sin (\frac{ω ξ}{2}) cos (\frac{ω ξ}{2}) (sin (ω) cos (\frac{ω ξ}{2}) - cos (ω) sin (\frac{ω ξ}{2}))$
	$=$	$(1 - θ) cos (\frac{ω ξ}{2}) cos (ω) + (1 + 2 {sin}^{2} (\frac{ω ξ}{2}) - θ + 2 {cos}^{2} (\frac{ω ξ}{2})) sin (\frac{ω ξ}{2}) sin (ω)$
	$=$	$(1 - θ) cos (\frac{ω ξ}{2}) cos (ω) + (3 - θ) sin (\frac{ω ξ}{2}) sin (ω) .$

Hence, the different versions of $θ$ in (10) verify the condition (15) for $m = 1 / ξ$ . Rearranging (15) as $tan (ω) = β_{0} / (2 m β_{1}) cot (ω / (2 m))$ , for $β_{0} > 0$ it is seen that exactly $m$ solutions exist in $(0, m π)$ since the function $ω \mapsto cot (ω / (2 m))$ is monotonically decreasing and positive. Finally, the vector norms are computed for $ω \neq 0$ . Abbreviating $ω / m =: Ω$ and using ${cos}^{2} (x) = (1 + cos (2 x)) / 2$ , we get

	$m \sum j = 1 {cos}^{2} ((j - \frac{1}{2}) \frac{ω}{m}) =$	$\frac{m}{2} + \frac{1}{2} m \sum j = 1 cos ((2 j - 1) Ω)$
	$=$	$\frac{m}{2} + \frac{1}{2} R m \sum j = 1 e^{i (2 j - 1) Ω} = \frac{}{m} 2 + \frac{1}{2} R \frac{e^{i 2 Ω m} - 1}{e^{i Ω} - e^{- i Ω}}$
	$=$	$\frac{m}{2} + \frac{1}{4} I \frac{e^{i 2 Ω m} - 1}{sin (Ω)} = \frac{m}{2} + \frac{1}{4} \frac{sin (2 ω)}{sin (Ω)},$

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

For later use, we introduce the diagonal matrix $Λ = diag (λ_{k})$

and the unitary matrix

V = (v^{[1]}, \dots, v^{[m]})

satisfying

M = V Λ V^{T}

Remark 2.1

The well-known frequencies for Dirichlet boundary conditions are $ω_{k} = (k - \frac{1}{2}) π$ and $ω_{k} = (k - 1) π$ , $k = 1, \dots, m$ for Neumann conditions. For general values $β_{0}, β_{1} > 0$ the equation (15) may be rewritten in fixed point form

ω = f_{k} (ω) := (k - 1) π + arctan (\frac{β_{0}}{2 m β_{1}} cot (\frac{ω}{2 m})), k = 1, \dots, m .

(17)

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 $2 m > β_{0} / β_{1} \geq 1$ , since $1 / | f_{k}^{'} (ω) | = (4 m^{2} β_{1} / β_{0}) {sin}^{2} (ω / (2 m)) + (β_{0} / β_{1}) {cos}^{2} (ω / (2 m))$ .

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 $y (t) = \sum_{k = 1, \dots, m} η_{k} (t) v^{[k]}$ leads to

m \sum k = 1 η_{k}^{'} (t) v^{[k]} = m \sum k = 1 λ_{k} η_{k} (t) v^{[k]} + γ e_{m} u (t) .

(18)

Since the matrix $M$ is symmetric, the inner product with $v^{[j]}$ yields the decoupled equations $η_{j}^{'} (t) = λ_{j} η_{j} (t) + γ v_{m}^{[j]} 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

y (t) =

m \sum k = 1 (e^{λ_{k} t} {v^{[k]}}^{T} ψ + γ v_{m}^{[k]} \int_{0}^{t} e^{λ_{k} (t - τ)} u (τ) d τ) v^{[k]} .

(19)

Remark 2.2

The presence of the terms $v_{m}^{[k]}$ 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 $u (t)$ inherited from the boundary condition (2) may be considered as a control to approach a given target profile $^y \in R^{m}$ at some given time $T > 0$ . In an optimal control context, controls are searched for minimizing an objective function like

C = \frac{1}{2} ∥ y (T) -^y ∥_{2}^{2} + \frac{α}{2} \int_{0}^{T} u (t)^{2} d t,

(20)

with $α > 0$ . Optima may be computed by using some multiplier function $p (t)$ for the ODE restriction (11) and considering the Lagrangian

	$L :=$	$- C + \int_{0}^{T} p^{T} (y^{'} - M y - γ e_{m} u) d t + p^{T} (0) (y (0) - ψ)$
	$=$	$- C - \int_{0}^{T} ({(p^{'})}^{T} y + p^{T} (M y + γ e_{m} u)) d t + p^{T} (T) y (T) - p^{T} (0) ψ .$		(21)

The partial derivatives of the Lagrangian $L$ with respect to $p (t)$ and $p (T)$ recover (11), (13) and the other ones are

$\partial_{y (t)} L =$	$- p^{'} - M^{T} p = - p^{'} - M p,$	(22)
$\partial_{y (T)} L =$	$^y - y (T) + p (T)$	(23)
$\partial_{u (t)} L =$	$- α u - γ e_{m}^{T} p .$	(24)

Hence, the Karush-Kuhn-Tucker conditions, $\partial_{(\cdot)} L = 0$ in (22)–(24), show that the control $u (t)$ may be eliminated by

u (t) = - \frac{γ}{α} e_{m}^{T} p (t) = - \frac{γ}{α} p_{m} (t),

(25)

and a necessary condition for the optimal solution is that it solves the following boundary value problem:

	$y^{'} =$	$M y - \frac{γ}{α} e_{m} e_{m}^{T} p, y (0) = ψ,$		(26)
	$p^{'} =$	$- M p, p (T) = y (T) -^y .$		(27)

The homogeneous differential equation (27) for $p$ has the simple solution

p (t) = e^{(T - t) M} p (T) = m \sum ℓ = 1 e^{λ_{ℓ} (T - t)} v^{[ℓ]} {v^{[ℓ]}}^{T} (y (T) -^y)

(28)

with the matrix exponential

e^{s M} =

V diag (e^{λ_{k} s}) V^{T}, 0 \leq s \leq T .

(29)

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

	$η_{k} (t) =$	$e^{λ_{k} t} η_{k} (0) - \frac{γ^{2}}{α} v_{m}^{[k]} \int_{0}^{t} e^{λ_{k} (t - τ)} p_{m} (τ) d τ$
	$=$	$e^{λ_{k} t} η_{k} (0) - \frac{γ^{2}}{α} v_{m}^{[k]} m \sum ℓ = 1 v_{m}^{[ℓ]} \int_{0}^{t} e^{λ_{k} t + λ_{ℓ} T - (λ_{k} + λ_{ℓ}) τ} d τ \cdot {v^{[ℓ]}}^{T} (y (T) -^y) .$		(30)

Considering the function $φ_{1} (z) = \int_{0}^{1} e^{z t} d t$ satisfying $φ_{1} (z) = (e^{z} - 1) / z$ for $z \neq 0$ and $φ_{1} (0) = 1$ , the integral may be written as

\int_{0}^{t} e^{λ_{k} t + λ_{ℓ} T - (λ_{k} + λ_{ℓ}) τ} d τ = t e^{λ_{ℓ} (T - t)} φ_{1} ((λ_{k} + λ_{ℓ}) t) .

(31)

The result (30) may be used in different ways. The first application is computing the solution for a given target profile $^y$ .

Lemma 2.3

Let $^y \in 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

(I + Q) η (T) = e^{T Λ} η (0) + Q V^{T}^y,

(32)

with the positive semi-definite matrix $Q = (q_{k ℓ})_{k, ℓ = 1}^{m}$ having the elements

q_{k ℓ} = \frac{γ^{2} T}{α} v_{m}^{[k]} φ_{1} ((λ_{k} + λ_{ℓ}) T) v_{m}^{[ℓ]}, k, ℓ = 1, \dots, m .

(33)

Proof: At the end point $T$ , the formula (30) simplifies to

η_{k} (T) =

e^{λ_{k} T} η_{k} (0) - \frac{γ^{2} T}{α} v_{m}^{[k]} m \sum ℓ = 1 v_{m}^{[ℓ]} φ_{1} ((λ_{k} + λ_{ℓ}) T) \cdot (η_{ℓ} (T) - {v^{[ℓ]}}^{T}^y) .

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

	$w^{T} Q w =$	$γ2Tα\emphm∑k,ℓ=1v[k]mwk∫10e(λk+λℓ)Tτdτ⋅v[ℓ]mwℓ$
	$=$	$\frac{γ^{2} T}{α} \int_{0}^{1} m \sum k, ℓ = 1 (e^{λ_{k} T τ} v_{m}^{[k]} w_{k}) (e^{λ_{ℓ} T τ} v_{m}^{[ℓ]} w_{ℓ}) d τ$
	$=$	$\frac{γ^{2} T}{α} \int_{0}^{1} (m \sum k = 1 e^{λ_{k} T τ} v_{m}^{[k]} w_{k})^{2} d τ \geq 0.$

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 $^y$ , 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]},

(34)

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

^y = y (T) - δ_{1} v^{[1]} - δ_{2} v^{[2]},

(35)

where, by (30), the coefficients of $y (T)$ are given by

η_{k} (T) =

e^{λ_{k} T} η_{k} (0) - \frac{γ^{2} T}{α} v_{m}^{[k]} 2 \sum ℓ = 1 γ_{ℓ} v_{m}^{[ℓ]} φ_{1} ((λ_{k} + λ_{ℓ}) T) .

(36)

We will use this construction in our numerical example.

Figure 1: Dirichlet heat problem with $m = 250, 500$ spatial points and given exact control $u (t)$ : $∥ y (T)$ $-$ $y_{h} (T) ∥_{\infty}$ (left) and $∥ p (0)$ $-$ $p_{h} (0) ∥_{\infty}$ (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:

	$min (y, u) C$	$:= \frac{1}{2} ∥ y (T) -^y ∥_{2}^{2} + \frac{α}{2} \int_{0}^{T} u (t)^{2} d t$
	$subject to y^{'} (t)$	$= M y (t) + γ e_{m} u (t), t \in (0, T],$
	$y (0)$	$= 1l,$

with $T = 1$ , $γ = 2 / ξ^{2}$ , $α = 1$ , $1l = (1, \dots, 1)^{T} \in R^{m}$ , state vector $y (t) \in R^{m}$ , and $M$ as defined in (12) with $θ = 3$ . We set $δ_{1} = δ_{2} = 1 / 75$ in (34) and compute the target profile $^y \in 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].

Figure 2: Dirichlet heat problem with $m = 250, 500$ spatial points, solved by gradient descent for $(y, p, u)$ : ${max}_{n, i} | u (t_{n i})$ $-$ $u_{h} (t_{n i}) |$

Two 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. The initial value for the multiplier is set to $p (T) = y_{h} (T) -^y$ 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_{n i} = t_{n} + c_{i} h$ , $i = 1, \dots, s,$ chosen by the time integrators on a time grid ${t_{0}, \dots, 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, \dots, 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.

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.

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

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

Strong Rates of Convergence for Space-Time Discretization of the Backward Stochastic Heat Equation, and of a Linear-Quadratic Control Problem for the Stochastic Heat Equation

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

Error Inhibiting Schemes for Initial Boundary Value Heat Equation

Discrete Adjoint Implicit Peer Methods in Optimal Control

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

Adjoint-based optimization of two-dimensional Stefan problems

1 Introduction

2 A discrete heat equation with boundary control

2.1 Finite difference discretization of the 1D heat equation

Lemma 2.1

Remark 2.1

2.2 Exact solution of the discrete heat equation

Lemma 2.2

Remark 2.2

2.3 Exact solution of an optimal control problem

Lemma 2.3

3 Test case: Dirichlet boundary control problem

References

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

Related Research

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

Strong Rates of Convergence for Space-Time Discretization of the Backward Stochastic Heat Equation, and of a Linear-Quadratic Control Problem for the Stochastic Heat Equation

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

Error Inhibiting Schemes for Initial Boundary Value Heat Equation

Discrete Adjoint Implicit Peer Methods in Optimal Control

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

Adjoint-based optimization of two-dimensional Stefan problems

1 Introduction

2 A discrete heat equation with boundary control

2.1 Finite difference discretization of the 1D heat equation

Lemma 2.1

Remark 2.1

2.2 Exact solution of the discrete heat equation

Lemma 2.2

Remark 2.2

2.3 Exact solution of an optimal control problem

Lemma 2.3

3 Test case: Dirichlet boundary control problem

References