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

module odeintfc$  % Enhanced ODE solver interface

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

% Use: odesolve(ode, y, x, conds, options)
%   or   dsolve(ode, y, x, conds, options)  (cf. Maple)

% The first argument must evaluate to an ODE or a list of ODEs.

% Each dependent variable y may be either an identifier or an
% operator, such as y(x), which is replaced by a new identifier
% internally.  If y(x) is specified and x is not specified then x is
% extracted from y(x) (cf. Maple).  Equations containing operators of
% the form y(x) with different arguments x are trapped, currently as
% an error, until differential-delay equation solving is implemented.

% If a dependent variable (y) does not depend on the independent
% variable (x) then y is automatically declared to depend on x, and a
% warning message to this effect is output.  Derivatives are not
% evaluated until this dependence has been enforced.  BUT NOTE THAT
% THIS DOES NOT WORK IF THE FIRST ARGUMENT IS AN ASSIGNMENT!  This is
% because the assignment is performed BEFORE the ode solver takes
% control.  This is something of an inconsistency in the current
% REDUCE algebraic processing model.

% All arguments after the first are optional but the order must be
% preserved.  If the first argument is a list of ODEs then y is
% expected to be a list of dependent variables.  If x is specified
% then y must also be specified (first).  An empty list can be used as
% a place-holder argument.  If x and/or y are missing then they are
% parsed out of the ODE.

% Thus, possible argument combinations, each of which may optionally
% be followed by conds, are: ode | ode, y | ode, y, x

% Currently, conditions can be specified only for a single ODE.

% If specified, conds must take the form of an unordered list of
% (unordered lists of) equations with either y, x, or a derivative of
% y on the left.  A single list of conditions need not be contained
% within an outer list.  Combinations of conditions are allowed.
% Conditions within one (inner) list all relate to the same x value.
% For example:

% Boundary conditions:
% {{y=y0, x=x0}, {y=y1, x=x1}, ...}

% Initial conditions:
% {x=x0, y=y0, df(y,x)=dy0, ...}

% Combined conditions:
% {{y=y0, x=x0}, {df(y,x)=dy1, x=x1}, {df(y,x)=dy2, y=y2, x=x2}, ...}

% Boundary conditions on the values of y at various values of x may
% also be specified by replacing the variables by equations with
% single values or matching lists of values on the right, of the form:

% y = y0, x = x0  |  y = {y0, y1, ...}, x = {x0, x2, ...}

% The final argument may be one of the identifiers

% implicit, explicit, laplace, numeric, series

% specifying an option.  The options "implicit" and "explicit" set the
% switches odesolve_implicit and odesolve_explicit locally.  The other
% options specify solution techniques -- they are not yet implemented.


% TO DO:
%   Improved condition code to handle eigenvalue-type BVPs.
%   Solve systems of odes, calling crack where appropriate

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

% User interface
% ==============

put('odesolve, 'psopfn, 'odesolve!-eval)$
put('dsolve, 'psopfn, 'odesolve!-eval)$ % alternative (cf. Maple)
listargp odesolve, dsolve$              % May have single list arg.

symbolic procedure odesolve!-eval args;
   %% Establish a suitable global environment:
   (if !*div or !*intstr or !*factor or not !*exp or not !*mcd then
      NoInt2Int reval result else NoInt2Int result) where result =
      begin scalar !*evallhseqp, !*multiplicities, !*div, !*intstr,
            !*exp, !*mcd, !*factor, !*ifactor, !*precise,
            !*nopowers, !*algint, !*echo;
         %% Turn echo off to stop Win32 PSL REDUCE (only)
         %% outputting its trigsimp lap code at the end of
         %% odesolve.tst.  (Don't ask!)
         !*evallhseqp := !*exp := !*mcd := t;
         return odesolve!-eval1 args
      end$

symbolic procedure odesolve(ode, y, x);
   %% Direct symbolic-mode interface equivalent to MAHM's original.
   %% Calls odesolve!-eval to ensure correct environment.
   odesolve!-eval{ode, y, x}$

global '(ODESolve!-tracing!-synonyms)$
ODESolve!-tracing!-synonyms := '(trode trace tracing)$

symbolic procedure odesolve!-eval1 args;
   %% args = (ode &optional y x conds)
   %% Parse variables from ode if necessary (like solve),
   %% automatically declare y to depend on x if necessary and
   %% optionally impose conditions on the general solution.
   %% Support for systems of odes partly implemented so far.
   ( begin scalar ode, system, y, x, yconds, xconds, conds, soln;
      if null args then RedErr
         "ODESolve requires at least one argument -- the ODE";

      begin scalar df_simpfn, !*uncached;  !*uncached := t;
         % Turn off simplification of df (in case y does not yet
         % depend on x) before evaluating args (which may be lists):
         df_simpfn := get('df, 'simpfn);
         put('df, 'simpfn, 'simpiden);
         args := errorset!*({'revlis, mkquote args}, t);
         put('df, 'simpfn, df_simpfn);
         if errorp args then error1();
         args := car args
      end;

      ode := car args;  args := cdr args;
      system := rlistp ode;
      %% Find dependent and independent variables:
      %% (rlistp is a smacro defined in rlisp.red)
      if args then <<  y := car args;
         if rlistp y then
            if null cdr y then          % empty list - ignore
               y := 'empty
            else if rlistp cadr y or eqnp cadr y then
               y := nil                 % condition rlist
            else if system then
               y := makelist            % rlist of dependent variables
                  for each yy in cdr y collect !*a2k yy
            else MsgPri("ODESolve: invalid second argument",
               y, nil, nil, t)
         else if system then TypErr(y, "dependent var list")
         else if eqnp y then
            if cadr y memq 'output . ODESolve!-tracing!-synonyms then % option
               y := nil
            else << yconds := caddr y;  y := !*a2k cadr y >>
         else if not smember(y:=!*a2k y, ode) then
            y := nil                    % option
      >>;
      if y then args := cdr args;
      if args and y then <<  x := car args;
         if rlistp x then
            if null cdr x then          % empty list - ignore
               x := 'empty
            else x := nil               % condition list
         else if eqnp x then
            if cadr x memq 'output . ODESolve!-tracing!-synonyms then % option
               x := nil
            else << xconds := caddr x;  x := !*a2k cadr x >>
         else if not smember(x:=!*a2k x, ode) then
            x := nil                    % option
      >>;
      if x then args := cdr args;
      if y eq 'empty then y := nil;
      if x eq 'empty then x := nil;
      %% If x not given and y an operator (list) then extract x:
      if null x and y then
         if rlistp y then
            begin scalar yy;  yy := cdr y;
               while yy and atom car yy do yy := cdr yy;
               if yy and cdar yy then x := cadar yy;
               if not idp x then x := nil
            end
         else if pairp y and cdr y then x := cadr y;
      %% Finally, attempt to parse variables from ode if necessary:
      if null y or null x then
         %%%%% NOTE: ODE ALREADY AEVAL'ED ABOVE
         %%%%% so some of the following is now redundant !!!!!
      begin scalar df_simpfn, k_list, df_list;
         % Turn off simplification of df (in case y does not yet
         % depend on x) before evaluating ode, which may be a list:
         df_simpfn := get('df, 'simpfn);
         put('df, 'simpfn, 'simpiden);
         k_list := errorset!*({'get_k_list, mkquote ode}, t);
         put('df, 'simpfn, df_simpfn);
         if errorp k_list then error1() else k_list := car k_list;
         df_list := get_op_knl('df, car k_list);
         for each knl in cdr k_list do
            df_list := union(df_list, get_op_knl('df, knl));
         %% df_list is set of derivatives in ode(s).
         if null df_list then RedErr
            "No derivatives found -- use solve instead.";
         %% df_list = ((df y x ...) ... (df z x ...) ... )
         if null y then <<
            y := cadar df_list . nil;
            for each el in cdr df_list do
               if not member(cadr el, y) then y := cadr el . y;
            %% y is a list at this point.
            if system then
               if length ode < length y then
                  RedErr "ODESolve: under-determined system of ODEs."
               else y := makelist y     % algebraic list of vars
            else
               if cdr y then
                  MsgPri("ODESolve -- too many dependent variables:",
                     makelist y, nil, nil, t)
               else y := car y;         % single var
            MsgPri("Dependent var(s) assumed to be", y, nil, nil, nil)
         >>;
         if null x then <<
            x := caddar df_list;
            MsgPri("Independent var assumed to be", x, nil, nil, nil)
         >>;
      end;

      %% Process the ode (re-simplifying derivatives):
      EnsureDependency(y, x);
      ode := aeval ode;
      %% !*eqn2a is defined in alg.red
      if system then
         if length ode > 2 then
            %% RedErr "Solving a system of ODEs is not yet supported."
            %% Skip conditions TEMPORARILY!
            return ODESolve!-Depend(
               makelist for each o in cdr ode collect !*eqn2a o, y, x, nil)
         else
            << ode := !*eqn2a cadr ode;  y := cadr y >>
      else ode := !*eqn2a ode;

      %% Process conditions (re-simplifying derivatives):
      if args then
         if rlistp(conds := aeval car args) then <<
            args := cdr args;
            conds := if not rlistp cadr conds then conds . nil
            else cdr conds
         >> else conds := nil;
      %% Now conds should be a lisp list of rlists (of equations).
      if yconds then
         yconds := if rlistp yconds then cdr yconds else yconds . nil;
      if xconds then
         xconds := if rlistp xconds then cdr xconds else xconds . nil;
      %% Concatenate separate x & y conds onto conds list:
      while yconds and xconds do <<
         conds := {'list, {'equal, x, car xconds},
            {'equal, y, car yconds}} . conds;
         yconds := cdr yconds;  xconds := cdr xconds
      >>;
      if yconds or xconds then
         RedErr "Different condition list lengths";
      if conds then
         %% Move this into odesolve!-with!-conds?
         conds := makelist odesolve!-sort!-conds(conds, y, x);

      %% Process remaining control option arguments:
      while args do
      begin scalar arg;  arg := car args;  args := cdr args;
         if eqnp arg then               % equation argument
            if cadr arg eq 'output then
               args := caddr arg . args
            else if cadr arg memq ODESolve!-tracing!-synonyms then
               !*trode := caddr arg
            else MsgPri("Invalid ODESolve option", arg,
               "ignored.", nil, nil)
                                        % keyword argument
         else if arg memq
            '(implicit explicit expand noint verbose
               basis noswap norecurse fast check) then
            set(mkid('!*odesolve_, arg), t)
         else if arg eq 'algint then on1 'algint
         else if arg eq 'full or !*odesolve_full then
            !*odesolve_expand := !*odesolve_explicit := t
         else if arg memq ODESolve!-tracing!-synonyms then !*trode := t
         else if arg memq '(laplace numeric series) then
            RedErr{"ODESolve option", arg, "not yet implemented."}
               %% Pass remaining args to routine called
         else RedErr{"Invalid ODESolve option", arg}
      end;

      if !*odesolve_verbose then algebraic <<
         write "ODE: ", num ode=0;
         write "Dependent variable: ", y, ";  independent variable: ", x;
         write "Conditions: ", symbolic(conds or "none");
      >>;

      %% Rationalize conflicting options:
      %% Conditions override basis
      if conds then !*odesolve_basis := nil;
%%       %% Basis overrides explicit
%%       if !*odesolve_basis then !*odesolve_explicit := nil;

      %% Finally, solve the ode!
      if not getd 'ODESolve!*0 then     % for testing
         return {'ODESolve, ode, y, x, conds};
%%       soln := if conds then
%%          odesolve!-with!-conds(ode, y, x, conds)
%%       else odesolve!-depend(ode, y, x);
      if null(soln := ODESolve!-Depend(ode, y, x, conds)) then
         return algebraic {num ode=0};
      %% Done as follows because it may be easier to solve after
      %% imposing conditions than before, and it would be necessary to
      %% remove root_of's before imposing conditions anyway.
      if !*odesolve_explicit and not ODESolve!-basisp soln then
         soln := ODESolve!-Make!-Explicit(soln, y, conds);
      if !*odesolve_expand then
         soln := expand_roots_of_unity soln;
      if !*odesolve_check then
         ODE!-Soln!-Check(if !*odesolve_noint then NoInt2Int soln else soln,
            ode, y, x, conds) where !*noint = t;
      return soln

   end ) where !*odesolve_implicit = !*odesolve_implicit,
               !*odesolve_explicit = !*odesolve_explicit,
               !*odesolve_expand = !*odesolve_expand,
               !*trode = !*trode,
               !*odesolve_noint = !*odesolve_noint,
               !*odesolve_verbose = !*odesolve_verbose,
               !*odesolve_basis = !*odesolve_basis,
               !*odesolve_noswap = !*odesolve_noswap,
               !*odesolve_norecurse = !*odesolve_norecurse,
               !*odesolve_fast = !*odesolve_fast,
               !*odesolve_check = !*odesolve_check$

symbolic procedure Odesolve!-Make!-Explicit(solns, y, conds);
   <<
      %% SHOULD PROBABLY CHECK THAT Y IS NOT INSIDE AN UNEVALUATED
      %% INTEGRAL BEFORE TRYING TO SOLVE FOR IT -- IT SEEMS TO UPSET
      %% SOLVE!
      solns := for each soln in cdr solns join
         if cadr soln eq y then {soln} else <<
            %% soln is an implicit solution of ode for y
            %% for each s in cdr reval aeval {'solve, soln, y} join
            %%    %% Make this test optional?
            %%    if eqcar(caddr s, 'root_of) or eval('and .
            %%       mapcar(cdr expand_roots_of_unity subeval{s, ode},
            %%          'zerop))
            %%    then {s}
            %%    else if !*trode then algebraic write "Solution ", s,
            %%    " discarded -- does not satisfy ODE";
            traceode
               "Solution before trying to solve for dependent variable is ",
               soln;
            cdr reval aeval {'solve, soln, y}
         >>;
      %% It is reasonable to return root_of's here.
      %% Solving can produce duplicates, so ...
      %% solns := union(solns, nil); % union still necessary?

      %% Check that each explicit solution still satisfies any
      %% conditions:
      if conds then
         for each cond in cdr conds do  % each cond is an rlist
         begin scalar xcond;
            xcond := cadr cond;
            cond := makelist for each c in cddr cond collect !*eqn2a c;
            solns := for each s in solns join
               if eqcar(caddr s, 'root_of) or
                  union(cdr %% trig_simplify
                     subeval{xcond, subeval{s, cond}}, nil) = {0}
               then {s}
               else algebraic traceode "Solution ", s,
               " discarded -- does not satisfy conditions";
         end;
      makelist solns
   >>$

% Should now be able to use the standard package `trigsimp' instead!
algebraic procedure trig_simplify u;
   u where tan_half_angle_rules$

algebraic(tan_half_angle_rules := {
   sin(~u) => 2tan(u/2)/(1+tan(u/2)^2),
   cos(~u) => (1-tan(u/2)^2)/(1+tan(u/2)^2) })$
%% Cannot include tan rule -- recursive!


symbolic procedure get_k_list ode;
   %% Return set of all top-level kernels in ode or rlist of odes.
   %% (Do not cache to ensure derivatives are [eventually] evaluated
   %% properly!)
   begin scalar k_list, !*uncached;  !*uncached := t;
      %% Do not make an assignment twice:
      if eqcar(ode, 'setk) then ode := caddr ode;
      if rlistp(ode := reval ode) then <<
         k_list := get_k_list1 cadr ode;
         for each el in cddr ode do
            k_list := union(k_list, get_k_list1 el)
      >>
      else k_list := get_k_list1 ode;
      return k_list
   end$

symbolic procedure get_k_list1 ode;
   union(kernels numr o, kernels denr o)
      where o = simp !*eqn2a ode$

symbolic procedure get_op_knl(op, knl);
   %% Return set of all operator kernels within knl with op as car.
   if pairp knl then
      if car knl eq op then knl . nil
      else
         ( if op_in_car then union(op_in_car, op_in_cdr) else op_in_cdr )
            where op_in_car = get_op_knl(op, car knl),
               op_in_cdr = get_op_knl(op, cdr knl)$

symbolic procedure EnsureDependency(y, x);
   for each yy in (if rlistp y then cdr y else y . nil) do
      if not depends(yy, x) then <<
         MsgPri("depend", yy, ",", x, nil);
         depend1(yy, x, t)
      >>$


symbolic procedure odesolve!-sort!-conds(conds, y, x);
   %% conds is a lisp list of rlists of condition equations.
   %% Return a canonical condition list.
   %% Collect conditions at the same value of x, check them for
   %% consistency and sort them by increasing order of derivative.
   begin scalar cond_alist;
      for each cond in conds do
      begin scalar x_cond, y_conds, x_alist;
         if not rlistp cond then TypErr(cond, "ode condition");

         %% Extract the x condition:
         y_conds := for each c in cdr cond join
            if not CondEq(c, y, x) then
               TypErr(c, "ode condition equation")
            else if cadr c eq x then << x_cond := c; nil >>
            else c . nil;
         if null x_cond then
            MsgPri(nil, x, "omitted from ode condition", cond, t);
         if null y_conds then
            MsgPri(nil, y, "omitted from ode condition", cond, t);

         %% Build the new condition alist, with the x condition as key:
         if (x_alist := assoc(x_cond, cond_alist)) then
            nconc(x_alist, y_conds)
         else cond_alist := (x_cond . y_conds) . cond_alist
      end;
      %% Now cond_alist is a list of lists of equations, each
      %% beginning with a unique x condition.

      %% Sort the lists and return a list of rlists:
      return for each cond in cond_alist collect makelist
         if null cddr cond then cond else car cond .
         begin scalar sorted, next_sorted, this, next, result;
            sorted := sort(cdr cond, 'lessp!-deriv!-ord);
            %% sorted is a list of equations.
            while sorted and (next_sorted := cdr sorted) do <<
               if cadr(this := car sorted) eq
                  cadr(next := car next_sorted) then
                  %% Two conds have same lhs, so ...
                  ( if caddr this neq caddr next then
                     MsgPri("Inconsistent conditions:",
                        {'list, this, next}, "at", car cond, t) )
                           % otherwise ignore second copy
               else result := this . result;
               sorted := next_sorted
            >>;
            return reversip(next . result)
         end
   end$

symbolic procedure CondEq(c, y, x);
   %% Return true if c is a valid condition equation for y(x).
   %% cf. eqexpr in alg.red
   eqexpr c and ( (c := cadr c) eq x or c eq y or
      (eqcar(c, 'df) and cadr c eq y and caddr c eq x
         %% Is the following test overkill?
         and (null cdddr c or fixp cadddr c)) )$

symbolic procedure lessp!-deriv!-ord(a, b);
   %% (y=y0) < (df(y,x)=y1) and df(y,x,m)=ym < df(y,x,n)=yn iff m < n
   %% But y might be a kernel rather than an identifier!
   if atom(a := cadr a) then            % a = (y=?)
      not atom cadr b                   % b = (df(y,x,...)=?)
   else if atom(b := cadr b) then       % b = (y=?)
      not atom cadr a                   % a = (df(y,x,...)=?)
   else if not(car a eq 'df) then       % a = (y(x)=?)
      car b eq 'df                      % b = (df(y(x),x,...)=?)
   else                                 % a = (df(y,x,...)=?), any y
      car b eq 'df and                  % b = (df(y,x,...)=?)
         if null(a := cdddr a) then     % a = (df(y,x)=?)
            cdddr b                     % b = (df(y,x,n)=?), 1 < n
         else                           % a = (df(y,x,m)=?), m > 1
            (b := cdddr b) and
               car a < car b$           % b = (df(y,x,n)=?), m < n


%%% THE FOLLOWING PROCEDURE SHOULD PROBABLY INCLUDE THE CODE TO MAKE
%%% SOLUTIONS EXPLICIT BEFORE RESTORING OPERATOR FORMS FOR Y.

symbolic procedure ODESolve!-Depend(ode, y, x, conds);
   %% Check variables and dependences before really calling odesolve.
   %% If y is an operator kernel then check whether ode is a
   %% differential-delay equation, and if not solve ode with y
   %% replaced by an identifier.
   ( begin scalar xeqt, ylist, sublist;
      y := if rlistp ode then cdr y else y . nil;
      %% Using `t' as a variable causes trouble when checking
      %% dependence of *SQ forms, which may contain `t' as their last
      %% element, so...
      if x eq t then <<
         xeqt := t;  x := gensym();
         for each yy in y do if idp yy then depend1(yy, x, t);
         %% Cannot simply use `sub' on independent variables in
         %% derivatives, so...
         ode := subst(x, t, reval ode); % reval subst?
         if conds then conds := subst(x, t, reval conds); % reval subst?
         sublist := (t.x) . sublist
      >>;
      for each yy in y do
         if idp yy and not(yy eq t) then <<
            %% Locally and quietly remove any spurious inverse
            %% implicit dependence of x on y:
            ylist := yy . ylist;
            depend1(x, yy, nil) where !*msg = nil;
         >> else                        % replace variable
         begin scalar yyy;
            yyy := gensym();  depend1(yyy, x, t);
            ylist := yyy . ylist;
            put(yyy, 'odesolve!-depvar, yy); % for later access
            sublist := (yy.yyy) . sublist;
            if xeqt then yy := subeval{{'equal,t,x}, yy};
            odesolve!-delay!-check(ode, yy);
            ode := subeval{{'equal,yy,yyy}, ode};
            if conds then conds := subeval{{'equal,yy,yyy}, conds}
         end;
      ylist := reverse ylist;
      ode := if rlistp ode then
         ODESolve!-System(cdr ode, ylist, x)
      else if conds then
         odesolve!-with!-conds(ode, car ylist, x, conds)
      else ODESolve!*0(ode, car ylist, x);
      if null ode then return;
      if sublist then
      begin scalar !*NoInt;
         %% Substitute into derivatives and integrals
         %% (and turn off integration for speed).
         ode := reval ode;              % necessary?
         for each s in sublist do
            ode := subst(car s, cdr s, ode);
         %% ode := reval ode;              % necessary?
      end;
      return ode
   end ) where depl!* = depl!*$

symbolic procedure ODESolve!-System(ode, y, x);
   %% TEMPORARY
   {'ODESolve!-System, makelist ode, makelist y, x}$
algebraic operator ODESolve!-System$

symbolic procedure odesolve!-delay!-check(ode, y);
   %% Check that ode is not a differential-delay equation in y,
   %% i.e. check that every occurrence of the operator y = y(x) has
   %% the same argument (without any shifts).  This could be used as a
   %% hook to call an appropriate solver -- if there were one!
   begin scalar odelist;
      odelist := if rlistp ode then cdr ode else ode . nil;
      for each ode in odelist do
         ( for each knl in kernels numr simp ode do
            for each yy in get_op_knl(y_op, knl) do
               if not(yy eq y) then
                  MsgPri("Arguments of", y_op, "differ --",
                     "solving delay equations is not implemented.", t) )
                        where y_op = car y
   end$

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

% Impose initial/boundary conditions
% ==================================

% A first attempt to impose initial/boundary conditions on the
% solution of a single ode returned by odesolve.

% Solving with conditions provides access to the general solution as
% the value of the global algebraic variable ode_solution.

% If the solution is explicit then the following code could be
% simplified and should be slightly more efficient, but is it worth
% testing for an explicit solution and adding the special case code?

algebraic procedure odesolve!-with!-conds(ode, y, x, conds);
   % conds must be a list of ordered lists of the form
   %    {x=x0, y=y0, df(y,x)=y1, df(y,x,2)=y2, ...}.
   % All conditions applied at the same value of x must be collected
   % into the same list.
   % More generality is allowed only by odesolve!-eval.
   % This code could perhaps be more efficient by building a list of
   % all required derivatives of the ode solution once and for all?
   begin scalar first!!arbconst, arbconsts;
      first!!arbconst := !!arbconst + 1;

      %% Find the general solution of the ode and assign it to the
      %% global algebraic variable ode_solution:
      %1.03% ode_solution := odesolve!-depend(ode, y, x);
      ode_solution := symbolic ODESolve!*0(ode, y, x);
      if not ode_solution then return;
      traceode "General solution is ", ode_solution;
      traceode "Applying conditions ", conds;

      arbconsts :=
         for i := first!!arbconst : !!arbconst collect arbconst i;

      return for each soln in ode_solution join
         odesolve!-with!-conds1(soln, y, x, conds, arbconsts)
   end$

algebraic procedure odesolve!-with!-conds1(soln, y, x, conds, arbconsts);
   begin scalar arbconsteqns;
      %% Impose the conditions (care is needed if the solution is
      %% implicit):
      arbconsteqns := for each cond in conds join
         begin scalar xcond, ycond, dfconds, arbconsteqns;
            xcond := first cond;  cond := rest cond;
            ycond := first cond;
            if lhs ycond = y then cond := rest cond else ycond := 0;
            %% Now cond contains only conditions on derivatives.
            arbconsteqns :=
               if ycond then         % Impose the condition on y:
                  {sub(xcond := {xcond, ycond}, soln)}
               else {};
            dfconds := {};
            %% Impose the conditions on df(y, x, n).  If the solution
            %% is implicit, then in general all lower derivatives will
            %% be introduced, so ...
            while cond neq {} do
            begin scalar dfcond, result;
               %% dfcond : next highest derivative
               %% result : of substituting for all derivatives
               %% dfconds : all derivatives so far including this one
               dfcond := first cond;  cond := rest cond;
               dfconds := dfcond . dfconds;
               %% All conditions on derivatives are handled before
               %% conditions on x and y to protect against
               %% substituting for x or y in df(y,x,...):
               result := sub(dfconds, map(y => lhs dfcond, soln));
               if not(result freeof df) then % see comment below
                  RedErr "Cannot apply conditions";
               arbconsteqns := sub(xcond, result) . arbconsteqns
            end;
            return arbconsteqns
         end;

      %% Solve for the arbitrary constants:
      arbconsts := solve(arbconsteqns, arbconsts);
      %% ***** SHOULD CHECK THAT THE SOLUTION HAS SUCCEEDED! *****
      %% and substitute each distinct arbconst solution set into the
      %% general ode solution:
      return for each cond in arbconsts collect
         if rhs soln = 0 then           % implicit solution
            num sub(cond, lhs soln) = 0
         else
            sub(cond, soln)
   end$

%% The above df error can happen only if the solution is implicit and
%% a derivative is missing from the sequence, which is unlikely.
%% Should try to recover by computing the value of the missing
%% derivative from the conditions on lower order derivatives, and
%% letting solve eliminate them.  Try this later IF it ever proves
%% necessary.

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

% Check the solution set
% ======================

% This code checks that the solution satisfies the ODE and that a
% general solution is general.

% A `solution' is either a basis solution:
%   {{b1(x), b2(x), ..., bn(x)}, PI(x)}; PI may be omitted if zero
% or a list of component solutions.

algebraic procedure ODE!-Soln!-Check(soln, ode, y, x, conds);
   %% SOLN is a LIST of solutions; ODE is a differential EXPRESSION; Y
   %% is the dependent variable; X is the independent variable; CONDS
   %% is true if conditions were specified.
   begin scalar n, !*allowdfint, !*expanddf;
      symbolic(!*allowdfint := !*expanddf := t);
      ode := num !*eqn2a ode;           % returns ode as expression
      %% Should compute order `on demand' in ODE!-Comp!-Soln!-Fails.
      n := ODE!-Order(ode, y);
      if symbolic ODESolve!-basisp soln then <<
         %% Basis solution (of linear ODE).
         %% Only arises as general solution.
         %% Remove contribution from PI if there is one:
         if arglength soln = 2 and second soln then
            ode := num sub(y = y + second soln, ode);
         if length(soln := first soln) neq n then
            write "ODESolve warning - ",
               "wrong number of functions in basis!";
         %% Test each basis function in turn:
         for each s in soln do
            if (s:=sub(y = s, ode)) and trigsimp s then
               write "ODESolve warning - ",
                  "basis function may not satisfy ODE: ", s
      >> else <<
         %% List of component solutions.
         %% Check generality:
         if not conds and ODESolve!-arbconsts soln < n then
            write "ODESolve warning - ",
               "too few arbitrary constants in general solution!";
         for each s in soln do
            if ODE!-Comp!-Soln!-Fails(s, ode, y, x, n) then
               write "ODESolve warning - ",
                  "component solution may not satisfy ODE: ", s;
      >>
   end$

% Each component solution may be
% explicit: y = f(x)
% implicit: f(x,y) = g(x,y); rhs MAY be 0
% unsolved: ode = 0, but CAN THIS CASE ARISE?
% parametric: {y = g(p), x = f(p), p}

algebraic procedure ODE!-Comp!-Soln!-Fails(soln, ode, y, x, n);
   %% SOLN is a SINGLE component solution; ODE is a differential
   %% EXPRESSION; Y is the dependent variable; X is the independent
   %% variable; N is the order of ODE.
   if symbolic eqnp soln then           % explicit, implicit or unsolved
      if lhs soln = y and rhs soln freeof y then
                                        % explicit: y = f(x)
         (if (ode := sub(soln, ode)) then trigsimp ode)
      else if rhs soln = 0 and lhs soln = ode then
         1                              % unsolved: ode = 0
      else                              % implicit: f(x,y) = 0
      begin scalar derivs, deriv;
         %% Construct in `derivs' a list of successive derivatives of
         %% the implicit solution f(x,y) up to the order of the ODE in
         %% decreasing order; each expression is linear in the highest
         %% derivative.
         derivs := {soln := num !*eqn2a soln};
         for i := 1 : n do
            derivs := (soln:=num df(soln,x)) . derivs;
         %% Substitute for each derivative in ODE in turn in
         %% decreasing order until the result is zero; if not the
         %% solution fails.
         while n > 0 and <<
            deriv := solve(first derivs, df(y,x,n)); % linear
            if deriv = {} then 0 else
               ode := num sub(first deriv, ode) >> do <<
               n := n - 1;
               derivs := rest derivs
               >>;
         if deriv = {} then <<
            write "ODESolve warning - cannot compute ", df(y,x,n);
            return 1
         >>;
         derivs := first derivs;
         ode := (ode where derivs => 0); % for tracing
         return ode                     % 0 for good solution
      end
   else if symbolic(rlistp soln and eqnp cadr soln) then
                                        % parametric: {y = g(p), x = f(p), p}
   begin scalar xx, yy, p, dp!/dx, deriv, derivs;
      yy := rhs first soln;             % Should not depend on ordering!
      xx := rhs second soln;            % Should not depend on ordering!
      p := third soln;                  % parameter
      %% Construct in `derivs' a list of successive derivatives of the
      %% parametric solution (yy,xx) up to the order of the ODE in
      %% decreasing order.
      dp!/dx := 1/df(xx,p);
      derivs := {deriv:=yy};
      for i := 1 : n do
         derivs := (deriv:=dp!/dx*df(deriv,p)) . derivs;
      %% Substitute for each derivative in ODE in turn in
      %% decreasing order until the result is zero; if not the
      %% solution fails.
      while n > 0 and (ode := num sub(df(y,x,n)=first derivs, ode)) do <<
         n := n - 1;
         derivs := rest derivs
      >>;
      return sub(y=yy, x=xx, ode)
   end
   else write "ODESolve warning - invalid solution type: ", soln$

%% Code to find the actual number of arbitrary constants in a solution:

fluid '(ODESolve!-arbconst!-args)$

symbolic operator ODESolve!-arbconsts$
symbolic procedure ODESolve!-arbconsts u;
   %% Return the number of distinct arbconsts in any sexpr u.
   begin scalar ODESolve!-arbconst!-args;
      ODESolve!-arbconsts1 u;
      return length ODESolve!-arbconst!-args
   end$

symbolic procedure ODESolve!-arbconsts1 u;
   %% Collect all the indices of arbconsts in u into a set in the
   %% fluid variable ODESolve!-arbconst!-args.
   if not atom u then
      if car u eq 'arbconst then
         (if not member(cadr u, ODESolve!-arbconst!-args) then
            ODESolve!-arbconst!-args := cadr u . ODESolve!-arbconst!-args)
      else
         << ODESolve!-arbconsts1 car u;  ODESolve!-arbconsts1 cdr u >>$

endmodule$

end$