R Under development (unstable) (2024-07-28 r86931 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #### Bessel functions for (very) large x/z and/or nu > > #### TODO: move some stuff from ./IJKY.R to here <<<<< > > ### The argument range is from David Scott, for K_nu(.) case: > > library(Bessel) > xs <- 100*(1:7) > nus <- 450 + 10*(0:9) > d.xn <- expand.grid(nu = nus, x = xs) > > M <- with(d.xn, + cbind(K = besselK(x,nu), + K_exp = besselK(x,nu, expon.scaled = TRUE), + K_nA.2 = besselK.nuAsym(x, nu, log = TRUE, k.max=2), + K_nA.3 = besselK.nuAsym(x, nu, log = TRUE, k.max=3), + K_nA.4 = besselK.nuAsym(x, nu, log = TRUE, k.max=4)) + ) > ## Transform into nicely labelled 3d array : > A <- M > datt <- attr(d.xn, "out.attrs") > dim(A) <- c(datt$dim, ncol(M)) > dimnames(A) <- c(datt$dimnames, list(colnames(M))) > A , , = K x nu x=100 x=200 x=300 x=400 x=500 nu=450 2.212120e+230 6.885989e+87 2.220613e-03 1.671939e-74 4.181206e-136 nu=460 9.596364e+239 4.074204e+94 3.731509e+02 2.859791e-70 1.459419e-132 nu=470 5.149411e+249 2.942827e+101 7.523288e+07 5.763696e-66 5.901696e-129 nu=480 3.403145e+259 2.585673e+108 1.814894e+13 1.365987e-61 2.760954e-125 nu=490 2.758469e+269 2.753982e+115 5.224585e+18 3.799358e-57 1.492104e-121 nu=500 2.731383e+279 3.543783e+122 1.790087e+24 1.237785e-52 9.301952e-118 nu=510 3.291180e+289 5.491328e+129 7.281307e+29 4.714314e-48 6.679872e-114 nu=520 4.808000e+299 1.021461e+137 3.507321e+35 2.095126e-43 5.517863e-110 nu=530 Inf 2.273910e+144 1.995796e+41 1.084462e-38 5.235745e-106 nu=540 Inf 6.040096e+151 1.338442e+47 6.525871e-34 5.698930e-102 x nu x=600 x=700 nu=450 2.530096e-192 3.263320e-245 nu=460 2.757857e-189 1.468056e-242 nu=470 3.431151e-186 7.441531e-240 nu=480 4.867215e-183 4.246967e-237 nu=490 7.863836e-180 2.726795e-234 nu=500 1.445585e-176 1.968080e-231 nu=510 3.020317e-173 1.595546e-228 nu=520 7.164844e-170 1.451815e-225 nu=530 1.927768e-166 1.481526e-222 nu=540 5.876892e-163 1.694194e-219 , , = K_exp x nu x=100 x=200 x=300 x=400 x=500 nu=450 5.946436e+273 4.975798e+174 4.313378e+127 8.729978e+99 5.868708e+81 nu=460 2.579615e+283 2.944009e+181 7.248182e+132 1.493231e+104 2.048429e+85 nu=470 1.384222e+293 2.126479e+188 1.461343e+138 3.009496e+108 8.283575e+88 nu=480 9.148051e+302 1.868400e+195 3.525298e+143 7.132459e+112 3.875254e+92 nu=490 Inf 1.990020e+202 1.014837e+149 1.983823e+117 2.094305e+96 nu=500 Inf 2.560728e+209 3.477111e+154 6.463058e+121 1.305615e+100 nu=510 Inf 3.968019e+216 1.414340e+160 2.461565e+126 9.375816e+103 nu=520 Inf 7.381054e+223 6.812712e+165 1.093964e+131 7.744830e+107 nu=530 Inf 1.643121e+231 3.876686e+171 5.662485e+135 7.348851e+111 nu=540 Inf 4.364557e+238 2.599826e+177 3.407464e+140 7.998974e+115 x nu x=600 x=700 nu=450 9.546102e+68 3.309763e+59 nu=460 1.040545e+72 1.488950e+62 nu=470 1.294580e+75 7.547439e+64 nu=480 1.836410e+78 4.307411e+67 nu=490 2.967041e+81 2.765603e+70 nu=500 5.454223e+84 1.996090e+73 nu=510 1.139572e+88 1.618253e+76 nu=520 2.703310e+91 1.472478e+79 nu=530 7.273509e+94 1.502611e+82 nu=540 2.217363e+98 1.718306e+85 , , = K_nA.2 x nu x=100 x=200 x=300 x=400 x=500 x=600 x=700 nu=450 530.3885 202.2544 -6.109972 -169.87731 -311.7210 -441.1681 -562.9506 nu=460 552.5792 217.8477 5.921983 -160.13021 -303.5632 -434.1741 -556.8417 nu=470 574.9826 233.6405 18.136099 -150.21904 -295.2582 -427.0479 -550.6133 nu=480 597.5942 249.6292 30.529633 -140.14581 -286.8076 -419.7906 -544.2665 nu=490 620.4101 265.8103 43.099907 -129.91252 -278.2126 -412.4030 -537.8018 nu=500 643.4260 282.1806 55.844306 -119.52110 -269.4748 -404.8865 -531.2201 nu=510 666.6383 298.7366 68.760278 -108.97348 -260.5956 -397.2419 -524.5222 nu=520 690.0432 315.4754 81.845331 -98.27155 -251.5764 -389.4703 -517.7088 nu=530 713.6371 332.3938 95.097032 -87.41715 -242.4185 -381.5728 -510.7808 nu=540 737.4166 349.4888 108.513006 -76.41212 -233.1234 -373.5503 -503.7389 , , = K_nA.3 x nu x=100 x=200 x=300 x=400 x=500 x=600 x=700 nu=450 530.3885 202.2544 -6.109972 -169.87731 -311.7210 -441.1681 -562.9506 nu=460 552.5792 217.8477 5.921983 -160.13021 -303.5632 -434.1741 -556.8417 nu=470 574.9826 233.6405 18.136099 -150.21904 -295.2582 -427.0479 -550.6133 nu=480 597.5942 249.6292 30.529633 -140.14581 -286.8076 -419.7906 -544.2665 nu=490 620.4101 265.8103 43.099907 -129.91252 -278.2126 -412.4030 -537.8018 nu=500 643.4260 282.1806 55.844306 -119.52110 -269.4748 -404.8865 -531.2201 nu=510 666.6383 298.7366 68.760278 -108.97348 -260.5956 -397.2419 -524.5222 nu=520 690.0432 315.4754 81.845331 -98.27155 -251.5764 -389.4703 -517.7088 nu=530 713.6371 332.3938 95.097032 -87.41715 -242.4185 -381.5728 -510.7808 nu=540 737.4166 349.4888 108.513006 -76.41212 -233.1234 -373.5503 -503.7389 , , = K_nA.4 x nu x=100 x=200 x=300 x=400 x=500 x=600 x=700 nu=450 530.3885 202.2544 -6.109972 -169.87731 -311.7210 -441.1681 -562.9506 nu=460 552.5792 217.8477 5.921983 -160.13021 -303.5632 -434.1741 -556.8417 nu=470 574.9826 233.6405 18.136099 -150.21904 -295.2582 -427.0479 -550.6133 nu=480 597.5942 249.6292 30.529633 -140.14581 -286.8076 -419.7906 -544.2665 nu=490 620.4101 265.8103 43.099907 -129.91252 -278.2126 -412.4030 -537.8018 nu=500 643.4260 282.1806 55.844306 -119.52110 -269.4748 -404.8865 -531.2201 nu=510 666.6383 298.7366 68.760278 -108.97348 -260.5956 -397.2419 -524.5222 nu=520 690.0432 315.4754 81.845331 -98.27155 -251.5764 -389.4703 -517.7088 nu=530 713.6371 332.3938 95.097032 -87.41715 -242.4185 -381.5728 -510.7808 nu=540 737.4166 349.4888 108.513006 -76.41212 -233.1234 -373.5503 -503.7389 > > ## Compare the different approximation levels k.max > stopifnot( + all.equal(M[,3], M[,4], tol=1e-12)# 2.826 e-13 + , + all.equal(M[,4], M[,5], tol=2e-15)# 5.357 e-16 + , + all.equal(M[,"K"], exp(M[,5]), tol= 1e-12)# on log.scale: 2e-16 ! + ) > > > cat('Time elapsed: ', proc.time(),'\n') # for ''statistical reasons'' Time elapsed: 0.62 0.14 0.71 NA NA > > proc.time() user system elapsed 0.62 0.14 0.71