share/contrib/diffequations/ode1_lie.mac

   1 /* ode1_lie.mac
   2
   3  Solution of first order ODEs using Lie symmetry methods
   4
   5
   6   References:
   7
   8   [1] E. S. Cheb-Terrab, A. D. Roche, Symmetries and First Order
   9       ODE Patterns, Computer Physics Communications 113 (1998), p 239.
  10       (http://lie.uwaterloo.ca/papers/ode_vii.pdf)
  11
  12   [2] E. S. Cheb-Terrab, T. Koloknikov,  First Order ODEs,
  13       Symmetries and Linear Transformations, European Journal of
  14       Applied Mathematics, Vol. 14, No. 2, pp. 231-246 (2003).
  15       (http://arxiv.org/abs/math-ph/0404014)
  16       (http://lie.uwaterloo.ca/papers/ode_iv.pdf)
  17
  18   [3] G W Bluman, S C Anco, Symmetry and Integration Methods for
  19       Differential Equations, Springer, (2002)
  20
  21
  22   Copyright (C) 2004 David Billinghurst
  23
  24   This program is free software; you can redistribute it and/or modify
  25   it under the terms of the GNU General Public License as published by
  26   the Free Software Foundation; either version 2 of the License, or
  27   (at your option) any later version.
  28
  29   This program is distributed in the hope that it will be useful,
  30   but WITHOUT ANY WARRANTY; without even the implied warranty of
  31   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  32   GNU General Public License for more details.
  33
  34   You should have received a copy of the GNU General Public License
  35   along with this program; if not, write to the Free Software
  36   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  37 */
  38
  39 put('ode1_lie,001,'version)$
  40
  41 /* Solve first order ODE y' = Phi(x,y) using Lie's
  42    method of symmetries */
  43 ode1_lie(phi,y,x) := block(
  44   [sym,xi,eta,mu,_r],
  45
  46   ode_disp("In ode1_lie"),
  47
  48   /* Check that phi is free of dy/dx */
  49   if not(freeof('diff(y,x),phi)) then (
  50     ode_disp("In ode1_lie: phi not free of dy/dx"),
  51     return(false)
  52   ),
  53   /* Look for symmetry generators [xi,eta] */
  54   sym:ode1_a(phi,y,x),
  55   if (sym=false) then (
  56     ode_disp("In ode1_lie: ode1_a returned false"),
  57     return(false)
  58   ),
  59   xi:sym[1], eta:sym[2],
  60   ode_disp2("In ode1_lie: xi = ",xi),
  61   ode_disp2("In ode1_lie: eta = ",eta),
  62
  63   /* Check symmetry generators */
  64   r:symtest(phi,xi,eta,y,x),
  65   ode_disp2("In ode1_lie: symtest = ",r),
  66
  67   /* Derive integrating factor from [xi,eta] */
  68   mu:lie_integrating_factor(phi,xi,eta),
  69
  70   /* Solve ODE using itegrating factor */
  71   lie_exact(phi,mu,y,x)
  72 )$
  73
  74 /* Determine if ODE  y' = Phi(x,y)
  75    has linear symmetry of form  xi = F(x), eta = P(x)*y + Q(x)
  76
  77    If so, we can transform into an ODE with symmetry
  78    xi = F(x), eta = H(x) for some F,H using lie_FxHx().
  79
  80    We then transform the symmetry back to the original
  81    coordinate system.  If x=t and y=y(u,t) then
  82
  83    F (d/dt) + H (d/du)
  84    = F (d/dx + (dy/dt) d/dy) + H (dy/du) d/dy
  85    = F d/dx + (F (dy/dt) + H (dy/du)) d/dy
  86
  87   so that xi  = F
  88           eta = F(dy/dt) + H(dy/du)
  89
  90   Return [xi,eta] if successful
  91          false    otherwise
  92 */
  93
  94 ode1_a(phi,y,x) := block(
  95   [ %A, %Ay, %Ayy, %Ayyy, %phi_yy, %phi_yyy ],
  96   ode_disp("In ode1_a ..."),
  97   ode_disp2("     phi = ",phi),
  98 /*
  99   Using [2], eq (13),(14)
 100 */
 101   /* Confirm that diff(phi,y,3)#0 */
 102   %phi_yy: ratsimp(diff(phi,y,2)),
 103   %phi_yyy: ratsimp(diff(%phi_yy,y)),
 104   if is(%phi_yyy=0) then
 105   (
 106     /* %phi_yyy=0 - Not covered in [2] */
 107     ode_disp("     %phi_yyy = 0"),
 108     return(lie_FxHx(phi,y,x))
 109   ),
 110   %A: ratsimp(%phi_yy/%phi_yyy),
 111   %Ay: ratsimp(diff(%A,y)),
 112   %Ax: ratsimp(diff(%A,x)),
 113
 114   if is(%Ay=0) then
 115     /* Case 1: diff(A,y)=0
 116        Change variables using y=A*u (15) */
 117     block(
 118       [u],
 119       ode_disp(" ..... Case 1 true"),
 120       ode_disp("       Change variables using y=A.u"),
 121       ode_disp2("       where A is ",%A),
 122       /* Here neweq is RHS only */
 123       neweq: ratsimp((subst(%A*u,y,phi)-%Ax*u)/%A),
 124       ode_disp2("       transformed equation is ",neweq),
 125
 126       /* Symmetry generators of u equation */
 127       ode_disp("       Finding symmetries of transformed equation"),
 128       gen:radcan(lie_FxHx(neweq,u,x)),
 129      ode_disp2("       [xi,eta] = ",gen),
 130
 131       if (gen=false) then
 132         false
 133       else block(
 134         /* transform [F,H] to [xi,eta] */
 135         [_F:gen[1],_H:gen[2],dydt:%Ax*y/%A,dudt:%A],
 136         [_F, ratsimp(_F*dydt+_H*dudt)]
 137       )
 138     )
 139   else if is((%Ayy:ratsimp(diff(%Ay,y)))=0) then
 140     /* Case 2: diff(A,y,2)=0 and diff(A,y)#0
 141        Change variables using u=ln(A) (17) */
 142     block(
 143       [neweq,sub,u],
 144       ode_disp(" ..... Case 2 true"),
 145       ode_disp2("       A is ",%A),
 146       ode_disp2("       ln(A) is ",log(%A)),
 147       sub: radcan(solve(u=log(%A),y)[1]), /* Robust enough? */
 148       ode_disp2("       ",sub),
 149       neweq: ratsimp(ev(phi,sub)-diff(rhs(sub),x))/diff(rhs(sub),u),
 150       neweq:radcan(neweq),
 151       ode_disp2("       neweq is ",neweq),
 152
 153       /* Symmetry generators of u equation */
 154       ode_disp("       Finding symmetries of transformed equation"),
 155       gen:lie_FxHx(neweq,u,x),
 156       if (gen=false) then
 157          false
 158       else block(
 159         /* transform [F,H] to [xi,eta] */
 160         [_F:gen[1],_H:gen[2],dydt,dudt,_xi,_eta],
 161         _xi:_F,
 162         dydt:diff(rhs(sub),x),
 163         dudt:diff(rhs(sub),u),
 164         _eta: _F*dydt+_H*dudt,
 165         /* Now eliminate u from eta */
 166         _eta:subst(log(%A),u,_eta),
 167         _eta:radcan(_eta),
 168         [ _xi, _eta ]
 169       )
 170     )
 171   else
 172     block(
 173       [%Axy,%eye,%Iy,%p,%px],
 174       %Axy: ratsimp(diff(%Ay,x)),
 175       %eye: ratsimp(%Axy/%Ayy),
 176       if %eye # 0 and linear2(y,%eye) then
 177         /* Case 3: diff(A,y,2)#0 and I=Axx/Ayy linear in y
 178            p=exp(integrate(Iy,x))
 179            Substitute y = u/p */
 180         block(
 181           [ _F, _H, _xi, _eta, u ],
 182           ode_disp(" ..... Case 3 true"),
 183           if get('contrib_ode,'verbose) then
 184             print("       A is ",%A),
 185           %Iy: diff(%eye,y),
 186           %Iy:ratsimp(%Iy),
 187           %Iy:trigsimp(%Iy),
 188           if not(freeof(y,%Iy)) then (
 189             print("I_y not free of y in ode1_a - impossible"),
 190             print("I_y = ",%Iy),
 191             return(false)
 192           ),
 193           %p: radcan(exp(integrate(%Iy,x))),
 194           ode_disp2("       p is ",%p),
 195           %px:diff(%p,x),
 196           neweq: ratsimp(%p*subst(u/%p,y,phi)+%px*u/%p),
 197           ode_disp2("       neweq is ",neweq),
 198
 199           /* Symmetry generators of u equation */
 200           ode_disp("       Finding symmetries of transformed equation"),
 201           gen:lie_FxHx(neweq,u,x),
 202           ode_disp2("       [xi,eta] = ",gen),
 203
 204           if (gen=false) then
 205             false
 206           else (
 207           /* Now we need to transform xi and eta back to original
 208              equation.  Have
 209                   F(x).(d/dx) + H(x).(d/du)
 210                   y = u/%p                    */
 211             _F:radcan(gen[1]),
 212             _H:radcan(gen[2]),
 213             _xi:_F,
 214             ode_disp2("  _xi = ",_xi),
 215             _eta: ratsimp(( _H - _F*%px*y )/%p),
 216             ode_disp2("  _eta = ",_eta),
 217             [ _xi, _eta ]
 218           )
 219         )
 220       else
 221         /* No such symmetry exists */
 222         false
 223     )
 224 )$
 225
 226
 227 /* Determine if ODE  y' = Phi(x,y)
 228    has linear symmetry of form  xi = F(x), eta = H(x)
 229
 230    This uses the notation from [1] E. S. Cheb-Terrab, A. D. Roche,
 231    Symmetries and First Order ODE Patterns
 232 */
 233 lie_FxHx(Phi,y,x) := block(
 234   [Phi_y, Phi_yy, Q, Q_y,_a,_b,_c,gen],
 235
 236   ode_disp2("In lie_FxHx with Phi:",Phi),
 237   Phi_y: ratsimp(diff(Phi,y)),
 238   Phi_yy: ratsimp(diff(Phi_y,y)),
 239   ode_disp2("Phi_yy = ",Phi_yy),
 240   if (Phi_yy=0) then (
 241     /* Phi must be linear in y, say Phi = _a*y + _b */
 242     ode_disp("    Phi_yy=0 so Phi must be linear in y"),
 243     _c: bothcoef(expand(Phi),y),
 244     _a: first(_c),
 245     _b: expand(second(_c)),
 246    ode_disp2("    _a: ",_a),
 247    ode_disp2("    _b: ",_b),
 248    if not(freeof(y,_a)) then error("lie_FxHx: _a not free of ",y),
 249    if not(freeof(y,_b)) then error("lie_FxHx: _b not free of ",y),
 250    if (_a=0) then (
 251      gen:lie_symgen_separable(_b,1,y,x)
 252    ) else (
 253      gen:lie_symgen_linear(_a,_b,y,x)
 254    ),
 255    ode_disp2("    xi: ",gen[1]),
 256    ode_disp2("    eta: ",gen[2]),
 257    return(gen)
 258   ),
 259   Q:ratsimp(Phi_y/Phi_yy),  /* (39) */
 260   Q_y:ratsimp(diff(Q,y)),
 261
 262   if (Q_y#0) then   /* Case 1: Q_y # 0 */
 263     block(
 264       [Phi_x,Upsilon,Upsilon_x,Upsilon_y,_t,_f,_h],
 265       Upsilon: ratsimp(diff(Q,x)/Q_y), /* (40) */
 266       Upsilon_y: ratsimp(diff(Upsilon,y)),
 267       if is(Upsilon_y#0) then return(false),
 268       Phi_x: diff(Phi,x),
 269       Upsilon_x: diff(Upsilon,x),
 270       _t: ratsimp((Upsilon*Phi_y-Upsilon_x-Phi_x)/(Phi+Upsilon)),
 271       if is(ratsimp(diff(_t,y))#0) then return(false),
 272       _f: ratsimp(exp(integrate(_t,x))), /* (43) */
 273       _h: ratsimp(-Upsilon*_f),
 274       [_f,_h]
 275     )
 276   else                 /* Case 2: Q_y = 0 */
 277     block(
 278       [_a,_b,_c,_f,_h],
 279       /* In this case Q is constant and Phi must have the form
 280              y' = Phi(x,y) = A(x) + B(x)*exp(y/Q)
 281          and F(x), H(x) are simple functions of A and B */
 282       _c: bothcoef(expand(Phi),exp(y/Q)),
 283       _a: second(_c),
 284       _b: first(_c),
 285       if (_b=0) then (
 286         print("Impossible case in lie_FxHx.  This cannot happen."),
 287         error("Unknown error in lie_FxHx") ),
 288       _f: exp(-integrate(_a/Q,x))/_b,
 289       _h: _a*_f,
 290       [_f,_h]
 291     )
 292 )$
 293
 294 /* Check if [ xi(x,y), eta(x,y) ] are symmetry generators
 295    of first order ode y'=phi(x,y)
 296
 297    Returns 0 when the symmetry generators are OK.  If the result is
 298    non-zero, then further simplification may reduce it to 0.
 299 */
 300 symtest(phi,xi,eta,y,x) :=
 301   ratsimp( diff(eta,x) + (diff(eta,y)-diff(xi,x))*phi
 302       - diff(xi,y)*phi^2 - xi*diff(phi,x) - eta*diff(phi,y) )$
 303
 304 /* Determine integrating factor of first order ode y'=phi(x,y)
 305    given symmetry generators [ xi(x,y), eta(x,y) ]
 306 */
 307 lie_integrating_factor(phi,xi,eta) :=
 308   1/ratsimp(eta-xi*phi)$
 309
 310 /* Solve DE given integrating factor */
 311 lie_exact(phi,mu,y,x) := block(
 312   [a,%c],
 313   a: integrate(mu*phi,x),
 314   method:'lie,
 315   intfactor: mu,
 316   ratsimp(-a+integrate(ratsimp(mu+diff(a,y)),y))=%c
 317 )$
 318
 319
 320
 321 /* from ode2.mac
 322    EXACT(M,N):=BLOCK([A,ynew,%c],
 323    INTFACTOR: SUBST(YOLD,YNEW,%q%),
 324    A: INTEGRATE(RATSIMP(M),X),
 325    METHOD: 'EXACT,
 326    RETURN(RATSIMP(A + INTEGRATE(RATSIMP(N-DIFF(A,Y)),Y)) = %C))$
 327 */
 328
 329 /* Find  symmetry generators [ xi, eta ] of separable
 330    first order ode 'diff(y,x) = phi(x,y) = f(x)*g(y)
 331
 332    Can find generator [ xi(x), 0 ] = [ 1/f, 0 ]
 333   */
 334 lie_symgen_separable(f,g,y,x) := block(return([1/f,0]))$
 335
 336 /* Find symmetry generators [ xi, eta ] of linear
 337    first order ode 'diff(y,x) = phi(x,y) = f(x)*y + g(x)
 338
 339    Can find generator [ xi(x), eta(x) ]
 340  */
 341 lie_symgen_linear(f,g,y,x) := block(
 342   [xi, eta, expintf],
 343   ode_disp("    In lie_symgen_linear"),
 344   xi: 1/f,
 345   ode_disp2("     xi: ",xi),
 346   expintf: radcan(exp(integrate(f,x))),
 347   ode_disp2("     expintf: ",expintf),
 348
 349   eta: expintf*(integrate((f*diff(g,x)-diff(f,x)*g)/(f^2*expintf),x)+%c),
 350   ode_disp2("     eta: ",eta),
 351   ode_disp("     and after eta:radcan(eta)"),
 352   eta:radcan(eta),
 353   ode_disp2("     eta: ",eta),
 354   [xi,eta]
 355 )$