RIA
laboria

laboratoire de recherche
en informatique
et automatique

ARCHIVES

identification
de la non-linéarité
d'une équation parabolique
quasilinéaire

G.Chavent,P.Lemonnier

66P

Rapport de Recherche n° 45
décembre 1973

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20
IDENTIFICATION DE LA NON-LINEARITÉ
D'UNE ÉQUATION PARABOLIQUE QUASILINEAIRE

G. Chavent - P. Lemonnier

- RESUME -

On étudie dans cet article le problème de la détermination de la non-linéarité dans une équation parabolique à partir de mesures ponctuelles sur la solution. Ceci correspond à la situation courante en physique, dans laquelle un paramètre dépend de l'état du système (par exemple un coefficient de conductibilité thermique dépend de la température).

Après avoir posé le problème comme un problème de contrôle dont on montre l'existence d'une solution, on calcule la dérivée du critère par rapport à la non-linéarité sans hypothèse sur la forme analytique de cette dernière.

Enfin on donne une application numérique de l'algorithme obtenu, qui montre que la situation considérée ici (estimation d'un paramètre dépendant de l'état) est plus favorable, du point de vue identification, que la situation où l'on se place habituellement (estimation d'un paramètre dépendant des variables d'espace et/ou de temps).

- ABSTRACT -

In this work, we study the problem of the determination of the non-linearity in a parabolic equation from measurements over its solution. This corresponds to an usual physical situation, in which the parameter depends on the state of the system (for instance, the heat conduction coefficient depends on the temperature).

We first state the problem as a control problem, for which we show the existence of a solution, and we calculate the derivative of the criterion with respect to the non-linearity without any assumption on the algebraic form of this latter.

We finally give a numerical application of the algorithm, which shows that the situation which is here taken in account (estimation of a parameter function depending on the state) is much better, from the identification point of view, than the usually studied situation for distributed system (estimation of a parameter function depending on space and/or time variables).
# PLAN

O - *Introduction*.

I - *Etude de l'équation parabolique quasilinéaire.*

II - *Dérivabilité de la solution en fonction du paramètre.*

1. Continuité de l'application $a \to y$.

2. Etude directe du système obtenu par différentiation formelle.

3. Dérivabilité de l'application $a \to y$.

III - *Etude du problème inverse.*

1. Position du problème, existence d'une solution.

2. Calcul du gradient de $J$ par rapport à $a$.

IV - *Discrétisation et Résolution numérique.*

1. Discrétisation de l'équation directe.

2. Discrétisation de l'équation adjointe.

3. Détermination du gradient continu.

4. Recherche de $a(y)$ sous forme libre.

5. Recherche de $a(y)$ sous forme paramétrique.

6. Introduction de contraintes sur $a(y)$.

7. Détermination du gradient discret.

V - *Résultats numériques.*

VI - *Conclusion.*
# INTRODUCTION.

Le progrès de l'analyse numérique et des mathématiques appliquées permet-
tent actuellement l'utilisation de modèles non linéaires dans de nombreux domaines
faisant intervenir des équations aux dérivées partielles (cf. par exemple [12] [16]
pour la chaleur, [1] [11] pour l'hydrologie, [10] [9] pour la physique).

Des résultats d'existence et d'unicité ont été donnés dans [14] pour de
nombreuses équations aux dérivées partielles comportant des non-linéarités de forme
analytique donnée. Le cas de non-linéarités de forme analytique non spécifiée a
été étudié dans des articles récents [15] [17] qui, admettant l'existence de la solu-
tion du système continu, donnent des résultats d'approximation de cette solution par
la méthode de Galerkin.

En ce qui concerne l'identification de tels systèmes, les études ont porté
jusqu'à présent sur l'estimation de paramètres fonctionnels dépendant d'une varia-
ble *indépendante* (tels l'espace ou le temps) (cf. [5] [6] et [8]), mais jamais, à notre
connaissance, sur la *non-linéarité* elle-même, à l'exception d'un article récent [2]
qui traite le cas très particulier d'une équation non linéaire se ramenant à une
équation linéaire par un changement de fonction.

Nous nous interessons ici au cas d'un système gouverné par une équation
simple (diffusion non linéaire), dans laquelle c'est la *non-linéarité* elle-même que
nous essayons de *déterminer* à partir de *mesures* sur la *solution* de l'équation.

En utilisant les techniques de [14], nous donnons des résultats d'existence
et d'unicité de la solution de l'équation parabolique non linéaire. Puis nous mon-
trons l'existence d'une solution au pro'lème inverse, et nous indiquons une méthode,
(généralisant [4] [5] [6] [7]) permettant une détermination effective de la non-linéarité
en tenant compte des impératifs physiques. Cette méthode utilise les concepts de la
théorie du contrôle dans les équations aux dérivées partielles (cf. [13] [5]) ainsi
qu'une technique originale de calcul du gradient par rapport à la fonction repré-
tant la non-linéarité. On peut ainsi s'affranchir de la nécessité de représenter
la non-linéarité par une formule algébrique, et la déterminer point par point.
Nous donnons enfin des résultats numériques illustrant les performances de l'algo-
rithme.
La méthode de calcul de la non-linéarité est très générale et peut s'appliquer à de nombreux autres exemples (nous l'appliquons actuellement à la détermination de perméabilités relatives dans des écoulements diphasiques en milieux poreux).

Plus précisement, nous nous sommes intéressés au problème suivant :
étant donné un domaine $\Omega$ de $\mathbb{R}^n$, de frontière $\Gamma$, on considère dans cet ouvert l'équation aux dérivées partielles :

$$ (1) \quad \left\{ \begin{array}{ll} \displaystyle \frac{\partial y}{\partial t} - \sum_{i=1}^{n} \frac{\partial}{\partial x_i} (a(y) \frac{\partial y}{\partial x_i}) = f(x,t) & \text{dans } \Omega \times \Gamma \cup \Gamma \times \partial \Omega \\ \\ \displaystyle \frac{\partial y}{\partial n} = 0 & \text{sur } \Gamma \times \partial \Omega \\ \\ y(x,0) = y_0(x) & \text{dans } \Omega \end{array} \right. $$

où $f, y_0$ sont des fonctions connues et où la fonction réelle d'une variable réelle $y \to a(y)$ est à identifier.

Pour cela, on mesure, par exemple à chaque instant, la valeur moyenne de la solution $y(x,t)$ dans $K$ sous-domaines $\Omega_1, \dots, \Omega_K$ de $\Omega$ ; ce qui donne $K$ fonctions $z_1(t), \dots, z_K(t)$, $t \in \Gamma \cup \partial \Omega$ connues.

On pose alors le problème de l'identification de la fonction $y \to a(y)$ sous la forme suivante :

Etant donnés deux nombres $\bar{\alpha}$ et $\bar{M}$ ($0 < \bar{\alpha} < \bar{M}$), trouver une fonction $a$ vérifiant

$$ (2) \quad 0 < \bar{\alpha} \le a(\bar{\zeta}) \le \bar{M} \quad \forall \bar{\zeta} \in \mathbb{R} $$

qui minimise

$$ (3) \quad J(a) = \sum_{k=1}^{K} \int_{0}^{T} \left( \frac{1}{|\Omega_k|} \int_{\Omega_k} (y(x,t)dx - z_k(t))^2 dt \right) $$

où $|\Omega_k|$ représente la surface du domaine $\Omega_k$.
L'équation (1) peut représenter une équation de la chaleur, les mesures $z_i(t)$ sont alors les températures mesurées en K points (par des thermocouples par exemple) et l'inconnue est alors la conductivité thermique $a$ en fonction de la température $y$, la capacité calorifique étant supposée constante et égale à 1 (cette restriction n'ayant rien de fondamental).

Pour déterminer numériquement $a$ en fonction de $y$, on minimisera $J(a)$ par une méthode de gradient, éventuellement avec contrainte. Un des buts de cet article est de montrer que ce gradient existe et comment le calculer.

Le plan de l'article est le suivant :

- Au paragraphe I, on étudiera les solutions de l'équation parabolique (1) pour $a$ fixé.

- Au paragraphe II, on montrera que, sous une hypothèse de régularité sur $y$, l'application $a \to y$ est dérivable.

- Au paragraphe III, on montrera que le problème d'identification de $a$, posé comme ci-dessus, admet au moins une solution, et on calculera $J'(a)$.

- Au paragraphe IV, on donnera la discrétisation utilisée pour les applications numériques.

- Enfin, le paragraphe V, présentera les résultats numériques et le paragraphe VI, les conclusions.
# I - ETUDE DE L'EQUATION PARABOLIQUE QUASILINEAIRE.

Nous utiliserons les méthodes développées dans [14].

Nous allons traiter un cas plus général que (1), englobant (1) et les équations semblables, mais présentant des conditions aux limites différentes.

Posons donc :

$$ (4) \qquad \left\{ \begin{array}{l} \Omega \subset \mathbb{R}^n \text{ ouvert, borné, de frontière } \Gamma \text{ régulière} \\ Q = \Omega x]_{\text{OT}[, \sum = \Gamma x]_{\text{OT}[.} \end{array} \right. $$

Soient V et H deux espaces de Hilbert, vérifiant :

$$ (5) \qquad H_0^1(\Omega) \subset V \subset H^1(\Omega) \qquad H = L^2(\Omega). $$

Selon le choix de V, on résoudra ainsi (1) ou son analogue avec conditions de Dirichlet, ou un problème intermédiaire.

On notera par :

$$ (6) \qquad \left\{ \begin{array}{ll} \| \| (( , )) & \text{la norme et le produit scalaire dans } V \\ | | ( , ) & \text{la norme et le produit scalaire dans } H \\ \| \|_* (( , ))_* & \text{la norme et le produit scalaire dans } V'. \end{array} \right. $$

On rappelle que :

$$ (7) \qquad ((u,v))_* = ((D^{-1}u, D^{-1}v)) \quad \forall u, v \in V' $$

où D est l'isomorphisme canonique de V sur V' défini par :

$$ (8) \qquad \forall v \in V, Dv \in V' \text{ est défini par : } (Dv,u)_{V'V} = ((v,u)) \quad \forall u \in V. $$

On identifiera enfin V à une partie de H et H à son dual H', et on aura la situation suivante :

$$ (9) \qquad \begin{array}{ccc} & \begin{array}{c} -1 \\ D \end{array} & \\ V \subset H \subset V' & \\ \begin{array}{c} \rightarrow \\ \downarrow \\ D \end{array} & \end{array} $$
(1)

Définissons maintenant un opérateur non linéaire $A$ de $V$ dans $V'$ par :

$$ (10) \qquad \left\{ \begin{array}{l} \forall u,v \in V : \\ a(u,v) = (A(u),v) = \int_{\Omega} \sum_{i=1}^{n} a(u) \frac{\partial u}{\partial x_i} \frac{\partial v}{\partial x_i} dx \end{array} \right. $$

où $a$ est une fonction de $\mathbb{R} \to \mathbb{R}$ vérifiant :

$$ (11) \qquad \mathbb{H} \alpha, M > 0 \text{ tels que : } 0 < \alpha \le a(\zeta) \le M \quad \text{pp sur } \mathbb{R}. $$

Prenons enfin :

$$ (12) \qquad \left\{ \begin{array}{l} f \in L^2(Q) = L^2(OT;H) \\ y \in L^2(\Omega) = H \end{array} \right. $$

On considère alors l'équation non linéaire suivante :

$$ (13) \qquad \left\{ \begin{array}{l} \left(\frac{dy}{dt}(t), v\right) + a(y(t),v) = (f(t),v) \quad \forall v \in V, \text{ pp sur } ]OT[ \\ y(0) = y_0 \end{array} \right. $$

Nous allons montrer le

**Théorème 1** : Avec les hypothèses et notations (4) à (12), le problème (13) admet une solution unique vérifiant :

$$ (14) \qquad y \in L^2(OT;V) \cap L^\infty(OT;H) $$

$$ (15) \qquad \frac{dy}{dt} \in L^2(OT;V') $$

**Démonstration** :

i) **existence** : soit $w_1, \dots, w_m$ une "base spéciale" de $V$ vérifiant :

$$ ((w_j,v))_V = \lambda_j(w_j,v) \quad \forall v \in V $$

<footnote>(1) On notera par a à la fois la fonction inconnue $\zeta \to a(\zeta)$ de $\mathbb{R} \to \mathbb{R}$ et l'application de $V \times V \to \mathbb{R}$ qui lui est associée par (10).</footnote>
Définissons ensuite le problème approché dans $V_m = [w_1, \dots, w_m]$ par

$$ (16) \quad \left\{ \begin{array}{l} \frac{dy}{dt} - w_j + a(y_m, w_j) = (f, w_j) \quad j=1,2,\dots,m, \quad t \in [0, T_f] \quad \text{OT} \\ y_m(0) = y_{om} \quad \longrightarrow \quad y_o \quad \text{dans } L^2(\Omega) \quad \text{lorsque } m \to \infty. \end{array} \right. $$

Ce système d'équations différentielles non linéaires admet une solution locale $y_m(t)$ pour $t \in [0, T_f[, \text{ avec } T_f \le T$.

**Estimation à priori I** : En remplaçant $w_j$ par $y_m(t)$ dans (16) et en intégrant de 0 à t on trouve facilement :

$$ (17) \quad \begin{aligned} \|y_m\|_{L^2(OT_f;V)} &\le C \\ \|y_m\|_{L^\infty(OT_f;H)} &\le C \end{aligned} $$

d'où il résulte en particulier que la solution $y_m$ de (1) existe sur tout l'intervalle $[OT]$, c.a.d : $T_f = T$

**Estimation à priori II** : Soit $P_m$ le projecteur de $H$ sur $V_m$ dans $H$. Graçons aux choix de la base $w_j$, $P_m$ coïncide avec le projecteur de $V$ sur $V_m$ dans $V$. Comme $P_m \in \mathcal{L}(V;V)$, avec $\|P_m\|_{\mathcal{L}(V;V)} = 1$, on a par transposition :

$$ (18) \quad P_m \in \mathcal{L}(V';V') \quad \text{avec} \quad \|P_m\|_{\mathcal{L}(V';V')} = 1 $$

Or la première équation de (16) s'écrit :

$$ \frac{dy_m}{dt} + P_m A(y_m) = (f, P_m v) \quad \forall v \in V. $$

Donc (16) se réécrit :

$$ (19) \quad \left\{ \begin{array}{l} \displaystyle P_m \frac{dy_m}{dt} + P_m A(y_m) = P_m f \quad \text{dans } V' \text{ pour presque tout } t \in [0, T] \quad \text{OT} \\ y_m(0) = y_{om} \\ \\ \text{Or } \frac{dy_m}{dt} (t) \in V_m. \text{ On en déduit facilement que } P_m \frac{dy_m}{dt} = \frac{dy}{dt} \in V_m \subset V'. \end{array} \right. $$
Donc (19) se réécrit :

$$
(20) \quad \begin{cases} \frac{dy_m}{dt} + P_m A(y_m) = P_m f \quad \text{dans } V' \text{ pour presque tout } t \in ]0,T[ \\ y_m(0) = y_{om}. \end{cases}
$$

De (17), (10), (11) et (12), on déduit facilement que $A(y_m)$ est dans
un borné de $L^2(OT;V')$. D'où, grâce à (20) :

$$
(21) \qquad \left\| \frac{dy_m}{dt} \right\|_{L^2(OT;V')} \le C
$$

Estimation à priori III : Soit $\beta$ la primitive de $a$ s'annulant à l'originale :

$$
\beta(\zeta) = \int_{0}^{\zeta} a(\eta) \, d\eta \qquad \forall \zeta \in \mathbb{R}.
$$

La fonction $\beta$ est donc *continue de $\mathbb{R} \to \mathbb{R}$*, et vérifie, grâce à (11) :

$$
(22) \quad \alpha |\zeta| \le |\beta(\zeta)| \le M |\zeta| \quad \forall \zeta \in \mathbb{R}.
$$

On déduit alors facilement de (17) (11) et (22) que :

$$
(23) \qquad \| \beta(y_m) \|_{L^2(OT;V)} \le C.
$$

De ces estimations à priori, on déduit qu'il existe une sous suite $y_\mu$ de
la suite $y_m$ telle que :

$$
(25) \quad y_{\mu} \to y \quad \text{dans } L^2(0,T;V) \text{ faible}
$$

$$
(26) \quad \frac{dy_\mu}{dt} \to \frac{dy}{dt} \quad \text{dans } L^2(OT;V') \text{ faible}
$$

$$
(27) \quad y_{\mu} \to y \quad \text{dans } L^2(Q) \text{ fort et presque partout}
$$

$$
(28) \quad \beta(y_\mu) \to \chi \quad \text{dans } L^2(OT;V) \text{ faible}
$$

De (27) (28) et du fait que la fonction $\beta$ est continue par rapport à $y$,
il résulte que :

$$
(29) \qquad \chi = \beta(y).
$$
et on peut passer à la limite dans (16)

ii) **Unicité** : Soient $y_1$ et $y_2$ deux solutions de (13). En posant :

$$ (30) \qquad v = y_1 - y_2 $$

il vient, d'après (13) :

$$ (31) \qquad \left\{ \begin{array}{l} \left( \frac{dw}{dt} , v \right) + \left( A(y_1(t)) - A(y_2(t)), v \right) = 0 \quad \forall v \in V \quad \text{pp sur } [0T] \\ w(0) = 0 \end{array} \right. $$

Prenons dans (31) $v=D^{-1}w(t)$. Il vient

$$ (32) \qquad \left\{ \begin{array}{l} \frac{1}{2} \frac{d}{dt} ((w(t),w(t)))_* + (A(y_1(t))-A(y_2(t)),D^{-1}w(t)) = 0 \\ w(0)=0 \end{array} \right. $$

Compte tenu de (10) et (8), (32) se réécrit :

$$ (33) \qquad \frac{1}{2} \frac{d}{dt} ((w(t),w(t)))_* + \int_{\Omega} (\beta(y_1)-\beta(y_2))(y_1-y_2)d\mathbf{x} = \int_{\Omega} (\beta(y_1)-\beta(y_2)) v \, d\mathbf{x} $$

d'où, grâce à (11) :

$$ \frac{1}{2} \frac{d}{dt} ((w(t),w(t)))_* + \alpha \int_{\Omega} |y_1 - y_2|^2 \, dx \le M \int_{\Omega} |y_1 - y_2| \, |v| \, dx $$

d'où :

$$ (34) \qquad \left\{ \begin{array}{l} \displaystyle \frac{1}{2} \frac{d}{dt} ((w(t),w(t)))_* + \frac{\alpha}{2} \int_{\Omega} |y_1 - y_2|^2 \, dx \le \frac{M^2}{2\alpha} |v|^2 L^2(\Omega) \le \frac{M^2}{2\alpha} ((w(t),w(t)))_* \\ w(0) = 0 \end{array} \right. $$

d'où le résultat en utilisant le lemme de Gronwall.

On a ensuite un théorème utilisant le principe du maximum.

**Théorème 2** : Avec les hypothèses et notations (4) à (12), ainsi que :

$$ (37) \qquad f \in L^{\infty}(Q) , y_o \in L^{\infty}(\Omega) $$

la solution de (13) vérifie (14) et (15) et :

$$ (38) \qquad y \in L^\infty(Q). $$
Démonstration : y étant la solution du théorème 1, posons :

$$ \dot{a}(x,t) = a(y(x,t)) $$

On voit alors que y est la solution faible classique d'un système parabolique linéaire avec conditions au bord homogènes, pour lequel le principe du maximum s'applique. ■
II - DERIVABILITE DE L'APPLICATION $a \to y$.

Posons :

$$
(39) \quad \mathcal{A} = \{ a \in C^1(\mathbb{R}) \mid a \text{ et a' soient bornés sur } \mathbb{R} \}
$$

$$
(40) \quad \mathcal{A}_c = \{ a \in \mathcal{A} \mid \exists \alpha, M > 0 \ \text{tq, } o < \alpha \le a(\zeta) \le M_{pp} \text{ sur } \mathbb{R} \}
$$

Alors, $f$ et $y_o$ étant donnés par (37), le système (13) admet une
solution unique (théo. 2) pour toute fonction $a \in \mathcal{A}_c$.

Nous voulons montrer que l'application $a \to y$ ainsi définie est
dérivable de $\mathcal{A}_c$ dans $L^2(Q)$. Nous aurons besoin pour cela de nous placer en un
point $a \in \mathcal{A}_c$ tel que la solution $y$ correspondante vérifie :

$$
(41) \quad \frac{\partial y}{\partial x_i} \in L^\infty(Q) \quad i=1,2,\dots,n
$$

Le problème de trouver les hypothèses de régularité à faire sur $a$, $f$ et $y_o$ pour avoir (41) est un problème ouvert.

1°) Continuité de l'application $a \to y$ :

**Théorème 3** : Avec les hypothèses et notations (4) à (12), (37), (39) (40) (41) l'application $a \to y$ définie par (13) est continue de $\mathcal{A}_c$ dans $L^p(Q)$ fort pour $p \ge 2$.

Si en plus $n \le 3$ on a les majorations

$$
(42) \quad \left\{
\begin{array}{ll}
\| \Delta y \|_{L^p(Q)} & \le \text{cste } \| \delta a \| \quad 2 \le p \le \frac{10}{3} \\
\\
\| \Delta y \|_{L^p(Q)} & \le \text{cste } \| \delta a \|^{3p} \quad \frac{10}{3} \le p < +\infty
\end{array}
\right.
$$

**Démonstration** : Soit $a \in \mathcal{A}_c$ fixé, et $\delta a \in \mathcal{A}$ un accroissement de $a$ tel que $a+\delta a \in \mathcal{A}_c$. Il en résulte pour $y$ un accroissement $\Delta y$. Posons :

$$
(43) \qquad \tilde{a} = a + \delta a \qquad \tilde{y} = y + \Delta y.
$$
En écrivant (13) pour les deux valeurs a et $\tilde{a}$ du paramètre, en prenant
$v = \Delta y(t)$, et en soustrayant les deux équations obtenues il vient :

$$
\frac{1}{2} \frac{d}{dt} |\Delta_y|^2 + \int_{\Omega} \tilde{a}(\tilde{y}) \sum_{i=1}^{n} \left( \frac{\partial \Delta y}{\partial x_i} \right)^2 dx + \int_{\Omega} (\delta a(\tilde{y}) + a(\tilde{y}) - a(y)) \sum_{i=1}^{n} \frac{\partial y}{\partial x_i} \frac{\partial \Delta y}{\partial x_i} dx = 0.
$$

D'où en utilisant le théorème des accroissements finis, pour a assez
petit :

$$
\frac{1}{2} \frac{d}{dt} |\Delta y|^2 + \frac{\alpha}{4} \sum_{i=1}^{n} \left| \frac{\partial \Delta y}{\partial x_i} \right|^2 \le \frac{1}{\alpha} \sup_i \left\| \frac{\partial y}{\partial x_i} \right\|_{L^\infty(Q)} x \left| \delta a(\tilde{y}) + a'(y + \theta \Delta y) \Delta y \right|^2
$$

d'où :

$$
(44) \quad \frac{1}{2} \frac{d}{dt} |\Delta y|^2 + \frac{\alpha}{4} \sum_{i=1}^{n} \left| \frac{\partial \Delta y}{\partial x_i} \right|^2 \le c_1 \| \delta a \|^2 + c_2 |\Delta y|^2
$$

d'où, en utilisant le lemme de Gronwall :

$$
(45) \quad \| \Delta y \|_{L^{\infty}(OT; L^2(\Omega))} \le \text{cste} \| \delta a \|
$$

$$
(46) \quad \| \Delta y \|_{L^2(OT;V)} \le \text{cste} \| \delta a \|.
$$

On a d'autre part, d'après le théorème 2 :

$$
(47) \quad \| \Delta y \|_{L^{\infty}(Q)} \le \text{cste.}
$$

Alors (45) et (47) montrent que $\Delta y \to 0$ dans $L^p(Q)$ fort, pour tout $p \ge 2$,
et la lère partie du théorème est démontrée.

Si maintenant la dimension d'espace n est inférieure ou égale à 3 on a,
d'après (46) et les inégalités de Sobolev (et grâce à l'hypothèse Ω borné si
η=2) :
$$
(48) \quad \| \Delta y \|_{L^2(\Omega; L^6(\Omega))} \le \text{cste} \| \delta a \| .
$$

En interpolant entre $L^2(\text{OT};L^6(\Omega))$ et $L^\infty(\text{OT};L^2(\Omega))$ on trouve

$$
L \quad \left(\frac{10}{3}\right) = \left[ L^2(\text{OT}; L^6(\Omega)), L^\infty(\text{OT}; L^2(\Omega)) \right] \quad 2/5
$$

Il résulte alors de (45) et (48) que :

$$
(49) \quad \| \Delta y \|_{10/3, L(Q)} \le \text{cste} \| \delta a \| .
$$

En interpolant entre $L^{10/3}(Q)$ et $L^\infty(Q)$, on trouve que, pour $p \ge \frac{10}{3}$ on a

$$
L^P(Q) = \left[ L^\infty(Q), L^{10/3}(\dot{Q}) \right] \frac{10}{3p}
$$

d'où

$$
\| \Delta y \|_{L^P(Q)} \le \| \Delta y \|_{L^\infty(Q)} \xrightarrow{1-\frac{10}{3p} \to} \| \Delta y \|_{L^{10/3}(Q)} \xrightarrow{\frac{10}{3p}}
$$

d'où compte tenu de (47) et (49) :

$$
(50) \quad \| \Delta y \|_{L^p(Q)} \le \text{cste} \| \delta a \|^{10/3p} , \quad p \ge \frac{10}{3}
$$

Remarque 1 : La lère partie du théorème 3 peut se démontrer sous les seules hypothèses $a \in L^\infty(R)$, $\delta a \to 0$ dans $L^\infty(R)$, et sans la condition (41) de régularité sur les dérivées de y en espace, en multipliant par $v = D^{-1}\Delta y$, on trouve alors
$|\Delta y|^2 \le d|\delta a|$, d'où le résultat annoncé, mais c'est insuffisant pour la suite

2°) Etude directe du système obtenu par différentiation formelle de (1).

Ce calcul nous amène à étudier directement les systèmes de la forme :
$$
(51) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u}{\partial t} - \sum_{i=1}^{n} \frac{\partial^2}{\partial x_i^2} (c(x,t)u) = h \quad \text{dans } Q \\[1em]
\displaystyle \frac{\partial}{\partial \bar{h}} (c(x,t)u) = 0 \quad \text{sur } \Sigma \\[1em]
\displaystyle u(x,0) = u_o(x) \quad \text{dans } \Omega.
\end{array}
\right.
$$

où c est une fonction de Ω x ] OT[ → R vérifiant :

$$
(52) \quad \mathfrak{H}\alpha, M \in \mathbb{R} \text{ tels que : } 0 < \alpha \le c(x,t) \le M \quad \text{pp sur } Q
$$

et où h vérifie

$$
(53) \qquad h \in L^2(Q).
$$

Nous allons montrer que (51) admet une solution "très faible" en prenant
V' comme espace pivot. Pour simplifier les notations, nous n'identifierons pas V'
à son dual V, ce qui nous permettra de garder H identifié à H'. Nous utiliserons
lorsque besoin sera, l'isomorphisme canonique D de V sur V' défini en (8).

Définissons sur H une forme bilinéaire continue par :

$$
(54) \quad \left\{
\begin{array}{@{}l@{}}
\forall \vec{u}, v \in H, \quad \forall t \in [0,T] : \\
\qquad c(t;u,v) = \int_{\Omega} c(x,t)u(x)(v(x) - D^{-1}v(x)) \,dx
\end{array}
\right.
$$

Donnons nous ensuite une famille de formes linéaires continues L(t) sur
H, telles que :

$$
(55) \qquad L \in L^2(\text{OT};H)
$$

et considérons le problème :

$$
(56) \quad \left\{
\begin{array}{l}
\displaystyle \frac{d}{dt} ((u,w))_* + c(t;u,w) = L(t) \quad w \quad \forall w \in H, \text{ pp sur } ]OT[ \\[1em]
\displaystyle u(0) = u_0 \in V'.
\end{array}
\right.
$$
On vérifie facilement que si, dans (51), toutes les fonctions sont régulières, alors la solution u de (51) vérifie bien (56) avec $L(t) w = ((h(t),w))_*$, et $V = H^1(\Omega)$.

On a alors le

**Théorème 4** : Avec les hypothèses et notations (52) à (55), le système (56) admet une solution unique u vérifiant

$$ (57) \qquad u \in L^2(\Omega) $$

$$ (58) \qquad \frac{d}{dt} D^{-1} u \in L^2(\Omega) $$

et dépendant continuellement des seconds membres.

**Démonstration** : C'est un résultat classique, car $c(t;u,u)$ vérifie :

$$ c(t;u,u) = \underbrace{\int_{\Omega} c(x,t) u^2(x) dx}_{A} - \underbrace{\int_{\Omega} D^{-1} u(x) dx}_{B} $$

où

$$ A \geq \alpha |u|_{L^2(\Omega)}^2 $$

$$ B \leq \frac{\alpha}{2} |u|_{L^2(\Omega)}^2 + K |D^{-1} u|_{L^2(\Omega)}^2 \leq \frac{\alpha}{2} |u|_{L^2(\Omega)}^2 + K \|u\|_*^2 $$

donc :

$$ c(t;u,u) + K \|u\|_*^2 \geq \frac{\alpha}{2} |u|^2. $$
3°) Dérivabilité de l'application a $\rightarrow y$ :

Plaçons nous en un point a $\in \mathcal{A}_c$. Alors, d'après le théorème 1, il existe une fonction unique $y \in L^2(\mathcal{O}T;\mathcal{V})$ telle que :

$$ (59) \quad \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} \int_{\Omega} y v \, dx + \int_{\Omega} \sum_{i=1}^{n} a(y) \frac{\partial y}{\partial x_i} \frac{\partial v}{\partial x_i} \, dx = \int_{\Omega} fv \, dx \quad \forall v \in V \\[1em] \qquad y(0) = y_0 \end{array} \right. $$

Prenons dans (59) $v=\mathcal{D}^{-1}w$. Alors (59) est vraie pour tout $w \in V'$, donc à fortiori pour tout $w \in H$. On a alors :

$$ \left| \begin{array}{l} \displaystyle \frac{d}{dt} ((y(t),w))_* + \int_{\Omega} \sum_{i=1}^{n} \frac{\partial}{\partial x_i} (\beta(y)) \frac{\partial v}{\partial x_i} \, dx = ((f(t),w))_* \quad \forall w \in H \\[1em] \qquad y(0) = y_0 \end{array} \right. $$

ce qui se réécrit, d'après (8) :

$$ (60) \quad \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} ((y(t),w))_* + \int_{\Omega} \beta(y) (w-\mathcal{D}^{-1}w) \, dx = ((f(t),w))_* \quad \forall w \in H \\[1em] \qquad y(0) = 0. \end{array} \right. $$

Donnons à a une variation $\delta a \in \mathcal{A}$ (telle que $a + \delta a \in \mathcal{A}_c$), et soit y la variation de y correspondante. La fonction $\Delta y$ vérifie, d'après (60), en posant $\tilde{y}=y+\Delta y$ :

$$ (61) \quad \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} ((\Delta y(t),w))_* + \int_{\Omega} [\beta(\tilde{y}) - \beta(y)] (w-\mathcal{D}^{-1}w) \, dx = - \int_{\Omega} \delta\beta(\tilde{y}) (w-\mathcal{D}^{-1}w) \, dx \\[1em] \qquad \Delta y(0) = 0. \end{array} \right. $$

Définissons ensuite la fonction $\delta y$ par :

$$ (62) \quad \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} ((\delta y(t),w))_* + \int_{\Omega} a(y(x,t)) \delta y(t) (w-\mathcal{D}^{-1}w) \, dx = - \int_{\Omega} \delta a(y) (w-\mathcal{D}^{-1}w) \, dx \quad \forall w \in H \\[1em] \qquad \delta y(0) = 0. \end{array} \right. $$
D'après le théorème 4, l'équation (62) admet une **solution unique** $\delta y \in L^2(Q)$
**dépendant de façon linéaire continue** de $\delta\beta(y)$ dans $L^2(Q)$, donc **de $\delta a$ dans $a$**.

Posons alors :

$$
(63) \qquad \Delta y = \delta y + \mathcal{R}
$$

Si nous montrons que la fonction $\mathcal{R}$ tend dans $L^2(Q)$ vers zéro plus
vite que $\|\delta a\|$, nous aurons bien montré la dérivabilité (Fréchet) souhaitée.
Or d'après (61) (62), la fonction $\mathcal{R}$ vérifie :

$$
(64) \qquad \left\{
\begin{aligned}
& \frac{d}{dt} \left( (\mathcal{R}(t), w) \right) * + \int_{\Omega} a(y(x,t)) \mathcal{R} (w - D^{-1}w) dx = \\
& \qquad - \int_{\Omega^2} \left[ \frac{1}{2} (\Delta y)^2 \right] a' (y + \theta \Delta y) + \Delta y \delta a (y + \theta' \Delta y) ] (w - D^{-1}w) dx \\
& \mathcal{R}(0) = 0.
\end{aligned}
\right.
$$

D'après le théorème 4, on voit que :

$$
(65) \qquad \begin{aligned} |\mathcal{R}| &\le c \left| \frac{1}{2} (\Delta y)^2 \right| a' (y + \theta \Delta y) + \Delta y \delta a (y + \theta' \Delta y) \right| \\ |\mathcal{R}| &\le c \left| (\Delta y)^2 \right| + \|\delta a\| \left| \Delta y \right|. \end{aligned}
$$

Or d'après le théorème 3 on a :

$$
\begin{gather*}
| \Delta y | \le c || \delta a || \\
\begin{split}
| (\Delta y)^2 | &= || \Delta y ||_4^2 \le c || \delta a || \\
&= c || \delta a ||^3
\end{split}
\end{gather*}
$$

En portant dans (65) il vient :

$$
|\mathcal{R}| \le c ||\delta a|| ( ||\delta a||^3 + ||\delta a|| )
$$

Ce qui prouve bien que $\mathcal{R}$ tend vers zéro dans $L^2(Q)$ plus vite que $\|\delta a\|$.

Nous avons donc démontré le
**Théorème 5** : Avec les hypothèses et notations (4) à (11), (37) (39) (40) (41), l'appli-
cation $a \to y$ définie par (13) est dérivable de $\mathcal{A}_c$ dans $L^2(Q)$, pour $n \le 3$.

La dérivée est l'application qui, à $\delta a \in \mathcal{A}$ fait correspondre $\delta y \in L^2(Q)$
solution de (62), c.a.d. solution très faible (dans le cas $V=H^1(\Omega)$), de :

$$
\text{(66)} \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial \delta y}{\partial t} - \sum_{i=1}^{n} \frac{\partial^2}{\partial x_i^2} (a(y(x,t)) \delta y) = \sum_{i=1}^{n} \frac{\partial}{\partial x_i} \left( \delta a(y(x,t)) \frac{\partial y}{\partial x_i} \right) \quad \text{pp ds Q} \\[2ex]
\displaystyle \frac{\partial}{\partial x} (a(y(x,t)) \delta y) = 0 \\[2ex]
\delta y(x,0) = 0.
\end{array}
\right.
\qquad \blacksquare
$$
### III - ÉTUDE DU PROBLEME INVERSE.

On se propose maintenant d'identifier la fonction $y \to a(y)$ sur le vu d'observations $z_i$ de la valeur moyenne de y sur les domaines $\Omega_i$. Nous prendrons donc :

$$ (67) \qquad z_i \in L^2(0,T) \quad i=1,2,...,K. $$

#### 1°) Position du problème, existence d'une solution

Nous poserons localement, dans ce paragraphe uniquement :

$$ (68) \qquad \mathcal{A} = L^\infty(\mathbb{R}) $$

et nous définissons, pour $\bar{\alpha}$ et $\bar{M} \in \mathbb{R}$ donnés, $0 < \bar{\alpha} < \bar{M}$, un ensemble de coefficients admissibles par :

$$ (69) \qquad \mathcal{A}_{ad} = \{ a \in \mathcal{A} \mid 0 < \bar{\alpha} \le a(\zeta) \le \bar{M} \text{ pp sur } \mathbb{R} \}. $$

Cet ensemble $\mathcal{A}_{ad}$ est une donnée du problème, et reflète les bornes inférieures et supérieures que l'utilisateur impose à son paramètre.

Le second membre $f$ et la valeur initiale $y_0$ étant donnés par (37), le théorème 2 montre que, pour tout $a \in \mathcal{A}_{ad}$, l'équation (13) admet une solution unique vérifiant :

$$ (70) \qquad y \in L^2(OT;V) \cap L^\infty(Q) \qquad \frac{dy}{dt} \in L^2(OT;V') $$

donc en particulier $y \in L^2(Q)$.

On peut donc définir une fonction coût $J$ par :

$$ (71) \qquad \forall a \in \mathcal{A}_{ad} : \quad J(a) = \sum_{i=1}^{K} \int_{O}^{T} \left( \frac{1}{|\Omega_i|} \int_{\Omega_i} y(x,t)dx - z_i(t) \right)^2 dt $$
où y est la solution de (13) calculée pour la valeur a du paramètre.

Le problème inverse se ramène donc à la recherche d'une fonction a ∈ a_{ad}
qui minimise la fonctionnelle J sur la convexe a_{ad}. On a pour cela le :

**Théorème 6** : Avec les hypothèses et notations (4) à (11), (37), (67) à (69), il existe au moins une fonction a ∈ a_{ad} qui réalise le minimum de J sur a_{ad}.

Démonstration :

Soit a_n une suite minimisante de J, et y_n la suite des états correspondants
on a pour tout n :

$$
(72) \quad \left\{
\begin{array}{l}
\displaystyle \left( \frac{dy_n}{dt} , v \right) + \int_{\Omega} \sum_{i=1}^{n} \frac{\partial}{\partial x_i} \beta_m(y_m) \frac{\partial v}{\partial x_i} dx = (f(t),v) \quad \forall v \in V, \text{ pp sur } \Omega \subset V \\
\\
y_m(o) = y_o.
\end{array}
\right.
$$

Grâce à (69) et (70), on peut passer à la limite dans (72) pour une sous-suite, encore notée $a_m$, $y_m$, à condition de montrer que $\beta_m(y_m) \to \beta(y)$, or :

$$
(73) \quad \left| \beta_m(y_m) - \beta(y) \right| \le \left| \beta_m(y_m) - \beta(y_m) \right| + \left| \beta(y_m) - \beta(y) \right|.
$$

Les fonctions $y \to \beta_m(y)$ sont **continues sur R**, **uniformément bornées** sur tout borné de R (cf. (22)), et **convergent partout vers β (car a_n → a dans L^∞(R) faible *)**, et théorème de Lebesgue), donc convergent vers **β uniformèment sur les bornés de R.** Comme $y_m$ est borné dans $L^\infty(Q)$ on a :

$$
(74) \quad \left| \beta_m(y_m) - \beta(y_m) \right| \underset{m \to \infty}{\longrightarrow} 0
$$

On a d'autre part, grâce à (11) et au théorème des accroissements finis :

$$
(75) \quad \left| \beta(y_m) - \beta(y) \right| < M \left| y_m - y \right| \to 0 \quad \text{pour } m \to \infty.
$$

De (73) (74) et (75) on déduit que $\beta_m(y_m) \to \beta(y)$ dans $L^2(Q)$ fort.
On a donc montré que $a_m \to a$ dans $L^\infty(R)$ faible * et $y_m \to y$ dans $L^2(Q)$ fort, et que $a$ et $y$ vérifient (13). Reste à montrer que $J(a_m) \to J(a)$, ce qui résulte de la convergence de $y_m \to y$ dans $L^2(Q)$ fort.

**Remarque 2** : Compte tenu de la non unicité de la solution du problème inverse, l'utilisateur risque de ne pas se satisfaire d'une fonction $a \in L^\infty(R)$ qui minimise $J$, car une telle fonction a des chances d'être assez peu régulière. On est donc amené à restreindre l'ensemble $\mathcal{A}_{\text{ad}}$ à des fonctions plus régulières, en imposant des conditions sur les dérivées de $a$, à préciser par l'utilisateur. Dans l'application numérique proposée plus loin, nous avons utilisé les conditions suivantes :

$$ \mathcal{A} = \{ a \mid a \in L^\infty(R) , a'' \in L^2(R) \} \\ \mathcal{A}_{\text{adm}} = \{ a \in \mathcal{A} \mid 0 < \bar{\alpha} \le a(y) \le M, \gamma_m(y) \le a''(y) \le \gamma_M(y) \} \text{ pp sur } \mathbb{R}. $$

où $\gamma_m$ et $\gamma_M$ sont deux fonctions de $L^2(R)$ données.

Le résultat du théorème 6 se généralise facilement à l'existence d'un minimum de $J$ sur ce convexe.

### 2°) Calcul du gradient de $J$ par rapport à $a$.

Pour avoir un algorithme numérique de minimisation de $J$, nous allons calculer $J'(a)$ et montrer que $J$ admet une dérivée partielle fonctionnelle par rapport à $a$.

Afin de pouvoir utiliser le théorème 5, nous prendrons donc de nouveau $\mathcal{A}$ et $\mathcal{A}_c$, comme en (39) (40). On a alors le

**Théorème 7** : Avec les hypothèses et notations (4) à (11), (37), (39) (40) et (41), la fonctionnelle $J$ définie en (71) est dérivable de $\mathcal{A}_c$ dans $\mathbb{R}$ pour $n \le 3$, et sa dérivée $J'(a)$ est donnée par :
$$
(76) \quad \left\{
\begin{array}{l}
\forall a \in \mathcal{A}_c, \quad \forall \delta a \in \mathcal{A}: \\
\displaystyle J'(a) \delta a = \int_Q \delta a(y(x,t)) \sum_{i=1}^n \frac{\partial y}{\partial x_i} \frac{\partial p}{\partial x_i} dx dt
\end{array}
\right.
$$

où y est la solution **du système direct** (10) et où p est **l'état adjoint**, donné par :

$$
(77) \quad p \in L^2(\text{OT};V) , \quad \frac{dp}{dt} \in L^2(Q) \quad \text{avec} \quad p(t) = D^{-1} \rho(t) \quad \text{pp sur } \partial \Omega
$$

où $\rho \in L^2(Q)$ est l'unique solution de :

$$
(78) \quad \left\{
\begin{array}{l}
- \frac{d}{dt} ((\rho(t,w))_*) + \int_{\Omega} a(y(x,t)) (\rho - D^{-1} \rho) w dx = -2 \sum_{i=1}^{K} \left[ \int_{\Omega_i} \frac{1}{|\Omega_i|} \left( \frac{1}{|\Omega_i|} \int_{\Omega_i} y dx - z_i(t) \right) w dx \right] \\
\qquad \rho(T) = 0 \quad \forall w \in H, \text{ pp sur } \partial \Omega .
\end{array}
\right.
$$

Dans le cas $V=H^1(\Omega)$, par exemple, $p$ est la solution forte, dans $L^2(\text{OT};H^2(\Omega))$
de (en notant par $\chi_i$ la fonction caractéristique de $\Omega_i$) :

$$
(79) \quad \left\{
\begin{array}{l}
- \frac{\partial p}{\partial t} - a(y(x,t)) \sum_{i=1}^{n} \frac{\partial^2 p}{\partial x_i^2} = -2 \sum_{i=1}^{K} \frac{\chi_i}{|\Omega_i|} \left( \frac{1}{|\Omega_i|} \int_{\Omega_i} y(x,t) dx - z_i(t) \right) \quad \text{dans } Q \\
\frac{\partial p}{\partial n} = 0 \quad \text{sur } \Sigma \\
p(x,T) = 0 \quad \text{dans } \Omega
\end{array}
\right.
$$

Démonstration :

D'après le théorème 4, l'équation (78) a une solution unique $\rho \in L^2(Q)$, et telle que $\frac{d}{dt} D^{-1} \rho \in L^2(Q)$, ce qui donne (77). Dans le cas où $V=H^1(\Omega)$, l'équivalence (78), (79) est facile à vérifier (en effet $D^{-1} \rho(t)$ est dans $H^2(\Omega)$ puisque $\rho(t)$ est dans $L^2(\Omega)$).

Enfin, grâce au théorème 5, la fonctionnelle J est dérivable, et la différentielle $\delta J$ de J correspondant à une variation $\delta a$ de $a$ s'écrit :

$$
\delta J = 2 \sum_{i=1}^{K} \int_{Q} \left| \frac{\chi_i}{|\Omega_i|} \left( \frac{1}{|\Omega_i|} \int_{\Omega_i} y(x,t) dx - z_i(t) \right) \right| \delta y(x,t) dx dt
$$
ce qui se réécrit, compte tenu de (78) :

$$ \delta J = - \int_{0}^{T} \left\{ \left(-\frac{d}{dt} D^{\rho-1} \rho, \delta y\right) + \int_{\Omega} a(y(x,t)\delta y (\rho - D^{-1}\rho)) dx \right\} dt $$

d'où, après intégration par parties en temps, compte tenu de (62) et (78)

$$ \delta J = - \int_{0}^{T} \left\{ \left( \frac{d}{dt} D^{-1} \delta y, \rho \right) + \int_{\Omega} a(y(x,t)) \delta y (\rho - D^{-1} \rho) dx \right\} dt $$

soit, d'après (62) :

$$ \delta J = \int_{Q} \delta\beta(y) (\rho - D^{-1}\rho) dx dt $$

ce qui donne (76) compte tenu de (77) et de (8) :

La formule (76) donnant J'(a) permet déjà de calculer facilement le gradient de J lorsque la fonction a est supposée dépendre d'un nombre fini de paramètres réels (cf. §IV-5). Cependant, cette formule peut être transformée pour faire apparaître une "dérivée partielle fonctionnelle" de J par rapport à a.

Posons pour cela :

$$ (80) \qquad \varphi(x,t) = \sum_{i=1}^{n} \frac{\partial y}{\partial x_i}(x,t) \frac{\partial p}{\partial x_i}(x,t) \qquad \varphi \in L^1(Q). $$

Pour tout $\zeta \in \mathbb{R}$ posons :

$$ (81) \qquad Q_{\zeta} = \{ (x,t) \in Q \mid y(x,t) \ge \zeta \} $$

$$ (82) \qquad \gamma(\zeta) = \int_{Q_{\zeta}} \varphi(x,t) \, dx \, dt $$

L'ensemble $Q_ζ$ est défini à un ensemble de mesure nulle près, mais la fonction $γ(ζ)$ est définie sans équivoque.

On a alors le
**Théorème 8** : Avec les hypothèses et notations du théorème 7, plus (80) à (82), la fonction $\gamma$ définie en (82) _est bornée et continue sur R._ Elle définit en particulier une distribution sur $R$, et on a

$$ (83) \qquad \left\{ \begin{array}{l} \forall \delta a \in \mathcal{A} : \\ J'(a) \delta a = \left< \frac{d\gamma}{d\zeta} , \delta a \right> . \end{array} \right. $$

**Démonstration :**

Pour tout $\zeta \in \mathbb{R}$ on a : $|\gamma(\zeta)| \le \|\varphi\|_{L^1(Q)}$, donc $\gamma$ est bornée sur $\mathbb{R}$.
Posons :

$$ Y(\zeta) = \begin{cases} 1 & \text{si } \zeta \ge 0 \\ 0 & \text{si } \zeta < 0 \end{cases} \quad \forall \zeta \in \mathbb{R}. $$

Alors en notant y _un_ représentant de la classe de fonctions y, (82) se réécrit :

$$ (84) \qquad \gamma(\zeta) = \int_Q Y(y(x,t)-\zeta) \varphi(x,t) \, dx \, dt. $$

Soit $\zeta_n \to \zeta, \zeta_n \in \mathbb{R}$. Montrons que $\gamma(\zeta_n) \to \gamma(\zeta)$. Il suffit pour cela, d'après le théorème de Lebesgue, de montrer que $Y(y(x,t)-\zeta_n) \varphi(x,t)$ converge presque partout sur $Q$ vers $Y(y(x,t)-\zeta) \varphi(x,t)$.

Or si nous posons :

$$ S_{\zeta} = \{ (x,t) \in Q \mid y(x,t) = \zeta \} $$

il est clair que $Y(y(x,t))-Q_n \to Y(y(x,t)-\zeta)$ partout sur $Q-S_\zeta$. Si $S_\zeta$ est de mesure nulle, la convergence a lieupresque partout, ce qu'il fallait démontrer.

Si $S_\zeta$ est de mesure non nulle, $\varphi$ est nulle sur $S_\zeta$ d'après (80) à un ensemble de mesure nulle près, donc $Y(y(x,t)-\zeta_n) \varphi(x,t)$ converge vers $Y(y(x,t)-\zeta) \varphi(x,t)$ sur $S_\zeta$, à un ensemble de mesure nulle près, ce qu'il fallait démontrer.
Pour montrer (83), divisons $\mathbb{R}$ en intervalle de longueur $h$, et posons :

$$y_i = ih \quad i \in \mathbb{Z}.$$

On a alors évidemment :

$$J'(a) \cdot \delta a = \sum_{i=-\infty}^{+\infty} \int_{Qy_{i+1}^-}^{Qy_i^-} \delta a(y) \varphi \, dx \, dt$$

d'où, d'après la formule des accroissements finis :

$$J'(a) \delta a = \sum_{i=-\infty}^{+\infty} \delta a(y_i) (\gamma(y_{i+1}) - \gamma(y_i)) + \mathbb{R}$$

$$\text{où :} \qquad \mathbb{R} = \sum_{i=-\infty}^{+\infty} \int_{Qy_{i+1}^-}^{Qy_i^-} \delta a'(y_i + \theta(y-y_i)) (y-y_i) \varphi \, dx \, dt$$

or :

$$| \mathbb{R} | \le \sum_{i=-\infty}^{+\infty} | \delta a | | h | \int_{Qy_{i+1}^-} | \varphi | \, dx \, dt = \| \delta a \| | h | \| \varphi \| L^1(Q)$$

ce qui prouve que $\mathbb{R} \to \mathbb{O}$ avec $h$, et que :

$$(85) \qquad J'(a) \cdot \delta a = \lim_{h \to 0} \sum_{i=-\infty}^{+\infty} \delta a(y_i) (\gamma(y_{i+1}) - \gamma(y_i))$$

ou encore que

$$\begin{align*}
J'(a) \cdot \delta a &= \lim_{h \to 0} \sum_{i=-\infty}^{+\infty} \frac{\delta a(y_i) - \delta a(y_{i-1})}{h} \times \gamma(y_i) \times h \\
&= \int_{\mathbb{R}} -\delta a'(y) \gamma(y) \, dy \quad \forall \delta a \in \mathcal{A}
\end{align*}$$

ce qui est bien (83).
**Remarque 3** : La démonstration de ce théorème est constructive, et c'est l'analogue discret de (85) que l'on utilisé au paragraphe V lors des applications numériques. ■

**Remarque 4** : La solution $y(x,t)$ étant dans $L^\infty(Q)$, on a :

$$ (86) \qquad y_{\min} \le y(x,t) \le y_{\max} \quad \text{pp sur } Q. $$

Ce qui, joint à (81) (82) montre que :

$$ (87) \qquad \left\{ \begin{array}{l} \gamma(\zeta) \equiv 0 \quad \forall \zeta \ge y_{\max} \\ \gamma(\zeta) \equiv \int_Q \varphi \, dx \, dt \quad \forall \zeta \le y_{\min} \end{array} \right. $$

On en déduit que

$$ (88) \qquad \operatorname{Supp} \frac{d\gamma}{d\zeta} \subset [y_{\min}, y_{\max}] $$

c.a.d. que seules influent sur J les valeurs de a aux points $y \in [y_{\min}, y_{\max}]$. ■
IV - DISCRETISATION ET RESOLUTION NUMERIQUE.

Nous avons traité le cas d'un problème à une dimension d'espace.

1°) **Discrétisation de l'équation directe (1) :**

a) **Discrétisation en espace**

Si $\Omega = ]0,1[$, considérons un réseau de pas h. Les points du réseau sont de la forme $(i-0,5)h, 1 \le i \le Nx$. Il y a donc $Nx$ points inérieurs, $h = \frac{1}{Nx}$. La discrétisation de l'équation directe (1) par un schéma aux différences finies donne le système suivant :

$$ (89) \qquad \begin{aligned} \frac{\partial y_i}{\partial t} - \frac{1}{h^2} \left( a(y_{i-\frac{1}{2}}) y_{i-1} - (a(y_{i-\frac{1}{2}}) + a(y_{i+\frac{1}{2}})) y_i \right) \\ \qquad \qquad + a(y_{i+\frac{1}{2}}) y_{i+1} = f_i \quad \text{pour } 1 \le i \le Nx \\ y_o = y_1 \cdot y_{Nx+1} = y_{Nx} \end{aligned} $$

Nous prendrons comme approximation de $a(y_{i+\frac{1}{2}})$ :

$$ (90) \qquad a(y_{i+\frac{1}{2}}) = a(\frac{1}{2} (y_i + y_{i+1})). $$

b) **Discrétisation en temps :**

L'équation étant non linéaire, nous avons utilisé le schéma suivant du type prédicteur-correcteur. L'intervalle en temps étant $[0,1]$, les points de discrétisation en temps sont de la forme $n \Delta t, 0 \le n \le Nt$ où $\Delta t$ est le pas en temps, $\Delta t = \frac{1}{Nt}$.

$$ (91) \qquad \frac{\hat{y}_i - y_i}{\Delta t} - \frac{\partial}{\partial x} \left( a(y_i^n) \frac{\partial}{\partial x} \left( \frac{y_i^n + \hat{y}_i}{2} \right) \right) = f_i^{n+\frac{1}{2}} $$

$$ (92) \qquad \frac{y_{i+1}^n - y_i^n}{\Delta t} - \frac{\partial}{\partial x} \left( a\left(\frac{y_i^n + \hat{y}_i}{2}\right) \frac{\partial}{\partial x} \left( \frac{y_{i+1}^n + \hat{y}_{i+1}}{2} \right) \right) = f_i^{n+\frac{1}{2}} $$

avec $y^o = y(x,o)$ donné.
Si $f(x,t)$ n'est pas connu au point $n+1/2$ on prendra

$$ (93) \qquad f^{n+\frac{1}{2}} = \frac{1}{2}(f^n + f^{n+1}). $$

**Remarque 5** : Nous avons aussi testé le schéma de type Crank-Nicholson :

$$ \frac{y^{n+1} - y^n}{\Delta t} - \frac{\partial}{\partial x} \left( a \left( \frac{y^n + y^{n+1}}{2} \right) \frac{\partial}{\partial x} \left( \frac{y^n + y^{n+1}}{2} \right) \right) = f^{n+\frac{1}{2}} $$

Pour n fixé, l'équation discrétisée en espace est non-linéaire. Elle est
de la forme :

$$ F_i(y_1, \dots, y_{N_x}) = 0 \quad \text{pour } 1 \le i \le N_x $$

Nous avons employé la méthode de sur-relaxation associée à la méthode de
Newton. L'algorithme est le suivant :

1 - On se donne une valeur initiale de y soit $y_i^0$, $1 \le i \le N_x$.

2 - i variant de 1 à Nₓ, on détermine par la méthode de Newton aᵢ tel que

$$ F_i(y_1, \dots, y_{i-1}, a_i, y_{i+1}, \dots, y_{N_x}) = 0 $$

$a_i$ est la limite de la suite :

$$ a_i = a_i^k - \frac{F_i(y_1, \dots, a_i, y_{i+1}, \dots, y_{N_x})}{\frac{\partial F_i}{\partial y_i}(y_1, \dots, a_i, \dots, y_{N_x}^m)} $$

avec $a_i^0 = y_i^m$.

Puis on définit :

$$ y_i^{m+1} = \omega a_i + (1-\omega)y_i^m $$
Cette méthode a donné des résultats satisfaisants lorsque $a(y)$ est défini sous forme paramétrique. Mais lorsque $a(y)$ est recherché sous forme libre, la méthode devient instable car la méthode de Newton peut diverger du fait des discontinuités de $a(y)$.

## 2°) **Discrétisation de l'équation adjointe (79):**

### a) **Discrétisation en espace :**

L'équation adjointe linéaire discrétisée par différences finies s'écrit :

$$ (94) \qquad -\frac{\partial p_i}{\partial t} - \frac{1}{h^2} a(y_i(p_{i-1} - 2p_i + p_{i+1})) = -2(y_i - zd_i) \quad 1 \le i \le Nx. $$
$$ p_o = p_l \cdot p_{Nx+1} = p_{Nx}. $$

### b) **Discrétisation en temps :**

Le schéma de Crank-Nicholson appliqué à l'équation adjointe permet d'obtenir :

$$ (95) \qquad \frac{p^{n+1} - p^n}{\Delta t} - \frac{1}{2} a \left( \frac{y^n + y^{n+1}}{2} \right) \left( \frac{\partial^2 p^n}{\partial x^2} + \frac{\partial^2 p^{n+1}}{\partial x^2} \right) \\ = - (y^n + y^{n+1} - z_d^n - z_d^{n+1}) \\ p_{Nt} = 0 $$
3°) Détermination du gradient "continu" défini par (76) :

Il y a deux façons de calculer un gradient par rapport au paramètre a,
selon qu'on considère les valeurs a_j comme associées à une fonction a(y) (nous par-
lerons alors du gradient "continu") ou bien comme formant un vecteur de R^{M+2} (nous par-
lerons alors de gradient "discret").

Le gradient "discret" sera déterminé au § V-6.

Nous allons étudier ici le gradient "continu", qui a servi à faire la ma-
jeure partie des essais numériques. Nous allons utiliser la formule continue (76),
dans laquelle nous allons discrétiser les fonctions.

On suppose que :

$$
\begin{align*}
& \text{y}_{\min} \le \text{y}(x,t) \le \text{y}_{\max} && \forall (x,t) \in \Omega \times [0,T] \\
& 0 < \text{a}_{\min} \le \text{a}(y) \le \text{a}_{\max} &&
\end{align*}
$$

Soit M le nombre de points de discrétisation de l'intervalle] y_min, y_max [
Soit Δ le pas de discrétisation

$$
(97) \qquad \Delta = \frac{y_{\max}^{-y_{\min}}}{M+1}
$$

La fonction a(y) à identifier est déterminée linéairement par morceaux :

$$
(98) \qquad y_j = y_{\min} + (j-1) \cdot \Delta \qquad 1 \le j \le M+2
$$

$$
(99) \qquad a(y) = \sum_{j=1}^{M+1} \psi_j(y) \chi_j(y)
$$

$$
\text{avec}
\begin{cases}
\psi_j(y) = a_j + \frac{1}{\Delta}(y-y_j)(a_{j+1}-a_j) \\
a_j = a(y_j) \\
\chi_j(y) = \begin{cases} 1 & \text{si } y \in [y_j, y_{j+1}] \\ 0 & \text{sinon} \end{cases}
\end{cases}
$$

(le dernier intervalle $[y_{M+1}, y_{M+2}]$ est pris fermé à droite).
Le gradient de la fonctionnelle des moindres carrés s'écrit :

$$
(100) \qquad J'(a) \cdot \delta a = \int_Q \delta a(y(x,t)) \frac{\partial y}{\partial x} \frac{\partial p}{\partial x} \, dx \, dt
$$

où y est la solution du système direct et p la solution du système adjoint.

Posons

$$
\varphi(x,t) = \frac{\partial y}{\partial x} \frac{\partial p}{\partial x}
$$

On a :

$$
(101) \qquad \delta \psi_j(y) = \delta a_j \left(1 - \frac{y-y_i}{\Delta}\right) + \delta a_{j+1} \left(\frac{y-y_i}{\Delta}\right).
$$

Le gradient s'écrit alors

$$
(102) \quad J'(a) \cdot \delta a = \sum_{i=1}^{N_x} \sum_{k=1}^{N_t} \int_{x_{i-1}}^{x_i} \int_{t_{k-1}}^{t_k} \sum_{j=1}^{M+1} \delta \psi_j(y) \chi_j(y) \varphi(x,t) \, dx \, dt
$$

soit

$$
J'(a) \cdot \delta a = \sum_{j=1}^{M+2} \delta a_j \text{ grad}_j
$$

où

$$
\text{grad}_j = \sum_{i=1}^{N_x} \sum_{k=1}^{N_t} \text{grad}_j^{ik}.
$$

La difficulté réside dans le calcul de la contribution apportée aux composantes du gradient par les termes $\text{grad}_j^{ik}$ :

$$
(103) \quad \text{grad}_j^{ik} = \int_{x_{i-1}}^{x_i} \int_{t_{k-1}}^{t_k} \varphi(x,t) \left[ \chi_j(y) \left(1 - \frac{y-y_i}{\Delta}\right) + \chi_{j-1}(y) \left(\frac{y-y_{j-1}}{\Delta}\right) \right] \, dx \, dt \\
\text{(avec } \chi_0 = \chi_{M+2} = 0\text{)}.
$$

La détermination des composantes $\text{grad}_j^{ik} = \frac{\partial J}{\partial a_j}$ se fait de la manière suivante :

* On effectue un balayage en temps et en espace (k=1,N_t et i=1,N_x). Pour (i,k) fixé, on détermine les contributions $\text{grad}_j^{i,k}$ apportées à chaque composante $\text{grad}_j$ du gradient par la maille (i,k), que l'on rajoute aux contributions des mailles précédentes.
Lorsque le balayage est terminé, on obtient ainsi, d'après (102), chacune des composantes du gradient.

Pour évaluer numériquement $\operatorname{grad}_j^{ik}$ à l'aide de (103), on divise la maille $[x_{i-1}, x_i] \times [t_{k-1}, t_k]$ du réseau de discrétisation en 4 triangles, et on suppose que $y(x,t)$ est linéaire sur chaque triangle. Si nous appelons $a_j^{ik}$ (resp. $a_{j-1}^{ik}$) la partie de la maille située entre les lignes de niveau $y_j$ et $y_{j+1}$ (resp. $y_{j-1}$ et $y_j$), l'expression (103) se réécrit :

$$ (103\text{bis}) \qquad \operatorname{grad}_j^{ik} = \int_{a_j^{ik}} \varphi(x,t) \left(1 - \frac{y-y_j}{\Delta}\right) dx \, dt + \int_{a_{j-1}^{ik}} \varphi(x,t) \frac{y-y_{j-1}}{\Delta} dx \, dt $$

<figure><img src="image_1.png" /><figcaption>(i-1,k-1)</figcaption></figure>

<caption>(i,k-1)</caption>

Le calcul des intégrales apparaissant dans (103bis) est fait de façon approchée en multipliant la valeur moyenne des intégrandes par l'aire du domaine $a_j^{ik}$ ou $a_{j-1}^{ik}$ correspondant.
### 4°) Recherche de a(y) sous forme libre.

L'algorithme consiste à minimiser la fonctionnelle des moindres carrés définie par (71) par la méthode du gradient avec projection lorsque la contrainte $a(y) \ge 0$ n'est pas satisfaite. Le schéma itératif est le suivant :

1. On se donne une valeur initiale arbitraire de la fonction $a(y)$ à identifier, soit $a^0(y)$. (1)

2. Connaissant $a^n(y)$ on résout le système direct (1). On obtient $y^n$.

3. Connaissant $y^n$ on résout le système adjoint (79). On obtient $p_n$.

4. On détermine $a^{n+1}(y)$ par l'analogue discret de

$$a^{n+1} = a^n - \rho \frac{\partial J^n}{\partial a} \quad \rho > 0$$

5. Test d'arrêt. Si négatif aller en 2.

L'indice n correspond à la n<sup>ème</sup> itération de la méthode de gradient.

### Remarque 6 :

La recherche du paramètre optimal $\rho$ à l'étape 4 est effectuée en testant des valeurs en progression géométrique de raison 2. On peut faire ensuite une interpolation parabolique pour obtenir la valeur de $\rho$ définitive. ■

### Remarque 7 :

La fonctionnelle $J(a)$ n'est pas convexe.

C'est pourquoi on ne peut pas prendre $\rho$ constant au cours des itérations de la méthode de gradient. La figure 5 présente la courbe de variation de $\rho$ en fonction des itérations de la méthode de gradient pour l'exemple traité. ■

### 5°) Recherche de a(y) sous forme paramétrique.

L'algorithme du paragraphe précédent correspond à une recherche sous forme libre de $a(y)$. On peut aussi envisager une recherche sous forme paramétrique.

Supposons que $a(y)$ soit de la forme :

<footnote>(1) $a(y)$ est déterminée en les points de discrétisation $y_j$ de $y$.</footnote>
$$ (103\text{ter}) \qquad a(y) = \sum_{i=0}^{p} \beta_i w_i(y) $$

où les $w_i$ sont des fonctions connues de y.

La valeur de $\beta_i^{n+1}$ au cours des itérations est définie par

$$ \beta_{i}^{n+1} = \beta_{i}^{n} - \rho \frac{\partial J}{\partial \beta_{i}} (\beta^{n}) $$

Mais compte tenu de la linéarité de $a$ par rapport aux paramètres $\beta$, l'expression précédente est équivalente à :

$$ a^{n+1}(y) = a^n(y) - \rho \sum_{i=0}^{p} \frac{\partial J}{\partial \beta_i} (\beta^n) w_i(y) $$

où $a^{n+1}$ et $a^n$ sont de la forme (103ter).

On constate que la programmation numérique de l'algorithme est la même que dans le cas d'une recherche libre du fait de la similarité des expressions définissant $a^{h+1}(y)$.

Nous avons considéré dans les exemples numériques la fonction $a$ définie par

$$ a(y) = \sum_{i=0}^{p} \beta_i y^i \quad (\text{polynômes}) $$

## 6°) Introduction de contraintes sur $a(y)$.

a) La fonction $a(y)$ à identifier doit vérifier la contrainte $a(y) \ge \alpha > 0$ afin d'assurer l'existence et l'unicité d'une solution du système direct. Il est possible qu'au cours des itérations de l'algorithme d'identification cette contrainte ne soit plus satisfaite. Si c'est le cas on projette alors la solution obtenue sur le convexe défini par la contrainte $a(y) \ge \alpha$. Lorsqu'on effectue une recherche sous forme libre, le convexe est défini par $a(y_i) \ge \alpha$, pour $i=1, M+2$ et la projection est immédiate. Dans le cas d'une recherche sous forme paramétrique le convexe est défini par

$$ \sum_{i=0}^{p} \beta_i w_i(y) \ge \alpha \quad \forall y \in [y_{\text{min}}, y_{\text{mdx}}] $$

la projection de $\beta = (\bar{\beta}_0, \dots, \bar{\beta}_p)$ sur ce convexe revient à résoudre le problème suivant :

$$ \min_{\beta} \sum_{i=0}^{p} |\beta_i - \bar{\beta}_i|^2 $$
avec $\sum_{i=0}^{p} f_i w_i (y_j) \ge \alpha$ pour $j=1, M+2$.

C'est un problème de minimisation avec des contraintes linéaires que nous avons résolu par la méthode du "gradient réduit".

b) L'introduction de contraintes supplémentaires sur $a(y)$ est envisageable, notamment lors d'une identification sous forme libre.

On peut en effet obtenir dans ce cas une solution présentant de fortes discontinuités. (cf. Paragraphe V.3). Pour réduire ce phénomène, il suffit d'imposer que la dérivée seconde de $a(y)$ soit bornée. L'ensemble des paramètres admissibles est alors un convexe défini par :

$$ (104) \qquad \begin{aligned} C_1 \Delta &\le a(y_{i+1}) - 2a(y_i) + a(y_{i-1}) \le C_2 \Delta, \quad 2 \le i \le M+1 \\ a_{\min} &\le a(y_i) \le a_{\max}, \quad 1 \le i \le M+2 \end{aligned} $$

$(C_1, C_2)$ sont les bornes inférieures et supérieures imposées à la division seconde de $a$. Les autres variables sont définies au paragraphe IV.3.

Le problème inverse se ramène alors à la minimisation de la fonctionnelle $J(a)$ définie en (71) sous les contraintes (104). Ce problème d'optimisation a été résolu par la méthode du gradient réduit, de Franck et Wolfe (cf. CEA [3] ) et du gradient avec projection. Le principe de cette dernière méthode est le suivant :

On définit

$$ a_k = P_c \left( a^n - \rho_k \frac{\partial J^n}{\partial a} \right) , \rho > 0 $$

où $P_c$ représente l'opérateur de projection sur le convexe $C$ défini par (104). On teste plusieurs valeurs de $\rho_k$ et on retient celle telle que $J(a_k)$ est minimale, soit $\rho_{k_o}$. On pose alors

$$ a^{n+1} = a_{k_o} $$

La projection sur le convexe $C$ est effectuée par la méthode du gradient réduit.
7°) Détermination du gradient "discret".

Nous avons jusqu'à présent calculé le gradient de la fonction coût J(a)
à partir des équations continues. Puis pour résoudre numériquement le problème d'iden-
tification, nous avons discrétisé l'équation d'état, l'équation adjointe et le gradient
obtenus. L'avantage de cette manière de procéder est la facilité d'écriture des équa-
tions en continu. L'inconvénient majeur dans le cas du problème traité est la comple-
xité du calcul de gradient lorsque la fonction a(y) est déterminée linéairement par
morceaux et non sous forme paramétrique (cf. IV-3). On peut éviter cette difficulté
en raisonnant directement sur les équations discrètes. On considère alors le système
direct et la fonctionnelle J(a) discrétisés et on calcule le gradient comme dans le
cas continu. Le système adjoint discret s'introduit alors naturellement. D'une manière
générale on n'a pas forcément équivalence entre le discrétisé du système adjoint et
l'adjoint discret (adjoint du système direct discrétisé).

Considérons le système direct discrétisé (1) :

$$
\begin{gather*}
(105) \quad \frac{y^{n+1} - y^n}{\Delta t} + A^{n+1} y^{n+1} = f^{n+1} \quad , n = 0, \dots, N_t \quad \text{avec} \\
A^n = A(a, y^n) \quad \text{définie par } (89) \ (90) \qquad A \in \mathbb{R}^{N_x \times N_x}
\end{gather*}
$$

La fonction coût (71) s'écrit sous forme discrète :

$$
(106) \quad J(a) = \sum_{i=1}^{K} \sum_{n=1}^{N_t} \Delta t (y_j^n(i) - z_i^n)^2
$$
où $j(i) = \text{indice du point d'observation } n^o = i$ (Il y a K points d'observations)

Le système direct (105) peut s'écrire :

$$
(107) \quad \Psi(a,y) = f \quad \text{avec} \quad f \in \mathbb{R}^q, q = N_t \times N_x
$$

(1) Nous prenons un schéma totalement implicite en temps pour faciliter les calculs qui permettent d'obtenir le système adjoint. Nous avions pris en fait un schéma prédicteur-correcteur défini en (91)(92). Mais cela ne modifie pas l'expression obtenue par le gradient..
La dérivée de la fonction coût J(a) définie par (106) est alors donnée par

$$
(108) \quad J'(a) \cdot \delta a = \left(\frac{\partial \Psi}{\partial a}(a,y) \cdot \delta a, p\right)_{\mathbb{R}^q}
$$

où p est solution du système adjoint défini par

$$
(109) \quad \left( \frac{\partial \Psi^*}{\partial y}(a,y) \cdot p, v \right)_{\mathbb{R}^q} = -2 \sum_{i=1}^{K} \sum_{n=1}^{N_t} \Delta t (y_j^{n(i)} - z_i^n) v_j^v(i) \quad \forall v \in \mathbb{R}^q
$$

En explicitant (109) nous obtenons :

$$
(110) \quad \sum_{n=1}^{N_t} (p^n - p^{n+1}) + \Delta t \left( A_n^* p^{n+1} + \Delta t B^n p^n \right)_{\mathbb{R}^q} = -2 \sum_{i=1}^{K} \sum_{n=1}^{N_t} \Delta t (y_j^n - z_i^n) v_j^v(i)
$$

$$
\text{avec } p^{N_t+1} = 0
$$

p est donc solution du système suivant :

$$
\begin{equation}
\begin{split}
& \left( \frac{p^n - p^{n+1}}{\Delta t} + A_n^* p^n + B^n p^n, v \right)_{\mathbb{R}^q} = -2 \sum_{i=1}^K (y_j^n - z_i^n) v_j(i) \\
& \qquad \qquad \qquad \qquad \qquad \qquad \text{N_t+1} \\
& \mathbf{p}^{N_t+1} = \mathbf{0} \qquad \forall v \in \mathbb{R}^q
\end{split}
\tag{111}
\end{equation}
$$

Le système (111) est un système rétrograde en temps.

La matrice Bⁿ est définie par :

$$
(112) \quad (B^n p^n)_i = - \frac{1}{2h^2} \left[ \begin{aligned}
& a'(y_{i-1/2}) (y_{i-1} - y_i) p_{i-1}^n + (a'(y_{i-1/2}) (y_{i-1} - y_i) \\
& + a'(y_{i+1/2}) (y_{i+1} - y_i) p_i^n + a'(y_{i+1/2}) (y_{i+1} - y_i) p_{i+1}^n
\end{aligned} \right]
$$

pour i = 1, Nx avec les conventions (89)(90)

La fonction a(y) étant supposée deux fois dérivable, écrivons la formule
de Taylor au voisinage de yᵢ :

$$
(113) \quad a(y_i) = a\left(\frac{y_i + y}{2} i + 1\right) + \frac{1}{2} a'\left(\frac{y_i + y}{2} i + 1\right) (y_i - y_{i+1}) + O(y_i - y_{i+1})^2
$$

Compte tenu de (112) et (113) on constate que le système (111) correspond
à la discrétisation du système adjoint continu considérée en (94)(95) lorsqu'on néglige les termes d'ordre O(h²) avec h = max_i=1,N_x |y_i - y_{i+1}|.

Explicitons à présent la dérivée de J(a) définie par (108) :

$$ (114) \quad J'(a) \cdot \delta a = \frac{\Delta t}{h^2} \sum_{n=1}^{N_t} \sum_{i=1}^{N_x} \left[ \delta a(y_{i-1/2}^n)(y_i^n - y_{i-1}^n) + \delta a(y_{i+1/2}^n)(y_i^n - y_{i+1}^n) \right] p_i^n $$

ce qui peut s'écrire compte tenu de (99) :

$$ (115) \quad J'(a) \cdot \delta a = \frac{\Delta t}{h^2} \sum_{n=1}^{N_t} \sum_{n=1}^{N_x} p_i^n \left[ \sum_{j=1}^{M+1} (\delta a_j + \frac{1}{\Delta}(y_{i-1/2} - y_j)(\delta a_{j+1} - \delta a_j))(y_i^n - y_{i-1}^n) x_j(y_{i-1/2} + \frac{1}{\Delta}(y_{i+1/2} - y_j)(\delta a_{j+1} - \delta a_j))(y_i - y_{i+1}) x_j(y_{i+1/2}) \right] $$

La détermination des composantes grad_j = ∂J/∂a_j se fait donc de la manière suivante :

- On effectue un balayage en temps et en espace (n=1,N_t et i=1,N_x).

- Pour (n,i) fixé, on détermine j tel que y_{i-1/2}^n ∈ [y_j, y_{j+1}]. On calcule alors les contributions grad_j^{n,i} apportées aux composantes grad_j du gradient par la maille (n,i), que l'on rajoute aux contributions des mailles précédentes :

$$ \text{grad}_j^{n,i} = (y_i^n - y_{i-1}^n) \left(1 - \frac{y_{i-1/2}^n - y_j^n}{\Delta}\right) p_i^n $$

$$ \text{grad}_{j+1}^{n,i} = \frac{1}{\Delta} (y_i^n - y_{i-1}^n) (y_{i-1/2} - y_j^n) p_i^n $$

- Puis, toujours pour (n,i) fixé, on détermine j tel que y_{i+1/2}^n ∈ [y_j, y_{j+1}]. et les nouvelles contributions grad_j^{n,i} :

$$ \text{grad}_j^{n,i} = (y_i^n - y_{i+1}^n) \left(1 - \frac{y_{i+1/2}^n - y_j^n}{\Delta}\right) p_i^n $$

$$ \text{grad}_{j+1}^{n,i} = \frac{1}{\Delta} (y_i^n - y_{i+1}^n) (y_{i+1/2} - y_j^n) p_i^n $$

Lorsque le balayage est terminé on obtient ainsi chacune des composantes du gradient sans avoir eu à évaluer des intégrales comme en (103bis).
Par contre, il faut bien prendre garde au fait que ce gradient contient les dérivées de la fonction coût discrétisée par rapport aux paramètres $a_j$ tels qu'ils interviennent effectivement dans l'équation discrète. La variable y ne prenant qu'un nombre fini ($N_x \times N_t$) de valeur $y_i^n$, il peut arriver que certains $a_j$ n'interviennent pas dans l'équation (105) : le gradient discret correspondant sera nul ! Il faudra donc n'utiliser ce gradient qu'avec des contraintes de régularités sur a, qui entraineront artificiellement les points "oubliés" par le gradient discret.

# V - <u>RESULTATS NUMERIQUES.</u>

## 1) <u>Exemple traité.</u>

Nous avons considéré un problème à <u>une dimension d'espace</u>, et nous avons essayer d'identifier la non-linéarité $a(y)$ à partir de mesures simulées de y.

Donc ici :

$$ \Omega = ]0,1[ \quad \text{et} \quad Q = ]0,1[ \times ]0,1[ $$

et l'équation parabolique s'écrit :

$$ (116) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial y}{\partial t} - \frac{\partial}{\partial x} (a(y) \frac{\partial y}{\partial x}) = \sum_{j=1}^{5} \delta(x-x_j) f_j(t) \quad \text{dans } Q \\[2ex] \displaystyle \frac{\partial y}{\partial x} (0,t) = \frac{\partial y}{\partial x} (1,t) = 0 \quad \forall t \in [0,1]. \end{array} \right. $$

où

$\delta(x,x_j)$ est la masse de Dirac au point $x_j \in [0,1]$. Les équations précédentes peuvent représenter l'évolution de la pression dans un gisement aquifère ou d'hydrocarbures. $f_j(t)$ est alors le débit soutiré à l'instant $t$ au puits de coordonnée $x_j$ et $y(x,t)$ est la pression du liquide. Le problème est alors d'identifier la fonction $a(y)$ sur le vu d'une observation de la pression.

Pour simuler l'observation nous avons choisi pour $a(y)$ la "valeur exacte" suivante :

$$ (117) \qquad a(y) = 0,021 - 0,28 y + 0,7 y^2 $$

puis nous avons intégré le système (116) afin d'obtenir les valeurs de la pression $y(x,t)$. Nous avons alors considéré deux types d'observation :
- une observation **répartie**: la pression est connue en tous les points en espace et en temps.

- une observation **ponctuelle**. La pression est supposée connue en tous les points en temps mais seulement aux cinq puits.
(nous avons aussi traité le cas d'une observation réduite à un seul puits).

Le deuxième type d'observation est plus réaliste du point de vue physique, la pression ne pouvant être mesurée qu'aux puits.

L'observation ainsi obtenue est alors considérée comme donnée et on cherche à retrouver la fonction a(y) définie par (116). Nous avons pris comme pas de discrétisation en espace, en temps et pour a(y) :

$$ 
\begin{array}{l@{\qquad}l}
h = 1/20, & N_x = 20 \\
h_t = 1/40, & N_t = 40 \\
M = 19, & 0.3 \le y(x,t) \le 2
\end{array}
$$

## 2) Recherche sous forme paramétrique (observation répartie ou ponctuelle).

On recherche a(y) sous la forme d'un polynôme du second degré suivant la technique du § IV-5.

### Observation répartie :

On suppose que y varie dans l'intervalle $[0,3;2]$ et a(y) dans l'intervalle $[0,001; 2,3]$. Nous avons pris comme valeur initiale la droite passant par les points extrémaux (cf. figure 1). Les résultats sont identiques à ceux obtenus dans le cas d'une observation ponctuelle (cf. figure 2). La différence réside dans la vitesse de convergence. Il suffit en effet de 20 itérations pour converger. Le temps d'exécution nécessaire est de 4 minutes sur CII 10070.

### Observation ponctuelle :

La pression étant observée seulement aux cinq puits, la convergence est alors plus lente. La figure 2 présente la convergence de a(y) suivant les itérations. Les figures 3 et 4 présentent les courbes des pressions aux cinq puits observées et calculées en intégrant le système avec la fonction a(y) obtenue après identification. La décroissance de la fonctionnelle des moindres carrés J(a) en fonction du nombre d'itérations est représentée par la figure 11. Le temps de calcul nécessaire à l'identification est de 7 minutes sur CII 10070. Il correspond à 32 itérations.
## Conclusion :

La recherche sous forme paramétrique donne des résultats satisfaisants. Lorsque la fonction $a(y)$ est trop irrégulière pour être approchée par un polynôme, une recherche sous forme paramétrique peut être intéressante pour fournir une valeur initiale à la recherche sous forme libre.

### 3) Recherche sous forme libre (observation répartie ou ponctuelle).

On utilise le gradient "continu" défini au paragraphe IV.3.

## Observation répartie :

La valeur initiale étant une droite (cf. figure 1), on recherche $a(y)$ sous forme libre. $a(y)$ est ainsi déterminée aux points de discrétisation de $y$: $y_i \le i \le M+2$. La convergence de $a(y)$ est présentée à la figure 6. Les courbes de pression aux puits obtenues avec les valeurs de $a(y)$ après identification sont semblables à celles des figures 3 et 4. Il faut 31 itérations et 8 minutes de CII 10070.

## Observation ponctuelle :

La valeur initiale est toujours la droite. On suppose connue la pression seulement aux cinq puits. La convergence de $a(y)$ est présentée à la figure 7. Il faut 30 itérations et 9 minutes de CII 10070. Le temps de calcul est supérieur à celui obtenu dans le cas d'une observation répartie car la recherche du paramètre optimal $\rho$ (cf. Remarque 1) a demandé plus de temps. On a tracé sur la figure 5 la courbe donnant la variation de la valeur, de $\rho$ suivant les itérations. La courbe de décroissance de la fonctionnelle $J(a)$ est présentée sur la figure 11.

## Commentaires :

On constate, lorsque la recherche de $a(y)$ est effectuée sous forme libre avec observation ponctuelle ou répartie, des discontinuités dans la courbe $a(y)$. (cf. figures 6 et 7 itération 30). La fonction $y(x,t)$ est discrétisée en 21 points dans l'intervalle $[0,3;2]$. La fonction $a(y)$ est donc déterminée en ces 21 points : $a_i, i=1,2,1$.

Les valeurs $a_2$ et $a_3$ (cf. figures 6 et 7) sont mal ajustées car elles ont peu d'influence dans l'équation directe. Elles sont donc restées à leurs valeurs initiales. La fonction $y(x,t)$ atteint en effet son minimum 0,376 en un seul point, aussi les valeurs $a_2$ et $a_3$ correspondant à $y=0,385$ et $y=0,470$ sont-elles difficiles à identifier.
La difficulté réside surtout dans l'ajustement de $a_{12}$. Elle est inhérente à la nature de $y(x,t)$: au voisinage de $y=1,23$ (correspondant à $a_{12}$), comme on peut le constater sur la figure 8 qui présente les lignes de niveau de $y(x,t)$, la fonction $y(x,t)$ est stationnaire et son gradient $\frac{\partial y}{\partial x}$ est donc nul ou voisin de 0. Par conséquent $\frac{\partial J}{\partial a_{12}}$ est voisin de 0 (cf. (100)) et la valeur $a_{12}$ est lente à s'ajuster.

4) Recherche sous forme libre avec contraintes (observation ponctuelle).

On utilise toujours le gradient "continu" du § IV-3. Pour obvier aux difficultés décrites au paragraphe précédent, nous avons introduit, en plus des contraintes $a(x) \ge 0$, des contraintes sur la dérivée seconde de la fonction à identifier, qui assurent une certaine régularité de la non linéarité et permettent éventuellement d'imposer le sens de sa concavité :

$$C_1 \le a''(y) \le C_2.$$

Nous avons testé pour résoudre le problème d'optimisation avec contraintes les trois algorithmes décrits au paragraphe IV.6, (gradient réduit, projeté et Frarck et Wolfe. C'est cette dernière méthode qui a donné les meilleurs résultats.

Nous avons fait plusieurs essais en faisant varier les valeurs de $C_1$, $C_2$ et la valeur initiale de l'algorithme $a^°(y)$. La valeur exacte de la dérivée seconde est : $a''(y) = 1,4$ pour l'exemple traité (cf. 117).

a) Valeur initiale : droite de la fig. 1 ; contraintes : $|a''(y)| \le 1,4$.

La convergence a lieu en deux itérations. Le choix de $C_1$ et $C_2$ est donc important pour la vitesse de convergence. Les résultats obtenus sont comparables à ceux des figures 2 (itération 32) et 3-4.

b) Valeur initiale ; droite de la fig. 1 ; contraintes : $|a''(y)| \le 5$.

La convergence de $a(y)$ est présentée à la figure 9. On pourra comparer la décroissance de la fonctionnelle $J(a)$ suivant les cas sur la figure 11 (recherche sous forme paramétrique, libre, libre avec contraintes $|a''| \le 5$). Le temps de calcul est de 8 minutes pour 22 itérations.

c) Valeur initiale : droite de la fig. 1 ; contraintes : $0 \le a''(y) \le 5$.

La convergence de $a(y)$ est présentée à la figure 10. On constate une convergence plus rapide. Le temps de calcul pour 19 itérations est de 8 minutes.
d) Valeur initiale : $a(y) \equiv 1$ ; contraintes : $0 \le a''(y) \le 3$.

Dans la mesure où les points extrémaux de la fonction $a(y)$ sont connus, il est préférable de prendre comme valeur initiale la droite passant par les deux points.
L'algorithme converge néanmoins quelle que soit la valeur initiale. Nous avons pris ici $a°(y)$ constant et égal à 1. Les résultats obtenus sont comparables à ceux des essais précédents. On peut difficilement améliorer la valeur de $a(y)$ obtneue à la 9ème itération (cf. figure 12). La valeur extrème correspondant à $y = 2$ est très mal ajustée car c'est la valeur initiale de la pression $y(x,t)$ donc aussi de l'observation $z_d(x,t)$ qui n'est atteinte qu'à l'instant initial $t = 0$. Le temps de calcul est de 4 minutes pour 10 itérations.

5) Recherche sous forme libre avec contraintes, observation ponctuelle réduite à un point, avec ou sans bruit, gradient "continu".

Devant les bons résultats obtenus au paragraphe précédent lors de la recherche de $a(y)$ sous forme libre avec contraintes, nous avons testé cette technique dans des situations plus difficiles : nous avons considéré ici que l'on n'observait la solution $y$, au cours du temps, qu'au puits n°1 et que cette mesure pouvait être bruitée. L'observation est donc la 1ère courbe de la figure n°3, éventuellement bruitée (cf. alors fig. 14). La valeur initiale $a°(y)$ est toujours celle de la figure 1.

a) Pas de bruit, contrainte $0 < a''(y) < 3$.

La figure 13 représente la convergence de $a(y)$. Le résultat est excellent : cela est dû en partie au choix de la contrainte $a''(y) \ge 0$.

b) Bruit gaussien d'écart type 0,05, contrainte $0 \le a''(y) \le 5$.

La figure 14 montre d'une part l'observation bruitée, et d'autre part la non linéarité $a(y)$ obtenue après 6 itérations de gradient, ainsi que la solution $y$ au puits n°1 correspondant à cette valeur de $a(y)$ (comparer avec la solution exacte non bruitée représentée sur la figure 3).
c) Bruit gaussien d'écart type 0,05, contrainte $|a''(y)| \le 5$.

Les résultats, représentés sur la figure 15 après 6 itérations de gradient, sont identiques à ceux de l'essai précédent en ce qui concerne la ressemblance observation-solution calculée, et sont un peu moins bons en ce qui concerne la non linéarité elle-même.

6) Recherche sous forme libre avec contraintes, observation répartie, sans bruit, gradient discret.

Nous avons fait des essais en utilisant le gradient discret défini au paragraphe IV.7. L'avantage, comme nous l'avons dit, réside dans la facilité d'écrire-ture du programme de calcul du gradient. En contre-partie, on constate numériquement un ralentissement de la vitesse de convergence et une moins bonne précision sur la détermination de la fonction a(y) solution du problème d'identification.

Nous présentons sur la figure 16, les résultats des deux essais effectués sur des exemples traités avec le gradient continu.

- Le cas d'une recherche sous forme libre avec contraintes $|a''(y)| \le 5$. (cf. 5.4.a).

- Le cas d'une recherche sous forme libre avec contraintes $0 \le a''(y) \le 3$ et observation ponctuelle réduite à un point (cf. 5.5.a).

On pourra comparer les résultats de la figure 16 avec ceux des figures 9 et 12. La décroissance de la fonctionnelle J(a) suivant le choix du gradient, dans le cadre de l'exemple avec contraintes $|a''(y)| \le 5$, est présentée sur la figure 11.
## VI - CONCLUSION.

D'un point de vue théorique, le problème de l'estimation de la fonction $a(y)$ a été résolu de façon satisfaisante (réserve faite de la nécessité pour les résultats de dérivabilité d'une hypothèse de régularité (41) sur la solution y). En particulier, on a un théorème d'existence de la solution du problème de contrôle associé, ce qui n'est pas le cas pour les coefficients dépendant de $x$ et/ou $t$. On n'a malheureusement pas de théorème d'unicité de la solution du problème inverse, mais les résultats numériques sont sur ce point, beaucoup plus encourageants. Enfin, la méthode est (au moins formellement) extensible à de nombreux autres systèmes.

D'un point de vue numérique, les résultats semblent très intéressants : l'approche classique consistant à représenter la non-linéarité par un polynôme $a$ donné, au moins pour les premières itérations de gradient, une convergence moins rapide que la recherche libre, et surtout que la recherche libre avec contraintes sur les dérivées secondes.

La possibilité d'imposer des contraintes sur la non-linéarité et ses dérivées première et seconde permet de prendre en compte le maximum d'information disponible (par exemple, les ingénieurs connaissent souvent le sens de la concavité de la non linéarité qu'ils veulent estimer). On a ainsi le maximum de chances d'obtenir une solution "physique", sinon proche de la valeur vraie.

D'autre part, l'observation de l'évolution de la fonction $a(y)$ au cours des itérations de gradient permet de se rendre compte de l'influence de $a(y)$ sur le système suivant les valeurs de $y$ , les valeurs prépondérantes s'ajustant plus rapidement.

Enfin, l'algorithme n'est pas trop sensible au choix de la valeur initiale, au nombre de points de mesure de l'observation, et se comporte bien en présence de bruit même important.

Dans l'exemple choisi on a constaté que les fonctions $y \to a(y)$ trouvées étaient proches (et souvent très proches) de leur valeur exacte et ceci même lorsqu'on n'utilisait qu'un seul point d'observation en espace : il doit pratiquement y avoir unicité de $y \to a(y)$ sur l'intervalle des valeurs effectivement atteintes par $y$ aux points de mesure.
Il semble donc que lorsque, en modélisant un système distri-
bué, on rencontre des paramètres qui varient expérimentalement d'un point à
un autre, on ait intérêt non pas à chercher ces paramètres comme fonctions
explicitement de l'espace et du temps, mais plutôt à chercher ces paramè-
tres comme fonctions de certaines variables d'état (à déterminer par des
considérations physiques) lorsque cela est possible : ceci diminue beaucoup
le nombre de degrés de liberté, et on a ainsi plus de chances d'obtenir une
solution unique au problème d'identification.
BIBLIOGRAPHIE

[1] J.C. BRUCH and A.M. ASCE : "Non linear Equation of Unsteady Groundwater flow". Jour. of the Hydraulics Division, 9589, March 1973, HY3.

[2] J.R. CANNON and P. DUCHATEAU : "Determining unknown coefficients in a nonlinear heat conduction problem". SIAM J. Appl. Math. Vol. 24, n°3. May 73.

[3] J. CEA : Optimisation. Dunod, Paris, 1971.

[4] G. CHAVENT : Sur une méthode de résolution du problème inverse dans les équations aux dérivées partielles paraboliques. Note CRAS, Paris, t. 260, Décembre 1969.

[5] G. CHAVENT : Analyse fonctionnelle et identification de coefficients répartis dans les équations aux dérivées partielles. Thèse, Paris, 1971.

[6] G. CHAVENT : Identification of Distributed Parameters. Third IFAC Symposium, The Hague/Delft, The Netherlands. Juin 1973.

[7] G. CHAVENT, M. DUPUY, P. LEMONNIER : History Matching by use of optimal control theory, SPE 4627 presented at the Annual Fall Meeting of the SPE of AIME, Las Vegas, Oct. 1973.

[8] W.H. CHEN and J.H. SEINFELD : "Estimation of spatially varying parameters in partial differential equations". Int. J. Control, 1972, Vol. 15, n°3, 487-495.

[9] M.I. DARBY and N. MORTON : "Numerical solution for the nonlinear penetration of magnetic flux into type II superconductors". J. Comp. Physics 13, 59-69 (1973).

[10] J.R. FREEMAN : "A method for computing the transition from ambipolar to free diffusion in decaying plasma". J. Comp. Physics 13, 59-69 (1973).
[11] W. GIESEL and M. RENGER : "Numerical Treatment of the Unsaturated Water Flow Equation : Comparison of experimental and computed results". Water Resources Research, Vol. 9, n°1, Febr. 1973.

[12] D.D. JUSEPH and T.S. LUNDGREN : "Quasilinear Dirichlet Problems driven by positive Sources". Arch. Rational Mech. Anal, Vol. 49.

[13] J.L. LIONS : Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles, Dunod, 1968.

[14] J.L. LIONS : Quelques méthodes de résolution de problèmes aux limites non linéaires, Dunod-Gauthier Villars, Paris 1969.

[15] H.D. MEYER : "The Numerical Solution of non-linear Parabolic problems by variational methods". SIAM J. Numer. Anal. Vol. 10, n°4, Sept. 1973.

[16] J. SALIN and J.G. SALIN : "Constant line sources of heat in infinite media, whose thermal resistivity are linear functions of the temperature". Int. J. Heat Mass Transfer, Vol. 16, pp 1193-1198.

[17] M.F. WHEELER : "A priori $L^2$ error estimates for Galerkin's approximations to parabolic partial differential equations". SIAM J. Numer. Anal. Vol. 10, n°4, Sept. 1973.
<figure><img src="image_2.png" /><figcaption>Figure 1 : Valeur initiale de a(y) .</figcaption></figure>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2 : Recherche de a(y) sous forme paramétrique (Observation ponct.)</figcaption></figure>
PUITS 1

<figure><img src="image_7.png" /></figure>

* PRESSION OBSERVEE

o PRESSION CALCULEE

PUITS 2

<figure><img src="image_8.png" /></figure>

* PRESSION OBSERVEE

o PRESSION CALCULEE

PUITS 3

<figure><img src="image_9.png" /></figure>

* PRESSION OBSERVEE

o PRESSION CALCULEE

<caption>Figure 3 : Courbes des pressions observées et calculées aux puits.</caption>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 4 : Courbes des pressions (suite)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 5 : Variations du paramètre optimal $\rho$ au cours des itérations.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

Figure 6 : Recherche de a(y) sous forme libre (Observation répartie)

Valeur initiale de la figure 1
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

Figure 7 : Recherche de a(y) sous forme libre (Observation ponctuelle)

Valeur initiale de la figure 1
<figure><img src="image_21.png" /></figure>

<table><caption>Figure 8 : Lignes de niveau de y(x,t)</caption><tr><td>VALEUR MAXIMUM DE LA FONCTION</td><td>0.200E01</td></tr><tr><td>VALEUR MINIMUM DE LA FONCTION</td><td>0.376E00</td></tr><tr><td>DIFFERENCE ENTRE 2 LIGNES CONSECUTIVES</td><td>0.850E-01</td></tr></table>
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>Figure 9 : Recherche de a(y) sous forme libre avec contraintes (|a''|< 5)</figcaption></figure>

Valeur initiale de la figure 1
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 10** : Recherche de a(y) sous forme libre avec contraintes ($0 \le a' \le 5$) : Valeur initiale de la figure 1</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 11 : Décroissance de la fonctionnelle J(a)</figcaption></figure>
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

Figure 12 : Recherche de a(y) sous forme libre avec contraintes (0<a"<3)

Valeur initiale constante
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /><figcaption>**Figure 13** : Recherche de la y) sous forme libre avec contraintes (0 ≤ a" ≤ 3)
Observation réduite à un point. Valeur initiale fig. 1</figcaption></figure>
<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /><figcaption>Figure 14 : Observation en un seul puits avec un bruit gaussien d'écart type 0,05 (à gauche). Valeur initiale de la figure 1, contraintes 0<a"(y)<5</figcaption></figure>
<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /><figcaption>Figure 15 : Observation en un seul puits avec un bruit gaussien d'écart type 0,05 (à gauche). Valeur initiale de la figure 1, contraintes |a''(y)|≤5.</figcaption></figure>
<figure><img src="image_43.png" /></figure>

Forme libre avec contraintes $|a'| \le 5$.

<figure><img src="image_44.png" /></figure>

Forme libre avec contraintes $0 < a' < 3$.

Observation réduite à un point.

<caption>**Figure 16** : Recherche de $a(y)$ avec le gradient "discret".</caption>