

1==1;

%tests of 1d fem
disp("equidistant mesh u=sin(t) c=0 a=0 b=pi");
N=20;M=6;
f=s=@sin;
err=zeros(M,3);
for k=1:M,
 x=linspace(0,pi,N)';
 [u,h,A,B]=FEM1dDirSol16(x,f);
 er=u-s(x);
 err(k,1)=eH1=sqrt(er'*A*er);
 err(k,2)=eL2=sqrt(er'*B*er);
 err(k,3)=einf=norm(er,'inf');
 N=2*N;
end
err
err(1:M-1,:)./err(2:M,:)


disp("equidistant mesh u=sin(2t) c=3 a=0 b=1");
N=20;M=6;
c=3;a=0;b=1;
s=@(x) sin(2*x); f=@(x) (4+c)*sin(2*x);
err=zeros(M,3);
for k=1:M,
 x=linspace(0,1,N)';
 [u,h,A,B]=FEM1dDirSol16(x,f,c,s(a),s(b));
 er=u-s(x);
 err(k,1)=eH1=sqrt(er'*A*er);
 err(k,2)=eL2=sqrt(er'*B*er);
 err(k,3)=einf=norm(er,'inf');
 N=2*N;
end
err
err(1:M-1,:)./err(2:M,:)



%mesh twice denser in the 2nd part of interval
disp("equidistant mesh on [0,1] twice denser on [1,2]  u=sin(2t)");
N=10;M=6;
c=0;a=0;b=2;
s=@(x) sin(2*x); f=@(x) (4+c)*sin(2*x);
err=zeros(M,3);
for k=1:M,
 x1=linspace(a,(a+b)/2,N)';x2=linspace((a+b)/2,b,2*N)';
 x=[x1;x2(2:end)];
 [u,h,A,B]=FEM1dDirSol16(x,f,c,s(a),s(b));
 er=u-s(x);
 err(k,1)=eH1=sqrt(er'*A*er);
 err(k,2)=eL2=sqrt(er'*B*er);
 err(k,3)=einf=norm(er,'inf');
 N=2*N;
end
err
err(1:M-1,:)./err(2:M,:)


%mesh scaled sin(x) for x equidistant in the 2nd part of interval
disp("irregular mesh");
N=10;M=6;
c=0;a=0;b=1;
s=@(x) sin(2*x); f=@(x) (4+c)*sin(2*x);
err=zeros(M,3);
for k=1:M,
 x1=linspace(a,b,N)';
 x=b*sin(x1)/sin(b);
 [u,h,A,B]=FEM1dDirSol16(x,f,c,s(a),s(b));
 er=u-s(x);
 err(k,1)=eH1=sqrt(er'*A*er);
 err(k,2)=eL2=sqrt(er'*B*er);
 err(k,3)=einf=norm(er,'inf');
 N=2*N;
end
err
err(1:M-1,:)./err(2:M,:)











