【问题标题】:Numerical Triple Integration in RR中的数值三重积分
【发布时间】:2017-11-17 00:31:09
【问题描述】:

是否可以在不使用 cubature 包的情况下在 R 中进行三重集成?

based on the answer in this post

InnerFunc = function(x) { x + 0.805 }
InnerIntegral = function(y) { sapply(y, 
    function(z) { integrate(InnerFunc, 15, z)$value }) }
integrate(InnerIntegral , 15, 50)
16826.4 with absolute error < 1.9e-10

例如编码this triple integral:

我试过了

InnerMostFunc = function(v) { v + y^2 }
InnerMostIntegral = function(w) { sapply(w, 
   function(x) { integrate(InnerMostFunc, 1, 2)$value }) }
InnerIntegral = function(y) { sapply(y, 
   function(z){integrate(InnerMostIntegral, 1, 2)$value }) }
integrate(InnerIntegral, 0, 1)

【问题讨论】:

  • 我无法读取您要计算的三重积分。能提供公式吗?也许使用您引用的帖子中的公式图像。
  • @G5W 我试着把它清理得更好!
  • 三重积分需要三个变量(和三个单独的范围)。您的集成只有一个变量,尽管我猜您打算在方形立方体上进行集成。
  • 所以你不能指望 R 提供那个符号结果,但是如果你为 k 指定一个值,你可以得到一个结果。
  • @42 - 是的,我试图在一个方形立方体上进行积分。不过,我一直在寻找具有重复范围的示例(例如我添加的最新示例)。你真的不能用重复范围做三重积分吗?

标签: r integral numerical-integration


【解决方案1】:

再往下是在 cmets 中要求引用的 previous post 的扩展名。但是这个顶部将展示如何计算更新帖子中给出的积分。

这比下面的积分类型更难写,因为函数必须完全嵌套。您不能像第二个示例那样将它们分开,因此单个表达式相当长且难以阅读。这是计算请求积分的代码。

integrate(Vectorize(function(x) { 
    integrate(Vectorize(function(y) { 
        integrate(function(z) { x^2 + y*z }, 1, 2)$value }), 1,2)$value }), 0,1)
2.583333 with absolute error < 2.9e-14

请注意,它会在 31/12 计算正确答案。引用的积分来源错误地给出了 31/4 的答案。

引用的上一篇文章的扩展

这是一个将上一篇文章中的过程扩展到三重积分的示例。我将计算的示例很容易进行分析,因此我们可以检查我们是否得到了正确的答案。我会用:

为了更容易理解,我将代码分成几个步骤。

InnerFunc = function(x) { x + 1 }
InnerIntegral = Vectorize(function(y) { integrate(InnerFunc, 0, y)$value})
Integral2     = Vectorize(function(z) { integrate(InnerIntegral , 0, z)$value})
(Tripleintegral = integrate(Integral2 , 0, 4))
21.33333 with absolute error < 2.4e-13

这扩展了前面的例子,但是没有在问题的新陈述中涵盖了积分的类型。

【讨论】:

  • 太棒了!并且没有变量的变化。我认为这不可能以如此紧凑的形式出现。但它确实需要明确的范围。
  • 让我“困惑”。如果最内层积分的上限是 5,你就不需要改变任何东西?
【解决方案2】:

对于这样的积分,cubature 包是要走的路。

> library(cubature)
> hcubature(function(v) v[1]^2+v[2]*v[3], c(0,1,1), c(1,2,2))
$integral
[1] 2.583333

【讨论】:

    猜你喜欢
    • 2013-12-12
    • 2018-01-17
    • 2018-10-29
    • 1970-01-01
    • 2017-04-06
    • 2017-08-28
    • 2021-11-09
    • 1970-01-01
    相关资源
    最近更新 更多