在 R 中建模无限系列

Question

我正在尝试编写代码以根据 R 中的 Theis 水文地质方程逼近以下无限泰勒级数。

我是函数式编程的新手，所以这是一个挑战！这是我的尝试：

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)
}

我将结果与此处提供的数据 table 的值进行了比较：https://pubs.usgs.gov/wsp/wsp1536-E/pdf/wsp_1536-E_b.pdf - 见下文（请原谅冗长的代码 - 事后看来应该制作一个 csv）：

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)

好消息是公式给出了重复次数 = 50、100、150 和 170 的相同结果（所以我刚刚给你上面的 100 版本）。坏消息是，虽然该函数在 u < ~10^-3 时表现良好，但它超出了 rails 并在 1 左右的数量级内给出负输出。这不会发生当我只是在一个单独的号码上调用该函数时。即：

> Wu(4.6)
[1] 0.001856671

2sf 的正确答案是什么。

谁能指出我做错了什么 and/or 提出一个更好的方法来编写这个等式？我认为问题与我的 for 循环有关 and/or 随着 u 变大，阶乘产生无限数的问题，但我完全不确定。

谢谢！

Answer 1

您可以编写一个函数，将一个值与重复次数一起输入并输出所需的值：

w=function(u,l){
  a=2:l
  -0.5772-log(u)+u+sum(u^(a)*rep(c(-1,1),length=l-1)/(a)/factorial(a))
}
transform(Wu_QC,new=Vectorize(w)(u,170))
         u    Wu_table          new
1  1.0e-15 3.39616e+01 3.396158e+01
2  4.1e-14 3.02480e+01 3.024800e+01
3  9.9e-13 2.70639e+01 2.706387e+01
4  7.0e-12 2.51079e+01 2.510791e+01
5  3.7e-11 2.34429e+01 2.344290e+01
6  2.3e-10 2.16157e+01 2.161574e+01
7  6.8e-09 1.82291e+01 1.822914e+01
8  5.7e-08 1.61030e+01 1.610301e+01
9  8.4e-07 1.34126e+01 1.341266e+01
10 6.3e-06 1.13978e+01 1.139777e+01
11 3.1e-05 9.80430e+00 9.804354e+00
12 7.4e-04 6.63240e+00 6.632400e+00
13 5.1e-03 4.70640e+00 4.706408e+00
14 2.9e-02 2.99200e+00 2.992051e+00
15 8.7e-01 2.74200e-01 2.741930e-01
16 4.6e+00 1.84100e-03 1.856671e-03
17 9.9e+00 4.63700e-06 2.030179e-05

随着数字越来越大，估计不太好，所以我们应该比170更进一步！但 R 不能那样做。也许您可以尝试其他平台。即Python

Answer 2

我想我自己可能已经解决了这个问题（尽管大量借鉴了 Onyambo 的回答！）这是我的代码：

well_func2 <- function (u, l = 100) {
    result <- numeric(length(u))
a <- 2:l
for(i in seq_along(u)){
result[i] <- -0.5772-log(u[i])+u[i]+sum(u[i]^(a)*rep(c(-1,1),length=l-1)/(a)/factorial(a))
}
return(result)
}

据我所知，这与 u <5 的表格结果很好地匹配（就像 Onyambo 的代码一样），并且它还给出了矢量与 single-value 输入相同的结果。

仍然需要更多的测试，并且可能有一种使用 map() 或类似方法而不是 for 循环来编码它的更简洁的方法，但我现在已经很高兴了。以为我会分享以防其他人遇到同样的问题。

Answer 3

正如您参考资料第 93 页所述，W 也称为指数积分。另见 here.

然后，例如，package expint 提供了一个计算 W(u) 的函数：

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

其中的结果与您提到的 table 完全相同。

在 R 中建模无限系列

modelling an infinite series in R

r

for-loop

taylor-series