share/contrib/diffequations/ode1_factor.mac

   1 /* Solution of first order ODEs by factoring
   2
   3   References:
   4
   5   Daniel Zwillinger, Handbook of Differential Equations, 3rd ed
   6   Academic Press, (1997), pp 265-266
   7
   8
   9   Copyright (C) 2004 David Billinghurst
  10
  11   This program is free software; you can redistribute it and/or modify
  12   it under the terms of the GNU General Public License as published by
  13   the Free Software Foundation; either version 2 of the License, or
  14   (at your option) any later version.
  15
  16   This program is distributed in the hope that it will be useful,
  17   but WITHOUT ANY WARRANTY; without even the implied warranty of
  18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19   GNU General Public License for more details.
  20
  21   You should have received a copy of the GNU General Public License
  22   along with this program; if not, write to the Free Software
  23   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  24 */
  25
  26 put('ode1_factor,001,'version)$
  27
  28 ode1_factor(eq,y,x):= block(
  29   [de,%p,ans:[],s,inflag:true,powflag:true,u,factors:[]],
  30
  31   ode_disp("    in ode1_factor"),
  32
  33   /* substitute %p=y' */
  34   de: expand(lhs(eq)-rhs(eq)),
  35   de: subst(%p,'diff(y,x),de),
  36
  37   /* At present restricted to simple cases such as p^2+5p+6=0
  38      What about repeated factors? */
  39
  40   /* if ( not(freeof(x,y,de)) ) then return(false), */
  41
  42   de:factor(de),
  43   /* fail if the ode wasn't simplified by factoring */
  44   if ( not(inpart(de,0)="*") or is(length(de)<2) ) then (
  45     ode_disp2("     cannot factor",de),
  46     return(false)
  47   ),
  48
  49   /* Count the factors that are DEs. There may be constant terms
  50      Some test cases failed when this was combined with the loop below */
  51   for u in de do (
  52     if not(freeof(%p,u)) then (
  53       factors:endcons(u,factors)
  54     )
  55   ),
  56
  57   if length(factors)<2 then (
  58     ode_disp2("     cannot factor",de),
  59     return(false)
  60   ),
  61   ode_disp2("    and after factoring is ",de),
  62
  63   for u in de do (
  64     /* Do not continue for higher order terms (including repeated factors) */
  65     if hipow(expand(u),%p)>1 then (
  66       ode_disp2("      unsuitable term ",u),
  67       powflag:false
  68     )
  69   ),
  70   if powflag=false then return(false),
  71
  72   /* For each factor, try and solve ODE */
  73   for u in factors do (
  74     ode_disp2("      Factor ",u),
  75     if not(freeof(%p,u)) then (
  76       s:ode2(subst('diff(y,x),%p,u)=0,y,x),
  77       ode_disp2("      has solution",s),
  78       if s#false then ans:endcons(s,ans)
  79     )
  80   ),
  81
  82   if ans#[] then (
  83     method:'factor,
  84     ans
  85   )
  86   else
  87     false
  88 )$