https://github.com/cran/pracma
Tip revision: 26e049d70b4a1c237987e260cba68f6a9413736c authored by Hans W. Borchers on 09 April 2019, 04:10:07 UTC
version 2.2.5
version 2.2.5
Tip revision: 26e049d
quadinf.R
##
## q u a d i n f . R Infinite Integrals
##
quadinf <- function(f, xa, xb, tol = 1e-12, ...) {
stopifnot(is.numeric(xa), length(xa) == 1,
is.numeric(xb), length(xb) == 1,
is.numeric(tol), length(tol) == 1)
fun <- match.fun(f)
f <- function(x) fun(x, ...)
if (xa == xb) {
return(list(Q = 0.0, reltol = 0.0))
} else if (xa > xb) {
stop("For the integration limits 'xa <= xb' is required.")
}
u <- -1; v <- 1
if (is.finite(xa) && is.finite(xb)) {
duv <- (xb-xa)/(v-u)
ff <- function(y) duv * f(xa + duv*(y-u))
} else if (is.finite(xa) && is.infinite(xb)) {
ff <- function(y) {
duv <- (v-u)/(1-y)^2
z <- (y-u)/(v-u)
duv * f(xa + z/(1-z))
}
} else if (is.infinite(xa) && is.finite(xb)) {
xa <- -xb
ff <- function(y) {
duv = (v-u)/(1-y)^2
z <- (y-u)/(v-u)
duv * f(-(xa + z/(1 - z)))
}
} else if (is.infinite(xa) && is.infinite(xb)) {
ff <- function(y) {
z <- pi * (y - u)/(v - u) - 0.5
duv <- (1/cos(z)^2)*pi / (v - u)
duv * f(tan(z))
}
} else {
stop("Other integration domains will be treated later.")
}
cc <- .quadinf_pre()
xx <- cc$nodes; ww <- cc$weights
h <- 1.0/2.0
x <- xx[[1]]; w <- ww[[1]]
s <- w[7] * ff(x[7])
for (j in 1:6) {
s <- s + w[j] * (ff(x[j]) + ff(-x[j]))
}
Q <- s * h
for (k in 2:7) {
x <- xx[[k]]; w <- ww[[k]]
s <- 0.0
for (j in 1:length(w)) {
s <- s + w[j] * (ff(x[j]) + ff(-x[j]))
}
h <- h/2.0
newQ <- s*h + Q/2.0
delta <- abs(newQ - Q)
Q <- newQ
if (delta < tol) break
}
return(list(Q = Q, relerr = delta, niter = k))
}
.quadinf_pre <- function() {
nodes <- list(
c(0.99999999999995703, 0.99999998887566488, 0.99997747719246155, 0.99751485645722437, 0.95136796407274693, 0.67427149224843574, 0),
c(0.99999999995285638, 0.99999920473711468, 0.9996882640283532, 0.98704056050737687, 0.85956905868989653, 0.37720973816403414),
c(0.99999999999823208, 0.99999999914270499, 0.99999989278161239, 0.99999531604122049, 0.9999093846951439, 0.99906519645578584,
0.99405550663140207, 0.97396686819567735, 0.91487926326457458, 0.78060743898320029, 0.53914670538796772, 0.19435700332493541),
c(0.99999999999970801, 0.99999999999039391, 0.99999999978973275, 0.99999999678719909, 0.99999996423908089, 0.99999969889415252,
0.99999801714059533, 0.99998948201481841, 0.99995387100562794, 0.99982882207287493, 0.99945143443527451, 0.99845420876769764,
0.99610866543750853, 0.99112699244169877, 0.98145482667733508, 0.96411216422354729, 0.93516085752198463, 0.88989140278426015,
0.8233170055064023, 0.73101803479256144, 0.61027365750063889, 0.46125354393958568, 0.28787993274271589, 0.09792388528783233),
c(0.99999999999988631, 0.99999999999927147, 0.99999999999582456, 0.99999999997845523, 0.99999999989927768, 0.99999999957078767,
0.99999999832336184, 0.99999999396413419, 0.99999997987450318, 0.99999993755407834, 0.9999998188937127, 0.99999950700571938,
0.99999873547186591, 0.9999969324491903, 0.99999293787666288, 0.99998451990227077, 0.99996759306794336, 0.99993501992508238,
0.99987486504878031, 0.99976797159956077, 0.9995847503515175, 0.99928111192179192, 0.9987935342988058, 0.99803333631543367,
0.99688031812819178, 0.9951760261553273, 0.99271699719682727, 0.98924843109013383, 0.98445883116743083, 0.97797623518666488,
0.96936673289691733, 0.95813602271021359, 0.94373478605275707, 0.92556863406861256, 0.90301328151357385, 0.87543539763040867,
0.84221924635075684, 0.80279874134324125, 0.75669390863372987, 0.70355000514714194, 0.64317675898520466, 0.5755844906351516,
0.50101338937930906, 0.41995211127844712, 0.33314226457763807, 0.24156631953888363, 0.14641798429058792, 0.049055967305077885),
c(0.99999999999992983, 0.99999999999981715, 0.99999999999953715, 0.99999999999886124, 0.99999999999727396, 0.9999999999936463,
0.99999999998556877, 0.99999999996803313, 0.99999999993088839, 0.99999999985405608, 0.9999999996987543, 0.99999999939177686,
0.9999999987979834, 0.99999999767323666, 0.99999999558563357, 0.99999999178645604, 0.9999999850030763, 0.99999997311323585,
0.99999995264266439, 0.99999991800479471, 0.99999986037121458, 0.9999997660233324, 0.9999996139885502, 0.9999993727073353,
0.99999899541068993, 0.99999841381096466, 0.99999752962380506, 0.9999962033471661, 0.99999423962761658, 0.9999913684483448,
0.99998722128200057, 0.99998130127012064, 0.99997294642523216, 0.99996128480785662, 0.99994518061445858, 0.99992317012928922,
0.99989338654759252, 0.99985347277311132, 0.99980048143113831, 0.99973076151980844, 0.9996398313456003, 0.99952223765121717,
0.99937140114093759, 0.99917944893488586, 0.99893703483351215, 0.99863314864067743, 0.99825491617199624, 0.99778739195890642,
0.99721334704346865, 0.99651305464025375, 0.99566407681695313, 0.99464105571251116, 0.99341551316926402, 0.99195566300267757,
0.99022624046752772, 0.98818835380074255, 0.98579936302528337, 0.9830127914811011, 0.97977827580061572, 0.97604156025657673,
0.97174454156548729, 0.96682537031235583, 0.9612186151511164, 0.9548554958050226, 0.94766419061515306, 0.93957022393327472,
0.93049693799715338, 0.92036605303195274, 0.90909831816302034, 0.89661425428007602, 0.88283498824466888, 0.86768317577564591,
0.85108400798784867, 0.83296629391941079, 0.81326360850297374, 0.79191549237614201, 0.76886868676824649, 0.74407838354734734,
0.71750946748732403, 0.68913772506166759, 0.65895099174335003, 0.62695020805104285, 0.59315035359195312, 0.55758122826077816,
0.52028805069123008, 0.48133184611690499, 0.44078959903390086, 0.39875415046723772, 0.3553338251650745, 0.31065178055284592,
0.26484507658344791, 0.218063473469712, 0.1704679723820105, 0.12222912220155764, 0.073525122985671293, 0.024539763574649157),
c(0.99999999999994504, 0.99999999999991063, 0.99999999999985567, 0.99999999999976874, 0.99999999999963207, 0.9999999999994188,
0.9999999999990884, 0.99999999999857991, 0.99999999999780287, 0.99999999999662348, 0.99999999999484512, 0.99999999999218125,
0.99999999998821665, 0.99999999998235345, 0.99999999997373656, 0.9999999999611503, 0.99999999994287725, 0.99999999991650601,
0.99999999987867016, 0.99999999982469845, 0.99999999974814635, 0.99999999964017294, 0.99999999948871821, 0.99999999927742189,
0.99999999898420988, 0.99999999857945798, 0.99999999802361939, 0.99999999726417366, 0.99999999623172697, 0.99999999483505064,
0.99999999295480446, 0.9999999904356327, 0.9999999870762647, 0.99999998261717338, 0.99999997672526708, 0.99999996897499022,
0.9999999588251014, 0.99999994559027294, 0.99999992840651408, 0.99999990618926848, 0.99999987758285502, 0.99999984089973593,
0.99999979404787598, 0.99999973444423285, 0.99999965891215925, 0.9999995635602319, 0.99999944363972881, 0.99999929337766846,
0.99999910578199569, 0.99999887241516183, 0.99999858313198442, 0.99999822577731134, 0.99999778583863519, 0.9999972460484261,
0.99999658593057072, 0.99999578128492839, 0.99999480360364879, 0.99999361941253884, 0.99999218953043356, 0.99999046823921378,
0.99998840235683317, 0.99998593020547477, 0.99998298046675649, 0.99997947091575146, 0.99997530702549198, 0.99997038043358921,
0.99996456726262606, 0.99995772628607793, 0.99994969693168922, 0.99994029711447929, 0.99992932089188413, 0.99991653593395113,
0.99990168080200281, 0.99988446202976611, 0.99986455100163496, 0.99984158062348205, 0.99981514178227304, 0.99978477959165113,
0.99974998942164883, 0.99971021271175098, 0.99966483256766459, 0.99961316914334686, 0.99955447481109672, 0.99948792912382323,
0.99941263357495247, 0.9993276061628299, 0.99923177576789879, 0.99912397635238981, 0.99900294099372877, 0.99886729576436373,
0.99871555347220875, 0.99854610727741266, 0.99835722420266371, 0.99814703855574838, 0.99791354528457699, 0.99765459328637784,
0.99736787869423538, 0.99705093816560908, 0.99670114219891259, 0.99631568850566088, 0.99589159546710015, 0.99542569570562278,
0.99491462980263878, 0.9943548401959208, 0.99374256529076155, 0.99307383382058434, 0.99234445949391814, 0.99155003596589164,
0.990685932173615, 0.989747288075987, 0.9887290108396023, 0.98762577151351061, 0.98643200223660843, 0.98514189402239793,
0.98374939516673121, 0.98224821032494103, 0.98063180030544861, 0.97889338262749392, 0.97702593289105788, 0.97502218700730625,
0.97287464433796378, 0.97057557179190035, 0.96811700892685626, 0.96549077410362039, 0.96268847173907823, 0.9597015007033366,
0.95652106390458091, 0.95313817910339371, 0.94954369099593627, 0.94572828460263214, 0.94168249999576925, 0.93739674839571907,
0.93286132966124002, 0.92806645119455511, 0.9230022482765543, 0.91765880584154924, 0.91202618169448679, 0.90609443116640409,
0.89985363319617007, 0.89329391781821088, 0.88640549502697863, 0.87917868497938934, 0.87160394948637765, 0.86367192473410515,
0.85537345516427143, 0.84669962843146362, 0.83764181134359994, 0.82819168667935705, 0.81834129076409778, 0.80808305167333849,
0.79740982792031223, 0.78631494747181951, 0.77479224692443527, 0.76283611066139234, 0.7504415097992404, 0.73760404072282515,
0.72431996299740742, 0.71058623643800578, 0.69640055710845916, 0.68176139201642727, 0.66666801226573758, 0.65112052442430413,
0.63511989986442174, 0.61866800183272797, 0.60176761000963341, 0.58442244232266394, 0.56663717378501877, 0.54841745213979232,
0.52976991010177454, 0.51070217400255813, 0.49122286866081144, 0.47134161831799842, 0.45106904350045196, 0.43041675369143706,
0.40939733572152948, 0.3880243378121177, 0.36631224923490407, 0.34427647557970487, 0.32193330965336914, 0.29929989806396046,
0.27639420357617861, 0.25323496335600021, 0.22984164325436074, 0.20623438831102875, 0.18243396969028913, 0.1584617282892995,
0.13433951528767221, 0.11008962993262801, 0.085734754877651045, 0.061297889413659976, 0.036802280950025079, 0.012271355118082201)
)
weights = list(
c(1.3581784274539089e-012, 2.1431204556943039e-007, 0.00026620051375271687, 0.018343166989927839, 0.23002239451478868,
0.96597657941230108, 1.5707963267948966),
c(1.1631165814255782e-009, 1.1983701363170719e-005, 0.0029025177479013132, 0.076385743570832304, 0.53107827542805397,
1.3896147592472563),
c(4.9378538776631926e-011, 1.8687282268736407e-008, 1.8263320593710658e-006, 6.2482559240744075e-005, 0.00094994680428346862,
0.0077426010260642402, 0.039175005493600777, 0.13742210773316771, 0.36046141846934365, 0.7374378483615478, 1.1934630258491568,
1.523283718634705),
c(8.6759314149796041e-012, 2.5216347918530147e-010, 4.8760060974240624e-009, 6.583518512718339e-008, 6.4777566035929716e-007,
4.8237182032615495e-006, 2.8110164327940134e-005, 0.00013205234125609973, 0.00051339382406790333, 0.0016908739981426396,
0.004816298143928463, 0.012083543599157953, 0.027133510013712, 0.055289683742240581, 0.10343215422333289, 0.17932441211072828,
0.29024067931245418, 0.44083323627385823, 0.63040513516474361, 0.85017285645662, 1.0816349854900702, 1.2974757504249779,
1.4660144267169657, 1.5587733555333301),
c(3.4841937670261058e-012, 2.0989335404511467e-011, 1.130605534749468e-010, 5.4828357797094976e-010, 2.409177325647594e-009,
9.649888896108962e-009, 3.543477717142195e-008, 1.199244278290277e-007, 3.7595411862360629e-007, 1.0968835125901263e-006,
2.9916615878138786e-006, 7.6595758525203149e-006, 1.8481813599879215e-005, 4.2183183841757599e-005, 9.1390817490710112e-005,
0.00018856442976700316, 0.00037166693621677759, 0.00070185951568424226, 0.0012733279447082382, 0.0022250827064786423,
0.003754250977431834, 0.0061300376320830297, 0.0097072237393916877, 0.01493783509605013, 0.022379471063648473, 0.032698732726609031,
0.046668208054846609, 0.065155533432536203, 0.089103139240941459, 0.11949741128869591, 0.15732620348436613, 0.20352399885860173,
0.25890463951405351, 0.32408253961152889, 0.39938474152571712, 0.48475809121475538, 0.57967810308778756, 0.68306851634426369,
0.79324270082051662, 0.90787937915489525, 1.0240449331118113, 1.1382722433763053, 1.2467012074518575, 1.3452788847662516,
1.4300083548722995, 1.4972262225410362, 1.543881116176959, 1.5677814313072218),
c(2.1835922099233607e-012, 5.518236946817488e-012, 1.3542512912336273e-011, 3.230446433325236e-011, 7.4967397573818219e-011,
1.6939457789411645e-010, 3.7299501843052787e-010, 8.0099784479729664e-010, 1.6788897682161906e-009, 3.437185674465009e-009,
6.8784610955899e-009, 1.3464645522302038e-008, 2.5799568229535891e-008, 4.8420950198072366e-008, 8.9071395140242379e-008,
1.6069394579076223e-007, 2.8449923659159806e-007, 4.9458288702754198e-007, 8.4473756384859861e-007, 1.4183067155493917e-006,
2.3421667208528095e-006, 3.8061983264644897e-006, 6.0899100320949032e-006, 9.598194128378471e-006, 1.4908514031870607e-005,
2.2832118109036146e-005, 3.4492124759343198e-005, 5.1421497447658797e-005, 7.5683996586201475e-005, 0.00011002112846666696,
0.00015802788400701192, 0.00022435965205008549, 0.00031497209186021199, 0.00043739495615911686, 0.00060103987991147413,
0.00081754101332469483, 0.0011011261134519382, 0.0014690143599429789, 0.0019418357759843675, 0.0025440657675291729,
0.00330446699403483, 0.0042565295990178572, 0.0054388997976239977, 0.0068957859690660034, 0.0086773307495391812, 0.010839937168255907,
0.01344653660528573, 0.016566786254247574, 0.020277183817500124, 0.024661087314753281, 0.029808628117310124, 0.035816505604196434,
0.042787652157725675, 0.050830757572570467, 0.060059642358636298, 0.070592469906866989, 0.082550788110701726, 0.096058391865189455,
0.11123999898874452, 0.12821973363120098, 0.14711941325785691, 0.16805663794826914, 0.19114268413342747, 0.21648020911729615,
0.24416077786983989, 0.2742622296890681, 0.3068459094179169, 0.34195379592301678, 0.37960556938665158, 0.41979566844501548,
0.46249039805536774, 0.50762515883190806, 0.55510187800363342, 0.6047867305784036, 0.6565082461316275, 0.71005590120546891,
0.76517929890895608, 0.82158803526696467, 0.87895234555278201, 0.93690461274566783, 0.99504180404613263, 1.0529288799552665,
1.1101031939653403, 1.1660798699324344, 1.2203581095793581, 1.2724283455378627, 1.3217801174437727, 1.3679105116808963,
1.4103329714462589, 1.4485862549613224, 1.482243297885538, 1.5109197230741696, 1.5342817381543032, 1.5520531698454121,
1.5640214037732321, 1.5700420292795931),
c(1.7237644036042717e-012, 2.7608587671398282e-012, 4.3888651899779303e-012, 6.9255280152681376e-012, 1.0849161834337119e-011,
1.6874519343915022e-011, 2.6061960502805292e-011, 3.9973389519930265e-011, 6.0893387064380662e-011, 9.2140564226518881e-011,
1.3850287525834143e-010, 2.0684203539029217e-010, 3.0692702078723327e-010, 4.525757610415382e-010, 6.6320863470162091e-010,
9.6594851858099294e-010, 1.3984414312565442e-009, 2.0126210218867171e-009, 2.8797013464471103e-009, 4.0967579139685196e-009,
5.7953495730966062e-009, 8.1527464828576535e-009, 1.1406465554850481e-008, 1.5872978090968231e-008, 2.1971648917089348e-008,
3.0255196464899464e-008, 4.1448233558505381e-008, 5.649576387167062e-008, 7.6623873974817583e-008, 1.0341528040745607e-007,
1.3890286996035951e-007, 1.8568491370025131e-007, 2.4706624511249697e-007, 3.2723037330517557e-007, 4.314482558716538e-007,
5.6633028402050766e-007, 7.401289349090861e-007, 9.6310052117659242e-007, 1.247935512117745e-006, 1.6102680094507651e-006,
2.0692761257364667e-006, 2.6483862254043407e-006, 3.3760952348055106e-006, 4.2869264939998223e-006, 5.4225358918240326e-006,
6.8329862774218501e-006, 8.5782093537079076e-006, 1.072967540685234e-005, 1.3372292284532372e-005, 1.6606555976508335e-005,
2.0550975944934244e-005, 2.5344798968850108e-005, 3.115105567744835e-005, 3.8159954120245754e-005, 4.6592644630494618e-005,
5.6705379853932874e-005, 6.8794093113496371e-005, 8.3199417240036467e-005, 0.00010031216460112419, 0.00012057928729056995,
0.00014451033429094585, 0.00017268441988592921, 0.00020575771467998817, 0.00024447146728689159, 0.00028966056108877171,
0.00034226260646297687, 0.00040332756454954089, 0.00047402789401816719, 0.00055566920742578549, 0.00064970141867429037,
0.00075773035782735918, 0.00088152982417296992, 0.001023054042974541, 0.0011844504858902771, 0.0013680730096097039,
0.001576495261910556, 0.0018125242991288641, 0.0020792143540086659, 0.0023798806881004382, 0.0027181134583502261,
0.0030977915233008206, 0.0035230961104429862, 0.0039985242627335309, 0.0045289019791562883, 0.0051193969614537812,
0.0057755308768065484, 0.0065031910442825787, 0.0073086414513132475, 0.0081985330052611986, 0.0091799129243109005,
0.010260233171410768, 0.011447357834799755, 0.012749569358731162, 0.01417557352833046, 0.015734503113059795, 0.01743592007397746,
0.019289816240845598, 0.02130661236612607, 0.023497155463988527, 0.025872714343617643, 0.028444973247334235, 0.03122602350533174,
0.034228353120175692, 0.037464834195628967, 0.040948708125867296, 0.04469356846276553, 0.048713341380701429, 0.053022263660287172,
0.05763485811465472, 0.062565906384455097, 0.067830419030657965, 0.073443602857638762, 0.079420825403008155, 0.085777576535268185,
0.092529427105778453, 0.099691984607788567, 0.10728084580255642, 0.11531154628093733, 0.12379950693841295, 0.13275997735244552,
0.14220797606340183, 0.15215822777419971, 0.16262509749938431, 0.17362252171163128, 0.18516393655277549, 0.19726220319743717,
0.20992953048020752, 0.22317739492218458, 0.23701645831941898, 0.25145648308451035, 0.26650624556313512, 0.28217344757959606,
0.29846462649944494, 0.31538506413267808, 0.3329386948377478, 0.35112801322442522, 0.36995398189211387, 0.38941593967921201,
0.4095115109381936, 0.43023651638978894, 0.45158488614754488, 0.47354857554061397, 0.49611748439731618, 0.51927938048424627,
0.54301982782484282, 0.5673221206467387, 0.59216722372820685, 0.6175337199299088, 0.6433977657082528, 0.66973305541029093,
0.69651079514654679, 0.72369968702683019, 0.75126592452435681, 0.77917319970479793, 0.80738272301876302, 0.83585325630826712,
0.86454115961966893, 0.89340045234719589, 0.92238288915245514, 0.9514380510163527, 0.98051345168085502, 1.0095546596294298,
1.0385054356373922, 1.0673078857975038, 1.0959026297929721, 1.1242289840506032, 1.1522251592625474, 1.1798284716173337,
1.206975566931308, 1.2336026567219467, 1.2596457651166706, 1.2850409853467271, 1.3097247444374396, 1.333634074575756,
1.3567068895156054, 1.3788822642731373, 1.4001007162694445, 1.4203044859996912, 1.4394378152464069, 1.4574472208125486,
1.4742817617280797, 1.4898932978832971, 1.5042367380636772, 1.5172702754050547, 1.5289556083545806, 1.5392581453118817,
1.5481471912355573, 1.5555961146316604, 1.5615824934918106, 1.5660882389174613, 1.5690996953516689, 1.570607716538275)
)
return(list(nodes = nodes, weights = weights))
}
# quadinf <- function(f, xa, xb, tol = .Machine$double.eps^0.5,
# method = NULL, ...) {
# stopifnot(is.numeric(xa), length(xa) == 1,
# is.numeric(xb), length(xb) == 1)
#
# fun <- match.fun(f)
# f <- function(x) fun(x, ...)
# g <- function(x) (1/x^2) * f(1/x)
#
# if (is.null(method)) {
# integ <- function(f, xa, xb)
# integrate(f, xa, xb, subdivisions = 512, rel.tol = tol)$value
#
# } else {
# methods <- c("Kronrod","Richardson","Clenshaw","Simpson","Romberg")
# method <- match.arg(method, methods)
# integ <- switch(method,
# "Kronrod" = function(fct, xmin, xmax) quadgk(fct, xmin, xmax, tol = tol),
# "Richardson" = function(fct, xmin, xmax) quadgr(fct, xmin, xmax, tol = tol)$value,
# "Clenshaw" = function(fct, xmin, xmax) quadcc(fct, xmin, xmax, tol = tol),
# "Romberg" = function(fct, xmin, xmax) romberg(fct, xmin, xmax, tol = tol)$value,
# "Simpson" = function(fct, xmin, xmax) simpadpt(fct, xmin, xmax, tol = tol)
# )
# }
#
# if (xa == xb) {
# Q <- 0
#
# } else if (xa > xb) {
# Q <- -1 * quadinf(f, xb, xa, tol = tol, method = method)
#
# } else if (is.finite(xa) && is.finite(xb)) {
# Q <- integ(f, xa, xb)
#
# } else if (xa == -Inf && xb == Inf) {
# Q <- integ(g, -1, 0) + integ(f, -1, 1) + integ(g, 0, 1)
#
# } else if (is.finite(xa) && xb == Inf) {
# if (xa > 0)
# Q <- integ(g, 0, 1/xa)
# else
# Q <- integ(f, xa, 1) + integ(g, 0, 1)
#
# } else if (xa == -Inf && is.finite(xb)) {
# if (xb < 0)
# Q <- integ(g, 1/xb, 0)
# else
# Q <- integ(g, -1, 0) + integ(f, -1, xb)
# }
#
# return(Q)
# }
