2014-10-03 221 views
1

我对python相当陌生,所以我希望我的措辞有意义。我目前正试图模拟一组要求积分乘积乘以浮点数的方程。我独自从获得的积分输出楠输出,我不知道为什么这是我的代码如下:python中的scipy quad集成问题

from __future__ import division 
import matplotlib.pylab as plt 
import numpy as np 
import scipy.special as sp 
import scipy as sci 
from scipy import integrate 
import math 
from sympy.functions import coth 
Tc = 9.25 
Tb = 7.2 
t = Tb/Tc 
Temperature = [] 
Temp1=[] 
Temp0=[] 
D=[] 
d = [] 
D1=[] 
d1 = [] 
n = 2*10**-6 
L = 100*10**-9 
W = 80*10**-9 
a = 3*10**-2 
s1 = W/ (2*n) 
y1 = (L+(W/2))/(2*n) 
x0 = 0.015 
r0 = 2*x0 
s2 = r0/n 
y0 = (x0/n)/1000000 
print x0, y0, y1 
A = ((W/n)**2) *(sp.kv(0, s1)+(math.pi/2)*sp.kv(1,s1)*coth(y1)) 
B = ((W/n)**2) *(sp.iv(0, s1)+(math.pi/2)*sp.iv(1,s1)*coth(y1)) 
print A, B 
def t1(t): 
    return (t**-1)*sp.kv(0, s2) 
def t2(t): 
    return (t**-1)*sp.iv(0, s2) 
print t2 
Fk2 =(math.pi**-2) * integrate.quad(t1, s1, s2, full_output=False)[0] 
FI2 =(math.pi**-2) * integrate.quad(t2, s1, s2, full_output=False)[0] 
print Fk2 , #FI2 
r1 = 0.0 
while r1 < y1: 
    #C0 = sp.kv(0,s2)*(1 + (A*FI2)-(B*Fk2))/A 
    #print C0 
    #D_ = 1 - B*Fk2 - A*Fk2*sp.iv(1, s1)/sp.kv(1, 1) 
    #print D_ 
    r1 += 0.0001 
    j = -1*r1 
    D.append(r1) 
    d.append(j) 
    #T = Tb + (Tc - Tb) * (sp.kv(0,s1) + (math.pi /2)* sp.kv(1, s1)*coth(r1))*(1- D_ * math.cosh(y1)) * (C0*A) 
    #Temp0.append(T) 
    #print Temp0, r1 

罪魁祸首似乎是在等式FI2 sp.iv修正贝塞尔函数(T2,S1)返回一个值楠,但其他公式结果FK2给出了0.1的一段时间,我发现了以下错误:

IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated. warnings.warn(msg, IntegrationWarning) 

但已经停止,现在我只得到0.0和楠。任何帮助真的很感激我很迷茫。

回答

1

您的s2值相当大(15000.0)。所以,当你在S2你零评估Bessel Function

>>> sp.kv(0, 15000.0) 
0.0 

所以你的函数T1总是返回零,使得积分为零。