标签: numerical-integration

使用Python进行三重积分的数值计算

我非常需要计算这个积分。几个月来,我一直在尝试这样做,使用 Python 中的 numpy 包,特别是集成.tplquad 函数。

from __future__ import division
from math import *
import numpy as np
import scipy.special as special
import scipy.integrate as integrate

a=1.e-19
b=1.e-09
zo=1.e7
H=1.e15

v=1.e18

def integrand(v,z,x,u):
    value=x**(-0.5)*special.kv(5/3,u)*(a*v*z/x-1/2.)*exp(-b*sqrt(z*v/x))
    return value

i=integrate.tplquad(lambda u,x,z: integrand(v,z,x,u),1.e7,1.e15,lambda z:0.,lambda z:np.inf, lambda x,z : x, lambda x,z : np.inf)

print i
Run Code Online (Sandbox Code Playgroud)

在上面的代码中,我尝试了 v=10^18 的值,以便对指数的参数进行归一化,并且不会得到太小或太大的系数。

但是,无论我插入的 v 值是多少,我总是得到

out: (0.0, 0.0)
Run Code Online (Sandbox Code Playgroud)

我不知道如何超越这个问题。

我也尝试将指数函数扩展为幂级数,但得到相同的结果。

现在,我知道积分对于所有 v 必须有一个有限的正值。 事实上,如果我可以为任何 v 计算它,我会很高兴。

如果有人遇到过类似的问题,如果他们能分享他们的智慧,我会很高兴。欢迎任何帮助。谢谢。

python numpy integral numerical-methods numerical-integration

1
推荐指数
1
解决办法
1351
查看次数

积分可能发散,或缓慢收敛

from scipy import integrate
import numpy as np
from mpmath import coth

DvDc = 6.5
dens = 5.65
Vs = 4.6e3
a = 3.3e-9
w0 = np.sqrt(2)*Vs/a
T = 50 
kb = 1.38064852e-5  # in eV`
j0 = (DvDc)**2 / ((2*np.pi)**2 *dens*Vs**5)

def func(x):
    return x*np.exp(-(x/w0)**2)*coth(x/(2*kb*T))
S = j0*integrate.quad(func, 0, np.inf)[0]
print(S)
Run Code Online (Sandbox Code Playgroud)

嗨,所以我没有使用数值积分的经验,但我想知道是否有人可以提供帮助。在这里,我定义了一些变量,并有一个要集成的函数。但是我得到了“积分可能是发散的,或者慢慢收敛的。” 尽管我在 coth 的分母中添加了内容,但几乎得到了相同的结果。

有谁知道我是否犯了一个错误,或者我只是不能做到这一点。

谢谢!

python numerical-integration

1
推荐指数
1
解决办法
2007
查看次数

复合辛普森的规则

我有复合Simpson规则的代码.但是,我一直在摆弄它已经有一段时间了,我似乎无法让它发挥作用.

我该如何修复这个算法?

function out = Sc2(func,a,b,N)
% Sc(func,a,b,N)
% This function calculates the integral of func on the interval [a,b]
% using the Composite Simpson's rule with N subintervals.
x=linspace(a,b,N+1);
% Partition [a,b] into N subintervals
fx=func(x);
h=(b-a)/(2*N);
%define for odd and even sums
sum_even = 0;
for i = 1:N-1
   x(i) = a + (2*i-2)*h;
   sum_even = sum_even + func(x(i));
end

sum_odd = 0;

for i = 1:N+1
   x(i) = a + (2*i-1)*h;
   sum_odd = sum_odd + func(x(i));
end …
Run Code Online (Sandbox Code Playgroud)

matlab numerical-analysis numerical-methods numerical-integration

0
推荐指数
1
解决办法
1947
查看次数

numpy 中的数字积分

我想做一些非常简单的事情,但我无法在numpy. 我想对一个由其值(而不是其公式!)给出的函数进行数值和连续积分。这意味着我只想要一个包含输入数组开头总和的数组。例子:

输入:

[ 4, 3, 5, 8 ]
Run Code Online (Sandbox Code Playgroud)

输出:

[ 4, 7, 12, 20 ]  # [ sum(i[0:1]), sum(i[0:2]), sum(i[0:3]), sum(i[0:4]) ]
Run Code Online (Sandbox Code Playgroud)

听起来很简单,所以我希望这一定很容易,numpy因为我目前无法找到一些功能。

我发现了类似scipy.integrate.quad()但似乎在给定范围(从 a 到 b)内积分并返回单个值的东西。我需要一个数组作为输出。

python arrays numpy numerical-integration

0
推荐指数
3
解决办法
1万
查看次数

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

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 …
Run Code Online (Sandbox Code Playgroud)

python numerical-integration

0
推荐指数
1
解决办法
265
查看次数