%***********************************************************************
%*****                                                             *****
%*****  Package  F I D E  --- T e s t Examples  Ver. 1.1  05/1991  *****
%*****                                                             *****
%***********************************************************************
%***********************************************************************
%*****                                                             *****
%*****    T e s t     Examples   ---   Module    E X P R E S       *****
%*****                                                             *****
%***********************************************************************

let cos th**2=1 - sin th**2,
    cos fi**2=1 - sin fi**2;
factor df;
on rat;
for all x,y let diff(x,y)=df(x,y);
depend u,r,th,fi;
depend v,r,th,fi;
depend f,r,th,fi;
depend w,r,th,fi;
% Spherical coordinate system
scalefactors 3,r*sin th*cos fi,r*sin th*sin fi,r*cos th,r,th,fi;
tensor a1,a2,a3,a4,a5;
vectors u,v;
dyads w;
a1:=grad f;
a2:=div u;
a3:=curl v;
a4:=lapl v;
a3:=2*a3+a4;
a5:=lapl f;
a1:=a1+div w;
a1:=u.dyad((a,0,1),(1,b,3),(0,c,d));
a2:=vect(a,b,c);
a1.a2;
%  Scalar product
u.v;
%  Vector product
u?v;
%  Dyadic
u&v;
%  Directional derivative
dirdf(u,v);
clear a1,a2,a3,a4,a5,u,v,w;
for all x,y clear diff(x,y);
clear cos th**2,
      cos fi**2;
remfac df;
off rat;
scalefactors 3,x,y,z,x,y,z;

%***********************************************************************
%*****                                                             *****
%*****      T e s t     Examples   ---   Module    I I M E T       *****
%*****                                                             *****
%***********************************************************************

% Example I.1 - 1-D Lagrangian Hydrodynamics

off exp;
factor diff;
on rat,eqfu;

% Declare which indexes will be given to coordinates
coordinates x,t into j,m;

% Declares uniform grid in x coordinate
grid uniform,x;

% Declares dependencies of functions on coordinates
dependence eta(t,x),v(t,x),eps(t,x),p(t,x);

% Declares p as known function
given p;

same eta,v,p;

iim a, eta,diff(eta,t)-eta*diff(v,x)=0,
    v,diff(v,t)+eta/ro*diff(p,x)=0,
    eps,diff(eps,t)+eta*p/ro*diff(v,x)=0;

clear a;
clearsame;
cleargiven;

%***********************************************************************

% Example  I.2 - How other functions (here sin, cos) can be used in
%                discretized terms

diffunc sin,cos;
difmatch all,diff(u*sin x,x),u=one,2,(u(i+1)*sin x(i+1)-u(i-1)
                   *sin x(i-1))/(dim1+dip1),
                 u=half,0,(u(i+1/2)*sin x(i+1/2)-u(i-1/2)*sin x(i-1/2))
                   /di;
difmatch all,cos x*diff(u,x,2),u=one,0,cos x i*(u(i+1)-2*u(i)+u(i-1))
                   /di^2,
                 u=half,3,(u(i+3/2)-u(i+1/2))/dip2/2 -
                          (u(i-1/2)-u(i-3/2))/dim2/2;
off exp;
coordinates x,t into j,m;
grid uniform,x,t;
dependence u(x,t),v(x,t);
iim a,u,diff(u,t)+diff(u,x)+cos x*diff(v,x,2)=0,
      v,diff(v,t)+diff(sin x*u,x)=0;
clear a;

%***********************************************************************

% Example I.3 - Schrodinger equation

factor diff;
coordinates t,x into m,j;
grid uniform,x,t;
dependence ur(x,t),ui(x,t);
same ui,ur;
iim a,ur,-diff(ui,t)+1/2*diff(ur,x,2)+(ur**2+ui**2)*ur=0,
       ui,diff(ur,t)+1/2*diff(ui,x,2)+(ur**2+ui**2)*ui=0;
clear a;
clearsame;

%***********************************************************************

% Example I.4 - Vector calculus in p.d.e. input
%               cooperation with expres module
%               2-D hydrodynamics

scalefactors 2,x,y,x,y;
vectors u;
off exp,twogrid;
on eqfu;
factor diff,ht,hx,hy;
coordinates x,y,t into j,i,m;
grid uniform,x,y,t;
dependence n(t,x,y),u(t,x,y),p(t,x,y);
iim a,n,diff(n,t)+u.grad n+n*div u=0,
     u,m*n*(diff(u,t)+u.grad u)+grad p=vect(0,0),
     p,3/2*(diff(p,t)+u.grad p)+5/2*p*div u=0;
clear a,u;

%***********************************************************************

% Example I.5 - 1-D hydrodynamics up to 3-rd moments (heat flow)

coordinates x,t into j,m;
grid uniform,x,t;
dependence n(x,t),u(x,t),tt(x,t),p(x,t),q(x,t);
iim a, n,diff(n,t)+u*diff(n,x)+diff(u,x)=0,
    u,n*m*(diff(u,t)+u*diff(u,x))+k*diff(n*tt,x)+diff(p,x)=0,
    tt,3/2*k*n*(diff(tt,t)+u*diff(tt,x))+n*k*tt*diff(u,x)+1/2*p
       *diff(u,x)+diff(q,x)=0,
    p,diff(p,t)+u*diff(p,x)+p*diff(u,x)+n*k*tt*diff(u,x)+2/5*diff(q,x)
       =0,
    q,diff(q,t)+u*diff(q,x)+q*diff(u,x)+5/2*n*k**2*tt/m*diff(tt,x)+n*k
       *tt*diff(p,x)-p*diff(p,x)=0;
clear a;

remfac diff,ht,hx,hy;
on exp;
off rat;

%***********************************************************************
%*****                                                             *****
%*****     T e s t     Examples   ---   Module    A P P R O X      *****
%*****                                                             *****
%***********************************************************************

% Example A.1

coordinates x,t into j,n;
maxorder t=2,x=3;
functions u,v;
approx( (u(n+1/2)-u(n-1/2))/ht=(v(n+1/2,j+1/2)-v(n+1/2,j-1/2)
                               +v(n-1/2,j+1/2)-v(n-1/2,j-1/2))/(2*hx) );

% Example A.2

maxorder t=3,x=3;
approx( (u(n+1)-u(n))/ht=(u(n+1,j+1/2)-u(n+1,j-1/2)
                               +u(n,j+1/2)-u(n,j-1/2))/(2*hx) );

% Example A.3

maxorder t=2,x=3;
center t=1/2;
approx( (u(n+1)-u(n))/ht=(v(n+1,j+1/2)-v(n+1,j-1/2)
                               +v(n,j+1/2)-v(n,j-1/2))/(2*hx) );

% Example A.4

approx( u(n+1)/ht=(v(n+1,j+1/2)-v(n+1,j-1/2)
                               +v(n,j+1/2)-v(n,j-1/2))/(2*hx) );

% Example A.5

maxorder t=3,x=3;
approx( (u(n+1)-u(n))/ht=(u(n+1,j+1/2)-u(n+1,j-1/2))/hx);

% Example A.6

approx( (u(n+1)-u(n))/ht=(u(n+1/2,j+1/2)-u(n+1/2,j-1/2))/hx);

% Example A.7;

approx((u(n+1)-u(n))/ht=(u(n+1/2,j+1)-2*u(n+1/2,j)+u(n+1/2,j-1))/hx**2);

%***********************************************************************
%*****                                                             *****
%*****     T e s t     Examples   ---   Module    C H A R P O L    *****
%*****                                                             *****
%***********************************************************************

% Example C.1

coordinates t,x into i,j;
grid uniform,t,x;
let cos ax**2=1-sin ax**2;
unfunc u,v;
matrix aa(1,2),bb(2,2);
aa(1,1):=(u(i+1)-u(i))/ht+(v(j+1)-v(j))/hx$
aa(1,2):=(v(i+1)-v(i))/ht+(u(j+1/2)-u(j-1/2))/hx$
bb:=ampmat aa;
bb:=denotemat bb;
factor lam;
pol:=charpol bb;
prdenot;
cleardenot;
clear aa,bb,pol;

%***********************************************************************

% Example C.2 : Reprint Vorozcov, Ganza, Mazurik: Simvolno-cislennyj
%               interfejs.  v zadacach ..., Novosibirsk 1986,  p.47.

unfunc u;
matrix aa(1,1),bb(1,1);
aa(1,1):=(u(i+1)-u(i))/ht+a*(u(j)-u(j-1))/hx$
bb:=ampmat aa;
bb:=denotemat bb;
pol:=charpol bb;
prdenot;
cleardenot;
clear aa,bb,pol;

%***********************************************************************

% Example C.3 : Reprint Vorozcov, Ganza, Mazurik: Simvolno-cislennyj
%               interfejs.  v zadacach ..., Novosibirsk 1986,  p.52.

coordinates t,x into m,j;
unfunc u,r;
matrix aa(1,2),bb(2,2);
aa(1,1):=(r(m+1)-r(m))/ht+u0*(r(m+1,j+1)-r(m+1,j-1))/2/hx
                         +r0*(u(m+1,j+1)-u(m+1,j-1))/2/hx$
aa(1,2):=(u(m+1)-u(m))/ht+u0*(u(m+1,j+1)-u(m+1,j-1))/2/hx
                         +c0**2/r0*(r(m,j+1)-u(m,j-1))/2/hx$
bb:=ampmat aa;
bb:=denotemat bb;
pol:=charpol bb;
prdenot;
cleardenot;
clear aa,bb,pol;

%***********************************************************************

% Example C.4 : Richtmyer, Morton: Difference methods for initial value
%               problems, &10.3.  p.262

coordinates t,x into n,j;
unfunc v,w;
matrix aa(1,2),bb(2,2);
aa(1,1):=(v(n+1)-v(n))/ht-c*(w(j+1/2)-w(j-1/2)+
                             w(n+1,j+1/2)-w(n+1,j-1/2))/(2*hx)$
aa(1,2):=(w(n+1,j-1/2)-w(n,j-1/2))/ht-c*(v(n+1,j)-v(n+1,j-1)+
                                         v(j)-v(j-1))/(2*hx)$
bb:=ampmat aa;
bb:=denotemat bb;
pol:=charpol bb;
prdenot;
cleardenot;
clear aa,bb,pol;

%***********************************************************************

% Example C.5: Mazurik: Algoritmy resenia zadaci..., Preprint no.24-85,
%               AN USSR SO, Inst. teor. i prikl. mechaniky, p.34

coordinates t,x,y into n,m,k;
grid uniform,t,x,y;
unfunc u1,u2,u3;
matrix aa(1,3),bb(3,3);
aa(1,1):=(u1(n+1)-u1(n))/ht+c/2*((-u1(m-1)+2*u1(m)-u1(m+1))/hx +
         (u2(m+1)-u2(m-1))/hx - (u1(k-1)-2*u1(k)+u1(k+1))/hy +
         (u3(k+1)-u3(k-1))/hy)$
aa(1,2):=(u2(n+1)-u2(n))/ht+c/2*((u1(m+1)-u1(m-1))/hx -
         (u2(m-1)-2*u2(m)+u2(m+1))/hx)$
aa(1,3):=(u3(n+1)-u3(n))/ht + c/2*((u1(k+1)-u1(k-1))/hy -
         (u3(k-1)-2*u3(k)+u3(k+1))/hy)$
off prfourmat;
bb:=ampmat aa;
pol:=charpol bb;
let
  cos ax=cos ax2**2-sin ax2**2,
  cos ay=cos ay2**2-sin ay2**2,
  sin ax=2*sin ax2*cos ax2,
  sin ay=2*sin ay2*cos ay2,
  cos ax2**2=1-sin ax2**2,
  cos ay2**2=1-sin ay2**2,
  sin ax2=s1,
  sin ay2=s2,
  hx=c*ht/cap1,
  hy=c*ht/cap2;
order s1,s2;
pol:=pol;
clear cos ax,cos ay,sin ax,sin ay,cos ax2**2,cos ay2**2,sin ax2,sin ay2,
      hx,hy;
pol:=complexpol pol;
pol1:=hurw pol;
denotid cp;
pol:=denotepol pol;
prdenot;
cleardenot;
clear aa,bb,pol,pol1;

%***********************************************************************

% Example C.6 : Lax-Wendrov (V. Ganzha)

coordinates t,x,y into n,m,k;
grid uniform,t,x,y;
let cos ax**2=1-sin ax**2,
    cos ay**2=1-sin ay**2;
unfunc u1,u2,u3,u4;
matrix aa(1,4),bb(4,4);
aa(1,1):=4*(u1(n+1)-u1(n))/ht+
         (w*(u1(m+2)-u1(m-2)+u1(m+1,k+1)+u1(m+1,k-1)-
         u1(m-1,k+1)-u1(m-1,k-1))+p*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+
         u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+
         v*(u1(m+1,k+1)+u1(m-1,k+1)-
         u1(m+1,k-1)-u1(m-1,k-1)+u1(k+2)-u1(k-2))+p*(u3(m+1,k+1)+
         u3(m-1,k+1)-u3(m+1,k-1)-u3(m-1,k-1)+u3(k+2)-u3(k-2)))/hx+ht*
         (2*w**2*(-u1(m+2)+2*u1(m)-u1(m-2))+4*w*p*(-u2(m+2)+2*u2(m)-
         u2(m-2))+2*(-u4(m+2)+2*u4(m)-u4(m-2))+2*v**2*(-u1(k+2)+
         2*u1(k)-u1(k-2))+4*v*p*(u3(k+2)+2*u3(k)-u3(k-2))+2*(-u4(k+2)+
         2*u4(k)-u4(k-2))+4*w*v*(-u1(m+1,k+1)+u1(m+1,k-1)+u1(m-1,k+1)-
         u1(m-1,k-1))+4*p*v*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)-
         u2(m-1,k-1))+4*w*p*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1)-
         u3(m-1,k-1)))/hx/hx$
aa(1,2):=4*p*(u2(n+1)-u2(n))/ht+
       (w*p*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+
         u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+u4(m+2)-u4(m-2)+
         u4(m+1,k+1)+
         u4(m+1,k-1)-u4(m-1,k+1)-u4(m-1,k-1)+
         p*v*(u2(m+1,k+1)+u2(m-1,k+1)+
         u2(k+2)-u2(k-2)-u2(m+1,k-1)-u2(m-1,k-1)))/hx+ht*(2*w**2*p*
         (-u2(m+2)+2*u2(m)-u2(m-2))+2*p*c**2*(-u2(m+2)+2*u2(m)-u2(m-2))
         +4*w*(-u4(m+2)+2*u4(m)-u4(m-2))+2*p*v**2*(-u2(k+2)+2*u2(k)-
         u2(k-2))+4*w*p*v*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)-
         u2(m-1,k-1))+2*p*c**2*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1)
         -u3(m-1,k-1))+4*v*(-u4(m+1,k+1)+u4(m+1,k-1)+u4(m-1,k+1)-
         u4(m-1,k-1)))/hx/hx$
aa(1,3):=4*p*(u3(n+1)-u3(n))/ht+(w*p*(u3(m+2)-u3(m-2)+u3(m+1,k+1)+
         u3(m+1,k-1)-u3(m-1,k+1)-u3(m-1,k-1))+u4(k+2)-u4(k-2)+
         u4(m+1,k+1)-u4(m+1,k-1)+u4(m-1,k+1)-u4(m-1,k-1)+
         p*v*(u3(m+1,k+1)+u3(m-1,k+1)+u3(k+2)-u3(k-2)-u3(m+1,k-1)-
         u3(m-1,k-1)))/hx+ht*(2*w**2*p*(-u3(m+2)+2*u3(m)-u3(m-2))+
         2*p*c**2*(-u3(k+2)+2*u3(k)-u3(k-2))+4*v*(-u4(k+2)+
         2*u4(k)-u4(k-2))+2*p*v**2*(-u3(k+2)+2*u3(k)-u3(k-2))+
         4*w*p*v*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1)-
         u3(m-1,k-1))+2*p*c**2*(-u2(m+1,k+1)+u2(m+1,k-1)+
         u2(m-1,k+1)-u2(m-1,k-1))+4*w*(u4(m+1,k+1)+u4(m+1,k-1)+
         u4(m-1,k+1)-u4(m-1,k-1)))/hx/hx$
aa(1,4):=4*(u4(n+1)-u4(n))/ht+(p*c**2*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+
         u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+w*(u4(m+2)-
         u4(m-2)+u4(m+1,k+1)+u4(m+1,k-1)-u4(m-1,k+1)-u4(m-1,k-1))+
         +p*c**2*(u3(m+1,k+1)+u3(m-1,k+1)-u3(m+1,k-1)-
         u3(m-1,k-1)+u3(k+2)-u3(k-2))+v*(u4(m+1,k+1)+u4(m-1,k+1)-
         u4(m+1,k-1)-u4(m-1,k-1)+u4(k+2)-u4(k-2)))/hx+ht*
         (2*w**2*(-u4(m+2)+2*u4(m)-u4(m-2))+4*w*p*c**2*(-u2(m+2)+
         2*u2(m)-u2(m-2))+2*c**2*(-u4(m+2)+2*u4(m)-u4(m-2))+
         4*p*v*c**2*(-u3(k+2)+2*u3(k)-u3(k-2))+2*c**2*(-u4(k+2)+
         2*u4(k)-u4(k-2))+2*v**2*(-u4(k+2)+2*u4(k)-u4(k-2))+
         4*p*v*c**2*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)-
         u2(m-1,k-1))+4*w*p*c**2*(-u3(m+1,k+1)+u3(m+1,k-1)+
         u3(m-1,k+1)-u3(m-1,k-1))+4*w*v*(-u4(m+1,k+1)+
         u4(m+1,k-1)+u4(m-1,k+1)-u4(m-1,k-1)))/hx/hx$
bb:=ampmat aa;
let
  sin(ax)=s1,
  cos(ax)=c1,
  sin(ay)=s2,
  cos(ay)=c2,
  w=k1*hx/ht,
  v=k2*hx/ht,
  c=k3*hx/ht,
  ht=r1*hx;
denotid a;
bb:=denotemat bb;
clear sin ax,cos ax,sin ay,cos ay,w,v,c,ht;
pol:=charpol bb;
denotid cp;
pol:=denotepol pol;
pol:=complexpol pol;
denotid rp;
pol:=denotepol pol;
prdenot;
cleardenot;
clear aa,bb,pol;

%***********************************************************************
%*****                                                             *****
%*****      T e s t     Examples   ---   Module    H U R W P       *****
%*****                                                             *****
%***********************************************************************

% Example H.1

x0:=lam-1;
x1:=lam-(ar+i*ai);
x2:=lam-(br+i*bi);
x3:=lam-(cr+i*ci);
hurwitzp x1;

% Example H.2

x:=hurw(x0*x1);
hurwitzp x;

% Example H.3

x:=(x1*x2);
hurwitzp x;

% Example H.4

x:=(x1*x2*x3);
hurwitzp x;
clear x,x0,x1,x2,x3;

%***********************************************************************
%*****                                                             *****
%*****    T e s t     Examples   ---   Module    L I N B A N D     *****
%*****                                                             *****
%***********************************************************************

on evallhseqp;  % So both sides of equations evaluate.

% Example L.1

operator v;
off echo;
gentran
<<literal
   tab!*,"DIMENSION U(200),V(200),ACOF(200,3),ARHS(200),XL(200,3)",cr!*$
  dx:=0.05$
  x:=0.1$
  for i:=1:101 do
    <<v(i):=x**2/2$
      x:=x+dx >> >>$
off period;
gentran
<<iacof:=200$
  iarhs:=200 >>$
on period;
genlinbandsol(1,1,{{u(1),u(1)=v(1)},{do,{k,2,100,1 },{u(k),u(k+1)-
   2*u(k)+u(k-1)=v(k+1)-2*v(k)+v(k-1)}},{u(101),u(101)=v(101)}})$
gentran
  <<amer:=0.0$
    arer:=0.0$
    for i:=1:101 do
      <<am:=abs(u(i)-v(i))$
        ar:=am/v(i)$
        literal tab!*,"IF(AM.GT.AMER) AMER=AM",cr!*$
        literal tab!*,"IF(AR.GT.ARER) ARER=AR",cr!* >>$
    literal tab!*,"WRITE(*,100)AMER,ARER",cr!*$
    literal tab!*,"STOP",cr!*$
    literal "100     FORMAT(' MAX. ABS. ERROR = ',E12.2,",
            "' MAX. REL. ERROR = ',E12.2)",cr!*$
    literal tab!*,"END",cr!* >>$
on echo;

%***********************************************************************

% Example L.2

on nag;
off echo;
gentran
<<literal
   tab!*,"DIMENSION U(200),V(200),ACOF(200,3),ARHS(200),XL(200,3)",cr!*$
  dx:=0.05$
  x:=0.1$
  for i:=1:101 do
    <<v(i):=x**2/2$
      x:=x+dx >> >>$
off period;
gentran
<<iacof:=200$
  iarhs:=200 >>$
on period;
genlinbandsol(1,1,{{u(1),u(1)=v(1)},{do,{k,2,100,1 },{u(k),u(k+1)-
   2*u(k)+u(k-1)=v(k+1)-2*v(k)+v(k-1)}},{u(101),u(101)=v(101)}})$
gentran
  <<amer:=0.0$
    arer:=0.0$
    for i:=1:101 do
      <<am:=abs(u(i)-v(i))$
        ar:=am/v(i)$
        literal tab!*,"IF(AM.GT.AMER) AMER=AM",cr!*$
        literal tab!*,"IF(AR.GT.ARER) ARER=AR",cr!* >>$
    literal tab!*,"WRITE(*,100)AMER,ARER",cr!*$
    literal tab!*,"STOP",cr!*$
    literal "100     FORMAT(' MAX. ABS. ERROR = ',E12.2,",
            "' MAX. REL. ERROR = ',E12.2)",cr!*$
    literal tab!*,"END",cr!* >>$
on echo;

%***********************************************************************

% Example L.3

on imsl;
off echo,nag;
gentran
<<literal
   tab!*,"DIMENSION U(200),V(200),ACOF(200,3),ARHS(200),XL(200,3)",cr!*$
  dx:=0.05$
  x:=0.1$
  for i:=1:101 do
    <<v(i):=x**2/2$
      x:=x+dx >> >>$
off period;
gentran
<<iacof:=200$
  iarhs:=200 >>$
on period;
genlinbandsol(1,1,{{u(1),u(1)=v(1)},{do,{k,2,100,1 },{u(k),u(k+1)-
   2*u(k)+u(k-1)=v(k+1)-2*v(k)+v(k-1)}},{u(101),u(101)=v(101)}})$
gentran
  <<amer:=0.0$
    arer:=0.0$
    for i:=1:101 do
      <<am:=abs(u(i)-v(i))$
        ar:=am/v(i)$
        literal tab!*,"IF(AM.GT.AMER) AMER=AM",cr!*$
        literal tab!*,"IF(AR.GT.ARER) ARER=AR",cr!* >>$
    literal tab!*,"WRITE(*,100)AMER,ARER",cr!*$
    literal tab!*,"STOP",cr!*$
    literal "100     FORMAT(' MAX. ABS. ERROR = ',E12.2,",
            "' MAX. REL. ERROR = ',E12.2)",cr!*$
    literal tab!*,"END",cr!* >>$
on echo;

%***********************************************************************

% Example L.4

on essl;
off echo,imsl;
gentran
<<literal
   tab!*,"DIMENSION U(200),V(200),ACOF(200,3),ARHS(200),XL(200,3)",cr!*$
  dx:=0.05$
  x:=0.1$
  for i:=1:101 do
    <<v(i):=x**2/2$
      x:=x+dx >> >>$
off period;
gentran
<<iacof:=200$
  iarhs:=200 >>$
on period;
genlinbandsol(1,1,{{u(1),u(1)=v(1)},{do,{k,2,100,1 },{u(k),u(k+1)-
   2*u(k)+u(k-1)=v(k+1)-2*v(k)+v(k-1)}},{u(101),u(101)=v(101)}})$
gentran
  <<amer:=0.0$
    arer:=0.0$
    for i:=1:101 do
      <<am:=abs(u(i)-v(i))$
        ar:=am/v(i)$
        literal tab!*,"IF(AM.GT.AMER) AMER=AM",cr!*$
        literal tab!*,"IF(AR.GT.ARER) ARER=AR",cr!* >>$
    literal tab!*,"WRITE(*,100)AMER,ARER",cr!*$
    literal tab!*,"STOP",cr!*$
    literal "100     FORMAT(' MAX. ABS. ERROR = ',E12.2,",
            "' MAX. REL. ERROR = ',E12.2)",cr!*$
    literal tab!*,"END",cr!* >>$
on echo;
off essl;

%***********************************************************************
%*****                                                             *****
%*****      T e s t     Complex Examples   ---   More Modules      *****
%*****                                                             *****
%***********************************************************************

% Example M.1

off exp;
coordinates t,x into n,j;
grid uniform,x,t;
dependence v(t,x),w(t,x);
%isgrid u(x..half),v(x..half);
%grideq u(x..one),v(x..one);
iim aa, v, diff(v,t)=c*diff(w,x),
        w, diff(w,t)=c*diff(v,x);
on exp;
center t=1/2;
functions v,w;
approx( aa(0,0)=aa(0,1));
center x=1/2;
approx( aa(1,0)=aa(1,1));
let cos ax**2=1-sin ax**2;
unfunc v,w;
matrix a(1,2),b(2,2),bt(2,2);
a(1,1):=aa(0,0);
a(1,2):=aa(1,0);
off prfourmat;
b:=ampmat a;
clear a,aa;
factor lam;
pol:=charpol b;
pol:=troot1 pol;
pol:=hurw num pol;
hurwitzp pol;
bt:=tcon b;
bt*b;
bt*b-b*bt;
clear aa,a,b,bt;

%***********************************************************************

% Example M.2 : Richtmyer, Morton: Difference methods for initial value
%               problems, &10.2.  p.261

coordinates t,x into n,j;
grid uniform,t,x;
let cos ax**2=1-sin ax**2;
unfunc v,w;
matrix a(1,2),b(2,2),bt(2,2);
a(1,1):=(v(n+1)-v(n))/ht-c*(w(j+1/2)-w(j-1/2))/hx$
a(1,2):=(w(n+1,j-1/2)-w(n,j-1/2))/ht-c*(v(n+1,j)-v(n+1,j-1))/hx$
off prfourmat;
b:=ampmat a;
clear a;
factor lam;
pol:=charpol b;
pol:=hurw num pol;
hurwitzp pol;
bt:=tcon b;
bt*b;
bt*b-b*bt;
clear a,b,bt;

%***********************************************************************

% Example M.3: Mazurik: Algoritmy resenia zadaci..., preprint no.24-85,
%               AN USSR SO, Inst. teor. i prikl. mechaniky, p.34

operator v1,v2;
matrix a(1,3),b(3,3),bt(3,3);
a(1,1):=(p(n+1)-p(n))/ht+c/2*((-p(m-1)+2*p(m)-p(m+1))/hx +
         (v1(m+1)-v1(m-1))/hx - (p(k-1)-2*p(k)+p(k+1))/hy +
         (v2(k+1)-v2(k-1))/hy)$
a(1,2):=(v1(n+1)-v1(n))/ht+c/2*((p(m+1)-p(m-1))/hx -
         (v1(m-1)-2*v1(m)+v1(m+1))/hx)$
a(1,3):=(v2(n+1)-v2(n))/ht + c/2*((p(k+1)-p(k-1))/hy -
         (v2(k-1)-2*v2(k)+v2(k+1))/hy)$
coordinates t,x,y into n,m,k;
functions p,v1,v2;
for k:=1:3 do approx(a(1,k)=0);
grid uniform,t,x,y;
unfunc p,v1,v2;
hy:=hx;
off prfourmat;
b:=ampmat a;
pol:=charpol b;
let
  cos ax=cos ax2**2-sin ax2**2,
  cos ay=cos ay2**2-sin ay2**2,
  sin ax=2*sin ax2*cos ax2,
  sin ay=2*sin ay2*cos ay2,
  cos ax2**2=1-sin ax2**2,
  cos ay2**2=1-sin ay2**2,
  sin ax2=s1,
  sin ay2=s2,
  hx=c*ht/cap;
factor lam;
order s1,s2;
pol:=troot1 pol;
clear cos ax,cos ay,sin ax,sin ay,cos ax2**2,cos ay2**2,sin ax2,sin ay2,
      hx,hy;
pol:=hurw num pol;
hurwitzp pol;
bt:=tcon b;
bt*b;
bt*b-b*bt;
clear a,b,bt,pol;

%***********************************************************************

end;
