File: symplectic_ode.texi

package info (click to toggle)
maxima-sage 5.45.1-2
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, sid, trixie
  • size: 113,788 kB
  • sloc: lisp: 440,833; fortran: 14,665; perl: 14,369; tcl: 10,997; sh: 4,475; makefile: 2,520; ansic: 447; python: 262; xml: 59; awk: 37; sed: 17
file content (157 lines) | stat: -rw-r--r-- 6,128 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
\input texinfo

@setfilename symplectic_ode.info
@settitle symplectic_ode

@ifinfo 
@macro var {expr}
<\expr\>
@end macro
@end ifinfo

@dircategory Mathematics/Maxima
@direntry
* symplectic_ode: (maxima/symplectic_ode).           Maxima share package symplectic_ode for solving Hamilton's equations.
@end direntry

@node Top, Introduction to symplectic_ode, (dir), (dir)
@top
@menu
* Introduction to symplectic_ode::
* Definitions for symplectic_ode::
* Function and variable index::
@end menu
@chapter symplectic_ode


@node Introduction to symplectic_ode, Definitions for symplectic_ode, Top, Top
@section Introduction to symplectic_ode

For a hamiltonian of the form F(p) + G(q), the function @code{symplectic_ode} numerically solves
Hamilton's equations of motion @math{p' = -dH/dq, q' = dH/dp}, where @math{H} is the 
hamiltonian. The method preserves the Poisson bracket of p and q. One time step is accomplished with the loop
@*
@sp 1
@verbatim
for k = 1 ... n 
   q <-  q + c(k)*diff(H,p)*dt ;update momentum p
   p <-  p - d(k)*diff(H,q)*dt ;use updated position q 
@end verbatim
@*
@sp 1
where @math{c1, c2, @dots{}, cn} and @math{d1, d2, @dots{}, dn} are constants and @math{H} is the hamiltonian. 

This code has built-in methods for the symplectic Euler, the Verlet method, and third and fourth order methods 
that are due to Ronald Ruth. For a complete description of these methods, see 
@uref{https://en.wikipedia.org/wiki/Symplectic_integrator}. Additionally, there is a fifth order method that is
due to Kostas Tselios and T. E. Simos.

The function @code{symplectic_ode} creates and compiles functions for updating @var{p} & @var{q}. 
The arguments to these functions are modedeclared to have type given by an optional argument 
(defaults to float).  Of course, hand coding of these functions could increase speed or accuracy, 
but the automatically generated functions are convenient.

Unlike adaptive methods such as RK45, @code{symplectic_ode} uses a fixed step size. Generally 
symplectic methods that use an adaptive step size lose their advantages over non-symplectic methods.

@node Definitions for symplectic_ode, Function and variable index, Introduction to symplectic_ode, Top
@section Definitions for symplectic_ode

@deffn{Function} poisson_bracket poisson_bracket(@var{f}, @var{g}, @var{p}, @var{q}) 
poisson_bracket(@var{f}, @var{g}, [@var{p1},@dots{}, @var{pn}], [@var{q1},@dots{}, @var{qn}]) 

Compute the Poisson bracket of the expressions @var{f} and @var{g} with respect to the canonical 
coordinates @var{p} and @var{q} (or @var{p1}, @var{p2},@dots{}, @var{pn} and @var{p1},  
@var{p2},@dots{}, @var{pn}).

@b{Examples:}
@example
(%i1) load("symplectic_ode")$
(%i2) poisson_bracket(p,p^2/2+q^4,p,q);
(%o2) -4*q^3

(%i3) poisson_bracket(q,p^2/2+q^4,p,q);
(%o3) p

(%i4) poisson_bracket(q1,p1^2/2+p2^2/2+q1^4+q2^4,[p1,p2],[q1,q2]);
(%o4) p1

(%i5) poisson_bracket(p1,p1^2/2+p2^2/2+q1^4+q2^4,[p1,p2],[q1,q2]);
(%o5) -4*q1^3

@end example
@end deffn


@deffn{Function} symplectic_ode symplectic_ode(ham,p,q,po,qo,dt,N)
     symplectic_ode(ham,p,q,po,qo,dt,N,method) 
     symplectic_ode(ham,p,q,po,qo,dt,N,method,type) 

Numerically solve Hamilton's equations of motion using a symplectic method. Specifically:
@*
@sp 1
@itemize @bullet

@item The hamiltonian is the Maxima expression @var{ham} that depends on the canonical coordinates 
@var{p} and @var{q}. The hamiltonian must be time independent. The method is symplectic when the 
hamiltonian is separable; that is when it has the form @code{f(p) + g(q)}.
   
@item The canonical coordinates are @var{p} and @var{q}. The arguments @var{p} and @var{q} should be 
symbols or equal length lists of symbols.

@item The arguments @var{po} and @var{q0} are the initial values of @var{p} and @var{q}, respectively. 
These should be expressions or equal length lists of expressions. Generally, the values of @var{po} and @var{q0}
should be numbers. When the optional argument @var{type} is float, the code attempts to convert the values
of @var{po} and @var{qo} into floating point numbers; when this isn't possible, the code signals an error.
  
@item @var{dt} is the fixed time step.

@item @var{N} is the number of time steps.

@item The optional argument @var{method} determines the integration method. It must be either 
symplectic_euler (default), verlet, symplectic_third_order, symplectic_fourth_order, or 
symplectic_fifth_order. For an explanation of these methods, see https://en.wikipedia.org/wiki/Symplectic_integrator.
     
 @item The optional argument @var{type} determines the value for mode_declare for various 
       automatically generated functions. The value @var{type} must be one of float (default), 
      rational, or any (no type). Since @var{float} is a Maxima option variable, the @var{type} 
      variable should be quoted, especially for type @var{float}.
  
  @end itemize
  @*
@sp 1
 
For both the scalar case (both @var{p} and @var{q} are mapatoms) and the nonscalar case 
(both @var{p} and @var{q} are lists of mapatoms), @code{symplectic_ode} returns a list 
of two lists. For the scalar case, the first list is a list of the values of @var{p} at 
the times @code{0, dt, 2*dt,@dots{}, N*dt} and similarly for the second list.
For a nonscalar case, the first list is a list of the form @math{[p1, p2,@dots{}, pn]} at 
the times @code{0, dt, 2*dt,@dots{}, N*dt}.

@b{Examples:}
@example

(%i2) load("symplectic_ode")$
(%i3) symplectic_ode(p^2/2 + q^4/4,p,q,1,0,1/10,2);
(%o3) [[1.0,1.0,0.9999],[0.0,0.1,0.19999]]

(%i4) symplectic_ode(p^2/2 + q^4/4,[p],[q],[1],[0],1/10,2);
(%o4) [[[1.0],[1.0],[0.9999]],[[0.0],[0.1],[0.19999]]]

(%i5) symplectic_ode(p^2/2 + q^4/4,p,q,1,0,1/10,2,verlet);
(%o5) [[1.0,0.9999875,0.9996500084374297],[0.0,0.099999375,0.1999812504218715]]

(%i6) symplectic_ode(p^2/2 + q^4/4,p,q,1.0b0,0.0b0, 0.1b0,2,verlet,'any);
(%o6) [[1.0b0,9.999875b-1,9.996500084374297b-1],[0.0b0,9.9999375b-2,1.999812504218715b-1]]

@end example


@end deffn

@node Function and variable index,  , Definitions for symplectic_ode, Top
@appendix Function and variable index
@printindex fn
@printindex vr

@bye