function[]=IIM_1D_Example1()
%Kevin Luna
%Group 3
%NCSU REU 2015
%Research Mentor: Dr.Li
%IIM_1D_Example.m is a matlab function that applies the IIM developed by Dr.Li with finite differences in 1D
% and ouputs plots of an exact solution, the approximate solution, a
% difference error plot, and a loglog plot of the error vs h
%The code as fits a regression line onto the solution error, and truncation
%errors at the irregular gridpoints j and j+1 of the IIM and prints out the
%slope of the loglog data. 


%The code applies the IIM to the problem:
%(\beta u_x)_x -\sigma u= f(x)+ c\delta(x-\alpha)
%u(a)=ua, u(b)=b , on the domain (a,b)



%Note that the  domain is split up into n subdivisons in the code

%A self adjoint ellptic equation discretization has been used at all the
%non-irregular grid points.

% Note that f,\sigma,u,u_x may all have finite jumps over [a,b] at \alpha

% The code computes the jumps at alpha for all functions.

%To use the function edit the IIM_1D_Example portion of the function to use
%new parameters and adjust the apprpriate subfunctions consistently in the section
%marked: INTERFACE BVP FUNCTIONS 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Interface BVP input information
global alpha  k1 k2 c bf W SJ % commonly used constants in code
close all 
a=0;        %left boundary
b=1;        %right boundary
ua=0;        %left boundary value 
n=20;        %number of subdivisons of (a,b)
alpha=a+(1/3); % location of point source
k1=120;         %constant in solution
k2=30;      % constant in solution
ub=cos(k2);      %right boundary value


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Grid Refinement analysis
maxind=15;       % Number of refinements
hvec=zeros(maxind,1);       % Preallocation vector for h values
e1inf=zeros(maxind,1);      % Preallocation vector for error  values
T1=zeros(maxind,1);     % Preallocation vector for truncation error values
T2=zeros(maxind,1);     % Preallocation vector for truncation values

A=uDerexactminus(alpha);%exact left deriavative
B=uDerexactplus(alpha);%exact right  deriavative

%preallocation vectors for derivatives at alpha and errors there
UalphaprimeR=zeros(maxind,1);
UalphaprimeRE=zeros(maxind,1);
UalphaprimeL=zeros(maxind,1);
UalphaprimeLE=zeros(maxind,1);
Ualphaprime3=zeros(maxind,1);
Ualphaprime3E=zeros(maxind,1);
DUalphaprime3=zeros(maxind,1);
DUalphaprime3E=zeros(maxind,1);

% grid refinement loop
for i=1:maxind
    [X,U,F,P,Q,Xhalf,gamma1,gamma2,Cj1,Cj2,j]=One_D_IIM_Solver(a,b,ua,ub,n); % calling IIM FD solver 
    u=zeros(n-1,1);% Preallocation vector for exact solution values
    for k=1:n-1
             u(k)=uexact(X(k));% filling u with exact solution values
    end
    
    
hvec(i)=(b-a)/n ;% computing h for current number of grid points
e1inf(i)=norm(u-U,inf); %computing error for current number of grid points
T1(i)=gamma1(1)*u(j-1)+gamma1(2)*u(j)+gamma1(3)*u(j+1)-Q(j)*u(j)-F(j); % Truncation error at X(j)  for current number of grid divisions
T2(i)=gamma2(1)*u(j)+gamma2(2)*u(j+1)+gamma2(3)*u(j+2)-Q(j+1)*u(j+1)-F(j+1);% Truncation error at X(j+1) for current number of grid divisions


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Using IIM to compute derivative finite difference coeffcients and
%%%%correction terms
B11d=[1 ,1, 1+(.5*(1/Bplus(alpha))*(X(j+1)-alpha)^2)*SJ];
B21d=[X(j-1)-alpha , X(j)- alpha, (Bminus(alpha)/Bplus(alpha))*(X(j+1)-alpha)+(Bderminus(alpha)/Bplus(alpha) - (Bminus(alpha)*Bderplus(alpha))/Bplus(alpha)^2)*((X(j+1)-alpha)^2)/2];
B31d=[.5*(alpha-X(j-1))^2, .5*(alpha-X(j))^2 , .5*(Bminus(alpha)/Bplus(alpha))*(X(j+1)-alpha)^2];
Bj1d=[B11d;B21d;B31d];
gammaRHS1d=[0;1;0];
gamma1d=Bj1d\gammaRHS1d; 

B12d=[1-(.5*(1/Bminus(alpha))*(X(j)-alpha)^2)*SJ  ,1 ,1];
B22d=[((Bplus(alpha)/Bminus(alpha))*(X(j)-alpha) + (Bderplus(alpha)/Bminus(alpha) - (Bplus(alpha)*Bderminus(alpha))/Bminus(alpha)^2)*((X(j)-alpha)^2)/2 ) , X(j+1)- alpha, (X(j+2)-alpha)];
B32d=[ .5*(Bplus(alpha)/Bminus(alpha))*(X(j)-alpha)^2, .5*(X(j+1)-alpha)^2 , .5*(X(j+2)-alpha)^2 ];
Bj2d=[B12d;B22d;B32d];
gammaRHS2d=[0;1;0];
gamma2d=Bj2d\gammaRHS2d; 

Cj1d=gamma1d(3)*(W+(X(j+1)-alpha)*(c/Bplus(alpha))-.5*((X(j+1)-alpha)^2)*(Bderplus(alpha)*c/(Bplus(alpha)^2)-(W*sigmaplus(alpha)/Bplus(alpha))- bf/Bplus(alpha)));
Cj2d=gamma2d(1)*(-W-(X(j)-alpha)*(c/Bminus(alpha))-.5*((X(j)-alpha)^2)*(-Bderminus(alpha)*c/(Bminus(alpha)^2)+(W*sigmaminus(alpha)/Bminus(alpha))+bf/Bminus(alpha)));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Computing two-sided derivatives
%left
Ualphaprime3(i)=gamma1d(1)*U(j-1)+gamma1d(2)*U(j)+ gamma1d(3)*U(j+1)-Cj1d;
Ualphaprime3E(i)=abs(A-Ualphaprime3(i));

%right
DUalphaprime3(i)=gamma2d(1)*U(j)+gamma2d(2)*U(j+1)+ gamma2d(3)*U(j+2)-Cj2d;
DUalphaprime3E(i)=abs(B-DUalphaprime3(i));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Computing one-sided derivatives

UalphaprimeR(i)= -(3/(2*hvec(i)))*U(j+1) +(2/hvec(i))*U(j+2) -(1/(2*hvec(i)))*U(j+3);
UalphaprimeRE(i)=abs(B-UalphaprimeR(i));
UalphaprimeL(i)= (3/(2*hvec(i)))*U(j) -(2/hvec(i))*U(j-1) +(1/(2*hvec(i)))*U(j-2);
UalphaprimeLE(i)=abs(A-UalphaprimeL(i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






n=2*n; %doubling number of grid divisions

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DUalphaprime3E=Ualphaprime3E./k2;% relative error
UalphaprimeRE=UalphaprimeLE./k2;% relative error


%%%OUTPUTS
% loglog plot of solution error
figure(1)
plot(log10(hvec),log10(e1inf),'or')
legend('Immersed Interface Method','Location','southeast')
xlabel('log_{10}(h)')
ylabel('log_{10}( ||U-u||_\infty )')
title('Solution Error Plot on a log-log Scale for the IIM')
axis([min(log10(hvec)) max(log10(e1inf(2:end))) min(log10(hvec)) max(log10(e1inf(2:end)))])
hold on
p1=polyfit(log10(hvec),log10(e1inf),1); %fit line onto loglog data
xh0=linspace(hvec(1),hvec(maxind),1000);
x0=log10(xh0);
y0=polyval(p1,x0);
plot(x0,y0)
axis('square')
print('soln1','-depsc')
hold off

% loglog plot of the left  derivative  error
figure(2)
plot(log10(hvec),log10(DUalphaprime3E),'or',log10(hvec),log10(UalphaprimeRE),'sq')
title('loglog Plot of U_x^+(\alpha) Error')
xlabel('log_{10}(h)')
ylabel('log_{10}( ||U_x^+(\alpha)-u_x^+(\alpha)||_\infty )')
legend('Two-Sided Method(IIM)','One-sided method','Location','southeast')
hold on
p4=polyfit(log10(hvec),log10(DUalphaprime3E),1);
p6=polyfit(log10(hvec),log10(UalphaprimeRE),1);%fit line onto loglog data
xh1=linspace(hvec(1),hvec(maxind),1000);
x1=log10(xh1);
y1=polyval(p4,x1);
plot(x1,y1)
xh2=linspace(hvec(1),hvec(maxind),1000);
x2=log10(xh2);
y2=polyval(p6,x2);
plot(x2,y2)
M=log10([DUalphaprime3E(2:end);UalphaprimeRE(2:end)]);
ch=max((M));
axis([min(log10(hvec)) ch min(log10(hvec)) ch])
axis('square')
print('alpha1','-depsc')
hold off


%difference error plot
figure(3)
plot(X,U-u)
title('Error (difference) Plot')
xlabel('x' )
ylabel('U-u(x)')
print('error1','-depsc')



% Exact and approximate solutions on same plot
figure(4)
plot(X,u,X,U,'or')
legend('Exact','Approximation')
xlabel('x' )
ylabel('u(x)')
title(' (\beta u_{x}(x))_x-\sigma(x)u(x)=f(x) + c\delta(x-\alpha) Approximation and Exact Solution with DBC')


p1=polyfit(log10(hvec),log10(e1inf),1); %fit line onto loglog data
fprintf('The slope of the  solution error plot on the  infinity norm plot is %.8f \n',p1(1))

p2=polyfit(log10(hvec),log10(T1),1); %fit line onto loglog data
fprintf('The slope of the truncation  error at j for the  loglog infinity norm data is %.8f \n',p2(1))


p3=polyfit(log10(hvec),log10(T2),1); %fit line onto loglog data
fprintf('The slope of the truncation  error at j+1 for the loglog infinity norm data is %.8f \n',p3(1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p4=polyfit(log10(hvec),log10(Ualphaprime3E),1);
fprintf('The slope of the left two sided U''(alpha) error plot is %.8f \n',p4(1))


p5=polyfit(log10(hvec),log10(DUalphaprime3E),1);%fit line onto loglog data
fprintf('The slope of the right two sided U''(alpha) error plot is %.8f \n',p5(1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p6=polyfit(log10(hvec),log10(UalphaprimeLE),1);%fit line onto loglog data
fprintf('The slope of the left  one sided U''(alpha) error plot is %.8f \n',p6(1))

p7=polyfit(log10(hvec),log10(UalphaprimeRE),1);%fit line onto loglog data
fprintf('The slope of the right  one sided U''(alpha) error plot is %.8f \n',p7(1))
end














function[X,U,F,P,Q,Xhalf,gamma1,gamma2,Cj1,Cj2,j]=One_D_IIM_Solver(a,b,ua,ub,n)
global alpha c bf W SJ % setting global varibles
c=ccalc(alpha); %[\beta u_x]_{x=\alpha} =c jump
bf=bfcalc(alpha); %[f] jump
W=bucalc(alpha);%  [u]=W jump
SJ=sigmaJump(alpha);% [\sigma] jump

h=(b-a)/n; % step size
hs=h*h;   % squared step size


%%%%%Preallocating memory for function evaluations
F=zeros(n-1,1);      %f preallocation
P=zeros(n,1);       % \beta preallocation
Q=zeros(n-1,1);      %\sigma preallocation
X=zeros(n-1,1);         %X_i preallocation
Xhalf=zeros(n,1);       % X_1/2 preallocation


Xhalf(1)=a+ h/2 ; % first X_1/2 value

%%%%%%%%%% Evaluating \beta, \sigma, \f , and filling X and X_1/2
for i=1:n-1
    X(i)=a+i*h;            %X values starting at a
    Xhalf(i+1)=X(i)+ h/2 ; %X_1/2 values. Loop starts at second value. Total of n entries.(1 more than X)
    F(i)=feval(@f1, X(i)); %f(x) values
    Q(i)=feval(@q,X(i));   %\sigma(x) values
    P(i)=feval(@B,Xhalf(i)); %\beta(x_1/2) values
end

P(n)=feval(@B,Xhalf(n)); % The  \beta(x_1/2) values have one more entry after n-1 entries
d1=zeros(n-1,1); % Preallocating memory for main diagonal of finite difference A matrix
d2=P(2:n); %left diagonals of finite difference A matrix
d3=P(1:n-1); %right diagonals of finite difference A matrix

for i=1:n-1
   d1(i)= -(P(i)+P(i+1)) -Q(i)*hs ;     %Main diagonal of finite difference A matrix
end
A=(1/hs).*spdiags([d2 d1 d3],[-1 0 1], n-1,n-1);  % Sparse tridigonal finite difference matrix with above diagonals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Computation for index just before or at \alpha 
Y=zeros(n-1,1); %Preallocation 
for i=1:n-1
    Y(i)=irreg_j_index1(X(i)); %vector of x values less than or equal to alpha
end
[M,j] = max(Y); % MATLAB function retrives the index of the largest x valueless than or equal to alpha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%IIM coeffcient computation at j
B11=[1 ,1, 1+ ((X(j+1)-alpha)^2)*.5*(1/Bplus(alpha))*SJ];
B21=[X(j-1)-alpha , X(j)- alpha, (Bminus(alpha)/Bplus(alpha))*(X(j+1)-alpha)+(Bderminus(alpha)/Bplus(alpha) - (Bminus(alpha)*Bderplus(alpha))/Bplus(alpha)^2)*((X(j+1)-alpha)^2)/2];
B31=[.5*(alpha-X(j-1))^2, .5*(alpha-X(j))^2 , .5*(Bminus(alpha)/Bplus(alpha))*(X(j+1)-alpha)^2];
Bj1=[B11;B21;B31];
gammaRHS1=[0;Bderminus(alpha);Bminus(alpha)];
gamma1=Bj1\gammaRHS1; 

%%%IIM coeffcient computation at j+1
B12=[1-((X(j)-alpha)^2)*.5*(1/Bminus(alpha))*SJ ,1 ,1];
B22=[((Bplus(alpha)/Bminus(alpha))*(X(j)-alpha) + (Bderplus(alpha)/Bminus(alpha) - (Bplus(alpha)*Bderminus(alpha))/Bminus(alpha)^2)*((X(j)-alpha)^2)/2 ) , X(j+1)- alpha, (X(j+2)-alpha)];
B32=[ .5*(Bplus(alpha)/Bminus(alpha))*(X(j)-alpha)^2, .5*(X(j+1)-alpha)^2 , .5*(X(j+2)-alpha)^2 ];
Bj2=[B12;B22;B32];
gammaRHS2=[0;Bderplus(alpha);Bplus(alpha)];
gamma2=Bj2\gammaRHS2; 

%%% correction term computation
Cj1=gamma1(3)*(W+(X(j+1)-alpha)*(c/Bplus(alpha))-.5*((X(j+1)-alpha)^2)*(Bderplus(alpha)*c/(Bplus(alpha)^2)-(W*sigmaplus(alpha)/Bplus(alpha)) -bf/Bplus(alpha)));
Cj2=gamma2(1)*(-W-(X(j )-alpha)*(c/Bminus(alpha))-.5*((X(j)-alpha)^2)*(-Bderminus(alpha)*c/(Bminus(alpha)^2)+(W*sigmaminus(alpha)/Bminus(alpha))+bf/Bminus(alpha)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Adding correction term to F
F(j)=F(j)+Cj1;
F(j+1)=F(j+1)+Cj2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Replacing the j and j+1 rows of the A matrix with the IIM coefficients
A(j,j-1)=gamma1(1);
A(j,j)=gamma1(2)-Q(j);
A(j,j+1)=gamma1(3);

A(j+1,j)=gamma2(1);
A(j+1,j+1)=gamma2(2)-Q(j+1);
A(j+1,j+2)=gamma2(3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F(1)=F(1)-(ua*P(1))/hs ; %  F(X_1) term
F(n-1)=F(n-1)-(ub*P(n))/hs; %F(X_n-1) F term
U=A\F ; % Solving  system for approximate solution u using "\"

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

end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%INTERFACE BVP FUNCTIONS 
% Note that "-" means values BEFORE \alpha(interface point)
%AND that "+" means values AFTER \alpha(interface point)
function[y]=Bminus(x) 
%  \beta- function
y=1+(cos(x)/2);
end


function[y]=Bplus(x)
%  \beta+ function
y=(1+x^2);
end

function[y]=Bderminus(x)
%  \beta_x- function
y=-sin(x)/2;
end


function[y]=Bderplus(x)
%   \beta_x+ function
y=2*x;
end

function[y]=uexactminus(x)
%exact u- solution
global k1 k2
 y=sin(k1*x);
end

function[y]=uexactplus(x)
% exact u+ solution

global k1 k2
 y=cos(k2*x);
end



function[y]=uDerexactminus(x)
%exact u_x-
global k1 k2
 y=k1*cos(k1*x);
end


function[y]=uDerexactplus(x)
%exact u_x+
global k1 k2
 y=-k2*sin(k2*x);
end


function[y]=fminus(x)
%f- RHS of ODE
global k1 k2
y=(k1^2)*(-((cos(x)/2)+1))*sin(k1*x)-(1/2)*(k1*sin(x)*cos(k1*x))-sigmaminus(x)*uexactminus(x);
end

function[y]=fplus(x)
%f+ RHS of ODE
global k1 k2
 y=k2*(-k2*(x^2+1)*cos(k2*x)-2*x*sin(k2*x))-sigmaplus(x)*uexactplus(x);
end

function[y]=sigmaminus(x)
% \sigma-
y=(1+x)^2;
end


function[y]=sigmaplus(x)
% \sigma+
y=log(2+x);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%











%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[Betader]=Bx(x)
% \beta(x)_x function
global alpha
if x <= alpha 
    Betader=Bderminus(x);    
else      
   Betader=Bderplus(x);
end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[y]=B(x)
% \beta(x) function
global alpha
if x<= alpha
    y=Bminus(x);
    else      
    y=Bplus(x);
end
    
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[RHS]=f1(x)
% f(x)function
global alpha k1 k2
if x <= alpha
    RHS=fminus(x); 
    
else 
        
   RHS=fplus(x);
   
end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[s]=q(x)
%\sigma(x) function
global alpha
if x<= alpha
    s=sigmaminus(x);
    else      
    s=sigmaplus(x);
end
    
end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y] = irreg_j_index1( x )
%outputs values of x <= \alpha and sets values of x > \alpha to 0
global alpha
if x <= alpha
    y=x;  
else   
  y=0 ;  
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[s]=sigmaJump(x)
%[\sigma] computation 
s=sigmaplus(x)-sigmaminus(x);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[v]=ccalc(x) 
%[\beta u_x] computation
v=Bplus(x)*uDerexactplus(x)-Bminus(x)*uDerexactminus(x);  
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[y]=bfcalc(x)
%[f] computation
y= fplus(x) - fminus(x) ; 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[y]=bucalc(x)
%[u] computation
y=uexactplus(x)-uexactminus(x);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[y]=uexact(x)
%Exact solution function
global alpha
if x <= alpha
    y=uexactminus(x); 
else 
    y=uexactplus(x); 
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%