% File: testchaos.m
%
% The goal is to apply Theorem 1. An upper bound for ||L^-1|| is needed.
% This bound is obtained with Lemma 2 and Lemma 3. In order to apply
% these lemmas we need the matrices An, Un and Rn. 
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First, interval variables are computed for f(u_i), f(v_i) and 
% Df(u_i), Df(v_i). These results are stored in the variables
% henonmapu.f, henonmapv.f, henonmapu.df and henonmapv.df.
% Au and Av contain point approximations of Df(u_i), Df(v_i) 
% (Au and Av are taken as the midpoints of the intervals for Df(u_i), 
% Df(v_i)).

  henonmapu=henon(uu,par);
  henonmapv=henon(vv,par);
  Au=henonmapu.df.mid;
  Av=henonmapv.df.mid;

%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we triangularize the matrices Au, Av. 
% QRu.o contains the (almost) orthogonal matrices U(u_i).
% QRu.o is an interval matrix.
% QRu.t contains the upper triangular matrices R(u_i).
% QRu.t is a point matrix.
% QRv.o and QRv.t contain U(v_i) and R(v_i).

  QRu=triang(Au,nn);
  QRv=triang(Av,nn);

%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We apply Lemma 3 to obtain a bound for ||R^-1||. 
% The variables au, bu and cu contain the absolute value of the nonzero 
% elements of R(u_i) (analogously for v).

  au=abs(squeeze(QRu.t(1,1,:)));
  bu=abs(squeeze(QRu.t(1,2,:)));
  cu=abs(squeeze(QRu.t(2,2,:)));
  av=abs(squeeze(QRv.t(1,1,:)));
  bv=abs(squeeze(QRv.t(1,2,:)));
  cv=abs(squeeze(QRv.t(2,2,:)));

% The numbers B, lambda and mu are computed. 

  B=intval(max([bu;bv]));
  lambda0=min([inf(prod(intval(au))^(1/nn)),inf(prod(intval(av))^(1/nn))]);
  lambda=intval(.8*lambda0+.2);
  mu0=max(sup(prod(intval(cu))^(1/nn)),sup(prod(intval(cv))^(1/nn)));
  mu=intval(.9*mu0+.1);

% The computation of P.
% First m_j^u, M_j^u, p_j^u and m_j^v, M_j^v, p_j^v are needed.

  mju=1; Mju=1; P1ju=1; pju=1; 
  mjv=1; Mjv=1; P1jv=1; pjv=1; 
  for j=1:nn,
    qu=au(j)/lambda;
    mju=min(inf(mju*qu),1);
    Mju=min(Mju,mju);
    P1ju=P1ju*qu;
    pju=min(inf(P1ju),pju);

    qv=av(j)/lambda;
    mjv=min(inf(mjv*qv),1);
    Mjv=min(Mjv,mjv);
    P1jv=P1jv*qv;
    pjv=min(inf(P1jv),pjv);
  end
  mn=intval(min(mju,mjv));
  Mn=min(Mju,Mjv);
  pn=intval(min(pju,pjv));
  P=intval(min(Mn,inf(mn*pn)));

% The computation of Q.
% First n_j^u, N_j^u, q_j^u and n_j^v, N_j^v, q_j^v are needed.

  nju=1; Nju=1; Q1ju=1; qju=1; 
  njv=1; Njv=1; Q1jv=1; qjv=1; 
  for j=1:nn,
    qu=cu(j)/mu;
    nju=max(sup(nju*qu),1);
    Nju=max(Nju,nju);
    Q1ju=Q1ju*qu;
    qju=max(sup(Q1ju),qju);

    qv=cv(j)/mu;
    njv=max(sup(njv*qv),1);
    Njv=max(Njv,njv);
    Q1jv=Q1jv*qv;
    qjv=max(sup(Q1jv),qjv);
  end
  nnn=intval(max(nju,njv));
  Nn=max(Nju,Njv);
  qn=intval(max(qju,qjv));
  Q=intval(max(Nn,sup(nnn*qn)));

% The bound Rinfnorm for ||R^-1|| is computed.

  alpha=1/(P*(lambda-1));
  gamma=1+Q*mu/(1-mu);
  alphaplusbeta=alpha*(1+B*gamma);
  Rinfnorm=sqrt(alphaplusbeta^2+gamma^2);

% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We apply Lemma 2 to obtain a bound for ||L^-1||. 
%
% Computation of delta2 and delta3.
% The matrices A(u_i), U(u_i), R(u_i) are already available in the
% variables Au, QRu.o and QRu.t (analogously for v_i).
%
% delta2 is computed as follows. The norm of 
% Du_i = A(u_i)*U(u_i)-U(u_i+1)*R(u_i) is computed as the square root
% of the sum of the squares of all elements (this is of the order of 
% rounding errors).
% For A(u_N)*U(u_N)-U(v_1)*R(u_N) and A(v_N)*U(v_N)-U(u_1)*R(v_N)
% improved estimates are computed. For a 2x2-matrix A=[a b; c d] define
% ss=a^2+b^2+c^2+d^2 and det=ad-bc. Then
% ||A|| = sqrt((ss+sqrt(ss^2-4*det^2))/2).
%
% For the computation of delta3 we use the simplified estimate of the
% norm of a matrix ||A|| <= ss (is of the order of rounding errors).

  delt2=0;
  delt3=0;
  for i=1:nn,
    Du=Au(:,:,i)*QRu.o(:,:,i)-QRu.o(:,:,i+1)*QRu.t(:,:,i);
    Dv=Av(:,:,i)*QRv.o(:,:,i)-QRv.o(:,:,i+1)*QRv.t(:,:,i);
    D3u=QRu.o(:,:,i)'*QRu.o(:,:,i)-eye(2);		% used for delta3
    D3v=QRv.o(:,:,i)'*QRv.o(:,:,i)-eye(2); 		% used for delta3
    d2u=sqrt(intval(sup(sum(sum(Du.*Du)))));
    d2v=sqrt(intval(sup(sum(sum(Dv.*Dv)))));
    d3u=sqrt(intval(sup(sum(sum(D3u.*D3u)))));		% used for delta3
    d3v=sqrt(intval(sup(sum(sum(D3v.*D3v)))));		% used for delta3
    d2=max(d2u.sup,d2v.sup);
    d3=max(d3u.sup,d3v.sup);				% used for delta3
    if delt2<d2, delt2=d2; end;
    if delt3<d3, delt3=d3; end;				% used for delta3
  end;
  Du=Au(:,:,nn)*QRu.o(:,:,nn)-QRv.o(:,:,1)*QRu.t(:,:,nn);
  Dv=Av(:,:,nn)*QRv.o(:,:,nn)-QRu.o(:,:,1)*QRv.t(:,:,nn);
  sumu=sum(sum(Du.*Du));
  sumv=sum(sum(Dv.*Dv));
  detu=Du(1,1)*Du(2,2)-Du(2,1)*Du(1,2);
  detv=Dv(1,1)*Dv(2,2)-Dv(2,1)*Dv(1,2);
  d2u=sqrt((sumu+sqrt(sumu^2-4*detu^2))/2);
  d2v=sqrt((sumv+sqrt(sumv^2-4*detv^2))/2);
  d2=max(d2u.sup,d2v.sup);
  if delt2<d2, delt2=d2; end;
  delta2=intval(delt2)
  delta3=intval(delt3)

% Condition (11) is verified. If the condition is satisfied then the
% bound of Lemma 2 for ||L^-1|| is computed (denoted by Linfnorm).

  if delta3>=1, 
    '**********************************'
    'delta3 >= 1, no estimates possible'
    '**********************************'
    return;
  end;
  cond=sqrt(1-delta3)-Rinfnorm*delta2;
  if cond>0,
    Linfnorm=Rinfnorm*sqrt(1-delta3)/cond;
  else
    '******************************************************************'
    '||R^-1||*delta2 is to large, no estimate of ||L^-1|| can be given.'
    '******************************************************************'
    return;
  end;


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Theorem 1 is applied. 
%
% The constants delta, delta1 and M are determined as small as possible 
% such that the estimates (2), (3) and (4) are satisfied. 
%
% First we compute delta.
% fy_i contains the interval values of f(u_i), f(v_i), f(u_N), f(v_N);
% y_i1 contains the point values of u_i+1, v_i+1, v_1, u_1 (We use the 
% convention  u_N+1=u_1,  v_N+1=v_1).
% These two variables are 2x(2N+2) matrices.
% The variable "diffnorm" is the vector of the norms of all possible 
% 2N+2 the differences y_i+1 - f(y_i).
% delta is the maximal element of the variable "diffnorm".

  fy_i=[henonmapu.f, henonmapv.f, henonmapu.f(:,nn), henonmapv.f(:,nn)];
  y_i1=[uu(:,[2:nn,1]), vv(:,[2:nn,1]), vv(:,1), uu(:,1)];
  diff=y_i1-fy_i;
  diffnorm=sqrt(diff(1,:).^2+diff(2,:).^2);
  delta=intval(max(diffnorm.sup));

% Computation of delta1.
% henonmapu.df.rad and henonmapv.df.rad are bounds for the matrices
% A_i-Df(y_i).
% Note that the second column of Df(y_i) is [1;0] and therefore it is
% known exactly. The norm of A_i-Df(y_i) may be computed as the norm
% of the first column vector.

  Adiff=[squeeze(henonmapu.df.rad(:,1,:)), squeeze(henonmapv.df.rad(:,1,:))];
  Adiff=intval(Adiff);
  normA=sqrt(Adiff(1,:).^2+Adiff(2,:).^2);
  delta1=intval(max(normA.sup));

% M = 2*a

  M=intval(sup(2*par(1)));

% compute Ltheorem = 1/(delta1+sqrt(2Mdelta)).

  Ltheorem = 1/(delta1+sqrt(2*M*delta));

% Condition (1) is verified. If it is satisfied the numbers r0, r1, 
% (rho=r1), are computed and a bound for ||L^-1|| is displayd.

  if Linfnorm>Ltheorem,
    '       **********************************************'
    '       Linfnorm too large, chaos cannot be guaranteed'
    '       **********************************************'
    return
  end
  '       **************************************'
  '       all conditions for chaos are satisfied'
  '       **************************************'
  diff=1/Linfnorm-delta1;
  hh=diff+sqrt(diff^2-2*delta*M);
  r0=2*delta/hh
  r1=hh/M
  Linfnorm
  Ltheorem

%
% written   05/20/99   D.Stoffer
%


