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

module odespcfn$  % Linear special function ODEs

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

% Old version temporarily preserved for testing other developments!

% A first attempt at pattern-matching solution of special (currently
% only second-order) linear odes.  Need to add more patterns.  But
% this approach may be too slow with more patterns anyway.

% If the specfn package is not loaded then we need this declaration:
algebraic operator Airy_Ai, Airy_Bi$

algebraic operator odesolve!-specfn!*$  % internal wrapper function

algebraic procedure ODESolve!-Specfn(odecoeffs1, driver1, x);
   %% Using monic coeffs for uniqueness.
   begin scalar ode, rules, soln;
      traceode1 "Looking for special-function solutions ...";
      ode := odesolve!-specfn!*(first odecoeffs1, second odecoeffs1);
      rules := {
         %% MUST use specific x, y (not ~x, ~y) for correct matching.
         %% ~~ does not seem to work in any of these rules.
         %% odesolve(df(y,x,2) - x*y, y, x);
         odesolve!-specfn!*(-x, 0) =>
            odesolve!-solns(Airy_Ai(x), Airy_Bi(x)),
         %% odesolve(df(y,x,2) - a3*x*y, y, x);
         odesolve!-specfn!*(-~a3*x, 0) =>
            odesolve!-solns(Airy_Ai(x), Airy_Bi(x), x=a3^(1/3)*x),
         %% odesolve(df(y,x,2) - (a3*x+a2b)*y, y, x);
         odesolve!-specfn!*(-(~a3*x+~a2b), 0) =>
            odesolve!-solns(Airy_Ai(x), Airy_Bi(x),
               x=a3^(1/3)*x+a2b/a3^(2/3)),

         %% The order of the following rules matters!

         %% odesolve(x^2*df(y,x,2) + x*df(y,x) - (x^2+n2)*y, y, x);
         odesolve!-specfn!*(-(1+~n2/x^2), 1/x)
            => odesolve!-solns(BesselI(n,x), BesselK(n,x),
               n = sqrt(n2)),
         %% odesolve(x^2*df(y,x,2) + x*df(y,x) + (x^2-n2)*y, y, x);
         odesolve!-specfn!*(1-~n2/x^2, 1/x)
            => odesolve!-solns(BesselJ(n,x), BesselY(n,x),
               n = sqrt(n2)),

         %% odesolve(x^2*df(y,x,2) + x*df(y,x) - (a2*x^2+n2)*y, y, x);
         odesolve!-specfn!*(-(~a2+~n2/x^2), 1/x)
            => odesolve!-solns(BesselI(n,a*x), BesselK(n,a*x),
               n = sqrt(n2), a = sqrt(a2)),
         %% odesolve(x^2*df(y,x,2) + x*df(y,x) + (a2*x^2-n2)*y, y, x);
         odesolve!-specfn!*(~a2-~n2/x^2, 1/x)
            => odesolve!-solns(BesselJ(n,a*x), BesselY(n,a*x),
               n = sqrt(n2), a = sqrt(a2)),

         %% odesolve(x*df(y,x,2) + df(y,x) - x*y, y, x);
%%          odesolve!-specfn!*(-1, 1/x)
%%             => odesolve!-solns(BesselI(0,x), BesselK(0,x)),
         %% odesolve(x*df(y,x,2) + df(y,x) - a2*x*y, y, x);
         odesolve!-specfn!*(-~a2, 1/x)
            => odesolve!-solns(BesselI(0,a*x), BesselK(0,a*x),
               a = sqrt(a2)),

         %% odesolve(x*df(y,x,2) + df(y,x) + x*y, y, x);
%%          odesolve!-specfn!*(1, 1/x)
%%             => odesolve!-solns(BesselJ(0,x), BesselY(0,x)),
         %% odesolve(x*df(y,x,2) + df(y,x) + a2*x*y, y, x);
         odesolve!-specfn!*(~a2, 1/x)
            => odesolve!-solns(BesselJ(0,a*x), BesselY(0,a*x),
               a = sqrt(a2))

                  }$
      soln := (ode where rules);        % `where' cannot produce a list!
      if soln neq ode then <<
         traceode
            "The reduced ODE can be solved in terms of special functions.";
         soln := part(soln, 1);
         %% if symbolic !*odesolve_load_specfn then load_package specfn;
         return if driver1 then
            %% BEWARE: This driver code is not well tested!
            %% traceode "But cannot currently handle the driver term! "
            { soln, ODESolve!-PI(soln, driver1, x) }
         else { soln }
      >>
   end$

algebraic operator ODESolve!-Solns!*$
listargp ODESolve!-Solns!*$

put('ODESolve!-Solns, 'psopfn, 'ODESolve!-Solns)$

symbolic procedure ODESolve!-Solns u; % (solns, subs)
   %% Avoid invalid lists on right of replacement rule, and build full
   %% optionally substituted basis data structure:
   begin scalar solns;
      %% u := revlis u;
      solns := {'list, car u, cadr u};  % algebraic list
      if (u := cddr u) then <<          % substitutions
         u := if cdr u then 'list . u else car u;
         solns := algebraic sub(u, solns)
      >>;
      return {'ODESolve!-Solns!*, solns}
   end$

endmodule$

end$