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

module odenon1$  % Special form nonlinear ODEs of order 1

% F.J.Wright@maths.qmw.ac.uk, Time-stamp: <11 August 2001>

% Original author: Malcolm MacCallum

% Basic layout is to test first for well-known special forms, namely:

%    separable
%    quasi-separable (separable after linear transformation)
%    (algebraically) homogeneous
%    quasi-homogeneous (homogeneous after linear transformation)
%    Bernoulli
%    Riccati
%    solvable for x or y, including Clairaut and Lagrange
%    exact  (FJW: general exact ode now handled elsewhere)

% and at a later stage of development to add more general methods,
% such as Prelle Singer

% Note that all cases of first order ODEs can be considered equivalent
% to searches for integrating factors.  (In particular Lie methods do
% not help here.)

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

%% algebraic procedure odesolve(ode, y, x);
%%    %% Temporary definition for test purposes.
%%    begin scalar !*precise, !*odesolve!-solvable!-xy, solution;
%%       ode := num !*eqn2a ode;           % returns ode as expression
%%       if (solution := ODESolve!-NonLinear1(ode, y, x)) then
%%          return solution
%%       else
%%          write "***** ODESolve cannot solve this ODE!"
%%    end$

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

global '(ODESolve_Before_Non1Grad_Hook ODESolve_After_Non1Grad_Hook)$

algebraic procedure ODESolve!-NonLinear1(ode, y, x);
   %% Top-level solver for non-linear first-order ODEs.
   begin scalar odecoeffs, gradient, solution, p, ode_p;
      traceode1 "Entering top-level non-linear first-order solver ...";
      %% traceode "This is a nonlinear ODE of order 1.";
      p := gensym();
      ode_p := num sub(df(y,x) = p, ode);
      %% If ode contains exponentials then the above sub can produce a
      %% denominator when there is none in ode!
      odecoeffs := coeff(ode_p, p);
      if length odecoeffs = 2 and not smember(p, odecoeffs)
      then <<                           % first DEGREE ODE
         gradient := -first odecoeffs / second odecoeffs;
         symbolic if (solution := or(
            ODESolve!-Run!-Hook(
               'ODESolve_Before_Non1Grad_Hook, {gradient, y, x}),
            ODESolve!-Separable(gradient, y, x),
            ODESolve!-QuasiSeparable(gradient, y, x),
            ODESolve!-Homogeneous(gradient, y, x),
            ODESolve!-QuasiHomog(gradient, y, x),
            ODESolve!-Bernoulli(gradient, y, x),
            ODESolve!-Riccati(gradient, y, x),
            ODESolve!-Run!-Hook(
               'ODESolve_After_Non1Grad_Hook, {gradient, y, x})))
         then return solution
      >>;
      %% If ode degree neq 1 or above solvers fail then ...

      %% A first-order ode is "solvable-for-y" if it can be put into
      %% the form y = f(x,p) where p = y'.  This form includes
      %% Clairaut and Lagrange equations as special cases.  The
      %% Lagrange form is y = xF(y') + G(y').  If F(p) = p then this
      %% is a Clairaut equation.  It is "solvable-for-x" if it can be
      %% put into the form x = f(y,p).
      if (solution := ODESolve!-Clairaut(ode, ode_p, p, y, x)) then
         return solution;
      %% Avoid infinite loops:
      symbolic if !*odesolve!-solvable!-xy then return;
      symbolic(!*odesolve!-solvable!-xy := t);
      %% "Solvable for y" includes Lagrange as a special case:
      symbolic return
         ODESolve!-Solvable!-y(ode_p, p, y, x) or
         ODESolve!-Solvable!-x(ode_p, p, y, x)
   end$

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

% Support routines

algebraic procedure ODENon!-Linear1(ode, y, x);
   %% Solve the linear ODE:  dy/dx = gradient(x,y) = P(x)*y + Q(x)
   %% Split into 2 procedures by FJW.
   begin scalar gradient;
      gradient := coeff(num ode, df(y,x));   % {low, high}
      gradient := -first gradient/second gradient;
      traceode!* "This is a first-order linear ODE solved by ";
      return if smember(y, gradient) then
      begin scalar P, Q;
         traceode "the integrating factor method.";
         P := lcof(num gradient,y)/den gradient;
         Q := gradient - P*y;
         return { y = ODENon!-Linear1PQ(P, Q, x) }
      end
      else <<
         traceode "quadrature.";
         % FJW: Optionally turn off final integration:
         { y = ODESolve!-Int(gradient, x) + newarbconst() }
      >>
   end$

algebraic procedure ODENon!-Linear1PQ(P, Q, x);
   %% Solve the linear ODE:  dy/dx = P(x)*y + Q(x)
   %% Called directly by ODESolve!-Bernoulli
   begin scalar intfactor, !*combinelogs;
      %% intfactor simplifies better if logs in the integral are
      %% combined:
      symbolic(!*combinelogs := t);
      intfactor := exp(int(-P, x));
      %% Optionally turn off final integration:
      return (newarbconst() + ODESolve!-Int(intfactor*Q,x))/intfactor
   end$

%% algebraic procedure unfactorize factorlist;
%%    %% Multiply out a factor list of the form
%%    %% { {factor, multiplicity}, ... }:
%%    for each fac in factorlist product first fac^second fac$

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

% Separable ODEs

algebraic procedure ODESolve!-Separable(gradient, y, x);
   %% The ODE has the form dy/dx = gradient = F(x,y).  If F(x,y) =
   %% f(x)g(y) then the ODE is separable, in which case return the
   %% solution; otherwise return nil.
   begin scalar f, g, !*redefmsg;
      traceode1 "Testing for a separable ODE ...";

      %% Handle implicit dependence on x (ignoring that via y):
      symbolic (<< depend1(y, x, nil);
         %% Hack sub to handle implicit dependences:
         copyd('ODESolve!-old!-subsublis, 'subsublis);
         copyd('subsublis, 'ODESolve!-subsublis);
         g := errorset!*(
            {'ODESolve!-Separable1, mkquote gradient, mkquote x}, nil);
         copyd('subsublis, 'ODESolve!-old!-subsublis);
         if errorp g then RedErr {"(in ODESolve!-Separable1)", emsg!*};
         g := car g;
         if depends(g, x) then g := nil;
      >> where depl!* = depl!*);
      if not g then return;

      %% Then F(alpha,y) = f(alpha)g(y), and F(x,y)/F(alpha,y) =
      %% f(x)/f(alpha) is free of y:
      if depends(f := gradient/g, y) then return;
      traceode "It is separable.";
      %% Any separation of F(x,y) will suffice!
      %% gradient := int(1/g, y) - int(f, x) + newarbconst();
      %% Try to optimize structure of solution to avoid a likely
      %% logarithm of a function of y:
      gradient := int(1/g, y);
      gradient := if part(gradient,0) = log then
         part(gradient,1) - newarbconst()*exp(int(f, x))
      else gradient - int(f, x) + newarbconst();
      return { num gradient = 0 }
   end$

algebraic procedure ODESolve!-Separable1(gradient, x);
   %% Find a small constant alpha such that F(alpha,y) exists and
   %% F(alpha,y) neq 0, and return F(alpha,y), where F = gradient:
   begin scalar numer, denom, alpha, d, n;
      numer := num gradient;
      denom := den gradient;
      alpha := 0;
      while (d:=sub(x=alpha,denom)) = 0 or
         (n:=sub(x=alpha,numer)) = 0 do
            alpha := alpha + 1;
      return n/d
   end$

symbolic procedure ODESolve!-subsublis(u,v);
   % NOTE: This definition assumes that with the exception of *SQ and
   % domain elements, expressions do not contain dotted pairs.

   % This is the standard `subsublis' plus support for one level of
   % indirect dependence (cf. freeof), in which case sub converts the
   % dependent variable into an independent variable with a different
   % but related name.

   % u is an alist of substitutions, v is an expression to be
   % substituted:
   begin scalar x;
      return if x := assoc(v,u) then cdr x
	      % allow for case of sub(sqrt 2=s2,atan sqrt 2).
	      else if eqcar(v,'sqrt)
		 and (x := assoc(list('expt,cadr v,'(quotient 1 2)),u))
	       then cdr x
	      else if atom v then
                 if (x := assoc(v,depl!*)) and % FJW
                    %% Run through variables on which v depends:
                    << while (x := cdr x) and not assoc(car x,u) do; x >> % FJW
                 then mkid(v, '!!) else % FJW
                 v
	      else if not idp car v
	       then for each j in v collect ODESolve!-subsublis(u,j)
	      else if x := get(car v,'subfunc) then apply2(x,u,v)
	      else if get(car v,'dname) then v
	      else if car v eq '!*sq then ODESolve!-subsublis(u,prepsq cadr v)
	      else for each j in v collect ODESolve!-subsublis(u,j)
   end$

algebraic procedure ODESolve!-QuasiSeparable(gradient, y, x);
   %% The ODE has the form dy/dx = gradient = F(x,y).  If F(x,y) =
   %% f(y+kx) then the ODE is quasi-separable, in which case return
   %% the solution; otherwise return nil.
   begin scalar k;
      traceode1 "Testing for a quasi-separable ODE ...";
      %% F(x,y) = f(y+kx) iff df(F,x)/df(F,y) = k = constant (using
      %% PARTIAL derivatives):
      k := (df(gradient,x)/df(gradient,y) where df(y,x) => 0);
      if depends(k, x) then return;     % sufficient since y = y(x)
      traceode "It is separable after letting ", y+k*x => y;
      %% Setting u = y+kx gives du/dx = dy/dx+k = f(u)+k:
      gradient := sub(x=0, gradient) + k;
      %% => int(1/(f(u)+k),u) = x + arbconst.
      gradient := sub(y=y+k*x, int(1/gradient, y)) - x + newarbconst();
      return { num gradient = 0 }
   end$

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

% Algebraically homogeneous ODEs

algebraic procedure ODESolve!-Homogeneous(gradient, y, x);
   %% The ODE has the form dy/dx = gradient = F(x,y).  If F(x,y) =
   %% f(y/x) then the ODE is algebraically homogeneous.
   %% Setting y = vx => v + x dv/dx = F(x,vx) = f(v)
   begin scalar v;  v := gensym();
      traceode1 "Testing for a homogeneous ODE ...";
      gradient := sub(y=v*x, gradient);
      if depends(gradient, x) then return;
      traceode
         "It is of algebraically homogeneous type ",
         "solved by a change of variables of the form `y = vx'.";
      %% Integrate to give  int(1/(f(v)-v),v) = int(1/x, x) + arbconst
      %% Hence  exp int(1/(f(v)-v),v) = arbconst * x
      gradient := exp(int(1/(gradient-v), v));
      gradient := (gradient where (sqrt ~a)*(sqrt ~b) => sqrt(a*b));
      gradient := sub(v=y/x, x/gradient) + newarbconst();
      return {num gradient = 0}
   end$

%% A quasi-homogeneous first-order nonlinear ODE has the form

%% dy/dx = F(..., ((a1*x + b1*y + c1)/(a2*x + b2*y + c2))^p, ...)

%% where F is an arbitrary composition of functions and p is an
%% arbitrary power.  F may have other arguments that do not depend on
%% x or y or that depend in the same way (e.g. F may be a sum or
%% product).  The argument of F that depends on x or y will be a
%% quotient of expanded polynomials if p is a positive integer.
%% Otherwise, it will be a power of a quotient (expt (quotient ... )),
%% in which case the `expt' can be treated as part of the composition
%% F, or a quotient of powers (quotient (expt ... ) (expt ... )) which
%% must be treated separately.

algebraic procedure ODESolve!-QuasiHomog(gradient, y, x);
   %% The ODE has the form dy/dx = gradient = F(x,y).  If F(x,y) =
   %% f((a1*x + b1*y + c1)/(a2*x + b2*y + c2)) where the function f
   %% may be arbitrary then the ODE is reducible to algebraically
   %% homogeneous.  Find the first linear-rational sub-expression, and
   %% use it to try to make the ODE homogeneous:
   begin scalar tmp, n, d, soln;
      traceode1 "Testing for a quasi-homogeneous ODE ...";
      %% First, find an "argument" that is a rational function, with
      %% numerator and denominator that both depend on x.
      if not(tmp := symbolic ODESolve!-QuasiHomog1(reval gradient, x))
      then return;
      n := num tmp;  d := den tmp;
      %% Now check that numerator and denominator have the same degree
      %% in x (and y):
      if (tmp := deg(n,x)) neq deg(d,x) then return;
      %% If that degree > 1 then extract the squarefree factor:
      if tmp = 1 then                % => deg(d,x) = 1
         %% ( if deg(n,y) neq 1 or deg(d,y) neq 1 then return )
      else <<
         %% Use partial squarefree factorization to find p'th root of
         %% numerator and denominator:
         %%   If f = poly = (a x + b y + c)^p
         %%   then f' = a p(a x + b y + c)^(p-1)
         %%   => gcd(f, f') = (a x + b y + c)^(p-1)
         %%   => f / gcd(f, f') = a x + b y + c.
         n := n / gcd(n, df(n, x));     % must be linear in x and y
         %% if deg(n,x) neq 1 or deg(n,y) neq 1 then return;
         d := d / gcd(d, df(d, x));     % must be linear in x and y
         %% if deg(d,x) neq 1 or deg(d,y) neq 1 then return
      >>;
      %% Check that numerator and denominator are really LINEAR
      %% POLYNOMIALS in x and y (where y depends on x):
      if length(tmp := coeff(n, y)) neq 2 or depends(tmp, y) then return;
      if depends(second tmp, x) or      % coeff of y^1
         length(tmp := coeff(first tmp, x)) neq 2 or % coeff of y^0
         depends(tmp, x) then return;
      if length(tmp := coeff(d, y)) neq 2 or depends(tmp, y) then return;
      if depends(second tmp, x) or      % coeff of y^1
         length(tmp := coeff(first tmp, x)) neq 2 or % coeff of y^0
         depends(tmp, x) then return;
      %% The degenerate case a1*b2=a2*b1 is now treated as quasi-separable.
      traceode "It is quasi-homogeneous if ",
         "the result of shifting the origin is homogeneous ...";
      soln := first solve({n,d},{x,y});
      n := rhs first soln;  d := rhs second soln;
      gradient := sub(x=x+n, y=y+d, gradient);
      %% ODE was quasi-homogeneous iff the new ODE is homogeneous:
      if (soln := ODESolve!-Homogeneous(gradient,y,x)) then
         return sub(x=x-n, y=y-d, soln);
      traceode "... which it is not!"
   end$

%%% The calls to `depends' below are inefficient!

symbolic procedure ODESolve!-QuasiHomog1(u, x);
   %% Assumes "algebraic" form!  Get the first argument of any
   %% composition of functions that is a quotient of polynomials or
   %% symbolic powers (expt forms) that both depend on `x'.
   if atom u then nil                   % not QH, else operator ...
   else if car u eq 'quotient and       % quotient, in which
      depends(cadr u, x) and            % numerator depends on x, and
      depends(caddr u, x) then          % denominator depends on x.
         if eqcar(cadr u, 'expt) then   % symbolic powers
            ( if eqcar(caddr u, 'expt) and (caddr cadr u eq caddr caddr u)
            then {'quotient, cadr cadr u, cadr caddr u} )
         else
            ( if eqcar(cadr u, 'plus) and eqcar(caddr u, 'plus)
            then u )                    % polynomials (?)
   else  % Process first x-dependent argument of operator u:
   begin
   a: if (u := cdr u) then
         if depends(car u, x) then return ODESolve!-QuasiHomog1(car u, x)
         else go to a
   end$

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

% Bernoulli ODEs

symbolic operator ODESolve!-Bernoulli$

symbolic procedure ODESolve!-Bernoulli(rhs, y, x);
   %% The ODE has the form df(y,x) = rhs.  If rhs has the Bernoulli
   %% form P(x)*y + Q(x)*y^n then extract P(x), Q(x), n and return the
   %% solution else return nil.
   ( begin scalar num_rhs, den_rhs, C1, C2, d, d1, d2, d3, P, Q, n;
      traceode1 "Testing for a Bernoulli ODE ...";
      %% Degrees will be constructed in true prefix form.
      %% Need sum of two terms, both with main var (essentially) y.
      depend1(x, y, nil) where !*msg = nil;
      << updkorder y;  rhs := simp rhs >> where kord!* = kord!*;
      num_rhs := numr rhs;  den_rhs := denr rhs;
      if domainp num_rhs or not(mvar num_rhs eq y) then return;

      %% Num may have the form y^d * (y^d1 C1(x) + y^d2 C2(x)) ...
      if null red num_rhs then <<
         d := ldeg num_rhs;  num_rhs := lc num_rhs;
         if domainp num_rhs then return
      >>;

      %% Now num must have the form y^d1 C1(x) + y^d2 C2(x),
      %% where d1 > d2 and d2 = 0 is allowed (if d <> 0 or d3 <> 0).
      if (C1 := get!!y!^n!*C(num_rhs, y)) then
         << d1 := car C1;  C1 := cdr C1 >>
      else return;
      num_rhs := red num_rhs;
      %% Allow d2 = 0 => num_rhs freeof y
      if not smember(y, num_rhs) then
         << d2 := 0;  C2 := num_rhs >>
      else if red num_rhs then return
      else if (C2 := get!!y!^n!*C(num_rhs, y)) then
         << d2 := car C2;  C2 := cdr C2 >>
      else return;

      %% Den must have the form C3(x) or y^d3 C3(x).
      %% In the latter case, combine the powers of y.
      if smember(y, den_rhs) then
         if null red den_rhs and
            (den_rhs := get!!y!^n!*C(den_rhs, y)) then <<
            d3 := car den_rhs;  den_rhs := cdr den_rhs;
            d1 := {'difference, d1, d3};
            d2 := {'difference, d2, d3}
         >> else return;

      %% Simplify the degrees of y and find which is 1:
      if d then << d1 := {'plus, d1, d};  d2 := {'plus, d2, d} >>;
      d1 := aeval d1;  d2 := aeval d2;
      if d1 = 1 then << P := C1; Q := C2; n := d2 >>
      else if d2 = 1 then << P := C2; Q := C1; n := d1 >>
      else return;
      %% A final check that P, Q, n are valid:
      if Bernoulli!-depend!-check(P, y) or
         Bernoulli!-depend!-check(Q, y) or
            Bernoulli!-depend!-check(den_rhs, y) or
               Bernoulli!-depend!-check(n, x) then return;
      %% ( Last test implies Bernoulli!-depend!-check(n, y). )
      %% For testing:
      %% return {'list, mk!*sq(P ./ den_rhs), mk!*sq(Q ./ den_rhs), n};
      P := mk!*sq(P ./ den_rhs);  Q := mk!*sq(Q ./ den_rhs);
      return ODESolve!-Bernoulli1(P, Q, y, x, n)
   end ) where depl!* = depl!*$

symbolic procedure Bernoulli!-depend!-check(f, xy);
   %% f is a standard form, xy is an identifier (kernel).
   if numr difff(f, xy) then <<         % neq 0 (nil)
      traceode "It is not of Bernoulli type because ...";
      MsgPri(nil, !*f2a f, "depends (possibly implicitly) on",
         get(xy, 'odesolve!-depvar) or xy, nil);
      %% (y might be gensym -- 'odesolve!-depvar set in odeintfc)
      t
   >>$

symbolic procedure get!!y!^n!*C(u, y);
   %% Assume that u is a standard form representation of y^n * C(x)
   %% with y the leading kernel.  Return (n . C) or nil
   begin scalar n, C;
      if mvar u eq y then <<            % ?? y^n * C(x), n nonnegint
         n := ldeg u;  C := lc u;
         return if not domainp C and smember(y, mvar C) then
            ( if (C := get!!y!^n!*C1(C, y)) then
               % y^n * y^n1 * C(x), n nonnegint
               {'plus, n, car C} . cdr C )
         else n . C
      >> else                           % (y^n1)^n * C(x), n nonnegint
         return get!!y!^n!*C1(u, y)
   end$

symbolic procedure get!!y!^n!*C1(u, y);
   % u = (y^n1)^n * C(x), n nonnegint
   begin scalar n, C;
      n := mvar u;
      if not(eqcar(n, 'expt) and cadr n eq y) then return;
      n := {'times, caddr n, ldeg u};
      C := lc u;
      return n . C
   end$

algebraic procedure ODESolve!-Bernoulli1(P, Q, y, x, n);
   begin scalar !*odesolve_noint;  % Force integration?
      traceode "It is of Bernoulli type.";
      n := 1 - n;
      return if symbolic !*odesolve_explicit then
         { y = ODENon!-Linear1PQ(n*P, n*Q, x)^(1/n)*
               newroot_of_unity(n) }    % explicit form
      else { y^n = ODENon!-Linear1PQ(n*P, n*Q, x) } % implicit form
   end$

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

% Riccati ODEs

algebraic procedure ODESolve!-Riccati(rhs, y, x);
   %% The ODE has the form df(y,x) = rhs.  If rhs has the Riccati form
   %% a(x)y^2 + b(x)y + c(x) then transform to a reduced linear
   %% second-order ODE and attempt to solve it.
   symbolic if not !*odesolve_fast then % heuristic solution
   algebraic begin scalar a, b, c, soln, !*ratarg;
      traceode1 "Testing for a Riccati ODE ...";
      %% rhs may have a denominator that depends on y, so ...
      symbolic(!*ratarg := t);
      c := coeff(rhs, y);
      if length c neq 3 or depends(c, y) then return;
      a := third c;  b := second c;  c := first c;
      c := a*c;  b := -(df(a,x)/a + b);
      traceode "It is of Riccati type ",
         "and transforms into the linear second-order ODE: ",
         num(df(y,x,2) + b*df(y,x) + c*y) = 0;
      soln := {c, b, 1};                % low .. high
      soln := ODESolve!-linear!-basis(soln, 0, y, x, 2,
         if c then 0 else if b then 1 else 2); % min_order
      if not soln then <<
         traceode "But ODESolve cannot solve it!";
         return
      >>;
      %% The solution of the linear second-order ode must have the
      %% form y = arbconst(1)*y1 + arbconst(2)*y2, in which only the
      %% ratio of the arbconst's is relevant here, so ...
      soln := first soln;
      soln := newarbconst()*first soln + second soln;
      return {y = sub(y = soln, -df(y,x)/(a*y))}
      %% BEWARE: above minus sign missing in Zwillinger, first edn.
   end$

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

% ODEs that are "solvable for y or x", including Clairaut and Lagrange
% as special cases.

algebraic procedure ODESolve!-Clairaut(ode, ode_p, p, y, x);
   %% Assuming that ode is first order, determine whether it is of
   %% Clairaut type f(xy'-y) = g(y'), and if so return its general
   %% solution together with a singular solution if one exists.
   begin scalar sing_soln;
      traceode1 "Testing for a Clairaut ODE ...";
      sing_soln := sub(y = x*p, ode_p);
      if depends(sing_soln, x) or depends(sing_soln, y) then
         return;                        % Not Clairaut
      traceode "It is of Clairaut type.";
      %% Look for a singular solution:
      sing_soln := solve(df(ode, x), df(y,x));
      sing_soln := for each soln in sing_soln join
         %% NB: Cannot use algebraic code to check soln because it may
         %% contain derivatives that evaluate to 0.
         if symbolic eqcar(caddr soln, 'root_of)
         then {} else {num sub(soln, ode) = 0};
      return (sub(p = newarbconst(), ode_p) = 0) . sing_soln
   end$

symbolic operator ODESolve!-Solvable!-y$
symbolic procedure ODESolve!-Solvable!-y(ode_p, p, y, x);
   %% ode_p has the form f(x,y,p) where p = y'.
   ( algebraic begin scalar c, lagrange, ode_y, ode_x;
      traceode1 "Testing for a ""solvable for y"" or Lagrange ODE ...";
      %% Is the ODE "solvable for y"?
      c := coeff(ode_p, y);
      if length c neq 2 or depends(c, y) then return;
      ode_y := -first c / second c;
      %% ode_y has the form f(x,p) where p = y' and ode is y = ode_y.
      %% If f(x,p) = xF(p) + G(p) then ode is a Lagrange equation.
      if not depends(den ode_y, x) and
         length(c:=coeff(num ode_y, x)) = 2 and
         not depends(c, x) then
         lagrange := 1;                 % Lagrange form
      symbolic depend1(p, x, t);
      ode_x := num(p - df(ode_y, x));
      if lagrange then <<
         %% ODE is a Lagrange equation, hence ode_x is a LINEAR ODE
         %% for x(p) that can be solved EXPLICITLY (using an
         %% integrating factor) for a single value of x.
         symbolic depend1(x, p, t);
         ode_x := num(ode_x where df(p,x) => 1/df(x,p));
         symbolic depend1(p, x, nil);
         traceode "It is of Lagrange type and reduces to this ",
            "subsidiary ODE for x(y'): ",
            ode_x = 0;
         ode_x := ODENon!-Linear1(ode_x, x, p)
      >> else if symbolic !*odesolve_fast then return
         traceode "Sub-solver terminated: fast mode, no heuristics!"
      else <<
         %% ode_x is an arbitrary first-order ODE for p(x), so ...
         traceode "It is ""solvable for y"" and reduces ",
            "to this subsidiary ODE for y'(x):";
         ode_x := ODESolve!-FirstOrder(ode_x, p, x);
         if not ode_x then <<
            traceode "But ODESolve cannot solve it!";
            return
         >>
      >>;
      traceode "The subsidiary solution is ", ode_x,
         " and the main ODE can be solved parametrically ",
         "in terms of the derivative.";
      return if symbolic(!*odesolve_implicit or !*odesolve_explicit) then
         %% Try to eliminate p between ode_y and ode_x else fail.
         %% Assume that the interface code will try to actually solve
         %% this for y:
         Odesolve!-Elim!-Param(ode_y, y, ode_x, p, y)
      else
         %% Return a parametric solution {y(p), x(p), p}:
         if lagrange then            % soln explicit for x
            for each soln in ode_x collect ODESolve!-Simp!-ArbParam
               sub(p=newarbparam(), {y = sub(soln, ode_y), soln, p})
         else
            for each soln in ode_x join <<
               %% Make solution as explicit as possible:
               soln := solve(soln, x);
               for each s in soln collect ODESolve!-Simp!-ArbParam
                  sub(p=newarbparam(),
                     if symbolic eqcar(caddr s, 'root_of) then
                        {y=ode_y, sub(part(rhs s, 2)=x, part(rhs s, 1)), p}
                     else {y=sub(s, ode_y), s, p})
            >>
   end ) where depl!* = depl!*$


symbolic operator ODESolve!-Solvable!-x$
symbolic procedure ODESolve!-Solvable!-x(ode_p, p, y, x);
   %% ode_p has the form f(x,y,p) where p = y'.
   not !*odesolve_fast and              % heuristic solution
   ( algebraic begin scalar c, ode_x, ode_y;
      traceode1 "Testing for a ""solvable for x"" ODE ...";
      %% Is the ODE "solvable for x"?
      c := coeff(ode_p, x);
      if length c neq 2 or depends(c, x) then return;
      ode_x := -first c / second c;
      %% ode_x has the form f(y,p) where p = y' and ode is x = ode_x.
      symbolic depend1(p, y, t);
      ode_y := num(1/p - df(ode_x, y));
      %% ode_y is an arbitrary first-order ODE for p(y), so ...
      traceode "It is ""solvable for x"" and reduces ",
         "to this subsidiary ODE for y'(y):";
      ode_y := ODESolve!-FirstOrder(ode_y, p, y);
      if not ode_y then <<
         traceode "But ODESolve cannot solve it! ";
         return
      >>;
      traceode "The subsidiary solution is ", ode_y,
         " and the main ODE can be solved parametrically ",
         "in terms of the derivative.";
      return if symbolic(!*odesolve_implicit or !*odesolve_explicit) then
         %% Try to eliminate p between ode_x and ode_y else fail.
         %% Assume that the interface code will try to actually solve
         %% this for y:
         Odesolve!-Elim!-Param(ode_x, x, ode_y, p, y)
      else
         for each soln in ode_y join <<
            %% Return a parametric solution {y(p), x(p), p}:
            %% Make solution as explicit as possible:
            soln := solve(soln, y);
            for each s in soln collect ODESolve!-Simp!-ArbParam
               sub(p=newarbparam(),
                  if symbolic eqcar(caddr s, 'root_of) then
                     {sub(part(rhs s, 2)=y, part(rhs s, 1)), x=ode_x, p}
                  else {s, x=sub(s, ode_x), p})
         >>
   end ) where depl!* = depl!*$

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

% The arbitrary parameters in solutions:

algebraic operator arbparam$
algebraic (!!arbparam := 0)$
algebraic procedure newarbparam();
   arbparam(!!arbparam:=!!arbparam+1)$

algebraic procedure Odesolve!-Elim!-Param(ode_y, y, ode_x, p, depvar);
   %% ode_y is an expression for y and ode_x is an rlist of odesolve
   %% solutions for x, both in terms of a parameter p.  Return a list
   %% of equations corresponding to the equations in ode_x but with p
   %% eliminated with ode_y, or nil if this is not possible.
   %% depvar is the true dependent variable.
   begin scalar soln, result;
      if polynomialp(ode_y := num(y - ode_y), p) then <<
         result := {};
         while ode_x neq {} and
            polynomialp(soln := num !*eqn2a first ode_x, p)
            do <<
               result := resultant(ode_y, soln, p) . result;
               ode_x := rest ode_x
            >>;
         if ode_x = {} then
            return Odesolve!-Tidy!-Implicit(result, y)
      >>;
      ode_y := solve(ode_y, p);
      %% solve here may return a one_of construct (zimmer (4) & (19)).
      %% I don't see why, but ...
      %% (expand_cases not defined until solve loaded.)
      ode_y := expand_cases ode_y;
      if not smember(root_of, ode_y) then <<
         result := for each soln in ode_y join
            for each s in ode_x join
               if rhs s = 0 then  % s is f(x,y) = 0
                  if (s:=sub(soln, num lhs s)) neq 0 then {num s}
                  else {}
               else {num(sub(soln, rhs s) - x)};  % s is x = f(x,y)
         return Odesolve!-Tidy!-Implicit(result, depvar)
      >>;
      traceode "But cannot eliminate parameter ",
         "to make solution explicit."
   end$

algebraic procedure Odesolve!-Tidy!-Implicit(solns, depvar);
   %% Remove repeated and irrelevant factors from implicit solutions.
   for each soln in solns join
      for each fac in factorize soln join
         if smember(depvar, fac:=first fac) then {fac = 0} else {}$

switch odesolve_simp_arbparam$          % DEFAULT OFF.  TEMPORARY?

symbolic operator ODESolve!-Simp!-ArbParam$
symbolic procedure ODESolve!-Simp!-ArbParam u;
   %% Simplify arbparam expressions within parametric solution u
   %% (cf. ODESolve!-Simp!-ArbConsts)
   begin scalar !*precise, x, y, ss_x, ss_y, arbexprns_x, arbexprns_y,
      arbexprns, param;
      if not(rlistp u and length u = 4) then
         TypErr(u, "parametric ODE solution");
      if not !*odesolve_simp_arbparam then return u; % TEMPORARY?
      x := lhs cadr u;  y := lhs caddr u;
      if not(ss_x := ODESolve!-Structr(caddr cadr u, x, y, 'arbparam))
      then return u;
      if not(ss_y := ODESolve!-Structr(caddr caddr u, x, y, 'arbparam))
      then return u;
      ss_x := cdr ss_x;  ss_y := cdr ss_y;
      arbexprns_x := for each s in cdr ss_x collect caddr s;
      arbexprns_y := for each s in cdr ss_y collect caddr s;
      if null(arbexprns := intersection(arbexprns_x, arbexprns_y))
      then return u;
      arbexprns_x := cdr ss_x;  ss_x := car ss_x;
      arbexprns_y := cdr ss_y;  ss_y := car ss_y;
      arbexprns_x := for each s in arbexprns_x join
         if member(caddr s, arbexprns) then {s}
         else << ss_x := subeval{s, ss_x}; nil >>;
      arbexprns_y := for each s in arbexprns_y join
         if member(caddr s, arbexprns) then {s}
         else << ss_y := subeval{s, ss_y}; nil >>;
      traceode "Simplifying the arbparam expressions in ", u,
         " by the rewrites ...";
      param := cadddr u;
      for each s in arbexprns_x do <<
         %% s has the form ansj = "expression in arbparam(n)"
         traceode rhs s => param;
         %% Remove other occurrences of arbparam(n):
         ss_x := algebraic sub(solve(s, param), ss_x);
         %% Finally rename ansj as arbparam(n):
         ss_x := subeval{{'equal, cadr s, param}, ss_x}
      >>;
      for each s in arbexprns_y do <<
         ss_y := algebraic sub(solve(s, param), ss_y);
         ss_y := subeval{{'equal, cadr s, param}, ss_y}
      >>;
      %% Try a further heuristic simplification:
      if smember(param, den ss_x) and smember(param, den ss_y) then
      begin scalar ss_x1, ss_y1;
         ss_x1 := algebraic sub(param = 1/param, ss_x);
         if smember(param, den ss_x1) then return;
         ss_y1 := algebraic sub(param = 1/param, ss_y);
         if smember(param, den ss_y1) then return;
         traceode "Simplifying further by the rewrite ",
            param => 1/param;
         ss_x := ss_x1;  ss_y := ss_y1
      end;
      return makelist {{'equal, x, ss_x}, {'equal, y, ss_y}, param}
   end$

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

symbolic operator Polynomialp$

%% symbolic procedure polynomialp(pol, var);
%%    %% Returns true if numerator of pol is polynomial in var.
%%    %% Assumes on exp, mcd.
%%    begin scalar kord!*;  kord!* := {var};
%%       pol := numr simp!* pol;
%%       while not domainp pol and mvar pol eq var and
%%          (domainp lc pol or not smember(var, mvar lc pol)) do
%%             pol := red pol;
%%       return not smember(var, pol)
%%    end$

symbolic procedure Polynomialp(pol, var);
   %% Returns true if numerator of pol is polynomial in var.
   %% Assumes on exp, mcd.
   Polynomial!-Form!-p(numr simp!* pol, !*a2k var)$

symbolic procedure Polynomial!-Form!-p(sf, y);
   %% A standard form `sf' is polynomial if each of its terms is
   %% polynomial:
   domainp sf or
      (Polynomial!-Term!-p(lt sf, y) and Polynomial!-Form!-p(red sf, y))$

symbolic procedure Polynomial!-Term!-p(st, y);
   %% A standard term `st' is polynomial if either (a) its
   %% leading power is polynomial and its coefficient is free of y, or (b)
   %% its leading power is free of y and its coefficient is polynomial:
   if tvar st eq y then not smember(y, tc st)
   else if not smember(y, tvar st) then Polynomial!-Form!-p(tc st, y)$

endmodule;

end;