share/contrib/symplectic_ode/symplectic_ode.texi

   1 \input texinfo
   2
   3 @setfilename symplectic_ode.info
   4 @settitle symplectic_ode
   5
   6 @ifinfo
   7 @macro var {expr}
   8 <\expr\>
   9 @end macro
  10 @end ifinfo
  11
  12 @dircategory Mathematics/Maxima
  13 @direntry
  14 * symplectic_ode: (maxima/symplectic_ode).           Maxima share package symplectic_ode for solving Hamilton's equations.
  15 @end direntry
  16
  17 @node Top, Introduction to symplectic_ode, (dir), (dir)
  18 @top
  19 @menu
  20 * Introduction to symplectic_ode::
  21 * Definitions for symplectic_ode::
  22 * Function and variable index::
  23 @end menu
  24 @chapter symplectic_ode
  25
  26
  27 @node Introduction to symplectic_ode, Definitions for symplectic_ode, Top, Top
  28 @section Introduction to symplectic_ode
  29
  30 For a hamiltonian of the form F(p) + G(q), the function @code{symplectic_ode} numerically solves
  31 Hamilton's equations of motion @math{p' = -dH/dq, q' = dH/dp}, where @math{H} is the
  32 hamiltonian. The method preserves the Poisson bracket of p and q. One time step is accomplished with the loop
  33 @*
  34 @sp 1
  35 @verbatim
  36 for k = 1 ... n
  37    q <-  q + c(k)*diff(H,p)*dt ;update momentum p
  38    p <-  p - d(k)*diff(H,q)*dt ;use updated position q
  39 @end verbatim
  40 @*
  41 @sp 1
  42 where @math{c1, c2, @dots{}, cn} and @math{d1, d2, @dots{}, dn} are constants and @math{H} is the hamiltonian.
  43
  44 This code has built-in methods for the symplectic Euler, the Verlet method, and third and fourth order methods
  45 that are due to Ronald Ruth. For a complete description of these methods, see
  46 @uref{https://en.wikipedia.org/wiki/Symplectic_integrator}. Additionally, there is a fifth order method that is
  47 due to Kostas Tselios and T. E. Simos.
  48
  49 The function @code{symplectic_ode} creates and compiles functions for updating @var{p} & @var{q}.
  50 The arguments to these functions are modedeclared to have type given by an optional argument
  51 (defaults to float).  Of course, hand coding of these functions could increase speed or accuracy,
  52 but the automatically generated functions are convenient.
  53
  54 Unlike adaptive methods such as RK45, @code{symplectic_ode} uses a fixed step size. Generally
  55 symplectic methods that use an adaptive step size lose their advantages over non-symplectic methods.
  56
  57 @node Definitions for symplectic_ode, Function and variable index, Introduction to symplectic_ode, Top
  58 @section Definitions for symplectic_ode
  59
  60 @deffn{Function} poisson_bracket poisson_bracket(@var{f}, @var{g}, @var{p}, @var{q})
  61 poisson_bracket(@var{f}, @var{g}, [@var{p1},@dots{}, @var{pn}], [@var{q1},@dots{}, @var{qn}])
  62
  63 Compute the Poisson bracket of the expressions @var{f} and @var{g} with respect to the canonical
  64 coordinates @var{p} and @var{q} (or @var{p1}, @var{p2},@dots{}, @var{pn} and @var{p1},
  65 @var{p2},@dots{}, @var{pn}).
  66
  67 @b{Examples:}
  68 @example
  69 (%i1) load("symplectic_ode")$
  70 (%i2) poisson_bracket(p,p^2/2+q^4,p,q);
  71 (%o2) -4*q^3
  72
  73 (%i3) poisson_bracket(q,p^2/2+q^4,p,q);
  74 (%o3) p
  75
  76 (%i4) poisson_bracket(q1,p1^2/2+p2^2/2+q1^4+q2^4,[p1,p2],[q1,q2]);
  77 (%o4) p1
  78
  79 (%i5) poisson_bracket(p1,p1^2/2+p2^2/2+q1^4+q2^4,[p1,p2],[q1,q2]);
  80 (%o5) -4*q1^3
  81
  82 @end example
  83 @end deffn
  84
  85
  86 @deffn{Function} symplectic_ode symplectic_ode(ham,p,q,po,qo,dt,N)
  87      symplectic_ode(ham,p,q,po,qo,dt,N,method)
  88      symplectic_ode(ham,p,q,po,qo,dt,N,method,type)
  89
  90 Numerically solve Hamilton's equations of motion using a symplectic method. Specifically:
  91 @*
  92 @sp 1
  93 @itemize @bullet
  94
  95 @item The hamiltonian is the Maxima expression @var{ham} that depends on the canonical coordinates
  96 @var{p} and @var{q}. The hamiltonian must be time independent. The method is symplectic when the
  97 hamiltonian is separable; that is when it has the form @code{f(p) + g(q)}.
  98
  99 @item The canonical coordinates are @var{p} and @var{q}. The arguments @var{p} and @var{q} should be
 100 symbols or equal length lists of symbols.
 101
 102 @item The arguments @var{po} and @var{q0} are the initial values of @var{p} and @var{q}, respectively.
 103 These should be expressions or equal length lists of expressions. Generally, the values of @var{po} and @var{q0}
 104 should be numbers. When the optional argument @var{type} is float, the code attempts to convert the values
 105 of @var{po} and @var{qo} into floating point numbers; when this isn't possible, the code signals an error.
 106
 107 @item @var{dt} is the fixed time step.
 108
 109 @item @var{N} is the number of time steps.
 110
 111 @item The optional argument @var{method} determines the integration method. It must be either
 112 symplectic_euler (default), verlet, symplectic_third_order, symplectic_fourth_order, or
 113 symplectic_fifth_order. For an explanation of these methods, see https://en.wikipedia.org/wiki/Symplectic_integrator.
 114
 115  @item The optional argument @var{type} determines the value for mode_declare for various
 116        automatically generated functions. The value @var{type} must be one of float (default),
 117       rational, or any (no type). Since @var{float} is a Maxima option variable, the @var{type}
 118       variable should be quoted, especially for type @var{float}.
 119
 120   @end itemize
 121   @*
 122 @sp 1
 123
 124 For both the scalar case (both @var{p} and @var{q} are mapatoms) and the nonscalar case
 125 (both @var{p} and @var{q} are lists of mapatoms), @code{symplectic_ode} returns a list
 126 of two lists. For the scalar case, the first list is a list of the values of @var{p} at
 127 the times @code{0, dt, 2*dt,@dots{}, N*dt} and similarly for the second list.
 128 For a nonscalar case, the first list is a list of the form @math{[p1, p2,@dots{}, pn]} at
 129 the times @code{0, dt, 2*dt,@dots{}, N*dt}.
 130
 131 @b{Examples:}
 132 @example
 133
 134 (%i2) load("symplectic_ode")$
 135 (%i3) symplectic_ode(p^2/2 + q^4/4,p,q,1,0,1/10,2);
 136 (%o3) [[1.0,1.0,0.9999],[0.0,0.1,0.19999]]
 137
 138 (%i4) symplectic_ode(p^2/2 + q^4/4,[p],[q],[1],[0],1/10,2);
 139 (%o4) [[[1.0],[1.0],[0.9999]],[[0.0],[0.1],[0.19999]]]
 140
 141 (%i5) symplectic_ode(p^2/2 + q^4/4,p,q,1,0,1/10,2,verlet);
 142 (%o5) [[1.0,0.9999875,0.9996500084374297],[0.0,0.099999375,0.1999812504218715]]
 143
 144 (%i6) symplectic_ode(p^2/2 + q^4/4,p,q,1.0b0,0.0b0, 0.1b0,2,verlet,'any);
 145 (%o6) [[1.0b0,9.999875b-1,9.996500084374297b-1],[0.0b0,9.9999375b-2,1.999812504218715b-1]]
 146
 147 @end example
 148
 149
 150 @end deffn
 151
 152 @node Function and variable index,  , Definitions for symplectic_ode, Top
 153 @appendix Function and variable index
 154 @printindex fn
 155 @printindex vr
 156
 157 @bye