tests/rtest_elliptic.mac

   1 /* Tests of Jacobi elliptic functions and elliptic integrals */
   2
   3 kill(all);
   4 done$
   5
   6 /* derivatives */
   7 diff(jacobi_sn(u,m),u);
   8 jacobi_cn(u,m)*jacobi_dn(u,m);
   9
  10 diff(jacobi_sn(u,m),m);
  11 jacobi_cn(u,m)*jacobi_dn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))
  12   /(2*m)
  13  +jacobi_cn(u,m)^2*jacobi_sn(u,m)/(2*(1-m));
  14
  15 diff(jacobi_cn(u,m),u);
  16 -jacobi_dn(u,m)*jacobi_sn(u,m);
  17
  18 diff(jacobi_cn(u,m),m);
  19 -(jacobi_dn(u,m)*jacobi_sn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))/(2*m))
  20   -(jacobi_cn(u,m)*jacobi_sn(u,m)^2/(2*(1-m)));
  21
  22 diff(jacobi_dn(u,m),u);
  23 -m*jacobi_cn(u,m)*jacobi_sn(u,m);
  24
  25 diff(jacobi_dn(u,m),m);
  26 -(jacobi_cn(u,m)*jacobi_sn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))/2)
  27   -(jacobi_dn(u,m)*'jacobi_sn(u,m)^2/(2*(1-m)));
  28
  29 diff(inverse_jacobi_sn(u,m),u);
  30 1/(sqrt(1-u^2)*sqrt(1-m*u^2));
  31
  32 diff(inverse_jacobi_sn(u,m),m);
  33 ((elliptic_e(asin(u),m)-(1-m)*elliptic_f(asin(u),m))/m-(u*sqrt(1-u^2)/sqrt(1-m*u^2)))/(1-m);
  34
  35 diff(inverse_jacobi_cn(u,m),u);
  36 -(1/(sqrt(1-u^2)*sqrt(m*u^2-m+1)));
  37
  38 diff(inverse_jacobi_cn(u,m),m);
  39 ((elliptic_e(asin(sqrt(1-u^2)),m)-(1-m)*elliptic_f(asin(sqrt(1-u^2)),m))/m
  40    -(sqrt(1-u^2)*abs(u)/sqrt(1-m*(1-u^2))))/(1-m);
  41
  42 diff(inverse_jacobi_dn(u,m),u);
  43 1/(sqrt(1-u^2)*sqrt(u^2+m-1));
  44
  45 diff(inverse_jacobi_dn(u,m),m);
  46 ((elliptic_e(asin(sqrt(1-u^2)/sqrt(m)),m)-(1-m)*elliptic_f(asin(sqrt(1-u^2)/sqrt(m)),m))/m
  47   -sqrt(1-(1-u^2)/m)*sqrt(1-u^2)/(sqrt(m)*abs(u)))/(1-m)
  48 -(sqrt(1-u^2)/(2*m^(3/2)*sqrt(1-(1-u^2)/m)*abs(u)));
  49
  50 diff(elliptic_e(phi,m),phi);
  51 sqrt(1-m*sin(phi)^2);
  52
  53 diff(elliptic_e(phi,m),m);
  54 (elliptic_e(phi,m)-elliptic_f(phi,m))/(2*m);
  55
  56 diff(elliptic_f(phi,m),phi);
  57 1/sqrt(1-m*sin(phi)^2);
  58
  59 diff(elliptic_f(phi,m),m);
  60 ((elliptic_e(phi,m)-(1-m)*elliptic_f(phi,m))/m
  61    -(cos(phi)*sin(phi)/sqrt(1-m*sin(phi)^2)))
  62  /(2*(1-m));
  63
  64 diff(elliptic_kc(m),m);
  65 (elliptic_ec(m)-(1-m)*elliptic_kc(m))/(2*(1-m)*m);
  66
  67 diff(elliptic_ec(m),m);
  68 (elliptic_ec(m)-elliptic_kc(m))/(2*m);
  69
  70 /* Integrals */
  71
  72 integrate(jacobi_sn(u,m),u); /* A&S 16.24.1 */
  73 log(jacobi_dn(u,m)-sqrt(m)*jacobi_cn(u,m))/sqrt(m);
  74
  75 integrate(jacobi_cn(u,m),u); /* A&S 16.24.2 */
  76 acos(jacobi_dn(u,m))/sqrt(m);
  77
  78 integrate(jacobi_dn(u,m),u); /* A&S 16.24.3 */
  79 asin(jacobi_sn(u,m));
  80
  81 integrate(jacobi_cd(u,m),u); /* A&S 16.24.4 */
  82 log(sqrt(m)*jacobi_sd(u,m)+jacobi_nd(u,m))/sqrt(m);
  83
  84 /* Use functions.wolfram.com 09.35.21.001.01, not A&S 16.24.5 */
  85 integrate(jacobi_sd(u,m),u);
  86 -(sqrt(1-m*jacobi_cd(u,m)^2)*jacobi_dn(u,m)*asin(sqrt(m)*jacobi_cd(u,m))
  87   /((1-m)*sqrt(m)));
  88
  89 /* Use functions.wolfram.com 09.32.21.0001.01, not A&S 16.24.6 */
  90 integrate(jacobi_nd(u,m),u);
  91 sqrt(1-jacobi_cd(u,m)^2)*acos(jacobi_cd(u,m))/((1-m)*jacobi_sd(u,m));
  92
  93 integrate(jacobi_dc(u,m),u); /* A&S 16.24.7 */
  94 log(jacobi_sc(u,m)+jacobi_nc(u,m));
  95
  96 integrate(jacobi_nc(u,m),u); /* A&S 16.24.8 */
  97 log(sqrt(1-m)*jacobi_sc(u,m)+jacobi_dc(u,m))/sqrt(1-m);
  98
  99 integrate(jacobi_sc(u,m),u); /* A&S 16.24.9 */
 100 log(sqrt(1-m)*jacobi_nc(u,m)+jacobi_dc(u,m))/sqrt(1-m);
 101
 102 integrate(jacobi_ns(u,m),u); /* A&S 16.24.10 */
 103 log(jacobi_ds(u,m)-jacobi_cs(u,m));
 104
 105 /* Use functions.wolfram.com 09.30.21.0001.01, not A&S 16.24.11 */
 106 integrate(jacobi_ds(u,m),u);
 107 log((1-jacobi_cn(u,m))/jacobi_sn(u,m));
 108
 109 integrate(jacobi_cs(u,m),u); /* A&S 16.24.12 */
 110 log(jacobi_ns(u,m)-jacobi_ds(u,m));
 111
 112 /* Check the integrals and derivatives by confirming
 113
 114        f(x,m)-diff(integral(x,m),x),x) = constant
 115
 116   Look at Taylor expansion about zero, rather than messing about with
 117   elliptic function.  This was sufficient to find several errors  */
 118 (te(f,n):=ratsimp(
 119            taylor(f(x,m)
 120            -diff(integrate(f(x,m),x),x),x,0,n)),
 121 ti(f,n):=ratsimp(
 122            taylor(integrate(f(x,m),x),x,0,n)
 123             -integrate(taylor(f(x,m),x,0,n-1),x)),
 124 td(f,n):=ratsimp(
 125            taylor(diff(f(x,m),x),x,0,n)
 126             -diff(taylor(f(x,m),x,0,n+1),x)),
 127 /* Compare analytic and numerical integral */
 128 ni(f,x1,x2,m):= block(
 129   [x,n,I,In,Ia,key:1,eps:1.0e-14],
 130   I:integrate(f(x,n),x),
 131   In:quad_qag(f(x,m),x,x1,x2,key),
 132   Ia:subst([x=float(x2),n=float(m)],I)-subst([x=float(x1),n=float(m)],I),
 133   if ( abs(Ia-In[1]) < eps ) then
 134     true
 135   else
 136     [Ia,In[1]]
 137 ),
 138 done);
 139 done;
 140
 141 te(jacobi_sn,4);
 142 0;
 143
 144 te(jacobi_cn,4);
 145 0;
 146
 147 te(jacobi_dn,4);
 148 0;
 149
 150 te(jacobi_cd,4);
 151 0;
 152
 153 assume(m>0,m<1); /* also ok for m<0 and m>0 */
 154 [m > 0, m < 1];
 155 te(jacobi_sd,4);
 156 0;
 157 forget(m>0,m<1);
 158 [m > 0, m < 1];
 159
 160 te(jacobi_nd,4);
 161 0;
 162
 163 te(jacobi_dc,4);
 164 0;
 165
 166 te(jacobi_nc,4);
 167 0;
 168
 169 te(jacobi_sc,4);
 170 0;
 171
 172 /* jacobi_ns, jacobi_ds and jacobi_cs are singular at x=0
 173    Compare numerical and analytic integrals for a single case.
 174 */
 175
 176 ni(jacobi_ns,1,2,1/2);
 177 true;
 178
 179 ni(jacobi_ds,1,2,1/2);
 180 true;
 181
 182 ni(jacobi_cs,1,2,1/2);
 183 true;
 184
 185 kill(te,ti,td,ni);
 186 done;
 187
 188 /* Slightly modified version of test_table taken from rtest_expintegral.mac */
 189
 190 (test_table(func,table,eps) :=
 191 block([badpoints : [],
 192        abserr    : 0,
 193        maxerr    : -1],
 194   for entry in table do
 195     block([z : entry[1],
 196            result, answer],
 197       z : expand(rectform(bfloat(entry[1]))),
 198       result : rectform(apply(func, z)),
 199       answer : expand(rectform(bfloat(entry[2]))),
 200       abserr : abs(result-answer),
 201       maxerr : max(maxerr,abserr),
 202       if abserr > eps then
 203         badpoints : cons ([z,result,answer,abserr],badpoints)
 204     ),
 205   if badpoints # [] then
 206     cons(maxerr,badpoints)
 207   else
 208     badpoints
 209 ),done);
 210 done;
 211
 212 /* These test values come from http://getnet.net/~cherry/testrf.mac */
 213
 214 (mrf:[[[1,2,0],1.3110287771461b0],
 215      [[%i,-%i,0],1.8540746773014b0],
 216      [[0.5,1,0],1.8540746773014b0],
 217      [[%i-1,%i,0],0.79612586584234b0-%i*(1.2138566698365b0)],
 218      [[2,3,4],0.58408284167715b0],
 219      [[%i,-%i,2],1.0441445654064b0],
 220      [[%i-1,%i,1-%i],0.93912050218619b0-%i*(0.53296252018635b0)]],
 221 done);
 222 done;
 223
 224 test_table('carlson_rf, mrf, 1.5b-13);
 225 [];
 226
 227 (mrc:[[[0,1/4],bfloat(%pi)],
 228       [[9/4,2],log(2b0)],
 229       [[0,%i],(1-%i)*(1.1107207345396b0)],
 230       [[-%i,%i],1.2260849569072b0-%i*(0.34471136988768b0)],
 231       [[1/4,-2],log(2b0)/3],
 232       [[%i,-1],0.77778596920447b0+%i*(0.19832484993429b0)],
 233       [[0,1/4],%pi],
 234       [[9/4,2],log(2)],
 235       [[2,1],-log(sqrt(2)-1)],
 236       [[-%i,%i],-log(sqrt(2)-1)/2+ %pi/4-%i*(log(sqrt(2)-1)/2+%pi/4)],
 237       [[1/4,-2],log(2)/3],
 238       [[%i,-1],sqrt(sqrt(2)/4-1/4)*atan(sqrt(sqrt(2)-1))-
 239                sqrt(sqrt(2)/16+1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1)+
 240                %i*(sqrt(sqrt(2)/4+1/4)*atan(sqrt(sqrt(2)-1))+sqrt(sqrt(2)/16-
 241                1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1))],
 242       [[0,1],%pi/2],
 243       [[%i,%i+1],%pi/4+%i*log(sqrt(2)-1)/2]],
 244 done);
 245 done;
 246
 247 test_table('carlson_rc, mrc, 2.5b-14);
 248 [];
 249
 250 (mrj:[[[0,1,2,3],0.77688623778582b0],
 251       [[2,3,4,5],0.14297579667157b0],
 252       [[2,3,4,-1+%i],0.13613945827771b0-%i*(0.38207561624427b0)],
 253       [[%i,-%i,0,2],1.6490011662711b0],
 254       [[-1+%i,-1-%i,1,2],0.9414835884122b0],
 255       [[%i,-%i,0,1-%i],1.8260115229009b0+%i*(1.2290661908643b0)],
 256       [[-1+%i,-1-%i,1,-3+%i],-0.61127970812028b0-%i*(1.0684038390007b0)],
 257       [[-1+%i,-2-%i,-%i,-1+%i],1.8249027393704b0-%i*(1.2218475784827b0)],
 258       [[2,3,4,-0.5],0.24723819703052b0],
 259       [[2,3,4,-5],-0.12711230042964b0]],
 260 done);
 261 done;
 262
 263 test_table('carlson_rj, mrj, 1b-13);
 264 [];
 265
 266 (mrd:[[[0,2,1],1.7972103521034b0],
 267       [[2,3,4],0.16510527294261b0],
 268       [[%i,-%i,2],0.6593385415422b0],
 269       [[0,%i,-%i],1.270819627191b0+%i*(2.7811120159521b0)],
 270       [[0,%i-1,%i],-1.8577235439239b0-%i*(0.96193450888839b0)],
 271       [[-2-%i,-%i,-1+%i],1.8249027393704b0-%i*(1.2218475784827b0)]],
 272 done);
 273 done;
 274
 275 test_table('carlson_rd, mrd, 6e-14);
 276 [];
 277
 278 /* Some tests of the Jacobian elliptic functions.
 279  * Just some tests at random points, to verify that we are accurate.
 280  * The reference values were obtained from Mathematica, but we could
 281  * also just compute the values using a much larger fpprec.
 282  */
 283
 284 test_table('jacobi_sn,
 285            [[[1b0+%i*1b0, .7b0], 1.134045971912365274954b0 + 0.352252346922494477621b0*%i],
 286             [[1b0+%i*1b0, 2b0], 0.98613318109123804740b0 + 0.09521910819727230780b0*%i],
 287             [[1b0+%i*1b0, 2b0+3b0*%i], 0.94467077879445294981b0 - 0.19586410083100945528b0* %i],
 288             [[1.785063352082689082581887b0 *%i + 8.890858759528509578661528b-1,
 289               9.434463451695984398149033b-1 - 1.476052973708684178844821b-1 * %i],
 290              1.345233534539675700312281b0 - 7.599023926615176284214056b-2 * %i]],
 291            1b-15);
 292 [];
 293
 294 test_table('jacobi_cn,
 295            [[[1b0+%i*1b0, .7b0], 0.571496591371764254029b0 - 0.698989917271916772991b0*%i ],
 296             [[1b0+%i*1b0, 2b0], 0.33759463268642431412b0 - 0.27814044708010806377b0*%i],
 297             [[1b0+%i*1b0, 2b0+3b0*%i], -0.52142079485827170824b0 - 0.35485177134179601850b0*%i],
 298             [[100b0, .7b0], 0.93004753815774770476196b0]],
 299            6b-14);
 300 [];
 301
 302 test_table('jacobi_dn,
 303            [[[1b0+%i*1b0, .7b0], 0.62297154372331777630564880787568b0 - 0.448863598031509643267241389621738b0 *%i ],
 304             [[1b0+%i*1b0, 2b0], 0.1913322443206041462495606602242b0 - 0.9815253294150083432282549919753b0 * %i],
 305             [[1b0+%i*1b0, 2b0+3b0*%i], 0.6147387452173944656984656771134b0 - 1.4819401302071697495918834416787b0 * %i],
 306             [[100b0, .7b0], 0.95157337933724324055428565654872978b0],
 307             [[1.785063352082689082581887b0 *%i + 8.890858759528509578661528b-1,
 308               9.434463451695984398149033b-1 - 1.476052973708684178844821b-1 * %i],
 309              -8.617730683333292717095686b-1 *%i - 2.663978258141280808361839b-1]],
 310            2b-15);
 311 [];
 312
 313 /* These routines for cn and dn work well for small (<= 1?) values of
 314  * u and m.  They have known issues for large real values of u.
 315  */
 316 (ascending_transform(u,m) :=
 317   block([root_m : expand(rectform(sqrt(m))), mu, root_mu1, v],
 318     mu : expand(rectform(4*root_m/(1+root_m)^2)),
 319     root_mu1 : expand(rectform((1-root_m)/(1+root_m))),
 320     v : expand(rectform(u/(1+root_mu1))),
 321     [v, mu, root_mu1]),
 322  elliptic_dn_ascending(u,m) :=
 323   if is(abs(m-1) < 4*10^(-fpprec)) then sech(u)
 324   else
 325     block([v, mu, root_mu1, new],
 326       [v, mu, root_mu1] : ascending_transform(u,m),
 327       new : elliptic_dn_ascending(v, mu),
 328       expand(rectform((1-root_mu1)/mu*(new^2 + root_mu1)/new))),
 329  elliptic_cn_ascending(u,m) :=
 330   if is(abs(m-1) < 4*10^(-fpprec)) then sech(u)
 331   else
 332     block([v, mu, root_mu1, new],
 333       [v, mu, root_mu1] : ascending_transform(u,m),
 334       new : elliptic_dn_ascending(v, mu),
 335       expand(rectform((1+root_mu1)/mu*(new^2-root_mu1)/new))),
 336  done);
 337 done;
 338
 339 /* Test with random values for the argument and parameter. */
 340 (test_random(n, eps, testf, truef) :=
 341   block([badpoints : [], maxerr : -1],
 342     for k : 1 thru n do
 343       block([z : bfloat(1-2*random(1d0) + %i * (1-2*random(1d0))),
 344              m : bfloat(1-2*random(1d0) + %i * (1-2*random(1d0))),
 345              ans, expected, abserr],
 346         ans : testf(z, m),
 347         expected : truef(z, m),
 348         abserr : abs(ans-expected),
 349         maxerr : max(maxerr, abserr),
 350         if abserr > eps then
 351           badpoints : cons([[z,m], ans, expected, abserr], badpoints)),
 352     if badpoints # [] then
 353       cons(maxerr, badpoints)
 354     else
 355       badpoints),
 356  done);
 357 done;
 358
 359 test_random(50, 2b-15, 'jacobi_dn, 'elliptic_dn_ascending);
 360 [];
 361
 362 test_random(50, 2b-15, 'jacobi_cn, 'elliptic_cn_ascending);
 363 [];
 364
 365 /* Test elliptic_f by using the fact that
 366
 367 inverse_jacobi_sn(x,m) = elliptic_f(asin(x), m)
 368
 369 */
 370 (test_ef(x,m) := jacobi_sn(elliptic_f(asin(x),m), m), id(x,m):=x, done);
 371 done;
 372
 373 test_random(100, 6b-13, 'test_ef, 'id);
 374 [];
 375
 376 /* Test elliptic_kc These values are from
 377  *
 378  * http://functions.wolfram.com/EllipticIntegrals/EllipticK/03/01/
 379  */
 380
 381 block([oldfpprec : fpprec, fpprec:100],
 382   test_table('elliptic_kc,
 383              [[[1/2], 8*%pi^(3/2)/gamma(-1/4)^2],
 384               [[17-12*sqrt(2)], 2*(2+sqrt(2))*%pi^(3/2)/gamma(-1/4)^2],
 385               [[-1], gamma(1/4)^2/4/sqrt(2*%pi)]],
 386              2b-100));
 387 [];
 388
 389 /* Some tests for specific values */
 390 inverse_jacobi_sn(1,m);
 391 elliptic_kc(m);
 392
 393 inverse_jacobi_sn(x,0);
 394 asin(x);
 395
 396 inverse_jacobi_sn(x,1);
 397 log(tan(asin(x)/2 + %pi/4));
 398
 399 inverse_jacobi_sd(1/sqrt(1-m),m);
 400 elliptic_kc(m);
 401
 402 inverse_jacobi_ds(sqrt(1-m),m);
 403 elliptic_kc(m);
 404
 405 /* elliptic_kc(1) is undefined */
 406 errcatch(elliptic_kc(1));
 407 [];
 408 errcatch(elliptic_kc(1.0));
 409 [];
 410 errcatch(elliptic_kc(1.0b0));
 411 [];
 412
 413 /* Test noun/verb results from elliptic functions */
 414 diff(inverse_jacobi_sn(x,0),x);
 415 1/sqrt(1-x^2);
 416
 417 diff(elliptic_pi(4/3,phi,0),phi);
 418 sec(phi)^2/(1-tan(phi)^2/3);
 419
 420 diff(elliptic_pi(3/4,phi,0),phi);
 421 sec(phi)^2/(1+tan(phi)^2/4);
 422
 423 diff(elliptic_pi(1,phi,0),phi);
 424 sec(phi)^2;
 425
 426 /* Special values */
 427 elliptic_ec(1);
 428 1;
 429
 430 elliptic_ec(0);
 431 %pi/2;
 432
 433 elliptic_ec(1/2);
 434 gamma(3/4)^2/(2*sqrt(%pi))+%pi^(3/2)/(4*gamma(3/4)^2);
 435
 436 elliptic_ec(-1);
 437 sqrt(2)*elliptic_ec(1/2);
 438
 439 /* Test periodicity of elliptic_e */
 440 elliptic_e(x, 1);
 441 2*round(x/%pi)+sin(x);
 442
 443 elliptic_e(3,1/3);
 444 elliptic_e(3-%pi,1/3)+2*elliptic_ec(1/3);
 445
 446 /* Bug #2629: elliptic_kc(3.0) not accurate */
 447 test_table('elliptic_kc, [[[3.0], elliptic_kc(3b0)]], 1b-15);
 448 [];
 449
 450 /* Bug #2630: inverse_jacobi_cn(-2.0, 3.0) generates an error */
 451 test_table('inverse_jacobi_cn, [[[-2.0, 3.0], 2.002154760912212-3.202503914656527*%i]],
 452   1b-15);
 453 [];
 454
 455 /* Bug #2615: Numeric evaluation of inverse Jacobi elliptic functions is wrong for some inputs
 456 */
 457 is(abs(jacobi_dn(inverse_jacobi_dn(-2.0,3.0), 3.0) + 2) < 1d-14);
 458 true;
 459