share/contrib/diffequations/ode1_clairaut.mac

   1 /* ode1_clairaut.mac
   2
   3   Solve Clairaut's equation f(x*y'-y)=g(y')
   4
   5   References:
   6
   7   Daniel Zwillinger, Handbook of Differential Equations, 3rd ed
   8   Academic Press, (1997), pp 216-218
   9
  10
  11   Copyright (C) 2004 David Billinghurst
  12
  13   This program is free software; you can redistribute it and/or modify
  14   it under the terms of the GNU General Public License as published by
  15   the Free Software Foundation; either version 2 of the License, or
  16   (at your option) any later version.
  17
  18   This program is distributed in the hope that it will be useful,
  19   but WITHOUT ANY WARRANTY; without even the implied warranty of
  20   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21   GNU General Public License for more details.
  22
  23   You should have received a copy of the GNU General Public License
  24   along with this program; if not, write to the Free Software
  25   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  26 */
  27
  28 put('ode1_clairaut,001,'version)$
  29
  30 /* Attempt to separate expression exp into f(x)+g(y)
  31    Returns [f(x),g(y)] is this is possible
  32            false otherwise
  33    Adapted from SEPARABLE() in ode2.mac */
  34
  35 plusseparable(eq,x,y) := block(
  36   [u,xpart:[],ypart:[],flag:false,inflag:true],
  37   eq:expand(eq),
  38   if atom(eq) or not(inpart(eq,0)="+") then eq: [eq],
  39   for u in eq do
  40     if freeof(x,u) then
  41       ypart: cons(u,ypart)
  42     else if freeof(y,u) then
  43       xpart: cons(u,xpart)
  44     else
  45       return(flag: true),
  46   if flag = true then return(false),
  47   if xpart = [] then xpart: 0 else xpart: apply("+",xpart),
  48   if ypart = [] then ypart: 0 else ypart: apply("+",ypart),
  49   return([xpart,ypart])
  50 )$
  51
  52 /* If we can write eqn as f(x*y'-y)=g(y') then the solution
  53    is given implicitly by f(x*%c-y)=g(%c)
  54
  55    A singular solution may exist
  56
  57  */
  58 ode1_clairaut(eq,y,x) := block(
  59   [ans,%c,de,%f,%g,%p,%r,_t,sep],
  60   /* substitute %p=y', %r=x*y'-y */
  61   de: expand(lhs(eq)-rhs(eq)),
  62   de: subst(%p,'diff(y,x),de),
  63   de: ratsimp(subst((%r+y)/%p,x,de)),
  64   if not(freeof(x,y,de)) then (
  65     ode_disp2("     Expression not free of x and y: ",de),
  66     return(false)
  67   ),
  68   sep:plusseparable(de,%r,%p),
  69   if is(sep=false) then
  70     /* Perhaps de is const+f(%p)*%r */
  71     sep:plusseparable(ratsimp(de/%r),%r,%p),
  72   if is(sep#false) then (
  73     %f: sep[1],
  74     %g:-sep[2]
  75   )
  76   else (
  77     ode_disp2("     Expression free of x and y but can't write as f(%r)=g(%p): ",de),
  78     /* Don't give up yet */
  79     _t:solve(de,%r),
  80     if _t=[] then return(false),
  81     ode_disp2("     Solving for %r: ",_t),
  82     %f:%r,
  83     %g:rhs(_t[1]),
  84     if not(freeof(%r,%g)) then return(false)
  85   ),
  86   /* Next line added to avoid testsuite failures 2003-01-02
  87      Guessed the fix and should analyse it further */
  88   if ( is(%f=0) ) then return(false),
  89   ode_disp2("     %f: ",%f),
  90   ode_disp2("     %g: ",%g),
  91   method: 'clairaut,
  92   ans:subst(x*%c-y,%r,%f)=subst(%c,%p,%g),
  93   /* Try and solve for y */
  94   _t:solve(ans,y),
  95   if length(_t)=0 then ans:[ans] else ans:_t,
  96   ans:append(ans,ode1_clairaut_singular(eq,y,x)),
  97   return(ans)
  98 )$
  99
 100 /* Is there a singular solution to Clairaut equation eq:f(x*p-y)-g(p)=0 ?
 101
 102      Differentiate eq wrt x to get 'diff(y,x,2)*eq2=0
 103      If possible, eliminate 'diff(y,x) between eq and eq2.
 104      else return parametric solution in variable %t
 105 */
 106
 107 ode1_clairaut_singular(eq,y,x) := block(
 108   [expr,eq2,%t,ans,u],
 109   eq:lhs(eq)-rhs(eq),
 110   expr:subst(u,y,eq),
 111   depends(u,x),
 112   expr:diff(expr,x),
 113   remove(u,dependency),
 114   expr:subst(y,u,expr),
 115   expr:expand(expr),
 116   eq2:ratsimp(ratcoeff(expr,'diff(y,x,2))),
 117   if not(freeof('diff(y,x,2),eq2)) then return([]),
 118   expr: ratsimp(expr-eq2*'diff(y,x,2)),
 119   if not(is(expr=0)) then return([]),
 120   /* Try and eliminate t */
 121   eq:subst(%t,'diff(y,x),eq),
 122   eq2:subst(%t,'diff(y,x),eq2),
 123   ans:eliminate([eq,eq2],[%t]),
 124   if ( not(freeof(%t,ans)) or (ans=1) ) then return([[eq=0,eq2=0]]),
 125   return(solve(ans,y))
 126 )$