reduce-historical/r35/lib/changevr.red

module changevar;   % Facility to perform CHANGE of independent VARs.

%*********************************************************************;
%                   -------------------------------                   ;
%                      C  H  A  N  G  E  V  A  R                      ;
%                   -------------------------------                   ;
%                                                                     ;
%   A REDUCE facility to perform CHANGE of independent VARiable(s)    ;
%   in differential equation (or a set of them).                      ;
%                                                                     ;
%                                                                     ;
%   Author : Gokturk Ucoluk                                           ;
%   Date   : Oct. 1989                                                ;
%   Place  : Middle East Tech. Univ., Physics Dept., Turkey.          ;
%   Email  : A07917 @ TRMETU.BITNET  or   A06794 @ TRMETU.BITNET      ;
%                                                                     ;
%   Version: 1.00                                                     ;
%                                                                     ;
%        ( *** Requires: REDUCE 3.0 or greater *** )                  ;
%        ( *** Requires: Matrix package to be present *** )           ;
%                                                                     ;
%   There exists a document written in LaTeX that explains the        ;
%   package in more detail.                                           ;
%                                                                     ;
%   Keywords : differential equation, change of variable, Jacobian    ;
%                                                                     ;

%*********************************************************************;

load!-package 'matrix;

fluid '(wtl!*);

global '(!*derexp !*dispjacobian);
switch derexp, dispjacobian ;  % on derexp       : Smart chain ruling
                               % on dispjacobian : Displays inverse
                               %                   Jacobian.
symbolic procedure simpchangevar v;
begin scalar u,f,j,dvar,flg;
     dvar := if pairp car v then cdar v else car v.nil;
     v := cdr v;                        % Dvar: list of depend. var
     u :=  if pairp car v then cdar v else car v.nil;
     v := cdr v;                        % u: list of new variables
     if eqcar(car v,'list) then v := append(cdar v,cdr v);
     while cdr v do
       << if caar v neq 'equal
            then rederr "improper new variable declaration";
          f := cdar v . f;  % f: list of entries (oldvar  func(newvrbs))
          v := cdr v >>;    %                           i     i
     v := reval car v;     % v: holds now last expression (maybe a list)
     if length u < length f then  rederr "Too few new variables"
      else if length u > length f then rederr "Too few old variables";
     % Now we form the Jacobian matrix ;
     j := for each entry in f collect
            for each newvrb in u collect
                reval list('df,cadr entry, newvrb);
     j := cdr aeval list('quotient,1,'mat.j);
          % j: holds inverse Jacobian.
     % We have to define the dependencies of old variables to new
     % variables.
     for each new in u do
       for each old in f do
          depend1(new, car old, t);
     % Below everything is perplexed :
     % The aim is to introduce  LET DF(new   ,old   ) = jacobian
     %                                   row    col            row,col
     % With the pairing trick below we do it in one step.
     %  new   : car row,   old   : caar col,   jacobian       : cdr col
     %     row                col                      row,col
     %
     for each row in pair(u,j) do
       for each col in pair(f,cdr row) do
        <<
          let2(list('df,car row,caar col), sqchk cdr col, nil, t);
          if !*dispjacobian and !*msg then
             mathprint list('equal,list('df,car row,caar col),
                            sqchk cdr col)
        >>;
     flg := !*derexp;
     !*derexp := t;

     v := changearg(dvar,u,v);

     for each  entry in f do
       v := subcare(car entry, cadr entry, v);
     % now here comes the striking point ... we evaluate the last
     % argument.
     v := simp!* v;
     % Now clean up the mess of LET;
     for each new in u do
       for each old in f do
          << let2(list('df,new,car old), nil, nil, nil);
             let2(list('df,new,car old), nil, t,   nil) >>;
     !*derexp := flg;
     return v;
end;

put('changevar,'simpfn,'simpchangevar);

symbolic procedure changearg(f,u,x);
if atom x then x
 else if car x memq f then car x . u
 else changearg(f,u,car x) . changearg(f,u,cdr x);


 symbolic procedure subcare(x,y,z);   % shall be used after changearg ;
 if null z then nil
  else if x = z then y
  else if atom z or get(car z,'subfunc) then z
  else (subcare(x,y,car z) . subcare(x,y,cdr z));


% Updated version of DIFFP.. chain rule handling is smarter. ;
% Example: If F is an operator and R has a implicit dependency on X,
%          declared by a preceding DEPEND R,X .. then the former version
%          of DIFFP, provided in REDUCE 3.3, was such that an algebraic
%          evaluation of    DF(F(R),X)    will evaluate to itself, that
%          means no change will happen. With the below given update this
%          is improved.   If the new provided flag DEREXP is OFF  then
%          the differentiation functions exactly like it was before,
%          but if DEREXP is ON then the chain rule is taken further to
%          yield the result:    DF(F(R),R)*DF(R,X)  .

remflag('(diffp),'lose);   % Since we want to reload it.

symbolic procedure diffp(u,v);
   %U is a standard power, V a kernel.
   % Value is the standard quotient derivative of U wrt V.
   begin scalar n,w,x,y,z,key; integer m;
        n := cdr u;     %integer power;
        u := car u;     %main variable;
        if u eq v and (w := 1 ./ 1) then go to e
         else if atom u then go to f
         %else if (x := assoc(u,dsubl!*)) and (x := atsoc(v,cdr x))
%               and (w := cdr x) then go to e   %deriv known;
             %DSUBL!* not used for now;
         else if (not atom car u and (w:= difff(u,v)))
                  or (car u eq '!*sq and (w:= diffsq(cadr u,v)))
          then go to c  %extended kernel found;
         else if (x:= get(car u,'dfn)) then nil
         else if car u eq 'plus and (w:=diffsq(simp u,v))
          then go to c
         else go to h;  %unknown derivative;
        y := x;
        z := cdr u;
    a:  w := diffsq(simp car z,v) . w;
        if caar w and null car y then go to h;  %unknown deriv;
        y := cdr y;
        z := cdr z;
        if z and y then go to a
         else if z or y then go to h;  %arguments do not match;
        y := reverse w;
        z := cdr u;
        w := nil ./ 1;
    b:  %computation of kernel derivative;
        if caar y
          then w := addsq(multsq(car y,simp subla(pair(caar x,z),
                                                   cdar x)),
                          w);
        x := cdr x;
        y := cdr y;
        if y then go to b;
    c:  %save calculated deriv in case it is used again;
        %if x := atsoc(u,dsubl!*) then go to d
        %else x := u . nil;
        %dsubl!* := x . dsubl!*;
  % d:   rplacd(x,xadd(v . w,cdr x,t));
    e:  %allowance for power;
        %first check to see if kernel has weight;
        if (x := atsoc(u,wtl!*))
          then w := multpq('k!* .** (-cdr x),w);
        m := n-1;
        % Evaluation is far more efficient if results are rationalized.
        return rationalizesq if n=1 then w
                else if flagp(dmode!*,'convert)
                     and null(n := int!-equiv!-chk
                                           apply1(get(dmode!*,'i2d),n))
                 then nil ./ 1
                else multsq(!*t2q((u .** m) .* n),w);
    f:  % Check for possible unused substitution rule.
        if not depends(u,v)
           and (not (x:= atsoc(u,powlis!*))
                 or not smember(v,simp cadddr x))
          then return nil ./ 1;
        w := list('df,u,v);
        go to j;
    h:  %final check for possible kernel deriv;
        y := nil;
        if car u eq 'df then key:=t;
        w := if key then 'df . cadr u . derad(v,cddr u)
              else list('df,u,v);
        y := cddr u;
        w := if (x := opmtch w) then simp x
              else if not depends(cadr w,caddr w) then nil ./ 1
              else if !*derexp then
                    begin
                        if atom cadr w then return mksq(w,1);
                        w := nil ./ 1;
                        for each m in cdr(if key then cadr u else u) do
                          w := addsq(multsq(
                                 if (x := opmtch (z :=
                                  'df . if key then (cadr u.derad(m,y))
                                         else list(u,m) )) then
                                  simp x else mksq(z,1),
                                  diffsq(simp m,v)),
                                  w);
                        return w
                    end
              else mksq(w,1);
        go to e;
    j:  w := if x := opmtch w then simp x else mksq(w,1);
        go to e
   end;

endmodule;

end;