function fdm
clear all;clc;close all;
% Finite Difference Method
%            "Q.Fang, T.Tsuchiya, T.Yamamoto, Finite difference, finite 
%            element and finite volume methods applied to two-point 
%            boundary value problems, Journal of Computational and Applied 
%            Mathematics 139 (2002) 9-19
%
%            p(x)=exp(1-x), f(x)=1+exp(1-x), u(0)=u(1)=0;
%            Exact solution: u(x)=x(1-exp(x-1))
%            Max norm: FDM 8.24e-5,8.31e-7,8.31e-9
%                      FEM 6.35e-5,6.36e-7,6.39e-9
%                      FVM 3.18e-5,3.18e-7,3.18e-9
p=@(x) exp(1-x);
f=@(x) 1+exp(1-x);
%f=@(x) exp(1-x).*(x.^(1.5).*(1-x).^(1.5).*(60-20*x)-15*x.^(0.5).*(1-x).^(0.5).*(1-2*x+2*x.^2));
ue=@(x) x.*(1-exp(x-1));%exact solution
%ue=@(x) 4*x.^2.5.*(1-x).^2.5;
maxnorm=[];L2norm=[];
for k=1:2
    ttt=cputime;
    if k==1,kind=ones(1,10);nodes=[1:10;2:11];nelem=10;
    elseif k==2,kind=ones(1,100);nodes=[1:100;2:101];nelem=100;
    else kind=ones(1,1000);nodes=[1:1000;2:1001];nelem=1000;
    end
    h=10^-k;%grid width
%for k=1:7
%    if k==1,kind=ones(1,2);nodes=[1:2;2:3];nelem=2;
%    elseif k==2,kind=ones(1,4);nodes=[1:4;2:5];nelem=4;
%    elseif k==3,kind=ones(1,8);nodes=[1:8;2:9];nelem=8;
%    elseif k==4,kind=ones(1,16);nodes=[1:16;2:17];nelem=16;
%    elseif k==5,kind=ones(1,32);nodes=[1:32;2:33];nelem=32;
%    elseif k==6,kind=ones(1,64);nodes=[1:64;2:65];nelem=64;
%    else kind=ones(1,128);nodes=[1:128;2:129];nelem=128;
%    end
%    h=2^-k;%grid width
    xa=0:h:1;n=length(xa)-2;
    H=1/h*eye(n);
    xg=0.5*(xa(1:end-1)+xa(2:end));
    F=f(xa(2:end-1)).';
    pxg=p(xg);a=1/h*pxg;
    A=diag(a(1:n)+a(2:end));
    for i=1:n-1,A(i:i+1,i:i+1)=A(i:i+1,i:i+1)+a(i+1)*[0 -1;-1 0];end
    u=(H*A)\F;u=[0;u;0].';
    maxnorm=[maxnorm max(abs(u-ue(xa)))];
    [enorm,sqnorm]=enorms(kind,nodes,xa,u,nelem);
    %plot(xa,u);hold;plot(xa,ue(xa),'r');pause
    L2norm=[L2norm sqnorm];close
    cputime-ttt
end
FDMmaxnorm=maxnorm
%FDML2norm=L2norm
%save journal1exp01.mat FDMmaxnorm FDML2norm -append;

%save fdmexp02.mat maxnorm L2norm;
function [enorm,sqnorm]=enorms(kind,nodes,x,u,nelem)
sqnorm=0;enorm=0;
for i=1:nelem
    n1=nodes(1,i);n=kind(i)+1;n2=nodes(n,i);wt=(x(n2)-x(n1))/50;
    xi=-1;
    for k=1:51
        [xx,uh,ux,sigh,sigx]=eval(i,xi,kind,nodes,x,u);
        [xk,xc,xb,xf]=getmat(xx);
        if k==1 || k==51,w=wt/2;else w=wt;end
        xi=xi+0.02;
        sqnorm=sqnorm+w*(ux-uh)^2;
        enorm=enorm+(sigx-sigh)^2*w/xk+(ux-uh)^2*w*xb;
    end
end
sqnorm=sqrt(sqnorm);enorm=sqrt(enorm);

function [psi,dpsi]=shape(xi,n)
switch n
case 2,
    % Linear base function
    psi(1)=(1-xi)/2;
    psi(2)=(1+xi)/2;
    dpsi(1)=-0.5;
    dpsi(2)=0.5;
    return
    case 3,
        % quadratic base function
        psi(1) = xi*(xi-1)*0.5;
        psi(2) = 1-xi*xi;
        psi(3) = xi*(xi+1)*0.5;
        dpsi(1) =  xi-0.5;
        dpsi(2) = -2*xi;
        dpsi(3) = xi + 0.5;
        return
        case 4,
            % cubic  base function
            psi(1) = 9*(1/9-xi*xi)*(xi-1)/16;
            psi(2) = 27*(1-xi*xi)*(1/3-xi)/16;
            psi(3) = 27*(1-xi*xi)*(1/3+xi)/16;
            psi(4) = -9*(1/9-xi*xi)*(1+xi)/16;
            dpsi(1) = -9*(3*xi*xi-2*xi-1/9)/16;
            dpsi(2) = 27*(3*xi*xi-2*xi/3-1)/16;
            dpsi(3) = 27*(-3*xi*xi-2*xi/3+1)/16;
            dpsi(4) = -9*(-3*xi*xi-2*xi+1/9)/16;
            return
end

function [xk,xc,xq,xf] = getmat(x)
%              -(k(x) u_x)_x + c(x) u_x + q(x) u(x) = f(x)                %
xk=exp(1-x);    xc=0;   xq=0;   xf=1+exp(1-x);

function [xx,uh,ux,sigh,sigx]=eval(nel,xi,kind,nodes,x,u)
n=kind(nel)+1;
[psi,dpsi]=shape(xi,n);
uh=0;sigh=0;
for i=1:n
    i1=nodes(i,nel);
    if i==1,x1=x(i1);end
    if i==n,x2=x(i1);end
    uh=uh+psi(i)*u(i1);
    sigh=sigh+dpsi(i)*u(i1);
end
dx=(x2-x1)*0.5;
xx=x1+(1+xi)*dx;
[xk,xc,xb,xf]=getmat(xx);
sigh=sigh*xk/dx;
[ux,sigx]=exact(xx);

function [u,sig]=exact(x)
u=(x.*(1-exp(x-1))).';sig=(1-exp(x-1)-x.*exp(x-1)).';
%u=(4*x.^2.5.*(1-x).^2.5).';sig=(10*x.^(3/2).*(1-x).^(5/2)-10*x.^(5/2).*(1-x).^(3/2)).';
%u=(x-sinh(x)/sinh(1)).';sig=(1-cosh(x)/sinh(1)).';