%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program solves Problem 7 in the homework assignment,
% Finite Differences and Finite Elements:
% Getting to Know You
%
% The problem is to solve the differential equation
%   -(a(x) u'(x))' + c(x) u(x) = f(x)  for x in (0,1)
%      u(0)=u(1)=0
% using various methods and various choices of a, c, and f,
% and compare the results.
%
% The function a should satisfy a(x) \ge a_0 > 0 on (0,1),
% and the function c should satisfy c(x) \ge 0 on (0,1).

% The methods are:
%
%  -- finite differences, with a first-order approximation to
%     u'  (finitediff1)
%  -- finite differences, with a second-order approximation to
%     u'  (finitediff2)
%  -- finite elements to generate a piecewise-linear 
%     approximation to u  (fe_linear)
%  -- finite elements to generate a piecewise-quadratic
%     approximation to u  (fe_quadratic)
%
% problem7.m  Dianne P. O'Leary       January 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Three global variables are used to choose among the
%  various definitions of a, c, and true solution u.
%  The choices we will use for these variables are stored in cases.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global acase ccase ucase

cases = [1 1 1
         1 2 1
         1 3 1
         2 1 1
         3 1 1
         1 1 2
         1 1 3 ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Three values of M, the number of mesh points, will be used:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mset = [11,101,1001];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We loop over the 7 test problems.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for problemno = 1:7,

  acase = cases(problemno,1);
  ccase = cases(problemno,2);
  ucase = cases(problemno,3);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use finite differences; one-sided difference for 1st derivative term.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  j = 0;
  for M = mset         
    [Afd1,gfd1,xmeshfd1,ufd1] = finitediff1(M,@a,@c,@f);
    utrue = trueu(xmeshfd1);
    j = j + 1;
    uerrorfd1(j) = norm(ufd1 - utrue,inf);
  end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use finite differences; central difference for 1st derivative term.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  j = 0;
  for M = mset         
    [Afd2,gfd2,xmeshfd2,ufd2] = finitediff2(M,@a,@c,@f);
    utrue = trueu(xmeshfd2);
    j = j + 1;
    uerrorfd2(j) = norm(ufd2 - utrue,inf);
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use piecewise linear finite elements.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  j = 0;
  for M=mset         
    [Afe1,gfe1,xmeshfe1,ufe1] = fe_linear(M,@a,@c,@f);
    utrue = trueu(xmeshfe1);
    j = j + 1;
    uerrorfe1(j) = norm(ufe1 - utrue,inf);
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use piecewise quadratic finite elements.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  j = 0;
  for M=mset         
    [Afe2,gfe2,xmeshfe2,ufeval2,ufe2] = fe_quadratic(M,@a,@c,@f);
    utrue = trueu(xmeshfe2);
    j = j + 1;
    uerrorfe2(j) = norm(ufe2 - utrue,inf);
  end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display the results for this problem
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  disp(sprintf(' '))
  disp(sprintf('                        PROBLEM %d',problemno))
  disp(sprintf('Using coefficient functions a(%d) and c(%d) with true solution u(%d)', ...
        acase,ccase,ucase)) 
  disp(sprintf('         Infinity norm of the error at the mesh points'))
  disp(sprintf('   for various methods and numbers of interior mesh points M'))
  disp(strcat(sprintf('                    M ='),sprintf('   %10d',mset-2)))
  disp(strcat(sprintf('1st order finite difference'), ...
              sprintf('   %10.4e',uerrorfd1)))
  disp(strcat(sprintf('2nd order finite difference'), ...
              sprintf('   %10.4e',uerrorfd2)))
  disp(strcat(sprintf('     Linear finite elements'), ...
              sprintf('   %10.4e',uerrorfe1)))
  disp(strcat(sprintf('  Quadratic finite elements'), ...
              sprintf('   %10.4e',uerrorfe2)))
  
end  % for problemno

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot results of Problem 7.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot(xmeshfd1,ufd1,'r-',xmeshfd1,utrue,'b','lineWidth',1.5)
xlabel('t')
ylabel('u')
legend('Computed solution','True solution')
print -depsc figtest.eps