RIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

laboratoire de recherche
en informatique
et automatique

A global shock fitting method.

W.L. Miranker, O. Pironneau

Rapport de Recherche n°123
Mai 1975
# A GLOBAL SHOCK FITTING METHOD

W.L. Miranker, O. Pironneau

Présenté par R. Glowinski

## RESUME :

Dans ce rapport, on montre comment l'onde de choc d'un problème évolutif peut être caractérisée comme la solution d'un problème d'optimisation.

Ce problème est résolu par une méthode de gradient. On obtient ainsi une méthode globale non dissipative pour le calcul des solutions de certaines équations aux dérivées partielles avec chocs.

## ABSTRACT :

We show how a shock can be characterized as the solution of an optimization problem. The optimization problem is solved directly by a gradient calculation and the method of steepest descent. The result is a global treatment of shocks eliminating dispersion.
# A GLOBAL SHOCK FITTING METHOD

by

W.L. Miranker and O.Pironneau

## ABSTRACT

We show how a shock can be characterized as the solution of an optimization problem. The optimization problem is solved directly by a gradient calculation and the method of steepest descent. The result is a global treatment of shocks eliminating dispersion.

## INTRODUCTION

In this paper we devise and study a method for shock fitting. We show how a shock can be characterized as the solution of an extremal problem in the calculus of variations and our method consists of an exploitation of this characterization.

Shock fitting is one of the methods used to calculate the weak solution of the shock wave problem. It avoids the introduction of dissipation characteristic of the class of methods which employ devices such as artificial viscosity.

Shock fitting consists of a process of equation solving to produce the shock profile step by step. The equations are obtained by coupling the Rankine-Hugoniot conditions and the difference equations for the state variables at each time step.

Our procedure gives a global characterization of the shock as the solution of an optimization problem. It replaces the usual step by step solution process by the global process of direct minimization of an appropriate functional. The minimization in turn is accomplished by a gradient calculation and the method of steepest descent. The optimization characterization which we employ is novel for shocks. It is a development of a method first used by one of us (cf. Pironneau [2]).
to characterize certain free boundary problems directly as
problems of optimal control. Shock fitting in the traditional
form has been employed by others as Quarles [5] and Thomas [7].

We begin with the treatment of a single equation of
evolution of order one in one dimension where the shocks springs
from a discontinuity in the initial conditions.

In §1 we give a statement of the shock wave problem and
the derivation of the optimization formulation. Then the gradient
calculation and a description of the minimization technique along
with a closed form example follow. In §2 we give a precise state-
ment of the descent algorithm, a theorem and several remarks
concerning the convergence of this algorithm. In §3 we deal with
computational aspects. We discretize the continuous characteri-
zation and derive from it a discrete form of the descent algo-
rithm. This is followed by an error analysis and then a compu-
tational example.

In §4 we discuss the extension of our method to systems
of equations in several space variables. We limit the discussion
to the case of a model problem.
# 1. STATEMENT OF THE METHOD.

Consider the initial value problem for Burger's equation

$$ (1.1) \qquad \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = g(x,t) $$
$$ u(x,0) = u_0(x) \quad x \in \mathbb{R}, 0 \le t \le T. $$

As is well known, the solution of this problem may develop shocks across which it is discontinuous. At such a shock, the Rankine-Hugoniot condition

$$ (1.2) \qquad \frac{[f(u)]}{[u]} = \dot{\xi}(t) \quad 0 \le t \le T $$

is satisfied. Here the shock is given by $x = \xi(t)$.

In addition the entropy condition

$$ (1.3) \qquad f'(u(\xi^+(t),t)) \le \dot{\xi}(t) \le f'(u(\xi^-(t),t)) $$

should be satisfied.

To study the problem of shock fitting, we assume that there is a single shock which results from a discontinuity in $u_0(x)$ at $t=0$.

We now introduce two auxiliary initial value problems $P^+$ and $P^-$ for Burger's equation with solution $u^+$ and $u^-$ respectively. For initial data, we take $u_0^+$ and $u_0^-$ in $C^1(-\infty, +\infty)$ with :

$$ u_0^- = u_0 \quad \text{for } x < 0 $$

$$ u_0^+ = u_0 \quad \text{for } x > 0. $$

We require that $u_0^+$ and $u_0^-$ be chosen so that the solutions $u^+$ and $u^-$, respectively, are free of shocks in $[0,T]$. We call $C^+(C^-)$ resp. the characteristic of $P^+$ ($P^-$ resp.) which emanates from the origin.
Let these characteristics be given by $x(t) = C^{\pm}(t)$ respectively.
Let :

$$
(1.4) \qquad \Omega = \{(x,t) | \theta^{+}(t) \le x \le \theta^{-}(t), 0 \le t \le T\}
$$

We note that $C^+$, $C^-$ and $\Omega$ are independent of the choice of
$u_o^+$ and $u_o^-$. Furthermore the values of $u^+$ and $u^-$ in $\Omega$ are inde-
pendent of this choice as well. We also note that the entropy condi-
tions implies that $\Omega$ is not empty. Let

$$
(1.5) \qquad J(\zeta(t),t) = \frac{\{f(u^+(\zeta(t),t)) - f(u^-(\zeta(t),t))\}}{\{u^+(\zeta(t),t) - u^-(\zeta(t),t)\}}.
$$

Now we consider the functional

$$
F(\zeta) = \int_{0}^{T} (\dot{\zeta}(t) - J(\zeta(t),t))^{2} dt .
$$

Note that $F \ge 0$ and that the Rankine-Hugoniot condition implies
that $F(\xi) = 0$. Therefore the shock $\xi$ is the solution of the follo-
wing extremal problem

$$
\min\{F(\zeta) | \zeta(0) = 0, \quad \zeta \in \Omega \cap C^1[0,T]\}
$$

We remark that the value of $F(\zeta)$ for any admissible $\zeta$ is
independent of the choice of $u_o^+$ and $u_o^-$ . Furthermore, if $u^+$, $u^-$, $f$
are continuously differentiable, (or $f \in C^2$ if the shock does not
persist for all $t \in [0,T]$) the differential equation :

$$
\dot{\zeta} = J(\zeta(t),t), \quad \zeta(0) = 0
$$

has a unique solution. Therefore the constraint $\zeta \in \Omega$ may be dropped
and $\zeta$ is also the unique solution of

$$
(1.6) \qquad \min_{\zeta \in C^1[0,T]} \{F(\zeta) \mid \zeta(0)=0\}.
$$
We proceed to solve the extremal problem by using techniques
of mathematical programming.

Gradient calculation.

We will now compute the Frechet derivative of F. Let ζ and
ζ + δζ be in C¹[0,T] with ζ(o) = δζ(o) = 0. Then

$$
\begin{align*}
F(\zeta + \delta\zeta) - F(\zeta) &= \int_{0}^{T} \left\{ [\dot{\zeta} + \delta\dot{\zeta} - J(\zeta + \delta\zeta, t)]^2 - [\dot{\zeta} - J(\zeta, t)]^2 \right\} dt \\
&= \int_{0}^{T} \left\{ [\dot{\zeta} + \delta\dot{\zeta} - J(\zeta, t) - \frac{\partial J(\zeta, t)}{\partial x}\delta\zeta]^2 + [\dot{\zeta} - J(\zeta, t)]^2 \right\} dt \\
&= 2 \int_{0}^{T} (\dot{\zeta} - J(\zeta, t))(\delta\dot{\zeta} - \frac{\partial J(\zeta, t)}{\partial x}\delta\zeta) dt + o(\delta\dot{\zeta}) .
\end{align*}
$$

Here we have used the fact that o(δζ) = o(δζ̇) and o(δζ̇δζ) = o(δζ̇),
since ζ(o) = 0. Now

$$
\int_{0}^{T} (\dot{\zeta} - J(\zeta, t)) \frac{\partial J}{\partial x}(\zeta, t) \delta\zeta dt = \int_{0}^{T} \delta\dot{\zeta} \int_{t}^{T} (\dot{\zeta} - J(\zeta, \tau)) \frac{\partial J}{\partial x}(\zeta, \tau) d\tau dt \text{ as an}
$$

integration by part and the fact that $\delta\zeta(0) = 0$ show. Thus :

$$
(1.7) \quad F(\zeta + \delta\zeta) - F(\zeta) = 2 \int_0^T \delta\dot{\zeta}(t) [\dot{\zeta}(t) - J(\zeta, t)] \left( \int_0^T (\zeta(\tau) - J(\zeta, \tau)) \frac{\partial J}{\partial x}(\zeta, \tau) d\tau \right) dt
$$

## Outline of the method.

To find $\dot{\xi}(t)$ we use a method of descent based on (1.6) and (1.7).
That is given a curve $\zeta(t)$ with $\zeta(o) = 0$ and an associated value for
$F(\zeta)$, we choose a successor $\zeta(t)+\delta\zeta(t)$ so that $F(\zeta+\delta\zeta)$ is diminished.
This can be done by choosing $\delta\zeta(o) = 0$ and

$$
(1.8) \qquad \delta\dot{\zeta}(t) = - \lambda \nabla F(\zeta)(t)
$$

where

$$
(1.9) \quad \nabla F(\zeta)(t) = 2[\dot{\zeta}(t) - J(\zeta(t), t)] - \int_{t}^{T} (\dot{\zeta}(\tau) - J(\zeta(\tau), \tau)) \frac{\partial J}{\partial x}(\zeta(\tau), \tau) d\tau
$$

and where $\lambda$ is a small positive constant. For this choice of $\delta\zeta$
we have :

$$ (1.10) \quad F(\zeta + \delta\zeta) - F(\zeta) = \lambda \int_0^T \nabla F(\zeta)(t)^2 dt + o(\lambda \nabla F). $$

For $\lambda \nabla F$ sufficiently small, the right member is $< 0$.

**Example.**

Let $f(u) = \frac{1}{2}u^2$, $g=0$ so that
$$ J(\zeta, t) = \frac{1}{2}[u^+(\zeta, t) + u^-(\zeta, t)]. \quad \text{Let} $$

$$ u_o(x) = \begin{cases} 1 & x \le 0 \\ 0 & x > 0 \end{cases} . \quad \text{Then} $$

$$ \Omega = \{(x,t) | o \le x \le t, o \le t \le T\}. $$

Furthermore

$$ \begin{gather*} u^{-}(x,t) \equiv 1, \quad \text{and} \quad u^{+}(x,t) \equiv 0 \quad \text{in } \Omega. \quad \text{Therefore} \quad J(\zeta,t) \equiv \frac{1}{2} \\ F(\zeta) = \int_{0}^{T} \left(\dot{\zeta}(t) - \frac{1}{2}\right)^{2} dt \quad \text{and} \quad \nabla F(\zeta)(t) = 2\dot{\zeta}(t). \end{gather*} $$

Let the initial choice for $\zeta(t)$ be $\zeta_o(t) = \beta_o t - \frac{\alpha_o}{2}t^2$; then the method generates the sequence $\{\zeta_i(t)\}_{i \ge 1}$ such that

$$ \zeta_{i+1}(t) - \zeta_i(t) = -2\lambda(\dot{\zeta}_i(t) - \frac{1}{2}) $$

according to (1.8). Therefore

$$ \zeta_{i+1}(t) = \beta_{i+1} t - \frac{\alpha_{i+1}}{2} t^2 $$

where

$$ \beta_{i+1} = \beta_i(1-2\lambda) + \lambda $$

$$ \alpha_{i+1} = \alpha_i(1-2\lambda) . $$
Thus

$$
\beta_i = \beta_o (1-2\lambda)^{i+1} + \sum_{j=0}^{i} \lambda (1-2\lambda)^j = \frac{1-(1-2\lambda)^{i+1}}{2} + \beta_o (1-2\lambda)^{i+1} \\
\alpha_i = \alpha_o (1-2\lambda)^i
$$

Thus if $0 < \lambda < \frac{1}{2}$ $$ \lim_{i \to \infty} \zeta_i(t) = \frac{t}{2} $$ .

Note that $\zeta_i$ converges geometrically towards its limit.

2. CONVERGENCE.

In this section, we give a precise statement of the shock
fitting algorithm and a theorem asserting the convergence of the
algorithm to the actual shock profile.

The algorithm is described as follows :

Algorithm 2.1.

1. Set i = 0 and choose θ ∈ (0,1)

2. Choose $\zeta_o \in C^1[0,T]$ with $\zeta_o(o) = 0$

3. Compute $u^+$ and $u^-$.

4. Compute $\nabla F$ according to (1.9) and (1.5).

5. Choose $\lambda = \theta^j$ where $j$ is the smallest natural number such that

$$
F(\zeta_i - \lambda \nabla F) - F(\zeta_i) \le - \frac{\lambda}{2} \|\nabla F\|_{L^2[0,T]}^2 \\
\text{is true.}
$$

6. Let

$$
\zeta_{i+1}(t) = \zeta_i(t) - \lambda \int_0^t \nabla F(\zeta)(\tau) d\tau
$$

7. Set $i = i+1$ and go to step 4.
Comments

1. $u^+$ and $u^-$ are computed only once.

2. choosing $\lambda$ according to step 5 is a procedure which is well known in descent methods. Such a choice assures the descent and provides a geometric rate of convergence of the iterates (cf. Polak [4]).

We now state and prove the following convergence theorem.

**THEOREM** Assume that $u^+, u^-$, $f$ are continuously differentiable (or $f \in C^2$ if the shock does not persist for all $t \in [0,T]$). Then Algorithm 2.1 generates a sequence all accumulation points ($L^\infty[0,T]$) of which coincides with $\xi(\cdot)$, the shock profile.

**Proof.** Algorithm 2.1 is a method of steepest descent applied to the problem :

$$ \min_{\zeta \in C[0,T]} F \left[ \int_{0}^{t} \dot{\zeta}(\tau) d\tau \right] $$

as is seen from (1.7). Equations (1.9) and (1.7) show that $F$ is continuously Frechet differentiable (with the $L^2[0,T]$ norm) with respect to $\dot{\zeta}$. Thus every accumulation point $\bar{\zeta}$ of the sequence $\{\zeta_i\}$ must satisfy $\nabla F(\bar{\zeta})(.) = 0$, i.e.

$$ (2.1) \qquad \dot{\bar{\zeta}} - J(\bar{\zeta}, t) = \int_{t}^{T} (\dot{\bar{\zeta}} - J(\zeta, \tau)) \frac{\partial J}{\partial x}(\bar{\zeta}, \tau) d\tau . $$

By hypothesis $f, u^+, u^-$ are continuously differentiable so that $\frac{\partial J}{\partial x}(\bar{\zeta}, t)$ is continuous. Thus (2.1) and the Bellman-Gronwall lemma imply that

$$ \dot{\bar{\zeta}} = J(\bar{\zeta}, t) . $$

This equation has a unique solution for $\bar{\zeta}(0) = 0$. This demonstrates the uniqueness of the accumulation points and the theorem.
We make the following remarks :

### Remark 2.1

The sequence produced by the algorithm could fail to have an accumulation point if its elements $\zeta_i$ tend to infinity. The form of F shows that the sequence $\{\zeta_i\}$ as well must tend to infinity in this case.

### Remark 2.2

We may choose the extremal problem by restricting $\zeta$ to be a sufficiently large compact set. The algorithm would then be altered to include this constraint, for example by means of a gradient projection or penalty technique. In this case the convergence to the solution is assured.

### Remark 2.3

Our numerical experience shows that $\zeta_i$ remains bounded and that the algorithm converges usually without activating the constraints described in the above remarks.

## 3. COMPUTATIONAL ASPECTS.

In this section we present a discretization of the algorithm introduced in §2 and a description of the properties of this discretization.

We begin by introducing the mesh $M = \{(x_j, t_n) | x_j = j\Delta x, j=0, \pm 1, \dots, t_n = (n-1)\Delta t, n=1, 2, \dots, N\}$.

We suppose that $u^+$ and $u^-$ are obtained by some numerical procedure. We remark that $u^+$ and $u^-$ are smooth by virtue of the way in which we composed the initial data. Thus these two functions are obtained easily on the mesh by any standard numerical scheme (such as the Lax-Wendroff scheme, cf. Richmeyer and Morton [6]). To implement the algorithm we first introduce a discrete functional $F_d$ approximating $F$. Then we formulate a descent algorithm for $F_d$.
Let $\zeta^n = \zeta(t_n)$ and let $j_n = \left[ \frac{\zeta^n}{\Delta x} + 0.5 \right]$
($j_n$, $n$ are the indices of the point in $M \cap \{t=t_n\}$ nearest to $(\zeta^n, t^n)$.

Let

$$
(3.1) \quad J_j^n \equiv J(x_j, t_n) = \frac{[f(u^+(x_j, t_n)) - f(u^-(x_j, t_n))]}{[u^+(x_j, t_n) - u^-(x_j, t_n)]} / \Delta t
$$

Then for $F_d$ we take

$$
(3.2) \quad F_d = \sum_{n=1}^{N} \left( \frac{\zeta^{n+1} - \zeta^n}{\Delta t} - J_{j_{n+1}}^{n+1} \right) \Delta t .
$$

Then taking variations we take

$$
\begin{align*}
(3.3) \quad \frac{1}{2} \delta F_d ={}& \sum_{n=1}^{N} \left( \frac{\zeta^{n+1} - \zeta^n}{\Delta t} - J_{j_{n+1}}^{n+1} \right) \left( \frac{\delta \zeta^{n+1} - \delta \zeta^n}{\Delta t} - \frac{n+1}{\Delta x} \right) \delta \zeta^{n+1} \Big|_{\Delta t = 0} \\
& - \sum_{n=2}^{N} \sum_{k=n}^{n} a_n \delta \zeta^{k+1} - \delta \zeta^k
\end{align*}
$$

where

$$
a_n = \left( \frac{\zeta^{n+1} - \zeta^n}{\Delta t} - J_{j_{n+1}}^{n+1} \right) \left( \frac{J_{j+1}^{n+1} - J_{j_{n+1}}^{n+1}}{\Delta x} \right) \Delta t
$$

and where we have used $\delta\zeta^1 = 0$.

Now

$$
\sum_{n=2}^{N} a_n \sum_{n=1}^{N} (\delta\zeta^{k+1} - \delta\zeta^k) = \sum_{n=1}^{N} \left( \sum_{k=n}^{N} a_k \right) (\delta\zeta^{n+1} - \delta\zeta^n)
$$

Then

$$
\begin{align*}
(3.4) \quad \frac{1}{2} \delta F_d &= \sum_{n=1}^{N} \left[ \frac{\zeta^{n+1} - \zeta^n}{\Delta t} J_{j_{n+1}}^{n+1} - \sum_{k=n}^{N} \left( \frac{\zeta^{n+1} - \zeta^n}{\Delta t} - J_{j_{n+1}}^{n+1} \right) \left( \frac{n+1}{\Delta x} \right) \right] \Delta t (\delta \zeta^{n+1} - \delta \zeta^n) \\
&= \sum_{n=1}^{N} \nabla F_d^n (\delta \zeta^{n+1} - \delta \zeta^n).
\end{align*}
$$
The discrete analogue of Algorithm 2.1 consists of making
the following replacements in that algorithm :

F by $F_d$, $\zeta(.)$ by $\{\zeta^n\}_{n=1}^N$, $u^+$ and $u^-$ by $u_j^{+n}$, $u_j^{-n}$ and
the integral by

$$
\sum_{k=1}^{n} \nabla_{F_d}^{k} \Delta t .
$$

We now turn to the error analysis.

## Error Analysis

Let us suppose that $u^+$ and $u^-$ are computed with an error $\epsilon$
and that the application of the discrete algorithm results in a
sequence converging to $\{\xi_d^n; n=1, \dots, N\}$. As in the continuous case
we have

$$
\frac{\xi_d^{n+1} - \xi_d^n}{\Delta t} = J_{n+1}^{n+1} .
$$

If $e^n = \xi_d^n - \xi(t_n)$ where $\xi(t)$ is the shock, we have

$$
\begin{align*}
\frac{e^{n+1} - e^n}{\Delta t} &= J_{n+1}^{n+1} - \frac{\xi^{n+1} - \xi^n}{\Delta t} = J_{n+1}^{n+1} - J(\xi^{n+1}, t_{n+1}) + \xi^{n+1} - \frac{\xi^{n+1} - \xi^n}{\Delta t} \\
&= \nabla_x J(\xi_d^{n+1} - \xi^{n+1} + \theta_{n+1} \Delta x) + O(\Delta t) + O(\epsilon)
\end{align*}
$$

where $\theta_{n+1} \in (-.5, +.5)$. Thus since $\nabla_x J$ is bounded and
$\Delta x = O(\Delta t)$, we have :

$$
e^{n+1} = (1 + \Delta t \nabla_x J) e^n + O(\Delta t^2) + O(\epsilon \Delta t) .
$$

Since $e^1 = 0$ a standard argument yields

$$
\|e^n\| \le O(\Delta t) + O(\epsilon) .
$$
A numerical example

As an illustration we treat the problem :

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0 \quad t > 0 \quad -\infty < x < +\infty $$

$$ u(x,0) = \begin{cases} .5 & x \le 0 \\ 1.5 & x > 0 \end{cases} $$

$$ \text{with } \xi = \frac{1}{2}(u(\xi^{-}(t),t) + u(\xi^{+}(t),t)) $$

as the Rankine-Hugoniot condition.

The functions $u^+(x,t)$ and $u^-(x,t)$ are calculated in the mesh by mean of the Lax-Wendroff scheme. Calling

$$ u_j^n = u(x_j, t_n) ; \quad F_j^n = \frac{1}{2}(u_j^n)^2 ; \quad A_{j+\frac{1}{2}}^n = \frac{1}{2}(u_{j+1}^n + u_j^n) , $$

this scheme is

$$ u_{j}^{n+1} = u_{j}^{n} - \frac{1}{2} \frac{\Delta t}{\Delta x} (F_{j+1}^{n} - F_{j-1}^{n}) + \frac{1}{2} \left(\frac{\Delta t}{\Delta x}\right)^2 \left[ A_{j+\frac{1}{2}}^{n} (F_{j+1}^{n} - F_{j}^{n}) - A_{j-\frac{1}{2}}^{n} (F_{j}^{n} - F_{j-1}^{n}) \right] $$

using the value $\Delta t=.05$, $\Delta x=.02$, $N = 20$ ($T=1$) and fixed descent step size $\lambda=.5$, we obtained the results summarized in the following table :

<table><thead><tr><th></th><th>ζ<sup>1</sup></th><th>ζ<sup>5</sup></th><th>ζ<sup>10</sup></th><th>ξ exact</th><th>ζ<sub>exact</sub>-ζ<sup>10</sup></th></tr></thead><tbody><tr><th>E<br>t</th><td>0.275</td><td>0.0045</td><td>6 x 10<sup>-6</sup></td><td></td><td></td></tr><tr><th>.1</th><td>0.0192</td><td>0.0271</td><td>0.0272</td><td>0.0259</td><td>.001</td></tr><tr><th>.2</th><td>0.0383</td><td>0.0559</td><td>0.0566</td><td>0.0535</td><td>.003</td></tr><tr><th>.3</th><td>0.058</td><td>0.0867</td><td>0.0879</td><td>0.0825</td><td>.005</td></tr><tr><th>.4</th><td>0.0776</td><td>0.119</td><td>0.122</td><td>0.113</td><td>.009</td></tr><tr><th>.5</th><td>0.0977</td><td>0.153</td><td>0.157</td><td>0.145</td><td>.01</td></tr><tr><th>.6</th><td>0.118</td><td>0.188</td><td>0.194</td><td>0.177</td><td>.02</td></tr><tr><th>.7</th><td>0.138</td><td>0.225</td><td>0.234</td><td>0.211</td><td>.02</td></tr><tr><th>.8</th><td>0.158</td><td>0.264</td><td>0.276</td><td>0.246</td><td>.03</td></tr><tr><th>.9</th><td>0.179</td><td>0.305</td><td>0.32</td><td>0.282</td><td>.04</td></tr><tr><th>1.0</th><td>0.2</td><td>0.347</td><td>0.367</td><td>0.318</td><td>.05</td></tr></tbody></table>

$$ \text{In this table } E = \|\nabla F_d\|_{L^2(N)}^2. $$

Thus we see that the minimum is rapidly approached by the numerical method. Note that the error is of the order of $\Delta t$ and that it increases with time only because the Lax-Wendroff scheme is dissipative.
The programming was done on a time-sharing system in APL :
the computing time taken by N iterations of the minimization
part is of the same order as the time taken to compute u+ and
u-.

4. EXTENSION TO MANY VARIABLES

The method may be extended to include systems of equations
in several independent variables. In the case of a single equation
with several independent variables, the technique is essentially
the same ; the shock being replaced by a surface discontinuity.
It is still possible to extend the discontinuous initial data
into C¹ data u₀⁺ and u₀⁻ on either side of the discontinuity.
Moreover the corresponding solutions are computable in the region
Ω once and for all, as above.

This desirable feature disappears in the case of systems
because in the systems case there are at least two families of
characteristics. Then through each point in the vicinity of the
shock, at least one characteristic will encounter the shock
before it reaches the initial manifold. To illustrate this diffi-
culty and possible ways to deal with it, we consider a model pro-
blem for a system of two equations in one space variable.

Model problem.

Consider the pair of differential equations :

$$
\begin{equation}
\begin{aligned}
(4.1) \quad & \frac{\partial r}{\partial t} + \frac{3r+s}{4} \frac{\partial r}{\partial x} = 0 \\
(4.2) \quad & \frac{\partial s}{\partial t} + \frac{3s+r}{4} \frac{\partial s}{\partial x} = 0
\end{aligned}
\quad x \in [a,b], t \in [0,T]
\end{equation}
$$

with the Rankine-Hugoniot condition

$$
(4.3) \quad [r] = f(\xi, s^+, s^-, r^+, r^-)
$$

$$
(4.4) \quad [s] = g(\xi, s^{+}, s^{-}, r^{+}, r^{-})
$$
and for the data, we take :

$$
\begin{align*}
(4.5) \quad & s(a,t) = s_a(t) && t \in [0,T] \\
(4.6) \quad & r(b,t) = r_b(t) && \\
(4.7) \quad & s(x,s) = s_o(x) \\
(4.8) \quad & r(x,o) = r_o(x) &&
\end{align*}
$$

Here we suppose that one of both of the initial data have
a discontinuity at a point in (a,b) which we fix at x=0.

There are many optimization problems whose solution
coincides with the solution of (4.1)-(4.8). Consider the
following one which employs three control functions ξ(t), u(t),
v(t).

$$
(4.9) \quad \text{Min}_{\xi, u, v} \int_0^T \left\{ \left[ [s] - g(\dot{s}, s^+, s^-, r^+, r^-) \right]^2 + \left[ [r] - f(\dot{r}, s^+, s^-, r^+, r^-) \right]^2 \right\} dt
$$

where $(s^-, r^-), (s^+, r^+)$ are defined as follows :

$$
- (s^{-}, r^{-}) \text{ is defined in the domain } a < x < \xi(t) \quad 0 < t < T
$$

as the solution of (4.1)(4.2) subject to the conditions (4.5),
(4.7), (4.8) and

$$
r^{-}(\xi(t), t) = u(t):
$$

$$
- (s^+, r^+) \text{ is defined in the domain } \xi(t) < x < b, \ 0 < t < T
$$

as the solution of (4.1), (4.2) subject to the conditions (4.6)-(4.8)
and

$$
s^{+}(\xi(t), t) = v(t).
$$

Clearly a solution ξ of (4.7) is the shock for (4.1)-(4.8).
Since (4.9) is a distributed parameter optimal control problem, its
solution may be found by techniques analogous to the ones used for
the single equations case treated above. The gradient of the cost
function may be computed by standard means (cf. Lions [1],
Pironneau [3]). However unlike the case of problems of sections
1-3, an iterative scheme which employs this gradient requires
the redetermination of the four functions $s^+, s^-, r^+$ and $r^-$ as
well as four associated adjoint functions at each step. We believe
that much of this work may be eliminated by a more judicious choice
of the optimization problem.

It is this wide range of variability in the choice of
optimisation problem to characterize the shock which strikes us
as a principle feature of our method. This wide choice suggests
that it will be possible to develope shock fitting as an important
method in shock calculation.

---
REFERENCES

[1] Lions J.L. Optimal control of systems governed by partial differential equations. Springer (1968).

[2] Pironneau O. Optimisation de structures : application à la mécanique des fluides. Proc. IFIP Conference IRIA, June 1974.

[3] Pironneau O. Optimum profiles in Stokes flow. J. Fluid Mech. 59, pp. 117-128 (1973).

[4] Polak E. Computational method in Optimization. Academic Press, (1971).

[5] Quarles D., Ph. D. Thesis, New-York University (1954).

[6] Richmeyer D and Morton K. Difference methods for initial value problems. Interscience (1957).

[7] Thomas Computation of one-dimensional compressible flow including shocks. Comm. Pure and Applied Math. Vol. 7, p. 195, (2954).

---