Artifact 7fe41733e5a9a816b7b5ea97ffeab72cfbef41d0b2af5490543ee1f61efa2798:

• Executable file r37/packages/odesolve/linearn.red — part of check-in [f2fda60abd] at 2011-09-02 18:13:33 on branch master — Some historical releases purely for archival purposes

  git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 2541) [annotate] [blame] [check-ins using] [more...]

module linearn;
% A solver for linear ordinary differential equations of order > 1
% ************************************************************

% Author: Malcolm MacCallum
% January 1988

% Basic ideas
%  first check for the simple case of linear equations with
%   constant coefficients. Use variation of parameters for the
%   particular integrals (see comments in module lccode),
%    [implemented up to here so far]
%   or if this fails use a Laplace transform ?
%  then look for methods for the variable coefficient case
%   starting with the homogeneous case (Euler equations)
%   [test implemented, solver still being written]
%  and I do not yet know what else....

algebraic procedure linearn(ode, y, x);
 begin scalar ans!*, rhs!*, auxeqn, lcoeff, arbconst!*;
  ans!* := 0;
  ode:=num ode;
  rhs!*:=-sub(y=0,ode);
  auxeqn := ode + rhs!*;
  lcoeff := lcof(auxeqn, df(y, x, odeorder));
  arbconst!* := lisp gensym();
  auxeqn := sub(y=e**(arbconst!**x), auxeqn/lcoeff)/e**(arbconst!**x);
  rhs!* := rhs!*/lcoeff;
   % Should form an equation in 3.3 for the solution (?), or is it
   % better to have an expression easily assigned to a variable ?
  if df(auxeqn, x) = 0 and 
     (ans!* := lccode(auxeqn, rhs!*, x, arbconst!*)) neq 0 
    then return {y=ans!*};
  odecoeffs := coeff(auxeqn, arbconst!*);
  if testeuler(odecoeffs, x)=0 and 
     (ans!* := solveeuler(odecoeffs, rhs!*, y, x, arbconst!*)) neq 0
    then return {y=ans!*};
  return odefailure(ode);
  % The kludgy way of testing if we got anywhere 
  % is to make it easy to add and debug extra cases
  % Note we assume AND only looks at the second part if first succeeds.
 end;

algebraic procedure testeuler(odecoeffs, x);
 begin scalar count!*, teuvar;
   count!* := 0;
a: teuvar := part(odecoeffs, count!*+1) * (x**(odeorder - count!*));
   if not freeof(teuvar, x) then return 1;
   count!* := count!* + 1;
   if not (count!* > odeorder) then go to a;
   traceode "This equation is of the homogeneous (Euler) type";
   return 0;
 end;

algebraic procedure solveeuler(coeffs, driver, y, x, arbconst!*);
 begin scalar sevar, rhs!*;
  sevar := part(coeffs,1)*x**odeorder;
  for i:=1:odeorder do
   sevar := sevar + 
       x**(odeorder-i)*part(coeffs,i+1)
            * (for j:=0:i-1 product (arbconst!*-j));
  rhs!* := sub(x=e**x,driver*x**odeorder);
  sevar := lccode1(sevar, rhs!*, x, arbconst!*);
  return if sevar=0 then sevar else sub(x=log x, sevar);
 end;

endmodule;      %linearn

end;