reduce-historical/r38/packages/odesolve/odepatch.red

module odepatch$  % Patches to standard REDUCE facilities

% F.J.Wright@Maths.QMW.ac.uk, Time-stamp: <18 September 2000>

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Integrator patches
% ==================

% Avoid trying to integrate an integral to speed up some odesolve
% calls.

% NB: subst copies, even if no substitution is made.  It is therefore
% very likely to destroy uniqueness of kernels!

%% load_package int$
%% apply1('load!-package, 'int)$           % not at compile time!
packages_to_load int$                   % not at compile time!

global '(ODESolveOldSimpInt)$
(if not(s eq 'NoIntInt_SimpInt) then ODESolveOldSimpInt := s) where
   s = get('int, 'simpfn)$              % to allow reloading
put('int, 'simpfn, 'NoIntInt_SimpInt)$

fluid '(!*NoInt !*Odesolve_NoInt)$

symbolic procedure NoIntInt_SimpInt u;
   %% Patch to avoid trying to re-integrate symbolic integrals in the
   %% integrand, because it can take forever and achieve nothing!
   if !*NoInt then
   begin scalar v, varstack!*;
      % Based on int/driver/simpint1.
      % Varstack* rebound, since FORMLNR use can create recursive
      % evaluations.  (E.g., with int(cos(x)/x**2,x)).
      u := 'int . u;
      return if (v := formlnr u) neq u then <<
         v := simp subst('int!*, 'int, v);
         remakesf numr v ./ remakesf denr v
      >> else !*kk2q u
   end
   else
   if atom u or null cdr u or cddr u and (null cdddr u or cddddr u)
     then rerror(int,1,"Improper number of arguments to INT")
    else if cddr u then simpdint u      % header from simpint
   else begin scalar car_u, result;
      %% put('int, 'simpfn, 'SimpNoInt);
      car_u := mk!*sq simp!* car u;
      %% car_u := subeval{{'equal, 'Int, 'NoInt}, car_u};
      car_u := subst('NoInt, 'Int, car_u);
      u := car_u . !*a2k cadr u . nil;
      %% Prevent SimpInt from resetting itself!
      put('int, 'simpfn, ODESolveOldSimpInt);    % assumed & RESET by simpint
      result := errorset!*({ODESolveOldSimpInt, mkquote u}, t);
      put('int, 'simpfn, 'NoIntInt_SimpInt); % reset INT interface
      if errorp result then error1();
      return NoInt2Int car result;
      %% Does this cause non-unique kernels?
   end$

algebraic operator NoInt$               % Inert integration operator

%% symbolic procedure SimpNoInt u;
%%    !*kk2q('NoInt . u)$                % remain symbolic

symbolic operator Odesolve!-Int$
symbolic procedure Odesolve!-Int(y, x);
   %% Used in SolveLinear1 on ode1 to control integration.
   if !*Odesolve_NoInt then formlnr {'NoInt, y, x}
   else mk!*sq NoIntInt_SimpInt{y, x}$  % aeval{'int, y, x}$

%% put('Odesolve!-Int, 'simpfn, 'Simp!-Odesolve!-Int)$
%% symbolic procedure Simp!-Odesolve!-Int u;
%%    %% Used in SolveLinear1 on ode1 to control integration.
%%    if !*Odesolve_NoInt then !*kk2q('NoInt . u)  % must eval u!!!
%%    else NoIntInt_SimpInt u$         % aeval{'int, y, x}$

symbolic procedure NoInt2Int u;
   %% Convert all NoInt's back to Int's, without algebraic evaluation.
   subst('Int, 'NoInt, u)$

switch NoIntInt$  !*NoIntInt := t$
put('NoIntInt, 'simpfg,
   '((nil (put 'int 'simpfn 'SimpInt) (rmsubs))
     (t (put 'int 'simpfn 'NoIntInt_SimpInt))))$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Differentiator patches
% ======================

% Differentiate integrals correctly!

% NB: `ON' is flagged ignore and so not compiled, so...
on1 'allowdfint$

% To replace the versions in `reduce/packages/poly/diff.red' once
% tested.

% deflist('((dfint ((t (progn (on1 'allowdfint) (rmsubs)))))), 'simpfg);

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; 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,'dfform) then return apply3(x,u,v,n)
	 else if x:= get(car u,dfn_prop u) 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 depends(cadddr x,v))
	   and null !*depend
	  then return nil ./ 1;
        % Derivative of a dependent identifier via the chain rule.
        % Suppose u(v) = u(a(v),b(v),...), i.e. given depend {u}, a,
        % b, {a, b}, v; then (essentially) depl!* = ((b v) (a v) (u b
        % a))
        if !*expanddf and not(v memq (x:=cdr atsoc(u, depl!*))) then <<
           w := nil ./ 1;
           for each a in x do
              w := addsq(w, multsq(simp{'df,u,a},simp{'df,a,v}));
           go to e
        >>;
	w := list('df,u,v);
	w := if x := opmtch w then simp x else mksq(w,1);
	go to e;
    h:  % Final check for possible kernel deriv.
	if car u eq 'df then <<         % multiple derivative
           if cadr u eq v then
              % (df (df v x y z ...) v) ==> 0 if commutedf
              if !*commutedf and null !*depend then return nil ./ 1
              else if !*simpnoncomdf and (w:=atsoc(v, depl!*))
                 and null cddr w % and (cadr w eq (x:=caddr u))
              then
                 % (df (df v x) v) ==> (df v x 2)/(df v x) etc.
                 % if single independent variable
                 <<
                    x := caddr u;
                    % w := simp {'quotient, {'df,u,x}, {'df,v,x}};
                    w := quotsq(simp{'df,u,x},simp{'df,v,x});
                    go to e
                 >>
           else if eqcar(cadr u, 'int) then
              % (df (df (int F x) A) v) ==> (df (df (int F x) v) A) ?
              % Commute the derivatives to differentiate the integral?
              if caddr cadr u eq v then
                 % Evaluating (df u v) where u = (df (int F v) A)
                 % Just return (df F A) - derivative absorbed
                 << w := 'df . cadr cadr u . cddr u;  go to j >>
              else if !*allowdfint and
                 % Evaluating (df u v) where u = (df (int F x) A)
                 % (If dfint is also on then this will not arise!)
                 % Commute only if the result simplifies:
                 not_df_p(w := diffsq(simp!* cadr cadr u, v))
              then <<
                 % Generally must re-evaluate the integral (carefully!)
                 w := 'df . reval{'int, mk!*sq w, caddr cadr u} . cddr u;
                 go to j >>;  % derivative absorbed
	   if (x := find_sub_df(w:= cadr u . merge!-ind!-vars(u,v),
					   get('df,'kvalue)))
			  then <<w := simp car x;
				 for each el in cdr x do
				    for i := 1:cdr el do
					w := diffsq(w,car el);
				 go to e>>
		       else w := 'df . w
        >> else if !*expanddf and not atom cadr u then <<
           % Derivative of an algebraic operator u(a(v),...) via the
           % chain rule: df(u(v),v) = u_1(a(v),b(v),...)*df(a,v) + ...
           x := intern compress nconc(explode car u, '(!! !! !_));
           y := cdr u;  w := nil ./ 1;  m := 0;
           for each a in y do
           begin scalar b;
              m:=m#+1;
              if numr(b:=simp{'df,a,v}) then <<
                 z := mkid(x, m);
                 put(z, 'simpfn, 'simpiden);
                 w := addsq(w, multsq(simp(z . y), b))
              >>
           end;
           go to e
        >> else w := {'df,u,v};
   j:   if (x := opmtch w) then w := simp x
	 else if not depends(u,v) and null !*depend then return nil ./ 1
	 else w := mksq(w,1);
      go to e
   end$

symbolic procedure dfform_int(u, v, n);
   % Simplify a SINGLE derivative of an integral.
   % u = '(int y x) [as main variable of SQ form]
   % v = kernel
   % n = integer power
   % Return SQ form of df(u**n, v) = n*u**(n-1)*df(u, v)
   % This routine is called by diffp via the hook
   % "if x := get(car u,'dfform) then return apply3(x,u,v,n)".
   % It does not necessarily need to use this hook, but it needs to be
   % called as an alternative to diffp so that the linearity of
   % differentiation has already been applied.
   begin scalar result, x, y, dx!/dv;
      y := simp!* cadr u;  % SQ form integrand
      x := caddr u;  % kernel
      result :=
      if v eq x then y
	 % Special case -- just differentiate the integral:
	 % df(int(y,x), x) -> y  replacing the let rule in INT.RED
      else if not !*intflag!* and       % not in the integrator
	 % If used in the integrator it can cause infinite loops,
	 % e.g. in df(int(int(f,x),y),x) and df(int(int(f,x),y),y)
	 !*allowdfint and               % must be on for dfint to work
         <<
            % Compute PARTIAL df(y, v), where y must depend on x, so
            % if x depends on v, temporarily replace x:
            result := if numr(dx!/dv:=diffp(x.**1,v)) then
               %% (Subst OK because all kernels.)
               subst(x, xx, diffsq(subst(xx, x, y), v)) where
                  xx = gensym()
            else diffsq(y, v);
            !*dfint or not_df_p result
         >>
      then
	 % Differentiate under the integral sign:
         % df(int(y,x), v) -> df(x,v)*y + int(df(y,v), x)
         addsq(
            multsq(dx!/dv, y),
            simp{'int, mk!*sq result, x})  % MUST re-simplify it!!!
	    % (Perhaps I should use prepsq -
            % kernels are normally true prefix?)
      else !*kk2q{'df, u, v};  % remain unchanged
      if n neq 1 then
	 result := multsq( (((u .** (n-1)) .* n) .+ nil) ./ 1, result);
      return result
   end$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Solve patches
% =============

% Needed for definition of `mkrootsoftag' function.
packages_to_load solve$                 % not at compile time!

endmodule$
end$

%                ***** NOT INSTALLED AT PRESENT *****

%% An algebraic solver appropriate to ODESolve, that never returns
%% implicit solutions and returns nil when it fails.  Also, it solves
%% quadratics in terms of plus_or_minus.

%% *** NB: This messes up `root_multiplicities'. ***

%% (Later also use the root_of_unity operator where appropriate.
%% Could do this by a wrapper around `solvehipow' in solve/solv1.)

% This switch controls `solve' globally once this file has been
% loaded:
switch plus_or_minus$                   % off by default

%% The integrator does not handle integrands containing the
%% `plus_or_minus' operator at all well.  This may require some
%% re-writing of the integrator (!).  Temporarily, just turn off the
%% use of the `plus_or_minus' operator when necessary.

% This switch controls some `solve' calls locally within `ODESolve':
switch odesolve_plus_or_minus$          % TEMPORARY
% off by default -- it's to odangerous at present!
% !*odesolve_plus_or_minus := t$          % TEMPORARY

symbolic operator AlgSolve$
symbolic procedure AlgSolve(u, v);
   %% Return either a list of EXPLICIT solutions of a single scalar
   %% expression `u' for variable `v' or nil.
   begin scalar soln, tail, !*plus_or_minus;
      !*plus_or_minus := t;
      tail := cdr(soln := solveeval1{u, v});
      while tail do
         if car tail eq v then tail := cdr tail
         else tail := soln := nil;
      return soln
   end$

algebraic procedure SolvePM(u, v);
   %% Solve a single scalar expression `u' for variable `v', using the
   %% `plus_or_minus' operator in the solution of quadratics.
   %% *** NB: This messes up `root_multiplicities'. ***
   symbolic(solveeval1{u, v}
      where !*plus_or_minus := !*odesolve_plus_or_minus)$

%% This is a modified version of a routine in solve/quartic

symbolic procedure solvequadratic(a2,a1,a0);
   % A2, a1 and a0 are standard quotients.
   % Solves a2*x**2+a1*x+a0=0 for x.
   % Returns a list of standard quotient solutions.
   % Modified to use root_val to compute numeric roots.  SLK.
   if !*rounded and numcoef a0 and numcoef a1 and numcoef a2
      then for each z in cdr root_val list mkpolyexp2(a2,a1,a0)
         collect simp!* z else
   begin scalar d;
      d := sqrtq subtrsq(quotsqf(exptsq(a1,2),4),multsq(a2,a0));
      a1 := quotsqf(negsq a1,2);
      return
         if !*plus_or_minus then        % FJW
            list(subs2!* quotsq(addsq(a1,
               multsq(!*kk2q newplus_or_minus(),d)),a2))
                  %% Must uniquefy here until newplus_or_minus does it
                  %% for itself!
         else
            list(subs2!* quotsq(addsq(a1,d),a2),
               subs2!* quotsq(subtrsq(a1,d),a2))
   end$

endmodule$
end$