% inverse state problem
% PDE model

clear
plotopts={'linewidth', 1.0};

dx=0.05;nx=21;x=0:dx:1;
b=[10;20];
tf=0.4;nt=11;dt=tf/(nt-1);

% two sensor locations
C=zeros(2,nx);C(1,5)=1;C(2,10)=1;

fo=dt/(b(1)*dx*dx);f1=2*fo*b(2)*dx;
A=(1+2*fo)*eye(nx);
for i=2:nx-1
   A(i,i+1)=-fo;A(i,i-1)=-fo;
end
A(1,2)=-2*fo;A(nx,nx-1)=-2*fo;
A(1,1)=A(1,1)+f1;
Tw=zeros(nx,1);

% the initial state 
% to be determined

for k=2:nt
   Tw(1)=Tw(1)+f1;
   Tnew=A\Tw;YM(:,k)=C*Tnew;
   Tw=Tnew;
end

T0=Tw;Y(:,1)=C*Tw;

% the output data with 
% the new initial state
% the input is switched to zero

f1=0; A(1,1)=1+2*fo;
tf=2;nt=51;dt=tf/(nt-1);t=0:dt:tf;
for k=2:nt
   Tnew=A\Tw;Y(:,k)=C*Tnew;
   Tw=Tnew;
end


% conjugate gradient algorithm

itm=100; f=zeros(nx,1);
mu = 0;

for it=1:itm
   
   % the output model
   % and the LS-criterion

   Tw=f;YM(:,1)=C*f;
   for k=2:nt
      Tnew=A\Tw;YM(:,k)=C*Tnew;
      Tw=Tnew;
   end
   e=Y-YM;
   lsi=norm(e)^2;
   ls(it)=lsi;
   
   % the adjoint variable
   % and the gradient
   Pw=zeros(nx,1);
   for k=nt:-1:2
      Pw=Pw+2*dt*C'*e(:,k-1);
      Pnew=A\Pw; Pw=Pnew;
   end
   g=b(1)*Pw-mu*f;  ng=norm(g)^2;
   
   % the descent direction
   if it==1
      w=-g;
   else
      w=-g+ng/ng1*w1;
   end
   
   % the descent step length
   % quadratic case
   dTw=w;dY(:,1)=C*dTw;
   for k=2:nt
      dTnew=A\dTw;dY(:,k)=C*dTnew;
      dTw=dTnew;
   end
   
   ro=0.5*g'*w/norm(dY)^2;
   
   % the new iterate
   
   f=f+ro*w;
   g1=g;ng1=ng;w1=w;
end

figure(1); clf;
semilogy(ls, plotopts{:})
grid on
xlabel('iteration'); ylabel('S');
figure(2); clf;
plot(x,T0,'k-',x,f,'ro', plotopts{:})
legend('Exact', 'estimate')
xlabel('x'); ylabel('Initial T_0(x)')

figure(3); clf;
plot(t,Y,'k-',t,YM,'ro', plotopts{:})
legend('Exact', 'Exact', 'estimate', 'estimate')
xlabel('time t'); ylabel('Temperature output Y(t)');
text(0.5,0.35,'x_1=0.2');
text(1, 0.17, 'x_2=0.5');
 
% fig=figure(3);
% print(fig,fullfile(pwd, 'Figure6_42.tif'),'-dtiff','-r600')
% fig=figure(2);
% print(fig,fullfile(pwd, 'Figure6_41.tif'),'-dtiff','-r600')
% fig=figure(1);
% print(fig,fullfile(pwd, 'Figure6_40.tif'),'-dtiff','-r600')