% Solve the wave equation with boundary feedback control
% u_tt=c^2 u_xx  
% u(0,t)=0, c^2u_x(L,t)=-k u_t(L,t)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;

fc=1; % index for movie frames

L=1; % the length of an interval
c = 1; % wave speed
% Solve the controlled wave equation 
% at the same time for three different
% control gains to compare:
gain=  1; % control gain
gain1= 0.5; % control gain
gain2 =5; % control gain

stepNo = 100000; % total time steps to run 
movieStep =500; % time interval to take movie frames
K=100;  % grid number for interval [0, L]
dx = L/K; % grid step 
x=0:dx:L; % grid
dt = 0.00005; % time step <dx^2/(2k)

u_0 = ones(1,K+1); % initial condition
u_1 = ones(1,K+1); % initial condition
% solution of  uncontrolled equation
compute_u = zeros(1,K+1); 
% solution of controlled equation with gain
compute_ctrlU = zeros(1,K+1); 
% solution of controlled equation with gain1
compute_ctrlU1 = zeros(1,K+1); 
% solution of controlled equation with gain2
compute_ctrlU2 = zeros(1,K+1); 

% energy
unctrl_E = zeros(1, stepNo);
ctrl_E = zeros(1, stepNo);
ctrl_E1 = zeros(1, stepNo);
ctrl_E2 = zeros(1, stepNo);
time_for_energy = zeros(1, stepNo);
 

% initial condition
for k=1:K+1
    u_0(k) = 10*(1-(k-1)*dx)^2 *((k-1)*dx);  
    u_1(k) = (1-(k-1)*dx)*(k-1)*dx;
end

tempU=u_0;
tempU_1=u_0+dt*u_1;
ctrlTempU = u_0;
ctrlTempU_1 =u_0+dt*u_1;

ctrlTempU1 = u_0;
ctrlTempU1_1 =u_0+dt*u_1;

ctrlTempU2 = u_0;
ctrlTempU2_1 =u_0+dt*u_1;
 
for i=1:stepNo
    
    % central difference to solve the wave equation
    compute_ctrlU(2:K) =2*ctrlTempU_1( 2:K)-ctrlTempU( 2:K)...
              +dt*dt* ( ctrlTempU(3:K+1)-2*ctrlTempU( 2:K) ...
              +ctrlTempU(1:K-1) )/(dx*dx);  
          
    compute_ctrlU1(2:K) =2*ctrlTempU1_1( 2:K)-ctrlTempU1( 2:K)...
              +dt*dt* ( ctrlTempU1(3:K+1)-2*ctrlTempU1( 2:K) ...
              +ctrlTempU1(1:K-1) )/(dx*dx);  
          
    compute_ctrlU2(2:K) =2*ctrlTempU2_1( 2:K)-ctrlTempU2( 2:K) ...
              +dt*dt* ( ctrlTempU2(3:K+1)-2*ctrlTempU2( 2:K) ...
              +ctrlTempU2(1:K-1) )/(dx*dx); 
          
    compute_u(2:K) = 2*tempU_1( 2:K)-tempU( 2:K)...
              +dt*dt* ( tempU(3:K+1)-2*tempU( 2:K) ...
              +tempU(1:K-1) )/(dx*dx); 
    
    % uncontrolled Dirchlet boundary condition
    compute_u(1)=0; 
    compute_u(K+1)=0;  
    % controlled Neumann boundary condition
    compute_ctrlU(1)=0; 
    compute_ctrlU(K+1) = 2*ctrlTempU_1( K+1)-ctrlTempU( K+1)...
             +dt*dt* ( (ctrlTempU(K+1) - ...
             gain*dx/dt*(ctrlTempU_1(K+1)-ctrlTempU(K+1))) ...
          -2*ctrlTempU( K+1) +ctrlTempU(K) )/(dx*dx);
    compute_ctrlU1(1)=0; 
    compute_ctrlU1(K+1) =2*ctrlTempU1_1( K+1)-ctrlTempU1( K+1)...
             +dt*dt* ( (ctrlTempU1(K+1) - ...
             gain1*dx/dt*(ctrlTempU1_1(K+1)-ctrlTempU1(K+1))) ...
          -2*ctrlTempU1( K+1) +ctrlTempU1(K) )/(dx*dx);
    compute_ctrlU2(1)=0; 
    compute_ctrlU2(K+1) =2*ctrlTempU2_1( K+1)-ctrlTempU2( K+1)...
             +dt*dt* ( (ctrlTempU2(K+1) - ...
             gain2*dx/dt*(ctrlTempU2_1(K+1)-ctrlTempU2(K+1))) ...
          -2*ctrlTempU2( K+1) +ctrlTempU2(K) )/(dx*dx);
    
    % plot the solution
    figplot = 1;
    if mod (i-1, movieStep) == 0
        i
        if figplot == 1
           figure(1);
           plot(x, tempU, 'k -', x, ctrlTempU, 'r --', ...
               'linewidth', 3);
           legend('Uncontrolled', 'Controlled');
           axis([0 L -2.5 2.5]);
           xlabel('x','fontsize', 20);
           ylabel('u','fontsize', 20);
           s1 = ['c=1, gain = ' num2str(gain) ', t=' ...
                   num2str((i-1)*dt)];
           title(texlabel(s1), 'fontsize', 16);
           grid on;
           set(gca,'fontsize', 16);
           pause(0.1);
           F(fc)=getframe(gcf); % take frames to make a movie
       end
           time(fc)=(i-1)*dt;
           % compute the controller
           ctrl(fc) =  - gain*(ctrlTempU_1(end) -ctrlTempU(end) )/dt; 
           % save solutions
           u(:, fc) =tempU';
           ctrlU(:, fc) = ctrlTempU'; 
           ctrlU1(:, fc) = ctrlTempU1';
           ctrlU2(:, fc) = ctrlTempU2';
           fc = fc+1;
    end  
    time_for_energy(i)=(i-1)*dt;
    % compute the energy
    unctrl_E(i) = dx*trapz(( ( tempU_1 -tempU) /dt).^2) /2 ...
         +c^2*dx*trapz(( ( tempU(2:end) -tempU(1:end-1) )/dx).^2 )/2;
    ctrl_E(i) = dx*trapz(( ( ctrlTempU_1 -ctrlTempU) /dt).^2) /2 ...
         +c^2*dx*trapz(( ( ctrlTempU(2:end) -ctrlTempU(1:end-1) )/dx).^2 )/2;
    ctrl_E1(i) = dx*trapz(( ( ctrlTempU1_1 -ctrlTempU1) /dt).^2) /2 ...
         +c^2*dx*trapz(( ( ctrlTempU1(2:end) -ctrlTempU1(1:end-1) )/dx).^2 )/2;
    ctrl_E2(i) = dx*trapz(( ( ctrlTempU2_1 -ctrlTempU2) /dt).^2) /2 ...
         +c^2*dx*trapz(( ( ctrlTempU2(2:end) -ctrlTempU2(1:end-1) )/dx).^2 )/2;
    % update solutions
    tempU=tempU_1;
    tempU_1=compute_u;
    ctrlTempU = ctrlTempU_1 ;
    ctrlTempU_1 = compute_ctrlU;
    ctrlTempU1 = ctrlTempU1_1 ;
    ctrlTempU1_1 = compute_ctrlU1;
    ctrlTempU2 = ctrlTempU2_1 ;
    ctrlTempU2_1 = compute_ctrlU2;
end

plotstep = 10;
% plot control
figure(1);
plot( time(1:plotstep:end),   ctrl(1:plotstep:end));
xlabel('t','fontsize', 20);
ylabel('x','fontsize', 20);
zlabel('Control','fontsize', 20);
set(gca,'fontsize', 20);
grid on

% plot the controlled solution 
figure(2);
mesh( time(1:plotstep:end),x, ctrlU(:,1:plotstep:end) );
xlabel('t','fontsize', 20);
ylabel('x','fontsize', 20);
zlabel('u','fontsize', 20);
set(gca,'fontsize', 20);
grid on

% plot the uncontrolled solution
figure(3);
mesh( time(1:plotstep:end), x, u(:,1:plotstep:end) );
xlabel('t','fontsize', 20);
ylabel('x','fontsize', 20);
zlabel('u','fontsize', 20);
set(gca,'fontsize', 20);
grid on

% plot energy
plotstep1=100;
figure(4);
plot(time_for_energy(1:plotstep1:end), unctrl_E(1:plotstep1:end), 'k--', ...
    time_for_energy(1:plotstep1:end), ctrl_E(1:plotstep1:end), 'r-', ...
    time_for_energy(1:plotstep1:end), ctrl_E1(1:plotstep1:end), 'g-.', ...
    time_for_energy(1:plotstep1:end), ctrl_E2(1:plotstep1:end), ... 
    'b:', 'linewidth', 3);
xlabel('t','fontsize', 20);
ylabel('E','fontsize', 20);
set(gca,'fontsize', 20);
str1=['Gain = 0'];
str2=['Gain = ' num2str(gain)];
str3=['Gain = ' num2str(gain1)];
str4=['Gain = ' num2str(gain2)];
legend(str1, str2, str3, str4);
grid on

% make an animated movie
moviename= ['wave-boundary-control-c-1-gain-'  num2str(gain) '.avi'];
movie2avi(F, moviename, 'fps', 12, 'colormap',  jet, 'quality', 100)