modelling an infinite series in R

Question

I'm trying to write a code to approximate the following infinite Taylor series from the Theis hydrogeological equation in R.

I'm pretty new to functional programming, so this was a challenge! This is my attempt:

Wu <- function(u, repeats = 100) {
result <- numeric(repeats) 
for (i in seq_along(result)){
result[i] <- -((-u)^i)/(i * factorial(i))
}
return(sum(result) - log(u)-0.5772)
}

I've compared the results with values from a data table available here: https://pubs.usgs.gov/wsp/wsp1536-E/pdf/wsp_1536-E_b.pdf - see below (excuse verbose code - should have made a csv, with hindsight):

Wu_QC <- data.frame(u = c(1.0*10^-15, 4.1*10^-14,9.9*10^-13, 7.0*10^-12, 3.7*10^-11, 
2.3*10^-10, 6.8*10^-9, 5.7*10^-8, 8.4*10^-7, 6.3*10^-6,
3.1*10^-5, 7.4*10^-4, 5.1*10^-3, 2.9*10^-2,8.7*10^-1,
4.6,9.90),
Wu_table = c(33.9616, 30.2480, 27.0639, 25.1079, 23.4429, 
21.6157, 18.2291, 16.1030, 13.4126, 11.3978,
9.8043,6.6324, 4.7064,2.9920,0.2742,
0.001841,0.000004637))

Wu_QC$rep_100 <-  Wu(Wu_QC$u,100)

The good news is the formula gives identical results for repeats = 50, 100, 150 and 170 (so I've just given you the 100 version above). The bad news is that, while the function performs well for u < ~10^-3, it goes off the rails and gives negative outputs for numbers within an order of magnitude or so of 1. This doesn't happen when I just call the function on an individual number. i.e:

> Wu(4.6)
[1] 0.001856671

Which is the correct answer to 2sf.

Can anyone spot what I've done wrong and/or suggest a better way to code this equation? I think the problem is something to do with my for loop and/or an issue with the factorials generating infinite numbers as u gets larger, but I'm not at all certain.

Thanks!

Answer 1

As it says on page 93 of your reference, W is also known as the exponential integral. See also here.

Then, e.g., the package expint provides a function to compute W(u):

library(expint)
expint(10^(-8))
# [1] 17.84347
expint(4.6)
# [1] 0.001841006

where the results are exactly as in your referred table.