% inverse parameter problem
% PDE model
clear global
clear
global b w  nx nt dx dt A Y YM TM C e


bzero = [20;5];  % [5; 35]  or  [20;5]
plotopts = {'linewidth', 1.0};

dx=0.05;nx=21;
tf=2;nt=51;dt=tf/(nt-1);t=0:dt:tf;

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

% the output data
Y=zeros(2,nt);
b=[10;20];bmin=[5; 5];
w=[0; 0];
phipde(0); Y=YM;
figure(1); clf;
plot(t,Y(1,:),t,Y(2,:),'r', plotopts{:} )
ylabel('Temperature output Y(t)');
xlabel('time')
text(0.5, 0.55, 'x_1 = 0.2' )
text(1.0, 0.22, 'x_2 = 0.5' )

% conjugate gradient algorithm

b=bzero;
bb(:,1)=b;itm=25;
iter=1:1:itm+1;
for it=1:itm
   
   % the output model
   % and the LS-criterion
   
   lsi=phipde(0);
   ls(it)=lsi;

   % the adjoint variable
   % and the gradient
   Pw=zeros(nx,1);g=[0; 0];
   for k=nt:-1:2
      Pw=Pw+2*dt*C'*e(:,k-1);
      Pnew=A\Pw; Pw=Pnew;
      dT=TM(:,k)-TM(:,k-1);
      g(1)=g(1)+dx*Pw'*dT;
      g(2)=g(2)-dt*Pw(1)*(1-TM(1,k));
   end
   ng=norm(g)^2;
   
   % the descent direction
   if it==1
      w=-g;
   else
      w=-g+ng/ng1*w1;
   end
 
  
   % the line search
   % and the new iterate

   rom=1e5;
   for j=1:2
      if w(j)<0
         rom=min(rom,(bmin(j)-b(j))/w(j));
      end
   end
   
   ro=fminbnd(@phipde,0,rom);
   fprintf('iter %2d  S=%g  ro=%f\n', it, lsi,ro)

   
   b=b+ro*w;bb(:,it+1)=b;
   g1=g;ng1=ng;w1=w;
end

figure(2); clf;
semilogy( ls, plotopts{:})
grid on
xlabel('iteration');
ylabel('S');
figure(3); clf;
plot(iter,bb(1,:),iter,bb(2,:),'r', 'linewidth', 1.0)
xlabel('iteration')
ylabel('\beta_{1,2}')
text(15, 21, '\beta_2' )
text(10, 11, '\beta_1' )

if( bzero(1) == 5 )
    fig=figure(1);
    print(fig,fullfile(pwd, 'Figure6_35.tif'),'-dtiff','-r600')
   fig=figure(3);
    print(fig,fullfile(pwd, 'Figure6_36.tif'),'-dtiff','-r600')
    fig=figure(2);
    print(fig,fullfile(pwd, 'Figure6_37.tif'),'-dtiff','-r600')
else
    fig=figure(3);
    print(fig,fullfile(pwd, 'Figure6_38.tif'),'-dtiff','-r600')
    fig=figure(2);
    print(fig,fullfile(pwd, 'Figure6_39.tif'),'-dtiff','-r600')
end




function s=phipde(r)
    global b w  nx nt dx dt A Y YM TM C e

    bw=b+r*w;
    fo=dt/(bw(1)*dx*dx);f1=2*fo*bw(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);TM=zeros(nx,nx);
    YM(:,1)=C*Tw;

    for k=2:nt
       Tw(1)=Tw(1)+f1;
       Tnew=A\Tw;YM(:,k)=C*Tnew;
       Tw=Tnew;TM(:,k)=Tw;
    end
    e=Y-YM;
    s=norm(e)^2;
end