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