%*****************************************************************************
%          4-e ORDE METHODE
%
% NAAM:      vierde.m
% DATUM:     mei 1995 / november 1997
%*****************************************************************************
% OMSCHRIJVING:
%   Deze module bepaalt een discrete oplossing van
%         dy/dx - k d2y/dx2 = RL
%   op het interval 0<x<1, met randvoorwaarden y(0)=0 en y(1)=1.
%   Er wordt gediscretiseerd op een gerekt rooster.
%   RL=0 of 1 [zodat y=0.5*{x+(1-exp(x/k))/(1-exp(1/k))} oplossing is].
%*****************************************************************************
% INVOER:   N         aantal roosterpunten
%           k         diffusiecoefficient
%
% UITVOER:  xrek(1:N)    roosterpunten
%           yexrek(1:N)  exacte oplossing (RL=0)
%           yrek(1:N)    discrete oplossing
%*****************************************************************************
clear              % remove all variables

% for communication with energy.m function
global h yexrek N

%
N0=input('voer in grofste N: ');
hoe=input('wilt u N met factor of met increment verhogen? (0/1):  ');
if hoe==0,
  increm=input('N met welke factor verhogen? :   ');
else
  increm=input('N met welk increment verhogen? :  ');
end
half=input('hoe vaak verhogen? :  ');  

k=input('voer in k: ');
rechts=input('wilt u een rechterlid? (0/1):  ');
u=1;
iew=input('wilt u eigenwaarden? (0/1):  ');
%iew=0;
%
% berekening van gerekte rooster; N/2 mazen links en rechts van x = 1 - Lk
% goede keuze L=5
L=5;
%
% generaliseer tot F*N mazen links en (1-F)*N mazen rechts
%
%expro=input('rooster abrupt / exponentieel? (0/1):  ');
%F=input('welk deel van punten links? :  ');
F=0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% begin loop met roosterverfijning

for teller=1:half+1,

%for expro=0:1,
  
for L=5:5:15,

if (L==5 | L==15),
    expro=0;
else
    expro=1;
end 

clear x xrek yrekex b hcrs hfin Nc Nf fac
clear Coef2L rl2L y2L Coef2A rl2A y2A upp low dia 
clear Lagrange rl4L y4LE lagm2 lagm1 lag0 lagp1 lagp2 
clear M4A rl4A y4A2 m4Am2 m4Am1 m4A0 m4Ap1 m4Ap2 

if hoe==0,
  N = round(N0 * increm^(teller-1));
else
  N = N0 + round((teller-1)*increm);
end
NN(teller)=N;

Nc=round(N*F);
Nf=N-Nc;
if Nf <= 3,
  disp(['maar ',num2str(Nf),' mazen in grenslaag!']);
end  
  
hcrs=(1.-L*k)/Nc;           % grove mazen links
hfin=L*k/Nf;                 % fijne mazen rechts
xrek(1:Nc)=hcrs*(1:1:Nc);                 % roosterpunten links
xrek(Nc+1:N)=xrek(Nc) + hfin*(1:1:Nf);   % roosterpunten rechts
if L*k > 0.5
   hfin=1/N;
   hcrs=1/N;
   xrek(1:N)=(1:1:N)/N;
end
for i=2:N
   h(i)=xrek(i)-xrek(i-1);
end
%     
if (expro == 1 & L*k <= 0.5),
% rekking zodanig dat zelfde verdeling als boven 
   fac=(L*k/(1-L*k))^(1/Nc);
   disp(['rekkingsfactor= ',num2str(fac)])
   mesh=(1-fac)/(1-fac^N);
   xrek(1)=mesh;
   for i=2:N
     mesh=fac*mesh;
     xrek(i)=xrek(i-1)+mesh;
     h(i)=mesh;
   end
%   xrek=atan((1:1:N-1)/(sqrt(k)*N)) / atan(1/sqrt(k));
%   xrek(N)=1;
end
%
h(1)=xrek(1);
for i=2:N,
  h(i)=xrek(i)-xrek(i-1);
end

% 
% exacte oplossing op gerekt rooster 
%
pe=1/k;
% aktie om underflow van exp(-pe) te voorkomen
if pe > 100,
   empe=0;
else
   empe=exp(-pe);
end
yexrek=1-(1-exp(pe*(xrek-1)))/(1-empe);
%
b = zeros(1,N-1);

if rechts ==1,
%% inhomogeen rechterlid
twopi=2*3.14159265359;
for i=1:N-1
    if xrek(i) <= 0.4 
      b(i)=2.-5.*xrek(i);
    else
      b(i)=0;
    end
    yexrek(i)=ns_exact(xrek(i),k);
%  if i==1,
%     xmid=0.5*xrek(2);
%  else
%     xmid=0.5*(xrek(i-1)+xrek(i+1));
%  end
%  xmid=xrek(i);
%  b(i)=u*twopi*cos(twopi*xmid) + k*(twopi^2)*sin(twopi*xmid);
%  yexrek(i)=yexrek(i)+sin(twopi*xrek(i));
%   b(i)=u*0.5;
%   yexrek(i)=0.5*yexrek(i)+0.5*xrek(i);
end  
  yexrek(N)=1;
end % inhomogeen 

 
%%%%%%%%%%%%%%  2e orde methoden  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% Tweede orde Lagrange discretisatie
%
for i=1:N-1
   upp(i) = (h(i)* u - 2*k) / h(i+1);
   low(i) = (-h(i+1) * u - 2*k) / h(i);
   dia(i) = -upp(i)-low(i);
   rl2L(i) = b(i) * (h(i)+h(i+1));   
end
%
Coef2L=diag(upp(1:N-2),1) + diag(dia(1:N-1),0) + diag(low(2:N-1),-1);
%
% randvoorwaarde invullen
%
rl2L(N-1)=rl2L(N-1)-upp(N-1);
rl2L=rl2L';
%
y2L = Coef2L\rl2L;
y2L(N) = 1;
if iew==1,
   [V2L,D2L]=eig(Coef2L);
end

subplot(221)
plot(xrek,y2L,xrek,yexrek)
set(gca,'XLim',[0 1])
l2rek=norm(yexrek'-y2L);
E2L=energy(y2L);
trek=sprintf('error=%6.1e',E2L);
text(0.1,0.7,trek)
title('2nd Lagrange');
     

%
% tweede orde (anti)symmetrische Coef2A met rechterlid rl2 vullen **********
%
for i=1:N-1
   upp(i) = u - 2*k/h(i+1);
   low(i) = -u - 2*k/h(i);
   dia(i) = -upp(i)-low(i);
   rl2A(i) = b(i) * (h(i)+h(i+1));
end
%
Coef2A=diag(upp(1:N-2),1) + diag(dia(1:N-1),0) + diag(low(2:N-1),-1);
% 
% randvoorwaarde invullen
%
rl2A(N-1)=rl2A(N-1)-upp(N-1);
rl2A=rl2A';
%
y2A = Coef2A\rl2A;
y2A(N) = 1;
if iew==1,
%   [V2A,D2A]=eig(Coef2A);
end
subplot(222)
plot(xrek,y2A,xrek,yexrek)
set(gca,'XLim',[0 1])
l2rek=norm(yexrek'-y2A);
E2A=energy(y2A);
trek=sprintf('error=%6.1e',E2A);
text(0.1,0.7,trek)
title('2nd order');


%%%%%%%%%%%   vierde orde methoden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x=xrek;
x(N+1)=2*x(N)-x(N-1);

for i=1:N-1

   if i==1,
      xm2=-2*x(1);
      xm1=-x(1);
   elseif i==2,
      xm2=-x(2);
      xm1=x(1)-x(2);
   else
      xm2=x(i-2)-x(i);
      xm1=x(i-1)-x(i);
   end
   xp2=x(i+2)-x(i);
   xp1=x(i+1)-x(i);

%  traditional Lagrange 4-th order

   noemer = (xm2-xm1)*xm2*(xm2-xp1)*(xm2-xp2);
   cm2 = -xm1*xp1*xp2 / noemer;
   dm2 = 2*(xp1*xp2+xm1*xp1+xm1*xp2) / noemer;
   lagm2(i) = u*cm2 - k*dm2;

   noemer = (xm1-xm2)*xm1*(xm1-xp1)*(xm1-xp2);
   cm1 = -xm2*xp1*xp2 / noemer;
   dm1 = 2*(xp1*xp2+xm2*xp1+xm2*xp2) / noemer;
   lagm1(i) = u*cm1 - k*dm1;

   noemer = (xp1-xm2)*(xp1-xm1)*xp1*(xp1-xp2);
   cp1 = -xm2*xm1*xp2 / noemer;
   dp1 = 2*(xm1*xp2+xm2*xm1+xm2*xp2) / noemer;
   lagp1(i) = u*cp1 - k*dp1;

   noemer = (xp2-xm2)*(xp2-xm1)*xp2*(xp2-xp1);
   cp2 = -xm2*xm1*xp1 / noemer;
   dp2 = 2*(xm1*xp1+xm2*xp1+xm2*xm1) / noemer;
   lagp2(i) = u*cp2 - k*dp2;

   lag0(i) = -lagm2(i) -lagm1(i) -lagp1(i) -lagp2(i);

%  modern (anti)symmetric 4th-order

   noemer4 = -xp2 + 8*xp1 -8*xm1 +xm2;
%   if i==1
%      bb(i)=(2+4*b(i)+b(i+1))/6;
%   elseif i < N-1
%      bb(i)=(b(i-1)+4*b(i)+b(i+1))/6;
%   else
%      bb(i)=0;
%   end
   b4A(i)=noemer4*b(i);
%  schaling irrelevant (aangezien noemer4 < 0 kan worden weglaten)
   noemer4=1.0;
   m4Am2(i)=(u-2*k/xm2)/noemer4;
   m4Am1(i)=-8*(u-2*k/xm1)/noemer4;
   m4Ap1(i)=8*(u-2*k/xp1)/noemer4;
   m4Ap2(i)=-(u-2*k/xp2)/noemer4;
   m4A0(i)=  -m4Am2(i) -m4Am1(i) -m4Ap1(i) -m4Ap2(i);

end

% traditionele 4e orde Lagrange discretisatie *******************************

% exacte randvoorwaarden 

ym1=1-(1-exp(pe*(-x(1)-1)))/(1-empe);
yNp1=1-(1-exp(pe*(x(N+1)-1)))/(1-empe);

if rechts == 1,
   ym1=ns_exact(-x(1),k);
   yNp1=ns_exact(x(N+1),k);
%  ym1=ym1+sin(-twopi*x(1));
%  yNp1=yNp1+sin(twopi*x(N+1));
%   ym1=0.5*ym1-0.5*x(1);
%   yNp1=0.5*yNp1+0.5*x(N+1);
end

rl4L=b;
rl4L(1)=rl4L(1)-lagm2(1)*ym1;
rl4L(N-2)=rl4L(N-2)-lagp2(N-2);
rl4L(N-1)=rl4L(N-1)-lagp2(N-1)*yNp1-lagp1(N-1);

Lagrange=diag(lagp2(1:N-3),2)+diag(lagp1(1:N-2),1)+diag(lag0(1:N-1),0)+ ...
             diag(lagm1(2:N-1),-1)+diag(lagm2(3:N-1),-2);

y4LE=Lagrange\rl4L';
y4LE(N)=1.0;
if iew==1,
   [V4LE,D4LE]=eig(Lagrange);
end

subplot(223)
plot(xrek,y4LE,xrek,yexrek)
set(gca,'XLim',[0 1])
l2rek=norm(yexrek'-y4LE);
E4LE=energy(y4LE);
trek=sprintf('error=%6.1e',E4LE);
text(0.1,0.7,trek)
title('4th Lag + exact BC');
	 
% moderne (anti)symmetrische methode ************************************

M4A=diag(m4Ap2(1:N-3),2)+diag(m4Ap1(1:N-2),1)+diag(m4A0(1:N-1),0)+ ...
         diag(m4Am1(2:N-1),-1)+diag(m4Am2(3:N-1),-2);

% laatste punten tweede orde 
     
M4A(1,1)=dia(1);
M4A(1,2)=upp(1);
M4A(1,3)=0;
M4A(N-1,N-3)=0;
M4A(N-1,N-2)=low(N-1);
M4A(N-1,N-1)=dia(N-1);

rl4A=b4A;
rl4A(N-2)=rl4A(N-2)-m4Ap2(N-2);
rl4A(1)=rl2A(1);
rl4A(N-1)=rl2A(N-1);

y4A2=M4A\rl4A';
y4A2(N)=1.0;

if iew==1,
%   [V4A2,D4A2]=eig(M4A);
end

subplot(224)
plot(xrek,y4A2,xrek,yexrek)
set(gca,'XLim',[0 1])
l2rek=norm(yexrek'-y4A2);
E4A2=energy(y4A2);
trek=sprintf('error=%6.1e',E4A2);
text(0.1,0.7,trek)
title('4th + 2nd order')


% overzicht ********************* 
 
ii=round(N/2);
for i=1:N-1
   if (x(i) <= 1-k & x(i+1) > 1-k),
       ii=i;
   end
end

ii=round(3*N/4);
ik=ii;

% store local errors

%local(teller,1+expro)=(y2L(ii)-yexrek(ii));
%local(teller,3+expro)=(y2A(ii)-yexrek(ii));
%local(teller,5+expro)=(y4LE(ii)-yexrek(ii));
%local(teller,7+expro)=(y4A2(ii)-yexrek(ii));

% eigenvalues

if iew==1,
  t2L=0;
  t2A=0;
  t4LE=0;
  t4A2=0;

  for i=1:N-1,
    d2L(i)=real(D2L(i,i));
%    d2A(i)=real(D2A(i,i));
    d4LE(i)=real(D4LE(i,i))/N;
%    d4A2(i)=real(D4A2(i,i));
    if d2L(i) < 0.0,
       t2L=t2L+1;
    end
%    if d2A(i) < 0.0,
%       t2A=t2A+1;
%    end
    if d4LE(i) < 0.0,
       t4LE=t4LE+1;
    end
%    if d4A2(i) < 0.0,
%       t4A2=t4A2+1;
%    end
  end
end


% store global errors and eigenvalues

if (L==5),

%  globaal5(teller,1+expro)=max(abs(y2L-yexrek'));
%  globaal5(teller,3+expro)=max(abs(y2A-yexrek'));
%  globaal5(teller,5+expro)=max(abs(y4LE-yexrek'));
%  globaal5(teller,7+expro)=max(abs(y4A2-yexrek'));

  globaal5(teller,1+expro)=energy(y2L);
  globaal5(teller,3+expro)=energy(y2A);
  globaal5(teller,5+expro)=energy(y4LE);
  globaal5(teller,7+expro)=energy(y4A2);

  if iew==1
    remin5(teller,1)=t2L;
    remin5(teller,2)=t2A;
    remin5(teller,3)=t4LE;
    remin5(teller,4)=t4A2;
  end

elseif L==10

%  globaal10(teller,1+expro)=max(abs(y2L-yexrek'));
%  globaal10(teller,3+expro)=max(abs(y2A-yexrek'));
%  globaal10(teller,5+expro)=max(abs(y4LE-yexrek'));
%  globaal10(teller,7+expro)=max(abs(y4A2-yexrek'));

  globaal10(teller,1+expro)=energy(y2L);
  globaal10(teller,3+expro)=energy(y2A);
  globaal10(teller,5+expro)=energy(y4LE);
  globaal10(teller,7+expro)=energy(y4A2);

  if iew==1,
    remin10(teller,1)=t2L;
    remin10(teller,2)=t2A;
    remin10(teller,3)=t4LE;
    remin10(teller,4)=t4A2;
  end

else

%  globaal15(teller,1+expro)=max(abs(y2L-yexrek'));
%  globaal15(teller,3+expro)=max(abs(y2A-yexrek'));
%  globaal15(teller,5+expro)=max(abs(y4LE-yexrek'));
%  globaal15(teller,7+expro)=max(abs(y4A2-yexrek'));

  globaal15(teller,1+expro)=energy(y2L);
  globaal15(teller,3+expro)=energy(y2A);
  globaal15(teller,5+expro)=energy(y4LE);
  globaal15(teller,7+expro)=energy(y4A2);

  if iew==1,
    remin15(teller,1)=t2L;
    remin15(teller,2)=t2A;
    remin15(teller,3)=t4LE;
    remin15(teller,4)=t4A2;
  end

end

% display errors

%disp(['  ']);
disp([' N','    E2L    ','   E2A    ','   E4LE   ','    E4A2']);

disp([num2str(N),'  ',num2str(E2L),'  ',num2str(E2A),'  ', ...
    num2str(E4LE),'  ',num2str(E4A2)])

disp([ num2str(yexrek(ik)),'  ', ...
       num2str(y2L(ik)-yexrek(ik)),'  ',num2str(y2A(ik)-yexrek(ik)),'  ', ...
       num2str(y4LE(ik)-yexrek(ik)),'  ',num2str(y4A2(ik)-yexrek(ik))]);
disp(['  ']);

%pause 
clf

yexrek=yexrek';
%
%semilogy(xrek,abs(y2L-yexrek),'r--',xrek,abs(y2A-yexrek),'g--', ...
%         xrek,abs(y4LE-yexrek),'r-', xrek,abs(y4A2-yexrek),'g-')
%set(gca,'XLim',[0 1],'YLim',[1e-6 1]);
%legend('2L','2A','4LE','4A2',2);
%pause 
%clf


% einde loop over middenpunt
end

% einde loop over soort rooster
%end

% einde loop met roosterverfijning
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% overzicht maaswijdte verfijning 

 subplot(231)
a = loglog(1./NN,abs(globaal5(:,1)),'r--', ...
        1./NN,abs(globaal5(:,3)),'g--', ...
        1./NN,abs(globaal5(:,5)),'r-', ...
	1./NN,abs(globaal5(:,7)),'g-'  );

  axis([3e-3,0.1,1e-7,1000]);
  set(gca,'Xtick',[0.01 0.1]);
%  legend('2L','2A','4LE','4A2',4)
  title('Abrupt: d=5k','Fontsize',12)
  xlabel('mean mesh size','Fontsize',12)
  ylabel('global error','Fontsize',12)
  set(a(2),'LineWidth',2);
  set(a(4),'Linewidth',2);


 subplot(232)
b = loglog(1./NN,abs(globaal15(:,1)),'r--', ... 
        1./NN,abs(globaal15(:,3)),'g--', ...
        1./NN,abs(globaal15(:,5)),'r-', ...
        1./NN,abs(globaal15(:,7)),'g-' );
    
  axis([3e-3,0.1,1e-7,1000]);
  set(gca,'Xtick',[0.01 0.1]);
  legend('2L','2A','4LE','4A2',2)
  title('Abrupt: d=15k','FontSize',12);    
  xlabel('mean mesh size','Fontsize',12)
%  ylabel('global error','Fontsize',12)
  set(b(2),'LineWidth',2);
  set(b(4),'Linewidth',2);


 subplot(233)
c = loglog(1./NN,abs(globaal10(:,2)),'r--', ...
         1./NN,abs(globaal10(:,4)),'g--', ...
	 1./NN,abs(globaal10(:,6)),'r-', ... 
         1./NN,abs(globaal10(:,8)),'g-'  );
			  
  axis([3e-3,0.1,1e-7,1000]);
  set(gca,'Xtick',[0.01 0.1]);
  %  legend('2L','2A','4LE','4A2',4)
  title('Exponential: d=10k','Fontsize',12)
  xlabel('mean mesh size','Fontsize',12)
%  ylabel('global error','Fontsize',12)
  set(c(2),'LineWidth',2);
  set(c(4),'Linewidth',2);
			      
			      
% subplot(224)
%d = loglog(1./NN,abs(globaal15(:,2)),'r--', ...
%        1./NN,abs(globaal15(:,4)),'g--', ... 
%        1./NN,abs(globaal15(:,6)),'r-', ...
%        1./NN,abs(globaal15(:,8)),'g-' );
			      
%%  axis([0.4e-3,.3,1e-12,10]);
%   legend('2L','2A','4LE','4A2',4)
%   title('Exponential: L=15k','FontSize',12);
%   xlabel('mean mesh size','Fontsize',12)
%   ylabel('error','Fontsize',12)
%   set(d(2),'LineWidth',2);
%   set(d(4),'Linewidth',2);

print -deps jem-1.eps

pause
clf


if iew==1,
   subplot(331)				      
   semilogx(1./NN,remin5(:,1),'r--', ...
            1./NN,remin5(:,3),'r-');
   axis([3e-3 0.1 0 15]);
   set(gca,'Xtick',[0.01 0.1]);
   title('Abrupt: d=5k','Fontsize',12)
   xlabel('mean mesh size', 'Fontsize',12);
   ylabel('# -ve eigenvalues','Fontsize',12);
   
   subplot(332)
   semilogx(1./NN,remin15(:,1),'r--', ...
         1./NN,remin15(:,3),'r-' );
   axis([3e-3 0.1 0 15]);
   set(gca,'Xtick',[0.01 0.1]);
   title('Abrupt: d=15k','FontSize',12);    
   xlabel('mean mesh size', 'Fontsize',12);
%   ylabel('# negative eigenvalues','Fontsize',12);
   legend('2L','4LE',2)
   
   subplot(333)
   semilogx(1./NN,remin10(:,1),'r--', ...
         1./NN,remin10(:,3),'r-' );
   axis([3e-3 0.1 0 15]);
   set(gca,'Xtick',[0.01 0.1]);
   title('Exponential: d=10k','Fontsize',12)
   xlabel('mean mesh size', 'Fontsize',12);
%   ylabel('number of negative eigenvalues','Fontsize',12);
end

print -deps jem-2.eps