我的积分一直在使用，或者在运行2小时或更长时间后给出了错误的答案

html5 • • 问答 

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from scipy import exp
import scipy as sp
import mpmath as mp
from sympy import polylog
from math import *

def F(s,y):
  return (16/pi) * (sqrt((1-s)/(s+y))) * (log(s+y/2+sqrt(s*(s+y)))-log(y/2))
def fd1(w):
  gw=np.exp((w-mu)/TTF)
  return 1/(gw+1)
def fd123(w,y):
  gwy=np.exp((w*(1+2*y)-mu)/TTF)
  return 1/(gwy+1)
def fd12(w,y,z):
  gwyz=np.exp((w*(1+y-z)-mu)/TTF)
  return 1/(gwyz+1)
def fd13(w,y,z):
  gwyz=np.exp((w*(1+y+z)-mu)/TTF)
  return 1/(gwyz+1)
def bounds_y():
  return [0,10]
def bounds_s(y):
  return [0,1]
def bounds_w(y,z):
  return [0,10]
def bounds_z(w,y,z):
  return [-y,y]
def integrand1(z,w,s,y):
  return 144 * w**(3) *  fd13(w,y,z) * fd12(w,y,z) * (1-fd1(w)) * (1-fd123(w,y)) * F(s,y)
mu=0.22
TTF=0.5;
integrate.nquad(integrand1,[bounds_z,bounds_w,bounds_s,bounds_y])

基本上我正在尝试计算取决于四个变量的四积分，但它一直花费太多时间并最终给出错误的答案。它还给了IntegrationWarning：积分可能是发散的，或者慢慢收敛的。我不知道我在哪里做错了什么。

回答

一大早，锻炼身体！

积分概述

我将使用一些屏幕截图，因为 stackoverflow 不支持方程。

这可以重新排列为

寻找被积函数，有一些希望z可以解析地解决积分。积分中最困难的部分是当s和y都接近于零时。表达式最右侧的对数项和左侧平方根中的分母都将接近于零。

数值复杂性实验

既然你提到它需要永远，我们在集成时发现了一些困难s，y第一个测试是对F(z, y)

integrate.nquad(F, [bounds_s, bounds_y])

 Wall time: 128 ms
 (15.340467397261975, 6.86478152545078e-07)

 Wall time: 128 ms
 (15.340467397261975, 6.86478152545078e-07)

它很容易积分，所以让我们来G(w, y)解决三重积分 

   Wall time: 1min 5s
  (5107.5853057484455, 0.00022839419762021862)

   Wall time: 1min 5s
  (5107.5853057484455, 0.00022839419762021862)

它更慢，精确到小数点后 7 位，并在 1 分钟内运行。因此，让我们尝试通过z分析来解决积分问题。

z 中的积分

我将利用同情这件作品

integrate.nquad(
   lambda w, s, y: w**2 * (1-fd1(w)) * (1-fd123(w,y)) * F(s,y) , 
   [bounds_w, bounds_s, bounds_y])

那给

在我们的例子b中总是非零。

最终计算

上面的表达式给了我们不定积分，现在计算这个表达式替换z=y和z=-y计算定积分的差。这给了我们更多的功能。

最后我们可以评估初始积分为

import sympy
from sympy.abc import a,b,z
Iz = sympy.integrate(1 / (1 + sympy.exp(2*a) + sympy.cosh(b*z)), z)
Iz

数值考虑

上面的集成将有问题，哪里w大，因为 ifa和b*ywill 很大，因此candtz将非常接近于 1。然后你可以应用近似值

考虑大阳性的情况a，b*y你有

因此，该术语log(c - tz)可以很好地近似为log(exp(-2*a) + 2*exp(-b*y))。您可以为符号 fora和的四种组合编写类似的表达式b*y。但是由于积分发生在非负w和y

更新的功能可以是

def Iwyz(w, y):
    a = (w*y + w - mu)/TTF
    b = 1/TTF;
    tz = np.tanh(b*y/2);
    ea = np.exp(-2*a)
    c = np.sqrt(1 + 2*ea);
    # Logarithm may be problematic non positive numbers are given
    # write the sum of logarithms as logarithm of product
    # since we know that the integral will be positive, the value
    # argument must be postive as well
    #
    # also takes advantage of tanh(-b*y/2) = -tanh(b*y/2) and reuse tz
    return (ea / c) * (np.log((-c-tz)*(c+tz)/((-c+tz)*(c-tz))));

您可能希望选择不同的条件来切​​换到渐近近似，您也可以选择不同的近似。

使用渐近近似的积分结果如果 eby + ea < 1e-6

3min 41s
(0.038781260929779036, 1.4898372223919173e-08)

为了评估近似值的影响，我将根据不同的条件重新计算应用近似值。

使用渐近近似的结果如果 eby + ea < 1e-9

Wall time: 3min 44s
(0.038781260929776414, 1.4898372223919173e-08)

结果匹配小数点后 13 位，这意味着积分误差高于渐近近似引入的误差。

以上是我的积分一直在使用，或者在运行2小时或更长时间后给出了错误的答案的全部内容。
THE END
分享
二维码
< <上一篇
下一篇>>