share/contrib/symplectic_ode/rtest_symplectic_ode.mac

   1 (load("symplectic_ode"),ty : elapsed_real_time(), 0);
   2 0$
   3
   4 (reset(), display2d : false, 0);
   5 0$
   6
   7 ([pp, qq] : symplectic_ode (p^2/2 + 42*q^4/4,[p], [q], [po],[qo], dt, 2, 'symplectic_euler, 'any),0);
   8 0$
   9
  10 expand( poisson_bracket(first(last(pp)), first(last(qq)),[po],[qo]));
  11 -1$
  12
  13 ([pp, qq] : symplectic_ode (p^2/2 + 42*q^4/4,p,q,po,qo, dt, 2, 'symplectic_euler, 'any),0);
  14 0$
  15
  16 expand( poisson_bracket(last(pp), last(qq),po,qo));
  17 -1$
  18
  19 ([pp, qq] : symplectic_ode(p1^2/2 + p2^2/2 + 42*q1^2 - q1*q2 + q2^4,[p1,p2],
  20      [q1,q2], [p10, p20],[q10, q20], dt, 2, 'symplectic_euler, 'any),0);
  21 0$
  22
  23 (pp : last(pp), qq : last(qq),0);
  24 0$
  25
  26 expand(outermap(lambda([f,g],  poisson_bracket(f,g, [p10,p20],[q10,q20])), pp,qq));
  27 [[-1,0],[0,-1]]$
  28
  29 ratsimp(expand(outermap(lambda([f,g],  poisson_bracket(f,g, [p10,p20],[q10,q20])), pp,pp)));
  30 [[0,0],[0,0]]$
  31
  32 ratsimp(expand(outermap(lambda([f,g],  poisson_bracket(f,g, [p10,p20],[q10,q20])), qq,qq)));
  33 [[0,0],[0,0]]$
  34
  35 /* check for correct order */
  36
  37 ([pp,qq] : symplectic_ode (p^2/2 + q^2/2,[p], [q],  [po],[qo], dt,  4, 'symplectic_euler, 'any),0);
  38 0$
  39
  40 taylor(first(last(pp)) - (po * cos(4 *dt) - qo * sin(4 *dt)),dt,0,2);
  41 2*po*dt^2$
  42
  43 taylor(first(last(qq)) - (qo * cos(4 *dt)  + po  * sin(4 *dt)),dt,0,2);
  44 -2*qo*dt^2$
  45
  46 /* Test verlet method */
  47
  48 ([pp, qq] : symplectic_ode(p^2/2 + q^4/4,[p],[q],[po],[qo],dt,1, 'verlet, 'any),0);
  49 0$
  50
  51 expand( poisson_bracket(first(last(pp)), first(last(qq)),[po],[qo]));
  52 -1$
  53
  54 ([pp, qq] : symplectic_ode(p1^2/2 + p2^2/2 + 42*q1^2 - q1*q2 + q2^4,[p1,p2],
  55      [q1,q2], [p10, p20],[q10, q20], dt, 2, 'verlet, 'any),0);
  56 0$
  57
  58 (pp : last(pp), qq : last(qq),0);
  59 0$
  60
  61 expand(outermap(lambda([f,g],  poisson_bracket(f,g, [p10,p20],[q10,q20])), pp,qq));
  62 [[-1,0],[0,-1]]$
  63
  64 expand(outermap(lambda([f,g],  poisson_bracket(f,g, [p10,p20],[q10,q20])), pp,pp));
  65 [[0,0],[0,0]]$
  66
  67 expand(outermap(lambda([f,g],  poisson_bracket(f,g, [p10,p20],[q10,q20])), qq,qq));
  68 [[0,0],[0,0]]$
  69
  70 ([pp,qq] : symplectic_ode (p^2/2 + q^2/2,[p], [q],  [po],[qo], dt, 4, 'verlet, 'any),0);
  71 0$
  72
  73 taylor(first(last(pp)) - (po * cos(4 *dt) - qo * sin(4 *dt)),dt,0,4);
  74 -(2*qo*dt^3)/3-(2*po*dt^4)/3$
  75
  76 taylor(first(last(qq)) - (qo * cos(4 *dt)  + po  * sin(4 *dt)),dt,0,4);
  77 -(po*dt^3)/3-(2*qo*dt^4)/3$
  78
  79 /* Test a third order method */
  80
  81 ([pp, qq] : symplectic_ode (p^2/2 + q^4/4,[p],[q],[po],[qo],dt,1, 'symplectic_third_order, 'any),0);
  82 0$
  83
  84 expand( poisson_bracket(first(last(pp)), first(last(qq)),[po],[qo]));
  85 -1$
  86
  87 ([pp, qq] : symplectic_ode (p1^2/2 + p2^2/2 + 42*q1^2 - q1*q2 + q2^4,[p1,p2],
  88      [q1,q2], [p10, p20],[q10, q20], dt, 2, 'verlet, 'any),0);
  89 0$
  90
  91 (pp : last(pp), qq : last(qq),0);
  92 0$
  93
  94 expand(outermap(lambda([f,g],  poisson_bracket(f,g,[p10,p20],[q10,q20])), pp,qq));
  95 [[-1,0],[0,-1]]$
  96
  97 expand(outermap(lambda([f,g],  poisson_bracket(f,g,[p10,p20],[q10,q20])), pp,pp));
  98 [[0,0],[0,0]]$
  99
 100 expand(outermap(lambda([f,g],  poisson_bracket(f,g,[p10,p20],[q10,q20])), qq,qq));
 101 [[0,0],[0,0]]$
 102
 103 ([pp,qq] : symplectic_ode (p^2/2 + q^2/2,[p], [q],[po],[qo], dt, 4, 'symplectic_third_order, 'any),0);
 104 0$
 105
 106 taylor(first(last(pp)) - (po * cos(4 *dt) - qo * sin(4 *dt)),dt,0,5);
 107 -(po*dt^4)/9-(qo*dt^5)/45$
 108
 109 taylor(first(last(qq)) - (qo * cos(4 *dt)  + po  * sin(4 *dt)),dt,0,5);
 110 (qo*dt^4)/9-(37*po*dt^5)/2160$
 111
 112 /* Test a fourth order method */
 113
 114 (algebraic : true,0);
 115 0$
 116
 117 ([pp, qq] : symplectic_ode (p^2/2,[p],[q],[po],[qo],dt,1, 'symplectic_fourth_order, 'any),0);
 118 0$
 119
 120 expand( poisson_bracket(first(last(pp)), first(last(qq)),[po],[qo]));
 121 -1$
 122
 123 ([pp, qq] : symplectic_ode (p1^2/2 + p2^2/2 + q1*q2,[p1,p2],
 124      [q1,q2], [p10, p20],[q10, q20], dt, 2, 'symplectic_fourth_order, 'any),0);
 125 0$
 126
 127 (pp : last(pp), qq : last(qq),0);
 128 0$
 129
 130 ratsimp(expand(outermap(lambda([f,g],  poisson_bracket(f,g,[p10,p20],[q10,q20])), pp,qq)));
 131 [[-1,0],[0,-1]]$
 132
 133 ratsimp(expand(outermap(lambda([f,g],  poisson_bracket(f,g,[p10,p20],[q10,q20])), pp,pp)));
 134 [[0,0],[0,0]]$
 135
 136 ratsimp(expand(outermap(lambda([f,g],  poisson_bracket(f,g,[p10,p20],[q10,q20])), qq,qq)));
 137 [[0,0],[0,0]]$
 138
 139 ([pp,qq] : symplectic_ode(p^2/2 + q^2/2,[p], [q],[po],[qo], dt, 4, 'symplectic_fourth_order, 'any),0);
 140 0$
 141
 142 taylor(first(last(pp)) - (po * cos(4 *dt) - qo * sin(4 *dt)),dt,0,5);
 143 ((3*(2^(1/3))^8-5*(2^(1/3))^7-5*(2^(1/3))^4+18)*qo*dt^5)/(225*(2^(1/3))^10-585*(2^(1/3))^8+1440)$
 144
 145 taylor(first(last(qq)) - (qo * cos(4 *dt)  + po  * sin(4 *dt)),dt,0,5);
 146 -((19*(2^(1/3))^8-155*(2^(1/3))^4-85*(2^(1/3))^2+325*2^(1/3)-6)*po*dt^5)/(225*(2^(1/3))^7-585*(2^(1/3))^5+720)$
 147
 148 /* type tests */
 149
 150 ([pp, qq] : symplectic_ode (p^2/2 + 42*q^4/4,[p], [q], [1],[0], 1/10, 57, 'symplectic_euler, 'float),0);
 151 0$
 152
 153 every('floatnump, xreduce('append, pp)) and every('floatnump, xreduce('append, qq));
 154 true$
 155
 156 ([pp, qq] : symplectic_ode (p^2/2 + 42*q^4/4, p, q, 1, 0, 1/10, 5, 'symplectic_euler, 'rational),0)$
 157 0$
 158
 159 every('ratnump, pp) and every('ratnump, qq);
 160 true$
 161
 162 ([pp, qq] : symplectic_ode (p^2/2 + 42*q^4/4, p, q, 1, 0, 1/10, 5, 'symplectic_euler, 'bfloat),0);
 163 0$
 164
 165 every('bfloatp, pp) and every('bfloatp, qq);
 166 true$
 167
 168 /*  poisson Brackets */
 169
 170  poisson_bracket(p,p,[p],[q]);
 171 0$
 172
 173  poisson_bracket(p,q,[p],[q]);
 174 -1$
 175
 176  poisson_bracket(p^2,q^2,[p],[q]);
 177 -4*p*q$
 178
 179  poisson_bracket(p^2,q^2,p,q);
 180 -4*p*q$
 181
 182 (ham : lambda([p,q], p^2/2 + q^4-q^2),dt : 1/10^4, N : 10^4, 0);
 183 0$
 184
 185 (x : symplectic_ode(ham(p,q), p,q,1,0,dt,N,'symplectic_euler, 'float),
 186  z : symplectic_ode(ham(p,q), p,q,1,0,dt, N, 'verlet, 'float),
 187  block([listarith : true], is(lmax(map('abs, xreduce('append,x-z))) < 2.0e-4)));
 188 true$
 189
 190 (z : symplectic_ode(ham(p,q),p,q,1,0,dt,N, 'symplectic_third_order, 'float),
 191  block([listarith : true], is(lmax(map('abs, xreduce('append,x-z))) < 2.0e-4)));
 192 true$
 193
 194 (z : symplectic_ode(ham(p,q), p,q,1,0,dt,N, 'symplectic_fourth_order, 'float),
 195  block([listarith : true], is(lmax(map('abs, xreduce('append,x-z))) < 2.0e-4)));
 196 true$
 197
 198 (z : symplectic_ode(ham(p,q), p,q,1,0,dt,N,'symplectic_fifth_order, 'float),
 199  block([listarith : true], is(lmax(map('abs, xreduce('append,x-z))) < 2.0e-4)));
 200 true$
 201
 202 /* How well does each method conserve energy? */
 203 (ham : lambda([p,q], p^2/2 + q^4-q^2),dt : 1/10, N : 10^4, 0);
 204 0$
 205
 206 (x : symplectic_ode(ham(p,q),p,q,1,0,dt,N,'symplectic_euler, 'float),
 207  x : map(ham,first(x),second(x)),
 208  is(lmax(x)-lmin(x) < 3.0e-1));
 209  true$
 210
 211  (x : symplectic_ode(ham(p,q), p,q,1,0,dt,N,'verlet, 'float),
 212   x : map(ham,first(x),second(x)),
 213   is(lmax(x)-lmin(x) < 2.0e-2));
 214  true$
 215
 216  (x : symplectic_ode(ham(p,q), p,q,1,0,dt,N,'symplectic_third_order, 'float),
 217   x : map(ham,first(x),second(x)),
 218   is(lmax(x)-lmin(x) < 2.0e-3));
 219  true$
 220
 221  (x : symplectic_ode(ham(p,q), p,q,1,0,dt,N,'symplectic_fourth_order, 'float),
 222   x : map(ham,first(x),second(x)),
 223   is(lmax(x)-lmin(x) < 5.0e-4));
 224  true$
 225
 226  (x : symplectic_ode(ham(p,q), p,q,1,0,dt,N,'symplectic_fifth_order, 'float),
 227   x : map(ham,first(x),second(x)),
 228   is(lmax(x)-lmin(x) < 2.0e-6));
 229  true$
 230
 231 /* bogus argument errors */
 232
 233 errcatch(symplectic_ode(p^2/2+(42*q^4)/4,[p],q,1,0, 1/10,2));
 234 []$
 235
 236 errcatch(symplectic_ode(p^2/2+(42*q^4)/4,p,[q],1,0, 1/10,2));
 237 []$
 238
 239 errcatch(symplectic_ode(p^2/2+(42*q^4)/4,p,q,[1],0, 1/10,2));
 240 []$
 241
 242 errcatch(symplectic_ode(p^2/2+(42*q^4)/4,p,q,1,[0], 1/10,2));
 243 []$
 244
 245 errcatch(symplectic_ode(p^2/2+(42*q^4)/4,p,q,1,0, [1/10],2));
 246 []$
 247
 248 errcatch(symplectic_ode(p^2/2+(42*q^4)/4,p,q,1,0,1/10,8.8));
 249 []$
 250
 251 errcatch( poisson_bracket(p^2,q^2,p,[q]));
 252 []$
 253
 254 errcatch( poisson_bracket(p^2,q^2,[p],q));
 255 []$
 256
 257 /* see mailing list "symplectic_ode and mapatom" 30 Nov 2020 */
 258 symplectic_ode (p^2/2 + q^2/2,p,q, %pi/7,%pi, 1/10, 2, 'symplectic_euler, 'any);
 259 [[%pi/7,(3*%pi)/70,-(403*%pi)/7000],[%pi,(703*%pi)/700,(69897*%pi)/70000]]$
 260
 261 errcatch(symplectic_ode (p^2/2 + q^2/2,[p],q, %pi/7,%pi, 1/10, 2, 'symplectic_euler, 'any));
 262 []$
 263
 264
 265 symplectic_ode (p^2/2 + q^2/2,[p],[q], [1],[1], 1/10, 2, 'symplectic_euler, 'any);
 266 [[[1],[9/10],[791/1000]],[[1],[109/100],[11691/10000]]]$
 267
 268 (print("Runtime was ", elapsed_real_time()-ty," seconds."),0);
 269 0$
 270
 271 /* make a mess? Clean it up.*/
 272 (remvalue(pp,qq,ty,ham,dt,N,x,z), reset(algebraic, display2d),0);
 273 0$
 274
 275
 276