## Lipschitz dependence of the coefficients on the resolvent and greedy approximation for scalar elliptic problems

Choulli M., Zuazua E.
Comptes Rendus Mathematique de l’Académie des Sciences, Volume 354, Issue 12, December 2016, Pages 1174–1187DOI: 10.1016/j.crma.2016.10.017

Abstract: We analyze the inverse problem of identifying the diffusivity coefficient of a scalar elliptic equation as a function of the resolvent operator. We prove that, within the class of measurable coefficients, bounded above and below by positive constants, the resolvent determines the diffusivity in an unique manner. Furthermore, we prove that the inverse mapping from resolvent to the coefficient is Lipschitz in suitable topologies. This result plays a key role when applying greedy algorithms to the approximation of parameter-dependent elliptic problems in an uniform and robust manner, independent of the given source terms. In one space dimension, the results can be improved using the explicit expression of solutions, which allows us to link distances between one resolvent and a linear combination of finitely many others and the corresponding distances on coefficients. These results are also extended to multi-dimensional elliptic equations with variable density coefficients. We also point out some possible extensions and open problems.

## Greedy Control

#### The problem

Control of a parameter dependent system in a robust manner.
Fix a control time , an arbitrary initial data , and a final target
Given we aim at determining a family of parameters in so that the corresponding controls are such that for every there exists steering the system to the state within the distance from the target .

#### The system

A finite dimensional linear control system:

(1)

• is a −matrix,
• is a control operator, ,
• is a parameter living in a compact set, of

Assumptions:

• the system is (uniform) controllable for all ,
• system dimension is large

#### The greedy approach

— a Banach space
— a compact subset.
The method approximates by a a series of finite dimensional linear spaces (a linear method).

A general greedy algorithm
The first step
Choose  such that
The general step
Having found , denote

Choose the next element

The algorithm stops
When  becomes less than the given tolerance .



The Kolmogorov width, measures optimal approximation of by a -dimensional subspace.

The greedy approximation rates have same decay as the Kolmogorov widths.

#### Greedy control

Each control can be uniquely determined by the relation

where is the unique minimiser of a quadratic functional associated to the adjoint problem.
This minimiser can be expressed as the solution of the linear system

where is the controllability Gramian

Perform a greedy algorithm to the manifold :

The (unknown) quantity to be maximised by the greedy algorithm is replaced by a surrogate:

The greedy control algorithm results in an optimal decay of the approximation rates.

#### Numerical examples

##### Wave equation

We consider the system (1) with:

The system corresponds to the discretisation of the wave equation problem with the control on the right boundary:

(2)

We take the following values:

The greedy control has been applied with and the uniform discretisation of in values.
The offline algorithm stopped after 24 iterations.
The 20-D controls manifold is well approximated by a 10-D subspace:

##### Heat equation

For , with given by:

and the control operator the system (1) corresponds to the space-discretisation of the heat equation problem with internal grid points and the control on the right boundary:

(3)

The parameter represents the diffusion coefficient and is supposed to range within the set .
The system satisfies the Kalman’s rank condition for any and any target time .
We aim to control the system from the initial state to zero in time .
The greedy control algorithm has been applied for the system of dimension with , and the uniform discretisation of in values.
The algorithm stops after only three iterations.

#### MATLAB CODE

Below you can see and download the MATLAB code executed to simulate and to obtain the data results related to the Greedy control problem. That code is divided into two main files, corresponding to the offline algorithm and the online algorithm respectively, and several additional functions used by both algorithms.

##### The Offline Algorithm

%% initiation

clear f A B y par_greedy sigma G G_nu_phi0 rhs phi0_sel nu_sel;
N=20; %system dimension
k=100; %discretisation constant of a parameter set

x0=zeros(N,1); %initiation of the initial condition

if 0 %wave eq. setting
cc=1; %example number
nt=10*N;  %discretisation of the time interval
a=1; b=10; %parameter interval
d_nu=(b-a)/k;
nu=a:d_nu:b;
n= size(nu,2)-1;
T=3; %target time
epsil=0.5; % greedy control error
%system matrix
M=zeros(N/2);%N should be even
I=eye(N/2);
A_pom=2*eye(N/2);
for i=1:1:N/2-1
A_pom(i,i+1)=-1;
A_pom(i+1,i)=-1;
end
A = zeros(N,N,n+1);
for l=1:1:n+1
A(:,:,l)=[M -I;  nu(l)*((N/2+1)^2)*A_pom  M];
end
%control operator
B=zeros(N-1,1);
B(N)=1*(N+1)^2;
%initial condition
for i=1:1:N/2
x0(i)=1*sin(pi*(i)/(N/2+1));
end
end

if 1 %heat eq. setting
cc=2; %example number
nt=30*N;  %discretisation of the time interval
a=1; b=2; %parameter interval
d_nu=(b-a)/k;
nu=a:d_nu:b;
T=0.1; %target time
epsil=0.0001; %greedy control error
%system matrix
A = zeros(N,N,n+1);
A_pom=2*eye(N);
for i=1:1:N-1
A_pom(i,i+1)=-1;
A_pom(i+1,i)=-1;
end
for l=1:1:n+1
A(:,:,l)=nu(l)*((N+1)^2)*A_pom;
end
%control operator
B=zeros(N-1,1);
B(N)=1*(N+1)^2;
%initial condition
for i=1:1:N
x0(i)=1*sin(pi*(i)/(N+1));
end
end

x1=0*ones(N,1); %target state

%% greedy control algorithm - offline part

k=N; %max. # of iterations
nu_iter=nu;

%Choosing the first parameter
y = zeros(k+1,1);
%calculation of free dynamics for all parameter values
for l=1:1:(k+1)
x=sys_sol(x0,T,nt,-A(:,:,l), zeros(N,1), zeros(1,nt+1));
rhs(:,l)=x1-x(:,nt+1);
y(l)=norm(rhs(:,l)); %error in the first step
if y(l)<epsil
nu_iter(l)=0;
end
end
%removing those parameter values for which we have already obtained
%approximation within a given error from further calculations
for l=size(nu_iter,2):-1:1
if  nu_iter(l)==0
nu_iter(l)=[];
end
end
%selection of the first parameter
par_greedy(1)=argmax(y);
nu_sel(1)=nu(par_greedy(1));

%Choosing subsequent parameters
m=1; sigma(1)=2;
while m<k+1 && sigma(max(1,m-1))>epsil % in each iteration one parameter is chosen
phi0_sel(:,m)=minFindAdjCG(x0,x1,T,nt,A(:,:,par_greedy(m)),B,0.00001,zeros(N,1)); %construction of  the minimisator corresponding to the parameter chosen in the previous iteration
nu_iter=nu;
clear proj

for l=1:1:k+1
%calculation of G_nu  for every
if l==par_greedy(m)
G_nu_phi0(:,m,l)=rhs(:,l);
else
u=B'*p;
x=sys_sol(zeros(N,1),T,nt,-A(:,:,l), B, u);
G_nu_phi0(:,m,l)=x(:,nt+1);
end;
%calculation of the error in the m-th step
proj(:,l)=projection(rhs(:,l),G_nu_phi0(:,:,l));
y(l)=norm(rhs(:,l)- proj(:,l));
if y(l)<epsil
nu_iter(l)=0;
end
end;
%removing those parameter values for which we have already obtained
%aprroximation within a given error from further calculations
for l=size(nu_iter,2):-1:1
if  nu_iter(l)==0
nu_iter(l)=[];
end
end
[par_greedy(m+1), sigma(m)]=argmax(y);
nu_sel(m+1)=nu(par_greedy(m+1)); %selection of the next parameter

m=m+1
end
m=m-1;%m represents the number of chosen parameters
par_greedy(m+1)=[];
nu_sel(m+1)=[];

savefile = 'data.mat';
save(savefile, 'cc','epsil', 'A' , 'B', 'x0', 'x1', 'T', 'nt','a', 'b', 'n', 'sigma','phi0_sel', 'par_greedy', 'proj', 'G_nu_phi0', 'nu_sel','nu','rhs' );

%% graphical part

figure;
%plot of greedy approximation rates
subplot(2,1,1);
pom=1:1:size(sigma,2);
p=plot(pom,sigma);
set(p,'Color','red')
title(['Greedy approximation rates    \newline N=',num2str(N),', n=',num2str(n),', T=',num2str(T),', nt=',num2str(nt),', \epsilon=',num2str(epsil)])
xlabel('Iteration #','Fontsize',14);

%plot of selected paramters
subplot(2,1,2);
plot(nu_sel,'-s');
title(['parameter choice from [',num2str(a),',',num2str(b),']']);
xlabel('Iteration #','Fontsize',14);


##### The Online Algorithm

%% initiation

clear A_nu B_nu Lambda_nu_phi0  elem j elem_sel;
k=size(nu_sel,2);
N=size(x0,1);
nu_online=input(['Choose a value of the parameter between ',num2str(a),' and ',num2str(b), ': '] );

if cc == 1 %System corresp. to the wave eq.
M=zeros(N/2);
I=eye(N/2);
A_pom=2*eye(N/2);
for i=1:1:N/2-1
A_pom(i,i+1)=-1;
A_pom(i+1,i)=-1;
end
A_nu=[M -I;  nu*((N/2+1)^2)*A_pom  M];
B=zeros(N-1,1); %control operator
B(N)=1*(N/2+1)^2;
end

if cc == 2 %System corresp. to the heat eq.
A_nu=2*eye(N);
for i=1:1:N-1
A_nu(i,i+1)=-1;
A_nu(i+1,i)=-1;
end
A_nu=nu_online*((N+1)^2)*A_nu;
B=zeros(N-1,1); %control operator
B(N)=1*(N+1)^2;
end

%% greedy control algorithm - online part

% calculating  (free dynamics)
if element(nu_online, nu)==1 %if the given param. value belongs to the disretised parameter set
[elem,j]=element(nu_online, nu);
rhs_nu=rhs(:,j); %free dynamics of the system is already calculated in the offline part
else
x=sys_sol(x0,T,nt,-A_nu, zeros(N,1), zeros(1,nt+1));%otherwise we calculate free dynamics of the system corresponding to given parameter value
rhs_nu=x1-x(:,nt+1);
end

j_sel=0;
if element(nu_online, nu_sel)==1 %if the given param. value belongs to the set of chosen  snapshots
[elem_sel,j_sel]=element(nu_online, nu_sel);
end

% calculating
for m=1:1:k
if m==j_sel
Lambda_nu_phi0(:,m)=rhs_nu;
else
u=B'*p;
x=sys_sol(zeros(N,1),T,nt,-A_nu, B, u);
Lambda_nu_phi0(:,m)=x(:,nt+1);
end
end

%projection
proj_nu=projection(rhs_nu,Lambda_nu_phi0);
alpha_nu=Lambda_nu_phi0\proj_nu;

%construction of the greedy control
phi0=zeros(N,1);
for i=1:1:k
phi0=phi0+alpha_nu(i)*phi0_sel(:,i);
end
u=B'*p;

%approximation error
x=sys_sol(x0,T,nt,-A_nu, B, u);
error=norm(x1-x(:,nt+1));
fprintf(['The error for nu=',num2str(nu_online),'\n |x(', num2str(T),')|=', num2str(error),'\n']);

%% graphical part
t = linspace(0,T,nt+1);

figure;
subplot(2,1,1);
set(gcf,'DefaultAxesColorOrder',[1 0 0; 0 1 0; 0 0 1;1 0 1;0 0 0]) %red green blue yellow black
lc=N-4; uc=N; %choose the components of the system to be plotted, lc - the lowest comp. to be plotted, uc the most upper one
for i=lc:1:uc
plot(t,x(i,:));
hold all;
end
title(['x_{',num2str(lc),'}(t) - x_{',num2str(uc),'}(t)'],'Fontsize',15);
xlabel('t','Fontsize',14);
hold off;
subplot(2,1,2);
set(gcf,'DefaultAxesColorOrder',[1 0 0]) %  red
plot(t,u);
title(['control for \nu=',num2str(nu_online)],'Fontsize',15);
xlabel('t','Fontsize',14);

if cc==1
%3d plot - wave eq.
figure;
y=0:1/(N/2+1):1;%discretisation of the domain [0,1]
sol=zeros(N/2+2,nt+1); %zero initiation
for i=2:N/2+1
sol(i,:)=x(i-1,:);
end
sol(N/2+2,:)=u;
mesh(y',t',sol','linewidth',2);
title(['v(t, x)'],'Fontsize',15);
xlabel('x','Fontname','Space variable','Fontsize',14,'Fontweight','b');
ylabel('t','Fontname','Time variable','Fontsize',14,'Fontweight','b');
end

if cc==2
%3d plot - heat eq.
figure;
y=0:1/(N+1):1;%discretisation of the domain [0,1]
sol=zeros(N+2,nt+1); %zero initiation
for i=2:N+1
sol(i,:)=x(i-1,:);
end
sol(N+2,:)=u;
mesh(y',t',sol','linewidth',2);
title(['v(t, x)'],'Fontsize',15);
xlabel('x','Fontname','Space variable','Fontsize',14,'Fontweight','b');
ylabel('t','Fontname','Time variable','Fontsize',14,'Fontweight','b');
end



%solves the adjoint problem p'=-A'p on [0,T] with datum p(T)=phi
%nt is the discretisation step of the time interval

p=zeros(size(phi,1),nt+1);
p(:,nt+1)=phi;
dt2=T/(2*nt);
R=speye(size(A))-dt2*A';
M=speye(size(A))+dt2*A';
for k=nt:-1:1
p(:,k)=R\(M*p(:,k+1));
end


##### Function argmax

%returns the (lowest) maximum position and the maximum value of the vector y
function[j,maks]=argmax(y)
j=1;
maks=y(1);
for i=2:1:length(y)
if y(i)> maks;
maks=y(i);
j=i;
end;
end
end




%finding minimisator associated to the adjoint problem of the orig. control problem x'+Ax=Bu
% by using conj. gradient method
%the minimisator satisfies the system  Gq =rhs, where G- Gramiam, rhs=xf-  e(-TA)xi

%xi -initial data,
%xf - final target
%T-target time
%nt -# of discretisation steps in time interval
%eps -precision
%qini - initial guess of the solution

N=size(xi,1); %dimension of the system
x=sys_sol(xi,T,nt,-A, zeros(N,1), zeros(1,nt+1)); %solves the original problem in order to find find e(-TA)xi
rhs=xf-x(:,nt+1);
nb=norm(rhs);
if nb==0, nb=1; end

%Finding G*qini
q=qini;
u_q=B'*p_q;
x_q=sys_sol(zeros(N,1),T,nt,-A, B, u_q); %x- solution of the original problem
Gq=x_q(:,nt+1);

M_inv=eye(N);%preconditioning matrix M_inv, symmetric, posit. def.
r=rhs-Gq; %residual r_0
z=M_inv*r; %if no preconditioning, z =r
h=z; % h stands for the next direction of the method
z_=z;
r_=r;
ng_=r'*r; %error norm
k=0;
nitermax=min(1000,size(xi,1)^2);
while k<nitermax && ng_>eps^2,
%Finding G*h
u_h=B'*p_h;
x_h=sys_sol(zeros(N,1),T,nt,-A, B, u_h);
Gh=x_h(:,nt+1);
alpha=(r'*z)/(h'*Gh);
q=q+alpha*h; %finding (k+1). aproximation
r=r-alpha*Gh;  %next residual r_(k+1)
z=M_inv*r;
beta=z'*r/(z_'*r_);
h=z+beta*h; % finding the next direction of the method
ng=r'*r;
ng_=ng;
z_=z; %store the residual for the next iteration
r_=r;
k=k+1;
end

if 1%check the error Gq-rhs
u_q=B'*p_q;
x_q=sys_sol(zeros(N,1),T,nt,-A, B, u_q);
Gq=x_q(:,nt+1);
err=norm(Gq-rhs); %it corresponds to deviation of the solution at time T from the target: |x(T)-xf|
fprintf(['minFindAdj:: number of iteration reached:', num2str(k),'\n error: ', num2str(err),'\n']);
end
if k>=nitermax,
fprintf(['minFindAdj:: maximal number of iteration reached, error: \n',num2str(ng)]);
end
end


##### Function projection

%provides the orthogonal projection of a vector v to the space spanned by
%columns of matrix B

function  [y]=projection(v,B)
g=orth(B);  %gram-schmidt orthogonalisation
m=size(g,1);
n=size(g,2);

y=zeros(m,1);
%calculation of the projection
for j=1:1:n
y=y+dot(v,g(:,j))*g(:,j);
end

end


##### Function sys_sol


%solves the system x'=Ax+bu on [0,T] with datum x0
%nt is the discretisation step of the time interval
function x=sys_sol(x0,T,nt,A,b,u)
x=zeros(size(x0,1),nt+1);
x(:,1)=x0;
dt2=T/(2*nt);
R=speye(size(A))-dt2*A;
M=speye(size(A))+dt2*A;
B=dt2*b;
for k=1:nt,
x(:,k+1)=R\(M*x(:,k)+B*(u(k)+u(k+1)));
end