ode_extra.mac

DEBUGFLAG : 1$
MAX_LENGTH_FOR_SIMPLIFICATION:30000$
/* ************************************************************************* */
/* ***** print expr only when flag<DEBUGFLAG                           ***** */
/* ************************************************************************* */
dprint(flag,[_expr])::= if flag <= DEBUGFLAG then buildq ([_expr], print (splice (_expr)))$
/* ************************************************************************* */


/*****************************************************************************************************/
/* ----- creates a list of user defined functions from expr. the functions should not have any properties
         with the possible exception of noun                                                   ----- */
/*****************************************************************************************************/
listUDF(_expr):= block( [ _counter:0,listUDF:[]],
        listUDFPriv (_expr, []),
        return(unique(listUDF))
)$

listUDFPriv(_expr, _opList) :=
 block ( [_x, _args, _newList],
        if atom(_expr) then ( 
           _opList 
        )
        else if udf(_expr) then (
          listUDF: cons(_expr,listUDF) 
        )
        else ( 
          _x:    op(_expr),
          _args: args(_expr),
          /* add the operators to the opList */
          _newList: cons(_x, _opList),
          for _arg in _args do _newList: listUDFPriv(_arg, _newList),
          _newList
        )
)$


listUDFGeneral(_expr):= block( [ _counter:0,listUDF:[]],
        listUDFGeneralPriv (_expr, []),
        return(unique(listUDF))
)$

listUDFGeneralPriv(_expr, _opList) :=
 block ( [_x, _args, _newList],
        if atom(_expr) then ( 
           _opList 
        )
        else if udfGeneral(_expr) then (
          listUDF: cons(_expr,listUDF) 
        )
        else (
          _x:    op(_expr),
          _args: args(_expr),
          /* add the operators to the opList */
          _newList: cons(_x, _opList),
          for _arg in _args do _newList: listUDFGeneralPriv(_arg, _newList),
          _newList
        )
)$

simplifyingp(_f):= symbolp(_f) and is(?get(_f,?operators)=false)$

/*udf(_f):=not(stringp(op(_f)) or symbolp(op(_f)) or (simplifyingp(op(_f))))$*/
/* a user defined function is: 
     not a a string, just like sin(x), but not like x 
     a symbolp, like sin(x), but not like x
     a function with no properties
*/


/*udf(_f):=not(stringp(op(_f))) and (symbolp(op(_f))) and (apply(properties,[op(_f)])=[] or apply(properties,[op(_f)])=[noun])$*/
/* returns true if it is a udf*/
udf(_f):=not(stringp(op(_f))) and (symbolp(op(_f))) and simplifyingp(op(_f)) and not ?fboundp(op(_f)) $

/* general functions, including sin(a*x+b*y) */
/* we do not want diff() stuff here, they are usually udf and we do not support unevaluated diffs for known functions */
udfGeneral(_f):=not(stringp(op(_f))) and (symbolp(op(_f))) and not ?fboundp(op(_f)) and not(op(_f)=op(diff(__f1(__x),__x)))$

listGDF(_expr,_y,_x):= block( [ _counter:0,listGDF:[]],

        
        listGDFPriv (_expr, _y,_x, []),
        return(unique(listGDF))
)$

listGDFPriv(_expr, _y,_x, _opList) :=
 block ( [_op, _args, _newList],
        if atom(_expr) then ( 
           _opList 
        )
        else if gdfGeneral1(_expr,_y,_x) then (
          listGDF: cons(_expr,listGDF) 
        )
        else (
          _op:    op(_expr),
          _args: args(_expr),
          /* add the operators to the opList */
          _newList: cons(_op, _opList),
          for _arg in _args do _newList: listGDFPriv(_arg,_y,_x, _newList),
          _newList
        )
)$

/* c^f(x,y)  or f(x,y)^c  or sqrt(f(x,y))*/
gdfGeneral1(_f,_y,_x):=(op(_f)="^" and (
                         ((atom(args(_f)[1]) and freeof(_x,_y,args(_f)[1])) and not(freeof(_x,args(_f)[2])) and not(freeof(_y,args(_f)[2])) ) 
                         or 
                         ((atom(args(_f)[2]) and freeof(_x,_y,args(_f)[2])) and not(freeof(_x,args(_f)[1])) and not(freeof(_y,args(_f)[1])) )
                         )) 
                      or (op(_f)=sqrt and not(freeof(_x,args(_f)[1])) and not(freeof(_y,args(_f)[1])) )
$


/* ************************************************************************** */
/* we introduce dependencies and gradefs for the ode                          */
/* ************************************************************************** */
explicit_form_to_dependencies_form1(_ode) :=block([_udfargs], 
/* ************************************************************************** */
  /* get a list of user defined functions */
  _dependencylist : copy(dependencies),
 
  _listudf : listUDF(rhs(_ode)),
  _listudf : sublist(_listudf,lambda([_i], not(length(args(_i))=1 and atom(args(_i)[1])))),
  _udf_op : map(op,_listudf),
  _udf_args : flatten(map(args,_listudf)), 

  /* get a list of generally defined functions */ 
  _listgdf : listUDFGeneral(rhs(_ode)),
  _listgdf : sublist(_listgdf,lambda([_i], not(length(args(_i))=1 and atom(args(_i)[1])))),
  _listgdf : append(_listgdf,listGDF(rhs(_ode),_y,_x)),
  /* remove everything that is already in listudf */
  _listgdf : unique(sublist(_listgdf,lambda([_i],not(member(_i,_listudf))))),
  _gdf_op : map(op,_listgdf),
  _gdf_args : unique(flatten(map(args,_listgdf))),

  /* for the udf, we need to introduce new variables      */
  /* note that we remove gradefs by removing dependencies */
  _varlist:makelist(concat(%g,_i),_i,1,length(_udf_args)),
  dprint(5,"varlist = ",_varlist), 
   depends(_varlist,[x,y]),
   _udfargs : _udf_args,  
   for _g in _varlist do (
      apply('gradef,[_g,_x,diff(first(_udfargs),_x)]),   
      apply('gradef,[_g,_y,diff(first(_udfargs),_y)]),
      _udfargs:rest(_udfargs)    
   ), 

  /* replace all udfs with their operators: f(ax+by)  -> f */
  _ode : subst(map("=",_listudf,_udf_op),_ode),

  /* only dependencies for the udf, not for the general free functions */
  map(depends,_udf_op,_varlist),   

  dprint(5,"canonical form of ode, replaced long arguments = ",_ode),
  dprint(5,"dependencies=",dependencies),
  dprint(5,"listudf = ",_listudf), 
  dprint(5,"udf_op = ",_udf_op), 
  dprint(5,"udf_args = ",_udf_args), 
  dprint(5,"listgdf = ",_listgdf), 
  dprint(5,"gdf_op = ",_gdf_op), 
  dprint(5,"gdf_args = ",_gdf_args), 
  /* ***** ************************************************* ***** */
  return(_ode)
)$
/* ************************************************************************** */





/*****************************************************************************************************/
/* ----- calculates the number of operators in the expression ----- */
/* ----- this is a simple measure of complexity               ----- */
/*****************************************************************************************************/
/* note: this is incredibly slow */
nrOps(_expression):=
 block( [ _counter:0],
        dprint(6,"expression:",_expression),
        nrOpsPriv (_expression, [])
       )$

nrOpsPriv(_expression, _opList) :=
 block ( [_x, _args, _newList],
        /* if expression is an atom, then we return opList */
        if atom(_expression)
           then _opList
           else (
              _x:    op(_expression),
              _args: args(_expression),
              /* add the operators to the opList */
              _newList: cons(_x, _opList),
              for _arg in _args do
                /* also expand all the subexpressions and count the operators */
                _newList: nrOpsPriv(_arg, _newList),
              _newList
             )
        )$






/* ************************************************************************** */

/* ************************************************************************** */
simplify(_S):=block([_N,_Nnew,_Snew,_ratvars,_lS,_isimaginary:false,_oldradexpand,_N0,_S0],
  dprint(5,"simplify::START, S=",grind(_S)),
  dprint(6,"dependencies = ",dependencies),

  /* some of the below (fullratsimp) has problems with pdiff, so if there are no integrals, it is save to switch temporarily to standard derivatives */
  
  if freeof('integrate,_S) then (use_pdiff:false, _S:convert_to_diff(_S)),

  /* simplify unfactored subexpressions like sqrt(f+g)/(fg+g^2) -> 1/sqrt(f+g) that are not factored */
  /* try to evaluate nouns */

  _N0 : slength(string(_S)),
  _N : _N0,
  _S0 : _S,
  _Snew : ev(_S,nouns),
  _Nnew:slength(string(_Snew)),
  if (_Nnew < 2.0*_N) then (_S : _Snew, _N : slength(string(_S))),
   

  dprint(6,"simplify:: S=",grind(_S)),
  if not freeof(%i,_S) then _isimaginary:true,

  if not freeof(%i,_S) then (
    dprint(2,"trying to get rid of imaginary numbers, check the result:",_S),
    /*_S : rectform(_S),*/
    /* for instance kamke 1.89*/
    _Snew : rootscontract(logarc(rootscontract(_S))),
    /* sometimes the imaginary number is only in the integration constant... */
    if (freeof(%c,subst(%i*%c = %k,_Snew))) then _Snew : subst(%i*%c=%c,_Snew),
    if freeof(%i,_Snew) then (_S : _Snew, _isimaginary:false, print("got rid of imaginary numbers...")),
    dprint(2,"trying to get rid of imaginary numbers, check the result:",_S)
  ),


  /* ***** we have a maximum limit for the expression, long expressions take forever to separate ***** */
  _lS : slength(string(_S)),  
  dprint(6,"   length S = ",_lS),  
  if (_lS>MAX_LENGTH_FOR_SIMPLIFICATION) then (
      dprint(1,"Warning: subexpression is too long to consider simplification. Length = ",_lS),
      dprint(1,"Returning original expression..."),
      dprint(1,"If you want to try anyway, put MAX_LENGTH_FOR_SIMPLIFICATION to a higher value. Current value: ",MAX_LENGTH_FOR_SIMPLIFICATION),
      return(_S)
  ),

  _ratvars : showratvars(_S),
  dprint(6,"ratvars = ",_ratvars),
  /* we use length of expression as a simple measure for complexity to check for expression swell*/
  _N : slength(string(_S)),
  _Snew : fullratsimp(_S),
  _Nnew : slength(string(_Snew)),
  dprint(6,"   N,Nnew = ",_N," ",_Nnew),  
  if _isimaginary=false and not freeof(%i,_Snew) then (
    dprint(2,"simplification introduced complex numers, ignoring fullratsimp-simplification")
  ) else if _Nnew < 2.0*_N then (_S : _Snew, _N : slength(string(_S))),



  dprint(6,"2. S=",grind(_S)),
  /* if we have trig functions, we only want to simplify if they depend on x or y */ 
  /*if sublist(_ratvars,lambda([_i],not freeof(sin,cos,tan,_i) and not freeof(_x,_y,_i))) # [] then (*/

  /* we should not do this when we have derivatives of unknown functions, trig functions do not like that... */
  /*the only ode in the kamke database that is trg + diff(f(x)) is kamke1.80 */
  if freeof(diff,_ratvars) and sublist(_ratvars,lambda([_i],not freeof(csc,sec,cot,sin,cos,tan,_i) )) # []  then (
    dprint(6,"3.0 S=",grind(_S)),
    _N : slength(string(_S)),

    _Strig : _S,
    _Snew : trigreduce(_S), 
    dprint(6,"3.1 trigreduce,S=",_Snew),
    _Nnew : slength(string(_Snew)),
    dprint(6,"   N,Nnew = ",_N," ",_Nnew),  
    if _Nnew < 25*_N then ( _Strig : _Snew /*,_N : slength(string(_Strig))*/ ),

    _Snew : trigsimp(_Strig),            /* sin^2+cos^2->1, sec -> 1/cos and tan(2x) -> sin(2x)/cos(2x)*/
    dprint(6,"3.2 trigsimp,S=",_Snew),
    _Nnew : slength(string(_Snew)),
    dprint(6,"   N,Nnew = ",_N," ",_Nnew),  
    if _Nnew < 25*_N then (_Strig : _Snew /*, _N : slength(string(_Strig))*/ ),

    _Snew : expand(trigexpand(_Strig)),  /* writes sin(2x) as 2sin(x)cos(x) and sin(a+b) as sin(a)cos(b)+cos(a)sin(b) and sin(3x) as 3cos^2sin - sin^2*/
    dprint(6,"3.3 trigexpand,S=",grind(_Snew)),
    _Nnew : slength(string(_Snew)),
    dprint(6,"   N,Nnew = ",_N," ",_Nnew),  
    if _Nnew < 25*_N then (_Strig : _Snew /*, _N : slength(string(_Strig))*/),

    _Snew : trigsimp(_Strig),            /* trigsimp again, this improves the probability of success of simplification! */
    dprint(6,"3.4 trigsimp,S=",grind(_Snew)),
    _Nnew : slength(string(_Snew)),
    dprint(6,"   N,Nnew = ",_N," ",_Nnew),  
    if _Nnew < 25*_N then (_Strig : _Snew /*, _N : slength(string(_Strig))*/),


    _Snew : trigreduce(_Strig),  /* this simplifies new sin^2,cos^2 to sin(2x),cos(2x) terms, for the backtransformation of the symmetries of kamke 1.80 for instance */
    dprint(6,"3.5 trigreduce,S=",grind(_Snew)),
    _Nnew : slength(string(_Snew)),
    dprint(6,"trig,   N,Nnew = ",_N," ",_Nnew),  
    if _Nnew < 25*_N then (_Strig : _Snew /*, _N : slength(string(_Strig))*/),


    if _isimaginary=false and not freeof(%i,_Snew) then (
      dprint(2,"simplification introduced complex numers, ignoring trig-simplification")
    ) else if _Nnew < 2.0*_N then (_S : _Strig, _N : slength(string(_S)))
    /*if _Nnew < 1.4*_N then (_S : _Snew, _N : slength(string(_S)))*/
  ),
  dprint(6,"3. S=",_S),

  /*if sublist(_ratvars,lambda([_i],not freeof(asinh,acosh,_i) and not freeof(_x,_y,_i))) # [] then (*/
  if sublist(_ratvars,lambda([_i],not freeof(asinh,acosh,_i) )) # [] then (
    _Snew : logarc(_S),
    _Nnew : slength(string(_Snew)),
    /*dprint(5,"trigh   N,Nnew = ",_N," ",_Nnew),  */
    if _Nnew < 2.0*_N then (_S : _Snew, _N : slength(string(_S)))
  ),

  /* important for kamke 68,69,70,71*/
  if sublist(_ratvars,lambda([_i],not freeof(sqrt,dispform(_i,all)) )) # [] then (
    _oldradexpand:radexpand,
    radexpand:all,
    _Snew : rootscontract(_S),
    radexpand:_oldradexpand,
    dprint(6,"sqrt: ",_Snew),
    dprint(6,"sqrt   N,Nnew = ",_N," ",_Nnew),  
    if _isimaginary=false and not freeof(%i,_Snew) then (
      dprint(2,"simplification introduced complex numers, ignoring sqrt-simplification")
    ) else 
    
    if _Nnew < 2.0*_N then (_S : _Snew, _N : slength(string(_S)))
  ),


  /*if sublist(_ratvars,lambda([_i],not freeof(%e,log,_i) and not freeof(_x,_y,_i))) # [] then (*/
  if sublist(_ratvars,lambda([_i],not freeof(%e,log,_i) )) # [] then (
    dprint(6,"log, exp"),
    _Snew : logcontract(radcan(_S)),
    dprint(6,"radcan: ",_Snew),
    
    _Nnew : slength(string(_Snew)),
    dprint(6,"log   N,Nnew = ",_N," ",_Nnew),  
    
    if _isimaginary=false and not freeof(%i,_Snew) then (
      dprint(2,"simplification introduced complex numers, ignoring radcan-simplification")
    ) else  
    if (_Nnew < 2.0*_N) then (_S : _Snew, _N : slength(string(_S)))
  ),

  /* should we always return a factored result? */

  dprint(5,"simplify::END, S=",grind(_S)),
  if freeof('integrate,_S) then (use_pdiff:true, _S:ev(_S,diff)),
  dprint(5,"returning from simplify"),
  
  /* the main criteria are that the simplified expression should have a smaller termcount */
  if (_N<2.0*_N0) then return(_S) else return(_S0)

)$
/* ************************************************************************** */

/* simplify the symmetries X=[xi,eta] */
simplifySymmetry(_X,_x,_y):=block([_xi,_eta,_gcd,_absP,_absQ],

    _xi : _X[1],
    _eta : _X[2],
    dprint(3,"simplifySymmetry: initial symmetry [xi,eta]=[",_xi,",",_eta,"]"),

    /* if the symmetries have the form [%i*xi,%i*eta] or even [a*xi,a*eta] then we can simplify to [xi,eta] */
    /* so we check if xi,eta have a common constant factor */
    _xi : simplify(_xi),
    _eta : simplify(_eta),
    dprint(4,"simplified [xi,eta]=[",_xi,",",_eta,"]"),
    _gcd : greatest_constant_divisor(_xi,_eta,[_x,_y]),
    dprint(4,"simplifying symmetries, found a common constant: ",_gcd),
    /* we need ratsimp because greatest_common_divisor can return the gcd of a different X*/
    if (_gcd#0) then (_xi:ratsimp(_xi/_gcd),_eta:ratsimp(_eta/_gcd)),
    dprint(4,"simplified symmetry (1) : [xi,eta]=[",_xi,",",_eta,"]"),

    /* ************************************** */
    /* often as symmetry with abs(x) has also a symmetry with x */
    /* we are simplifying abs here. abs leads to huge delays in obtaining the solution. */
    /* when the gcd is of the form P/Q, we assume that either P or Q has abs(_x)*/
    _gcd : gcd(_xi,_eta),
    dprint(4,"simplifying symmetries, found a common constant: ",_gcd),
    if not freeof(abs,_gcd) then (
        _absP : num(_gcd),
        _absQ : denom(_gcd),
        if not atom(_absP) then (
          if op(_absP)=abs then (dprint(4,"case 1"),_xi : simplify( _xi*_absP/args(_absP)[1]), _eta : simplify(_eta*_absP/args(_absP)[1])))
        else 
        if not atom(_absQ) then (
          if op(_absQ)=abs then (dprint(4,"case 2"),_xi : simplify( _xi*args(_absQ)[1]/_absQ), _eta : simplify(_eta*args(_absQ)[1]/_absQ)))
        else
        if freeof(abs(_x),_absP) then (dprint(4,"case 3"),_xi:simplify(_xi*_x/abs(_x)),_eta:_eta*_x/abs(_x))
        else
        if freeof(abs(_x),_absQ) then (dprint(4,"case 4"),_xi:simplify(_xi*abs(_x)/_x),_eta:_eta*abs(_x)/_x),
      dprint(4,"simplified symmetry (2) [xi,eta]=[",_xi,",",_eta,"]")
     ),
     /* second check if the first one fails */
     if not freeof(abs(_x),_xi) and not freeof(abs(_x),_eta) then (
       [_xi,_eta] : subst(abs(_x)=_x,[_xi,_eta])
       /*isit : is(checkSymmetries([_xi,_eta],_ode_depform,_y,_x)=0), */
       /*if isit then _X:[_xi,_eta]*/
     ),
     /* second check if the first one fails */
     if not freeof(abs(_y),_X[1]) and not freeof(abs(_y),_X[2]) then (
       [_xi,_eta] : subst(abs(_y)=_y,[_xi,_eta])
       /*isit : is(checkSymmetries([_xi,_eta],_ode_depform,_y,_x)=0), */
       /*if isit then _X:[_xi,_eta]*/
     ),

  dprint(3,"simplifySymmetry: simplified symmetry [xi,eta]=[",_xi,",",_eta,"]"),
  return([_xi,_eta])
)$