%
% Delay DE well-conditioned according to eigenvalues,
% seen to be ill-conditioned by pseudospectra
%
close all
clear all
n = 10;
% Differentiation matrix on tau = cos(pi*(0:n-1)/(n-1))
D2=gallery('chebspec',n);
e = eig(D2)-1;
figure
% Eigenvalues exactly -1, multiple
plot(e,'ko')
hold on
nu=linspace(-pi/2,pi/2,407);
% Envelope where Re(W(k,lambda)) < 0
plot(-nu.*sin(nu),nu.*cos(nu),'k')
axis([-1.75,0.25,-1,1])
axis('square')
% sin(pi*tau) initial vector, constant history
hist = dde23( @(t,y,z) zeros(size(y)), [], ...
              sin(pi*cos(pi*(0:n-1)/(n-1))), [-1,0]);
opts=ddeset('abstol',1.0e-12,'reltol',1.0e-12);
% y'(t) = (D2-I)*y(t-1)
chebspecdde = @(t,y,z) D2*z-z;
% dde23 takes 58,664 steps to get to 4...
% ddesd takes 4,676
Tf = 8;
sol=ddesd(chebspecdde,1,hist,[0,Tf],opts);
figure
semilogy(sol.x,abs(sol.y),'k.')
axis([0,Tf,1.0e-3,1.0e9])
% Residual is less than 1.0e-4 (1.0e-7*||y||)
h=diff(sol.x);
samp = [sol.x(1:end-1)+h/3;sol.x(1:end-1)+2*h/3];
t=samp(:);
k=find(t>0);
tt=t(k);
[y,yp]=deval(sol,tt);
z=deval(sol,tt-1);
r=zeros(size(z));
[dummy,m]=size(z);
nrms=zeros(size(z));
for i=1:m,
    r(:,i)=yp(:,i)-chebspecdde(tt(i),y(:,i),z(:,i));
    nrms(:,i)=norm(y(:,i),2)*ones(n,1);
end
figure
semilogy(tt,abs(r)./nrms,'k.')
set(gca,'fontsize',16);
xlabel('t','fontsize',16);
ylabel('relative residual','fontsize',16)



