doc/info/contrib_ode.texi

   1 @menu
   2 * Introduction to contrib_ode::
   3 * Functions and Variables for contrib_ode::
   4 * Possible improvements to contrib_ode::
   5 * Test cases for contrib_ode::
   6 * References for contrib_ode::
   7 @end menu
   8
   9 @node Introduction to contrib_ode, Functions and Variables for contrib_ode, Package contrib_ode, Package contrib_ode
  10
  11 @section Introduction to contrib_ode
  12
  13 Maxima's ordinary differential equation (ODE) solver @mref{ode2} solves
  14 elementary linear ODEs of first and second order.  The function
  15 @mref{contrib_ode} extends @code{ode2} with additional methods for linear
  16 and non-linear first order ODEs and linear homogeneous second order ODEs.
  17 The code is still under development and the calling sequence may change
  18 in future releases.  Once the code has stabilized it may be
  19 moved from the contrib directory and integrated into Maxima.
  20
  21 This package must be loaded with the command @code{load("contrib_ode")}
  22 before use.
  23
  24 The calling convention for @code{contrib_ode} is identical to @code{ode2}.
  25 It takes
  26 three arguments: an ODE (only the left hand side need be given if the
  27 right hand side is 0), the dependent variable, and the independent
  28 variable.  When successful, it returns a list of solutions.
  29
  30 The form of the solution differs from @code{ode2}.
  31 As non-linear equations can have multiple solutions,
  32 @code{contrib_ode} returns a list of solutions.  Each  solution can
  33 have a number of forms:
  34 @itemize @bullet
  35 @item
  36 an explicit solution for the dependent variable,
  37
  38 @item
  39 an implicit solution for the dependent variable,
  40
  41 @item
  42 a parametric solution in terms of variable @code{%t}, or
  43
  44 @item
  45 a transformation into another ODE in variable @code{%u}.
  46
  47 @end itemize
  48
  49 @code{contrib_ode} uses the global variables @mrefcomma{%c}
  50 @mrefcomma{%k1} @mrefcomma{%k2} @mref{method} and @mref{yp}
  51 similarly to @code{ode2}.
  52
  53 If @code{contrib_ode}
  54 cannot obtain a solution for whatever reason, it returns @code{false}, after
  55 perhaps printing out an error message.
  56
  57 It is necessary to return a list of solutions, as even first order
  58 non-linear ODEs can have multiple solutions.  For example:
  59
  60 @c ===beg===
  61 @c load("contrib_ode")$
  62 @c eqn:x*'diff(y,x)^2-(1+x*y)*'diff(y,x)+y=0;
  63 @c contrib_ode(eqn,y,x);
  64 @c method;
  65 @c ===end===
  66 @example
  67 (%i1) load("contrib_ode")$
  68 @group
  69 (%i2) eqn:x*'diff(y,x)^2-(1+x*y)*'diff(y,x)+y=0;
  70                     dy 2             dy
  71 (%o2)            x (--)  - (1 + x y) -- + y = 0
  72                     dx               dx
  73 @end group
  74 @group
  75 (%i3) contrib_ode(eqn,y,x);
  76                     dy 2             dy
  77 (%t3)            x (--)  - (1 + x y) -- + y = 0
  78                     dx               dx
  79
  80               first order equation not linear in y'
  81
  82                                              x
  83 (%o3)             [y = log(x) + %c, y = %c %e ]
  84 @end group
  85 @group
  86 (%i4) method;
  87 (%o4)                        factor
  88 @end group
  89 @end example
  90
  91 Nonlinear ODEs can have singular solutions without constants of
  92 integration, as in the second solution of the following example:
  93
  94 @c ===beg===
  95 @c load("contrib_ode")$
  96 @c eqn:'diff(y,x)^2+x*'diff(y,x)-y=0;
  97 @c contrib_ode(eqn,y,x);
  98 @c method;
  99 @c ===end===
 100 @example
 101 (%i1) load("contrib_ode")$
 102 @group
 103 (%i2) eqn:'diff(y,x)^2+x*'diff(y,x)-y=0;
 104                        dy 2     dy
 105 (%o2)                 (--)  + x -- - y = 0
 106                        dx       dx
 107 @end group
 108 (%i3) contrib_ode(eqn,y,x);
 109                        dy 2     dy
 110 (%t3)                 (--)  + x -- - y = 0
 111                        dx       dx
 112
 113               first order equation not linear in y'
 114
 115                                            2
 116                                  2        x
 117 (%o3)              [y = %c x + %c , y = - --]
 118                                           4
 119 @group
 120 (%i4) method;
 121 (%o4)                       clairaut
 122 @end group
 123 @end example
 124
 125
 126 The following ODE has two parametric solutions in terms of the dummy
 127 variable @code{%t}.  In this case the parametric solutions can be manipulated
 128 to give explicit solutions.
 129
 130 @c ===beg===
 131 @c load("contrib_ode")$
 132 @c eqn:'diff(y,x)=(x+y)^2;
 133 @c contrib_ode(eqn,y,x);
 134 @c method;
 135 @c ===end===
 136 @example
 137 (%i1) load("contrib_ode")$
 138 @group
 139 (%i2) eqn:'diff(y,x)=(x+y)^2;
 140                           dy          2
 141 (%o2)                     -- = (x + y)
 142                           dx
 143 @end group
 144 @group
 145 (%i3) contrib_ode(eqn,y,x);
 146 (%o3) [[x = %c - atan(sqrt(%t)), y = (- x) - sqrt(%t)],
 147                      [x = atan(sqrt(%t)) + %c, y = sqrt(%t) - x]]
 148 @end group
 149 @group
 150 (%i4) method;
 151 (%o4)                       lagrange
 152 @end group
 153 @end example
 154
 155 The following example (Kamke 1.112) demonstrates an implicit solution.
 156
 157 @c ===beg===
 158 @c load("contrib_ode")$
 159 @c assume(x>0,y>0);
 160 @c eqn:x*'diff(y,x)-x*sqrt(y^2+x^2)-y;
 161 @c contrib_ode(eqn,y,x);
 162 @c method;
 163 @c ===end===
 164 @example
 165 (%i1) load("contrib_ode")$
 166 @group
 167 (%i2) assume(x>0,y>0);
 168 (%o2)                    [x > 0, y > 0]
 169 @end group
 170 @group
 171 (%i3) eqn:x*'diff(y,x)-x*sqrt(y^2+x^2)-y;
 172                      dy           2    2
 173 (%o3)              x -- - x sqrt(y  + x ) - y
 174                      dx
 175 @end group
 176 @group
 177 (%i4) contrib_ode(eqn,y,x);
 178                                   y
 179 (%o4)                  [x - asinh(-) = %c]
 180                                   x
 181 @end group
 182 @group
 183 (%i5) method;
 184 (%o5)                          lie
 185 @end group
 186 @end example
 187
 188
 189
 190 The following Riccati equation is transformed into a linear
 191 second order ODE in the variable @code{%u}.  Maxima is unable to
 192 solve the new ODE, so it is returned unevaluated.
 193 @c ===beg===
 194 @c load("contrib_ode")$
 195 @c eqn:x^2*'diff(y,x)=a+b*x^n+c*x^2*y^2;
 196 @c contrib_ode(eqn,y,x);
 197 @c method;
 198 @c ===end===
 199 @example
 200 (%i1) load("contrib_ode")$
 201 @group
 202 (%i2) eqn:x^2*'diff(y,x)=a+b*x^n+c*x^2*y^2;
 203                     2 dy      2  2      n
 204 (%o2)              x  -- = c x  y  + b x  + a
 205                       dx
 206 @end group
 207 @group
 208 (%i3) contrib_ode(eqn,y,x);
 209                d%u
 210                ---                            2
 211                dx        2  a       n - 2    d %u
 212 (%o3)  [[y = - ----, %u c  (-- + b x     ) + ---- c = 0]]
 213                %u c          2                 2
 214                             x                dx
 215 @end group
 216 @group
 217 (%i4) method;
 218 (%o4)                        riccati
 219 @end group
 220 @end example
 221
 222
 223 For first order ODEs @code{contrib_ode} calls @code{ode2}.  It then tries the
 224 following methods: factorization, Clairaut, Lagrange, Riccati,
 225 Abel and Lie symmetry methods.  The Lie method is not attempted
 226 on Abel equations if the Abel method fails, but it is tried
 227 if the Riccati method returns an unsolved second order ODE.
 228
 229 For second order ODEs @code{contrib_ode} calls @code{ode2} then @mref{odelin}.
 230
 231 Extensive debugging traces and messages are displayed if the command
 232 @code{put('contrib_ode,true,'verbose)} is executed.
 233
 234 @opencatbox{Categories:}
 235 @category{Differential equations}
 236 @category{Share packages}
 237 @category{Package contrib_ode}
 238 @closecatbox
 239
 240
 241 @node Functions and Variables for contrib_ode, Possible improvements to contrib_ode, Introduction to contrib_ode, Package contrib_ode
 242 @section Functions and Variables for contrib_ode
 243
 244 @anchor{contrib_ode}
 245 @deffn {Function} contrib_ode (@var{eqn}, @var{y}, @var{x})
 246
 247 Returns a list of solutions of the ODE @var{eqn} with
 248 independent variable @var{x} and dependent variable @var{y}.
 249
 250
 251 @opencatbox{Categories:}
 252 @category{Package contrib_ode}
 253 @closecatbox
 254
 255 @end deffn
 256
 257 @anchor{odelin}
 258 @deffn {Function} odelin (@var{eqn}, @var{y}, @var{x})
 259
 260 @code{odelin} solves linear homogeneous ODEs of first and
 261 second order with
 262 independent variable @var{x} and dependent variable @var{y}.
 263 It returns a fundamental solution set of the ODE.
 264
 265 For second order ODEs, @code{odelin} uses a method, due to Bronstein
 266 and Lafaille, that searches for solutions in terms of given
 267 special functions.
 268
 269 @c ===beg===
 270 @c load("contrib_ode")$
 271 @c odelin(x*(x+1)*'diff(y,x,2)+(x+5)*'diff(y,x,1)+(-4)*y,y,x);
 272 @c ===end===
 273 @example
 274 (%i1) load("contrib_ode")$
 275 @group
 276 (%i2) odelin(x*(x+1)*'diff(y,x,2)+(x+5)*'diff(y,x,1)+(-4)*y,y,x);
 277        gauss_a(- 6, - 2, - 3, - x)  gauss_b(- 6, - 2, - 3, - x)
 278 (%o2) @{---------------------------, ---------------------------@}
 279                     4                            4
 280                    x                            x
 281 @end group
 282 @end example
 283
 284 @opencatbox{Categories:}
 285 @category{Package contrib_ode}
 286 @closecatbox
 287
 288 @end deffn
 289
 290 @anchor{ode_check}
 291 @deffn {Function} ode_check (@var{eqn}, @var{soln})
 292
 293 Returns the value of ODE @var{eqn} after substituting a
 294 possible solution @var{soln}.  The value is equivalent to
 295 zero if @var{soln} is a solution of @var{eqn}.
 296
 297 @c ===beg===
 298 @c load("contrib_ode")$
 299 @c eqn:'diff(y,x,2)+(a*x+b)*y;
 300 @c ans:[y = bessel_y(1/3,2*(a*x+b)^(3/2)/(3*a))*%k2*sqrt(a*x+b)
 301 @c          +bessel_j(1/3,2*(a*x+b)^(3/2)/(3*a))*%k1*sqrt(a*x+b)];
 302 @c ode_check(eqn,ans[1]);
 303 @c ===end===
 304 @example
 305 (%i1) load("contrib_ode")$
 306 @group
 307 (%i2) eqn:'diff(y,x,2)+(a*x+b)*y;
 308                          2
 309                         d y
 310 (%o2)                   --- + (b + a x) y
 311                           2
 312                         dx
 313 @end group
 314 (%i3) ans:[y = bessel_y(1/3,2*(a*x+b)^(3/2)/(3*a))*%k2*sqrt(a*x+b)
 315          +bessel_j(1/3,2*(a*x+b)^(3/2)/(3*a))*%k1*sqrt(a*x+b)];
 316                                   3/2
 317                     1  2 (b + a x)
 318 (%o3) [y = bessel_y(-, --------------) %k2 sqrt(a x + b)
 319                     3       3 a
 320                                           3/2
 321                             1  2 (b + a x)
 322                  + bessel_j(-, --------------) %k1 sqrt(a x + b)]
 323                             3       3 a
 324 @group
 325 (%i4) ode_check(eqn,ans[1]);
 326 (%o4)                           0
 327 @end group
 328 @end example
 329
 330 @opencatbox{Categories:}
 331 @category{Package contrib_ode}
 332 @closecatbox
 333
 334 @end deffn
 335
 336
 337 @anchor{gauss_a}
 338 @deffn {Function} gauss_a (@var{a}, @var{b}, @var{c}, @var{x})
 339
 340 @code{gauss_a(a,b,c,x)} and @code{gauss_b(a,b,c,x)} are 2F1
 341 hypergeometric functions.  They represent any two independent
 342 solutions of the hypergeometric differential equation
 343 @code{x*(1-x) diff(y,x,2) + [c-(a+b+1)x] diff(y,x) - a*b*y = 0} (A&S 15.5.1).
 344
 345 The only use of these functions is in solutions of ODEs returned by
 346 @mref{odelin} and @mref{contrib_ode}.  The definition and use of these
 347 functions may change in future releases of Maxima.
 348
 349 See also @mrefcomma{gauss_b} @mref{dgauss_a} and @mrefdot{gauss_b}
 350
 351 @opencatbox{Categories:}
 352 @category{Package contrib_ode}
 353 @closecatbox
 354
 355 @end deffn
 356
 357 @anchor{gauss_b}
 358 @deffn {Function} gauss_b (@var{a}, @var{b}, @var{c}, @var{x})
 359 See @mrefdot{gauss_a}
 360
 361 @opencatbox{Categories:}
 362 @category{Package contrib_ode}
 363 @closecatbox
 364
 365 @end deffn
 366
 367 @anchor{dgauss_a}
 368 @deffn {Function} dgauss_a (@var{a}, @var{b}, @var{c}, @var{x})
 369 The derivative with respect to @var{x}
 370 of @mref{gauss_a}@code{(@var{a}, @var{b}, @var{c}, @var{x})}.
 371
 372 @opencatbox{Categories:}
 373 @category{Package contrib_ode}
 374 @closecatbox
 375
 376 @end deffn
 377
 378 @anchor{dgauss_b}
 379 @deffn {Function} dgauss_b (@var{a}, @var{b}, @var{c}, @var{x})
 380 The derivative with respect to @var{x}
 381 of @mref{gauss_b}@code{(@var{a}, @var{b}, @var{c}, @var{x})}.
 382
 383 @opencatbox{Categories:}
 384 @category{Package contrib_ode}
 385 @closecatbox
 386
 387 @end deffn
 388
 389
 390 @anchor{kummer_m}
 391 @deffn {Function} kummer_m (@var{a}, @var{b}, @var{x})
 392
 393 Kummer's M function, as defined in Abramowitz and Stegun,
 394 @i{Handbook of Mathematical Functions}, Section 13.1.2.
 395
 396 The only use of this function is in solutions of ODEs returned by
 397 @mref{odelin} and @mref{contrib_ode}.  The definition and use of this
 398 function may change in future releases of Maxima.
 399
 400 See also @mrefcomma{kummer_u} @mrefcomma{dkummer_m} and @mrefdot{dkummer_u}
 401
 402 @opencatbox{Categories:}
 403 @category{Package contrib_ode}
 404 @closecatbox
 405
 406 @end deffn
 407
 408 @anchor{kummer_u}
 409 @deffn {Function} kummer_u (@var{a}, @var{b}, @var{x})
 410
 411 Kummer's U function, as defined in Abramowitz and Stegun,
 412 @i{Handbook of Mathematical Functions}, Section 13.1.3.
 413
 414 See @mrefdot{kummer_m}
 415
 416 @opencatbox{Categories:}
 417 @category{Package contrib_ode}
 418 @closecatbox
 419
 420 @end deffn
 421
 422 @anchor{dkummer_m}
 423 @deffn {Function} dkummer_m (@var{a}, @var{b}, @var{x})
 424 The derivative with respect to @var{x}
 425 of @mref{kummer_m}@code{(@var{a}, @var{b}, @var{x})}.
 426
 427 @opencatbox{Categories:}
 428 @category{Package contrib_ode}
 429 @closecatbox
 430
 431 @end deffn
 432
 433 @anchor{dkummer_u}
 434 @deffn {Function} dkummer_u (@var{a}, @var{b}, @var{x})
 435 The derivative with respect to @var{x}
 436 of @mref{kummer_u}@code{(@var{a}, @var{b}, @var{x})}.
 437
 438 @opencatbox{Categories:}
 439 @category{Package contrib_ode}
 440 @closecatbox
 441
 442 @end deffn
 443
 444 @anchor{bessel_simplify}
 445 @deffn {Function} bessel_simplify (@var{expr})
 446 Simplifies expressions containing Bessel functions @mrefcomma{bessel_j}
 447 @mrefcomma{bessel_y} @mrefcomma{bessel_i} @mrefcomma{bessel_k}
 448 @mrefcomma{hankel_1} @mrefcomma{hankel_2} @mref{struve_h}
 449 and @mrefdot{struve_l}
 450 Recurrence relations (DLMF §10.6(i))(A&S 9.1.27)
 451 are used to replace functions of highest order n
 452 by functions of order n-1 and n-2.
 453
 454 This process is repeated until all the orders
 455 differ by less than 2.
 456
 457 @c ===beg===
 458 @c load("contrib_ode")$
 459 @c bessel_simplify(4*bessel_j(n,x^2)*(x^2-n^2/x^2)
 460 @c   +x*((bessel_j(n-2,x^2)-bessel_j(n,x^2))*x
 461 @c   -(bessel_j(n,x^2)-bessel_j(n+2,x^2))*x)
 462 @c   -2*bessel_j(n+1,x^2)+2*bessel_j(n-1,x^2));
 463 @c bessel_simplify( -2*bessel_j(1,z)*z^3 - 10*bessel_j(2,z)*z^2
 464 @c  + 15*%pi*bessel_j(1,z)*struve_h(3,z)*z - 15*%pi*struve_h(1,z)
 465 @c    *bessel_j(3,z)*z - 15*%pi*bessel_j(0,z)*struve_h(2,z)*z
 466 @c  + 15*%pi*struve_h(0,z)*bessel_j(2,z)*z - 30*%pi*bessel_j(1,z)
 467 @c    *struve_h(2,z) + 30*%pi*struve_h(1,z)*bessel_j(2,z));
 468 @c ===end===
 469 @example
 470 (%i1) load("contrib_ode")$
 471 @group
 472 (%i2) bessel_simplify(4*bessel_j(n,x^2)*(x^2-n^2/x^2)
 473   +x*((bessel_j(n-2,x^2)-bessel_j(n,x^2))*x
 474   -(bessel_j(n,x^2)-bessel_j(n+2,x^2))*x)
 475   -2*bessel_j(n+1,x^2)+2*bessel_j(n-1,x^2));
 476 (%o2)                           0
 477 @end group
 478 @group
 479 (%i3) bessel_simplify( -2*bessel_j(1,z)*z^3 - 10*bessel_j(2,z)*z^2
 480  + 15*%pi*bessel_j(1,z)*struve_h(3,z)*z - 15*%pi*struve_h(1,z)
 481    *bessel_j(3,z)*z - 15*%pi*bessel_j(0,z)*struve_h(2,z)*z
 482  + 15*%pi*struve_h(0,z)*bessel_j(2,z)*z - 30*%pi*bessel_j(1,z)
 483    *struve_h(2,z) + 30*%pi*struve_h(1,z)*bessel_j(2,z));
 484 (%o3)                           0
 485 @end group
 486 @end example
 487
 488 @opencatbox{Categories:}
 489 @category{Package contrib_ode}
 490 @category{Bessel functions}
 491 @category{Special functions}
 492 @closecatbox
 493
 494 @end deffn
 495
 496 @anchor{expintegral_e_simplify}
 497 @deffn {Function} expintegral_e_simplify (@var{expr})
 498 Simplify expressions containing exponential integral @mref{expintegral_e}
 499 using the recurrence (A&S 5.1.14).
 500
 501    expintegral_e(n+1,z)
 502         = (1/n) * (exp(-z)-z*expintegral_e(n,z))      n = 1,2,3 ....
 503
 504 @opencatbox{Categories:}
 505 @category{Package contrib_ode}
 506 @category{Exponential Integrals}
 507 @category{Special functions}
 508 @closecatbox
 509
 510 @end deffn
 511
 512
 513 @node Possible improvements to contrib_ode, Test cases for contrib_ode, Functions and Variables for contrib_ode, Package contrib_ode
 514 @section Possible improvements to contrib_ode
 515
 516
 517 These routines are work in progress.  I still need to:
 518
 519 @itemize @bullet
 520
 521 @item
 522 Extend the FACTOR method @code{ode1_factor} to work for multiple roots.
 523
 524 @item
 525 Extend the FACTOR method @code{ode1_factor} to attempt to solve higher
 526   order factors.  At present it only attempts to solve linear factors.
 527
 528 @item
 529 Fix the LAGRANGE routine @code{ode1_lagrange} to prefer real roots over
 530   complex roots.
 531
 532 @item
 533 Add additional methods for Riccati equations.
 534
 535 @item
 536 Improve the detection of Abel equations of second kind.  The existing
 537   pattern matching is weak.
 538
 539 @item
 540 Work on the Lie symmetry group routine @code{ode1_lie}.  There are quite a
 541   few problems with it: some parts are unimplemented; some test cases
 542   seem to run forever; other test cases crash; yet others return very
 543   complex "solutions".  I wonder if it really ready for release yet.
 544
 545 @item
 546 Add more test cases.
 547
 548 @end itemize
 549
 550 @node Test cases for contrib_ode, References for contrib_ode, Possible improvements to contrib_ode, Package contrib_ode
 551 @section Test cases for contrib_ode
 552
 553
 554 The routines have been tested on a approximately one thousand  test cases
 555 from Murphy,
 556 Kamke, Zwillinger and elsewhere.  These are included in the tests subdirectory.
 557
 558 @itemize @bullet
 559 @item
 560 The Clairaut routine @code{ode1_clairaut} finds all known solutions,
 561   including singular solutions, of the Clairaut equations in Murphy and
 562   Kamke.
 563
 564 @item
 565 The other routines often return a single solution when multiple
 566   solutions exist.
 567
 568 @item
 569 Some of the "solutions" from @code{ode1_lie} are overly complex and
 570   impossible to check.
 571
 572 @item
 573 There are some crashes.
 574
 575 @end itemize
 576
 577 @node References for contrib_ode, ,Test cases for contrib_ode, Package contrib_ode
 578 @section References for contrib_ode
 579
 580
 581 @enumerate
 582 @item
 583 E. Kamke, Differentialgleichungen L@"osungsmethoden und L@"osungen, Vol 1,
 584     Geest & Portig, Leipzig, 1961
 585
 586 @item
 587 G. M. Murphy, Ordinary Differential Equations and Their Solutions,
 588     Van Nostrand, New York, 1960
 589
 590 @item
 591 D. Zwillinger, Handbook of Differential Equations, 3rd edition,
 592     Academic Press, 1998
 593
 594 @item
 595 F. Schwarz, Symmetry Analysis of Abel's Equation, Studies in
 596     Applied Mathematics, 100:269-294 (1998)
 597
 598 @item
 599 F. Schwarz, Algorithmic Solution of Abel's Equation,
 600     Computing 61, 39-49 (1998)
 601
 602 @item
 603 E. S. Cheb-Terrab, A. D. Roche, Symmetries and First Order
 604     ODE Patterns, Computer Physics Communications 113 (1998), p 239.
 605 @c dead link
 606     (@url{http://lie.uwaterloo.ca/papers/ode_vii.pdf})
 607
 608 @item
 609 E. S. Cheb-Terrab, T. Kolokolnikov,  First Order ODEs,
 610     Symmetries and Linear Transformations, European Journal of
 611     Applied Mathematics, Vol. 14, No. 2, pp. 231-246 (2003).
 612 @c dead link
 613     (@url{http://arxiv.org/abs/math-ph/0007023},@*
 614     @url{http://lie.uwaterloo.ca/papers/ode_iv.pdf})
 615
 616
 617 @item
 618 G. W. Bluman, S. C. Anco, Symmetry and Integration Methods for
 619     Differential Equations, Springer, (2002)
 620
 621 @item
 622 M. Bronstein, S. Lafaille,
 623 Solutions of linear ordinary differential equations in terms
 624 of special functions,
 625 Proceedings of ISSAC 2002, Lille, ACM Press, 23-28.
 626 (@url{http://www-sop.inria.fr/cafe/Manuel.Bronstein/publications/issac2002.pdf})
 627
 628
 629 @end enumerate