

%test convergence order fot -u''+cu=f on [a,b]
%u(a)=al u(b) be for known solution

%1/  u=sin(x) f=sin(x) 

a=0;
b=1;
c=1;
s=@ sin;
f=@(x) c*s(x)+s(x);
al=s(a);be=s(b)
N=20;
M=6;
for k=1:M,
 [u,x]=FDMDsolver(f,a,b,al,be,c,N);
 h=x(2)-x(1);
 plot(x,u,"+",x,ss=s(x));pause(1);
 erinf(k)=norm(er=u-ss,'inf');
 el2(k)=norm(er,2)*sqrt(h);
 N*=2;
end

er=[erinf',el2']
%ratio of erinf(h)/erinf(h/2)
rer=er(1:M-1,:)./er(2:M,:)
%order
log2(rer)


%test what happens if the discrete boundary is on the original bnd
%let x(N)<b and b-x(N)< d*h 0<d<1
%let h=(b-a)/(N+d) d random
%and in our FDM scheme take the last equation x(N)=be=s(b)

a=0;
b=1;
c=1;
s=@ sin;
f=@(x) c*s(x)+s(x);
al=s(a);be=s(b)
N=20;
M=6;
d=0.2+0.6*rand;
for k=1:M,
 h=(b-a)/(N+d); 
 xN=a+N*h;
 [u,x,A,F]=FDMDsolver(f,a,xN,al,be,c,N);
 xx=linspace(a,b,200);
 ss=s(x);
 plot(x,u,"+",xx,s(xx));pause(1);
 plot(x,u-ss,";error;");pause(1);
 erinf(k)=norm(er=u-ss,'inf');
 el2(k)=norm(er,2)*sqrt(h);
 N*=2;
end

er=[erinf',el2']
%ratio of erinf(h)/erinf(h/2)
rer=er(1:M-1,:)./er(2:M,:)
%order
log2(rer)

 %now repeat the tests replacing the last equation of FDM scheme with
 %the Collatz approximation:
 %-u(N-1)*d +u(N)(1+d)=be  here b-x(N)=d*h
 
a=0;
b=1;
c=1;
s=@ sin;
f=@(x) c*s(x)+s(x);
al=s(a);be=s(b)
N=20;
M=6;
d=0.2+0.6*rand;
for k=1:M,
 h=(b-a)/(N+d); 
 xN=a+N*h;
 [u,x,A,F]=FDMDsolver(f,a,xN,al,be,c,N);%just to get A and F
 %we modify only the last equation!!!!! (for d=0 there is no modification!)
 
 A(N+1,N+1)=1+d;A(N+1,N)=-d;
 u=A\F;
 xx=linspace(a,b,200);
 ss=s(x);
 plot(x,u,"+",xx,s(xx));pause(1);
 plot(x,u-ss,"+;Collatz error;");pause(1);
 erinf(k)=norm(er=u-ss,'inf');
 el2(k)=norm(er,2)*sqrt(h);
 N*=2;
end

er=[erinf',el2']
%ratio of erinf(h)/erinf(h/2)
rer=er(1:M-1,:)./er(2:M,:)
%order
log2(rer)

%test mixed boundary condition i.e.


%test convergence order fot -u''+cu=f on [a,b]
%u(a)=al u'(b)= be for known solution
%we approximate the Neumann bnd at b by the backward difference 
%(u(N)-u(N-1))/h=be
%1/  u=sin(x) f=sin(x) 

a=0;
b=1;
c=1;
s=@ sin;
f=@(x) c*s(x)+s(x);
al=s(a);
be=cos(b) %d/dx s(b)=cos(b) since s(x)=sin(x)
N=20;
M=6;
for k=1:M,
 [u,x,A,F]=FDMDsolver(f,a,b,al,be,c,N);%just to get A and F
 %we modify only the last equation!!!!! 
 h=x(2)-x(1);
 A(N+1,N+1)=1/h;A(N+1,N)=-1/h;%backward difference
 u=A\F;
 xx=linspace(a,b,200);
 ss=s(x);
 plot(x,u,"+",xx,s(xx));pause(2);
 %plot(x,u-ss,"+;error;");pause(1);
 erinf(k)=norm(er=u-ss,'inf');
 el2(k)=norm(er,2)*sqrt(h);
 N*=2;
end

er=[erinf',el2']
%ratio of erinf(h)/erinf(h/2)
rer=er(1:M-1,:)./er(2:M,:)
%order
log2(rer)
%we got the 1st order convergence

%add a ghost point and use central difference in order to get 2nd order convegence 
%i.e. assume that the ODE is satsfied at x=b and add a 3 point stencil equation for x(n) and 
%central difference equation at x(N) i.e.
%(-x(N-1)+2*x(N)-x(N+1))/(h*h)=f(b)
%(x(N+1)-x(N-1))/(2*h)=be


a=0;
b=1;
c=1;
s=@ sin;
f=@(x) c*s(x)+s(x);
al=s(a);
be=cos(b) %d/dx s(b)=cos(b) since s(x)=sin(x)
N=20;
M=6;
for k=1:M,
 h=(b-a)/N; 
 [u,x,A,F]=FDMDsolver(f,a,b+h,al,be,c,N+1);%just to get A and F
 %we modify only the last equation but be careful the system is larger !!!!! 
 A(N+2,N+2)=1/(2*h); A(N+2,N)=-1/(2*h);%backward difference
 u=A\F;
 x=x(1:N+1);u=u(1:N+1);
 xx=linspace(a,b,200);
 ss=s(x);
 plot(x,u,"+",xx,s(xx));pause(2);
 %plot(x,u-ss,"+;error;");pause(1);
 erinf(k)=norm(er=u-ss,'inf');
 el2(k)=norm(er,2)*sqrt(h);
 N*=2;
end

er=[erinf',el2']
%ratio of erinf(h)/erinf(h/2)
rer=er(1:M-1,:)./er(2:M,:)
%order
log2(rer)
%we got the second order convergence


%in order to get 2nd order convegence 
%assume that the ODE is satisfied at x=b and 
%add a 3 point stencil equation for x(n) and 
%use the equation to get 2nd order scheme
%we have (u(b)-u(b-h))/h-be= -0.5*u''(b)h + O(h^2)
%we use the eqution to replace u''(b) 
%we have -u''+cu=f or u''(b)=cu(b)-f(b) thus we get 
%BDu(b)-be + u''(b)h= BDu(b)-be + 0.5*(cu(b)-f(b))h=O(h^2) thus we repace the last 
%equation in FDM scheme (u(N)-u(N-1))/h=be by 
%(u(N)-u(N-1))/h +0.5*cu(N)*h=be+0.5*h*f(b) or 
%(1/h+0.5*c*h)*u(N)-(1/h)u(N-1)=be+0.5*h*f(b) 


a=0;
b=1;
c=1;
s=@ sin;
f=@(x) c*s(x)+s(x);
al=s(a);
be=cos(b) %d/dx s(b)=cos(b) since s(x)=sin(x)
N=20;
M=6;
for k=1:M,
 h=(b-a)/N; 
 [u,x,A,F]=FDMDsolver(f,a,b,al,be,c,N);%just to get A and F
 %we modify only the last equation but be careful the system is larger !!!!! 
 A(N+1,N+1)=1/h+0.5*c*h;
 A(N+1,N)=-1/(h);%backward difference +modification
 F(N+1)+=0.5*h*f(b);
 u=A\F;
 x=x(1:N+1);u=u(1:N+1);
 xx=linspace(a,b,200);
 ss=s(x);
 plot(x,u,"+",xx,s(xx));pause(2);
 %plot(x,u-ss,"+;error;");pause(1);
 erinf(k)=norm(er=u-ss,'inf');
 el2(k)=norm(er,2)*sqrt(h);
 N*=2;
end

er=[erinf',el2']
%ratio of erinf(h)/erinf(h/2)
rer=er(1:M-1,:)./er(2:M,:)
%order
log2(rer)
%we got the second order convergence
