share/contrib/diffequations/ode1_lagrange.mac

   1 /* ode1_lagrange.mac
   2
   3    Solution of Lagrange equation y = x f(y') + g(y')
   4
   5    References:
   6
   7    Daniel Zwillinger, Handbook of Differential Equations, 3rd ed
   8    Academic Press, (1997), pp 328-331
   9
  10   Copyright (C) 2004 David Billinghurst
  11
  12   This program is free software; you can redistribute it and/or modify
  13   it under the terms of the GNU General Public License as published by
  14   the Free Software Foundation; either version 2 of the License, or
  15   (at your option) any later version.
  16
  17   This program is distributed in the hope that it will be useful,
  18   but WITHOUT ANY WARRANTY; without even the implied warranty of
  19   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20   GNU General Public License for more details.
  21
  22   You should have received a copy of the GNU General Public License
  23   along with this program; if not, write to the Free Software
  24   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  25
  26 */
  27
  28 put('ode1_lagrange,001,'version)$
  29
  30 ode1_lagrange(eq,y,x):= block(
  31   [eqn,ans1,ans:[],e,t,f,g,t2,s,p,q,%t],
  32
  33   ode_disp("In ode1_lagrange"),
  34   /* Express equation as y = x*f(p)+g(p), where p=y' */
  35   eqn:subst(p,'diff(y,x),eq),
  36   if freeof(y,eqn) then return(false),
  37   eqn:solve(eqn,y),
  38   if eqn=[] then return(false),
  39
  40   /* There may be multiple solutions */
  41   for e in eqn do (
  42     ans1:block(
  43       if lhs(e)#y then return(false),
  44       if not(freeof(y,rhs(e))) then return(false),
  45       t:rhs(e),
  46       f:ratcoeff(t,x),
  47       g:t-f*x,
  48       if not(freeof(x,y,f)) then return(false),
  49       if not(freeof(x,y,g)) then return(false),
  50       if f=p then return(false),
  51
  52       /* Transform ode */
  53       t2:'diff(x,p)=x*diff(f,p)/(p-f)+diff(g,p)/(p-f),
  54       ode_disp2("ode1_lagrange: Transform successful and new equation is ",t2),
  55       s:ode2(t2,x,p),
  56       if (s=false) then return(false),
  57       ode_disp("Transformed equation solved"),
  58       method:'lagrange,
  59
  60       t:y=x*f+g,
  61       /* Simple test to prevent failures */
  62       ode_disp("Need to eliminate p from"),
  63       ode_disp2("  Equation t: ",t),
  64       ode_disp2("  Solution s: ",s),
  65       /* Eliminate is only for polynomials.
  66          Punt that it returns 1 when called inappropriately */
  67       q:first(eliminate([t,s],[p])),
  68       ode_disp("  and after eliminating p the result is"),
  69       ode_disp2("  q: ",q),
  70       if ( freeof(p,q) and (q#1) ) then (
  71         q=0
  72       )
  73       else (
  74         /* Return the parametric equation in %t=p */
  75         [subst(%t,p,s),subst(%t,p,t)]
  76       )
  77     ),
  78     if ans1#false then ans:endcons(ans1,ans)
  79   ),
  80   if ans=[] then
  81     false
  82   else
  83     ans
  84 )$