share/contrib/diffequations/contrib_ode.mac

   1 /* contrib_ode.mac
   2
   3   Driver for contributed ODE routines
   4
   5
   6   Copyright (C) 2004,2006,2007 David Billinghurst
   7
   8   This program is free software; you can redistribute it and/or modify
   9   it under the terms of the GNU General Public License as published by
  10   the Free Software Foundation; either version 2 of the License, or
  11   (at your option) any later version.
  12
  13   This program is distributed in the hope that it will be useful,
  14   but WITHOUT ANY WARRANTY; without even the implied warranty of
  15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16   GNU General Public License for more details.
  17
  18   You should have received a copy of the GNU General Public License
  19   along with this program; if not, write to the Free Software
  20   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  21 */
  22
  23 if get('ode1_nonlinear,'version)=false then
  24   load('ode1_nonlinear)$
  25
  26 if get('odelin,'version)=false then
  27   load('odelin)$
  28
  29 if get('kovacic,'version)=false then
  30   load('kovacicODE)$
  31
  32 /* Set default debug output for kovacicODE */
  33 DEBUGFLAG:0;
  34
  35 load('opsubst)$
  36
  37 contrib_ode(eqn,y,x) := block(
  38   [soln,degree:derivdegree(eqn,y,x)],
  39
  40   ode_disp("-> ode_contrib"),
  41
  42   /* Try ode2 if first or second order*/
  43   if (degree < 1) then (
  44     ode_disp("   Degree of equation less than 1"),
  45     false
  46   )
  47   else if (degree > 2) then (
  48     ode_disp("   Degree of equation greater then 2"),
  49     false
  50   )
  51   else if (degree=1) then (
  52     ode_disp("   First order equation"),
  53     ode_disp("   -> ode2"),
  54     soln:ode2(eqn,y,x),
  55     if (soln#false) then (
  56       ode_disp("    Successful"),
  57       return([soln])
  58     )
  59     else (
  60       return(ode1_nonlinear(eqn,y,x))
  61     )
  62   )
  63   else if (degree=2) then (
  64     return(ode2_lin(eqn,y,x))
  65   )
  66   else (
  67     error("contrib_ode: This case cannot occur"),
  68     false
  69   )
  70 )$
  71
  72 /* Try methods for linear 2nd order ODEs.
  73    eqn is a 2nd order ODE.
  74    FIXME: Remove multiple checks for linearity
  75 */
  76 ode2_lin(eqn,y,x) := block(
  77   [de,a1,a2,a3,a4,soln,%k1,%k2],
  78   ode_disp("   Second order equation"),
  79   de: desimp(lhs(eqn)-rhs(eqn)),
  80   a1: coeff(de,'diff(y,x,2)),
  81   a2: coeff(de,'diff(y,x)),
  82   a3: coeff(de,y),
  83   a4: expand(de - a1*'diff(y,x,2) - a2*'diff(y,x) - a3*y),
  84   if (freeof(y,[a1,a2,a3,a4])) then (
  85     ode_disp("     Linear 2nd order ODE"),
  86     /* Try the ode2 routines first */
  87     ode_disp("    -> ode2"),
  88     soln:ode2(eqn,y,x),
  89     /* ode2 can return [] on failure */
  90     if (soln#false and soln#[]) then (
  91       ode_disp("     Successful"),
  92       [soln]
  93     ) else if is(a4=0) then (
  94       /* Now try the routines in odelin */
  95       ode_disp("    -> odelin"),
  96       soln:odelin(eqn,y,x),
  97       if (soln#false) then (
  98         ode_disp("     Successful"),
  99         method:'odelin,
 100         ode_disp(soln),
 101         /* odelin returns a fundamental soln set */
 102         soln:listify(soln),
 103         [y=%k1*soln[1]+%k2*soln[2]]
 104       ) else (
 105         false
 106       )
 107     ) else (
 108       /* Try kovacic algorithm */
 109       ode_disp("    -> kovacicODE"),
 110       soln:kovacicODE(de=0,y,x),
 111       if (soln#false) then (
 112         ode_disp("     Successful"),
 113         method:'kovacic,
 114         ode_disp(soln),
 115         soln
 116       ) else (
 117         false
 118       )
 119     )
 120   ) else (
 121     ode_disp("     Non-linear 2nd order ODE"),
 122     false
 123   )
 124 )$
 125
 126 ode_disp(msg) := block(
 127   if get('contrib_ode,'verbose) then print(msg)
 128 )$
 129
 130 ode_disp2(msg,e) := block(
 131   if get('contrib_ode,'verbose) then print(msg,e)
 132 )$
 133
 134 /* recurse through expression eq looking for a derivative
 135    Return the derivative if found, false otherwise
 136  */
 137 ode_deriv(eq) := block(
 138   [u,v:false],
 139   eq:expand(eq),
 140   if mapatom(eq) then return(false),
 141   if operatorp(eq,nounify('diff)) then return(eq),
 142   for u in eq do
 143     (if (v:ode_deriv(u))#false then return(true)),
 144   v
 145 )$
 146
 147 /* Check the solution of an ode. */
 148 ode_check(e_,a_) := block(
 149   [deriv,x_,y_,diff_e,e,u,v],
 150
 151   deriv:ode_deriv(e_),
 152   if deriv=false then (
 153     print("Not a differential equation"),
 154     return(false)
 155   ),
 156   x_:part(deriv,2),  /* Independent variable */
 157   y_:part(deriv,1),  /* Dependent variable */
 158
 159   if not(mapatom(x_)) then (
 160     print("Independent variable ",x_," not an atom"),
 161     return(false)
 162   ),
 163   if not(mapatom(y_)) then (
 164     print("Dependent variable ",y_," not an atom"),
 165     return(false)
 166   ),
 167
 168   /* Is it a simple solution */
 169   if lhs(a_)=y_ then (
 170     return(ode_check_tidy(ratsimp(ev(lhs(e_)-rhs(e_),a_,diff))))
 171   ),
 172
 173   /* Is a_ a parametric solution */
 174   if length(a_)=2 and lhs(a_[1])=x_ and lhs(a_[2])=y_ then (
 175     return(ode_check_tidy(fullratsimp(ode_check_parametric(e_,a_,x_,y_))))
 176   ),
 177
 178   /* Is a_ a parametric solution in disguise?
 179        For example, Kamke 1.555 returns
 180        [(-y)+%t*x+sqrt(%t^2+1)=0, (sqrt(%t^2+1)*x+%t)/sqrt(%t^2+1)=0]
 181    */
 182   if length(a_)=2 and not freeof(%t,a_) then (
 183     s:solve(a_,[x_,y_]),
 184     /* Check to see if the solution s is a valid
 185        FIXME: Have assumed a single solution
 186      */
 187     if s#[]
 188         and freeof(x_,y_,rhs(s[1][1]))
 189         and freeof(x_,y_,rhs(s[1][2])) then (
 190       return(ode_check_tidy(fullratsimp(ode_check_parametric(e_,s[1],x_,y_))))
 191     )
 192   ),
 193
 194   /* Must be an implicit solution */
 195   diff_e:diff(subst(v(x),y_,a_),x_),
 196   diff_e:subst(y_,v(x),diff_e),
 197   u:solve(diff_e,'diff(y_,x_)),
 198   if u=[] then (
 199     print("Problem - dy/dx is ",u),
 200     return(false)
 201   ),
 202   u:first(u),
 203   if ( lhs(u)#'diff(y_,x_) or not(freeof('diff(y_,x_),rhs(u))) ) then (
 204     print("Problem - dy/dx is ",u),
 205     return(false)
 206   ),
 207   ode_check_tidy(fullratsimp(ev(lhs(e_)-rhs(e_),u)))
 208 )$
 209
 210 /* Attempt to simplify the result from ode_check() */
 211 ode_check_tidy(e) := block(local(f),
 212   for f in ['bessel_simplify, 'expintegral_e_simplify,'ode_trig_simp]
 213   do
 214     (if (f(e)=0) then return(e:0)),
 215   return(expand(e))
 216 )$
 217
 218 /*  This Bessel function simplification code can be improved,
 219     but is does work well on the contrib_ode testsuite.
 220 */
 221 bessel_simplify(e) := block(local(f),
 222  for f in ['besksimp, 'besisimp, 'besysimp, 'besjsimp,
 223            'hankel1simp, 'hankel2simp, 'struvehsimp, 'struvelsimp]
 224  do
 225   e:f(e),
 226  return(e))$
 227
 228 freeof_bessel(e):=freeof(bessel_j,bessel_y,bessel_i,bessel_k,e);
 229
 230 /* Simplify expressions containing bessel_F, where F in {j,y,i,k}.
 231    Repeatedly find the highest order term F and substitute Fsub for F
 232
 233    When F is a Bessel function of order n and argument x, then
 234    Fsub is contains Bessel functions of order n-1 and n-2.
 235
 236    We separate the set of arguments into classes that have identical
 237    second arguments and orders that differ by integers, then work on each
 238    class until all the orders differ by less than 2.
 239 */
 240 besFsimp(ex,F,Fsub) := block([arglist,ords,ord2,max,min,s,t,u,m,n,x,e:ex],
 241   if freeof(F,ex) then return(ex),
 242   argset:setify(gatherargs(ex,F)),
 243   /*  partition arguments into sets with same second argument and
 244       orders that differ by an integer */
 245   argset:equiv_classes(argset,lambda([x,y],
 246     (second(x)=second(y)) and integerp(ratsimp(first(x)-first(y))))),
 247
 248   /* Iterate over each equivalence class */
 249
 250   for s in argset do (
 251     x:second(first(s)),  /* The common argument for s */
 252     ords:map('first,s),  /* List of orders */
 253     t:first(ords),          /* One element of the list of orders */
 254     ord2:map(lambda([m],ratsimp(m-t)),ords), /* List of differences */
 255     max:lmax(ord2),      /* maximum order is t+max */
 256     min:lmin(ord2),      /* minimum order is t+min */
 257     for m: max step -1 unless ratsimp(m-min)<2 do (
 258       n:subset(ords,lambda([u],is(ratsimp(u-t-m)=0))),
 259       n:if emptyp(n) then ratsimp(t+m) else first(n),
 260       e:subst(apply(Fsub,[n,x]),apply(F,[n,x]),e),
 261       e:ratsimp(e)
 262     )
 263   ),
 264   ratsimp(e)
 265 )$
 266
 267 /* Recurrence for Bessel and Hankel functions A&S 9.1.27, DLMF 10.6.1 */
 268
 269 besjsimp(ex):=besFsimp(ex,bessel_j,
 270   lambda([n,x], 2*(n-1)*bessel_j(n-1,x)/x-bessel_j(n-2,x)));
 271
 272 besysimp(ex):=besFsimp(ex,bessel_y,
 273   lambda([n,x], 2*(n-1)*bessel_y(n-1,x)/x-bessel_y(n-2,x)));
 274
 275 besisimp(ex):=besFsimp(ex,bessel_i,
 276   lambda([n,x],-2*(n-1)*bessel_i(n-1,x)/x+bessel_i(n-2,x)));
 277
 278 hankel1simp(ex):=besFsimp(ex,hankel_1,
 279   lambda([n,x], 2*(n-1)*hankel_1(n-1,x)/x-hankel_1(n-2,x)));
 280
 281 hankel2simp(ex):=besFsimp(ex,hankel_2,
 282   lambda([n,x], 2*(n-1)*hankel_2(n-1,x)/x-hankel_2(n-2,x)));
 283
 284
 285 /* bessel_k requires special treatment due to the simplification
 286    bessel_k(-n,x) => bessel_k(n,x), which upsets the recurrence
 287    relationship.  */
 288 besksimp(ex):=block(
 289   [bk:gensym(),e],
 290   e:opsubst(bk,bessel_k,ex),
 291   e:besksimp_1(e,bk),
 292   /* Recurrence relation for bessel_k */
 293   e:besFsimp(e,bk,lambda([n,x], 2*(n-1)*bk(n-1,x)/x+bk(n-2,x))),
 294   opsubst(bessel_k,bk,e)
 295 );
 296
 297 /* This is a helper function for besksimp.
 298   - the function bk has been substituted for bessel_k in ex
 299   - collect all the arguments of each occurrence of bk in e
 300   - partition arguments into sets with same argument and where
 301     either the orders m and n differ by integer or m and (-n)
 302     differ by integer (so m+n is an integer).
 303   - for the elements where  m+n is an integer, substitute
 304     bk(-n,x) for bk(n,x) in e
 305   - return e
 306 */
 307 besksimp_1(e,bk):=block(
 308   [s,argset],
 309   /* Break into equivalence classes */
 310   argset:equiv_classes(setify(gatherargs(e,bk)),
 311     lambda([x,y],
 312       (second(x)=second(y))
 313       and
 314            (integerp(ratsimp(first(x)-first(y)))
 315         or  integerp(ratsimp(first(x)+first(y)))))),
 316   /* For each equivalnce class */
 317   for s in argset do block(
 318     [x,ords,t,n],
 319     x:second(first(s)),  /* The common argument for s */
 320     ords:map('first,s),  /* Set of orders */
 321     /* t is an element of ords, If ords are numbers then t is the largest */
 322     t:first(ords), if every(numberp,ords) then t:lmax(ords),
 323     /* Find the subset that do not differ by an integer from t
 324        and change sign of the order for those orders */
 325     for n in subset(ords,lambda([i],integerp(ratsimp(t+i)))) do (
 326       e:subst(bk(-n,x),bk(n,x),e) )
 327   ),
 328   e
 329 );
 330
 331 /* Recurrence for struve_h DLMF 11.4.23 */
 332 struvehsimp(ex):=besFsimp(ex,struve_h,
 333   lambda([n,x],2^(1-n)*x^(n-1)/(sqrt(%pi)*gamma(n+1/2))
 334               +2*(n-1)*struve_h(n-1,x)/x-struve_h(n-2,x)))$
 335
 336 /* Recurrence for struve_l DLMF 11.4.25 */
 337 struvelsimp(ex):=besFsimp(ex,struve_l,
 338   lambda([n,x],-2^(1-n)*x^(n-1)/(sqrt(%pi)*gamma(n+1/2))
 339                -2*(n-1)*struve_l(n-1,x)/x+struve_l(n-2,x)))$
 340
 341 /* Simplify expressions containing exponential integral expintegral_e
 342    using the recurrence (A&S 5.1.14).
 343
 344    expintegral_e(n+1,z)
 345         = (1/n) * (exp(-z)-z*expintegral_e(n,z))      n = 1,2,3 ....
 346
 347    functions.wolfram.com indicates that it also applies for non-integer n,
 348    and this is supported by some numerical tests.
 349
 350    This is a support routine for checking solutions of differential equations
 351    in the contrib_ode testsuite.  Assumes numberp(n) but this could be relaxed,
 352    as was done for Bessel functions.  Not required for existing testsuite.
 353 */
 354
 355 expintegral_e_simplify(ex) := block( [argset],
 356   if freeof(expintegral_e,ex) then return(ex),
 357   /* Get the set of arguments. Partition into sets with orders that
 358      differ by an integer and with common second arguments.  */
 359   argset:setify(gatherargs(ex,expintegral_e)),
 360   argset:equiv_classes(argset,
 361            lambda([x,y], integerp(x[1]-y[1]) and (second(x)=second(y)))),
 362
 363   /* Iterate over each equivalence class */
 364   for s in argset do block( [z, ords],
 365     z:second(first(s)),  /* The common argument for s */
 366     ords:map('first,s),  /* List of orders */
 367     for n: lmax(ords)-1 step -1 thru lmin(ords) do (
 368       ex:subst((1/n)*(exp(-z)-z*expintegral_e(n,z)), expintegral_e(n+1,z), ex)
 369     )
 370   ),
 371   /* Note: When ratsimp was in the loop, then z within the expression
 372      was sometimes altered and subst didn't work */
 373   ratsimp(ex)
 374 );
 375
 376 ode_freeof_trig(e) :=
 377   freeof(cos,sin,tan,cot,sec,csc,cosh,sinh,tanh,coth,sech,csch,e)$
 378
 379 ode_trig_simp(e) := block(
 380   if ode_freeof_trig(e) then
 381     e
 382   else
 383     radcan(trigrat(trigsimp(e)))
 384 )$
 385
 386
 387 /* Check parametric solution of first order ode e_
 388    Solution is of form  [x = X(%t), y = Y(%t) ]
 389
 390    May have x=X(%t,y) or y=Y(%t,x) but not both
 391 */
 392 ode_check_parametric(e_, a_, x, y) := block(
 393   [X,Y,dydx,s],
 394   X:rhs(a_[1]),
 395   Y:rhs(a_[2]),
 396   if not freeof(y,X) then X:subst(Y,y,X),
 397   if not freeof(x,Y) then Y:subst(X,x,Y),
 398   dydx:diff(Y,%t)/diff(X,%t),
 399   s:subst(dydx,'diff(y,x),lhs(e_)-rhs(e_)),
 400   s:subst(X,x,s),
 401   s:subst(Y,y,s),
 402   return(s)
 403 )$