% EXAMPLE3_53 Computations for Example 3_53.
%
fprintf('\n')
f1='%10.0f %12.3e %12.3e %12.3e %12.3e\n';
disp('         n       E2          Einf0        Einf1        Einf2')
N=5; s=2; same=1;
%
% Sobolev least squares
%
s=2; Ns=N*(s+1);
t=linspace(0,2,N)'; nd=N*ones(1,s+1);
xw=ones(N,2*s+2)/N;
for sig=1:s+1
  xw(:,sig)=t;
end
xw0=xw(:,s+1:2*s+2);
a0=sum(t)/N;
[B,normsq]=stieltjes_sob(3*N-1,s,nd,xw,a0,same);
f=[exp(t.^2).*erfc(t) 2*t.*exp(t.^2).*erfc(t)-2/sqrt(pi)...
   2*(1+2*t.^2).*exp(t.^2).*erfc(t)-4*t/sqrt(pi)];
for n=0:3*N-1
  [phat,c]=least_squares_sob(n,f,xw0,B);
end
x=linspace(0,2)';
fx=[exp(x.^2).*erfc(x) 2*x.*exp(x.^2).*erfc(x)-2/sqrt(pi)...
    2*(1+2*x.^2).*exp(x.^2).*erfc(x)-4*x/sqrt(pi)];
E2=zeros(3*N,1); 
Einf0=zeros(3*N,1); Einf1=zeros(3*N,1); Einf2=zeros(3*N,1);
for n=1:3*N
  for sig=1:s+1
    sigN=sig:s+1:Ns;
    E2(n)=E2(n)+sum(xw0(:,sig+1)'.*...
      (phat(n,sigN)-f(:,sig)').^2);
  end
  E2(n)=sqrt(E2(n));
  p=clenshaw_sob(n-1,3*N-1,s,x,B,c');
  Einf0(n)=max(abs(p(:,1)-fx(:,1)));
  Einf1(n)=max(abs(p(:,2)-fx(:,2)));
  Einf2(n)=max(abs(p(:,3)-fx(:,3)));
  fprintf(f1,n-1,E2(n),Einf0(n),Einf1(n),Einf2(n))
end
%
% Ordinary least squares
%
fprintf('\n')
s=0; Ns=N*(s+1);
t=linspace(0,2,N)'; nd=N*ones(1,s+1);
xw=ones(N,2*s+2)/N;
for sig=1:s+1
  xw(:,sig)=t;
end
xw0=xw;
a0=sum(t)/N;
[B,normsq]=stieltjes_sob(N-1,s,nd,xw,a0,same);
f=[exp(t.^2).*erfc(t) 2*t.*exp(t.^2).*erfc(t)-2/sqrt(pi)...
   2*(1+2*t.^2).*exp(t.^2).*erfc(t)-4*t/sqrt(pi)];
for n=0:N-1
  [phat,c]=least_squares_sob(n,f,xw0,B);
end
x=linspace(0,2)';
fx=[exp(x.^2).*erfc(x) 2*x.*exp(x.^2).*erfc(x)-2/sqrt(pi)...
    2*(1+2*x.^2).*exp(x.^2).*erfc(x)-4*x/sqrt(pi)];
E2=zeros(N,1); 
Einf0=zeros(N,1); Einf1=zeros(N,1); Einf2=zeros(N,1);
for n=1:N
  for sig=1:s+1
    sigN=sig:s+1:Ns;
    E2(n)=E2(n)+sum(xw0(:,sig+1)'.*...
      (phat(n,sigN)-f(:,sig)').^2);
  end
  E2(n)=sqrt(E2(n));
  p=clenshaw_sob(n-1,N-1,s+2,x,B,c');
  Einf0(n)=max(abs(p(:,1)-fx(:,1)));
  Einf1(n)=max(abs(p(:,2)-fx(:,2)));
  Einf2(n)=max(abs(p(:,3)-fx(:,3)));
  fprintf(f1,n-1,E2(n),Einf0(n),Einf1(n),Einf2(n))
end
