

%we test local truncation error of the FDM scheme
%we want to solve -u''+cu=f in (a,b) with u()=al u(b)=be
%by a FDM method on uniform mesh (x(k)=a+(k-1)*h) k=1,N+1
%u(1)=al bnd cpndition at a
% (1/(h*h)(-u(k-1)+2u(k)-u(k+1)) +c*u(k)=f(x(k)) k=2,..,N  
%u(N+1)=be
%||Lh rh(u) -fh|| =O(h^p)  ||lh rh(u)-gf||
%where Lh is the FDM operator,lh boundary FDM operator;
% and rh restriction operator
%fh, gh the respective rhs of the discrete problems 
%In our case we have to compute 
% e(h)=discrete max norm of 
%Lh rh(u) -fh=(1/(h*h)(-u(k-1)+2u(k)-u(k+1)) - f(x(k)) k=2,..N
% (boundary conditions are satisfied with the error zero)
%We assume that we know the solution,
%we can use the so called halving method: e(h)/e(0.5h) \approx 2^p

%consider -u''=sin(x) u(0)=0 u(b)=sin(b) (the solution is sin function)

a=0;b=pi;
sol=@sin;
f=@sin;
N=10;
M=7;
e=zeros(M,1);
for k=1:M,
 h=(b-a)/N;
 xh=linspace(a,b,N+1);xh=xh';
 rhu=sol(xh);
 Lhrhu=(-rhu(1:N-1)+2*rhu(2:N)-rhu(3:N+1))/(h*h);
 rhf=f(xh(2:N));
 e(k)=norm(Lhrhu-rhf,"inf"); %discrete max norm
 N*=2;
 [Lhrhu,rhf]
end
ratio=e(1:M-1)./e(2:M) 
log2(ratio)
 
 
 
 