2017-09-07 73 views
-2
import matplotlib.pyplot as plt 
import numpy as np 
import math 
from scipy import * 
from scipy.integrate import quad, dblquad, tplquad 
q=range(1,6) 
L=range(1,6) 
sigmak=range(1,6) 
x_lower = -3000 
x_upper = 3000 
y_lower = -3000 
y_upper = 3000 #Integrate range 
def final(a,b): #final(a,b)=0 to be plotted on a-b plane 
    m=a 
    n=b 
    def f3(x,y): 
     mass=0 
     for i in range(len(q)): 
      mass+=(L[i]*exp(-(x*x+y*y/(q[i]*q[i]))/(2*sigmak[i]*sigmak[i])))/(2*3.1415926*q[i]*sigmak[i]*sigmak[i]) 
     return mass*(m-x)/((x-m)**2+(y-n)**2) 
    val=dblquad(f3,x_lower, x_upper, lambda x : y_lower, lambda x: y_upper) 
    return val[0] 
y,x=np.ogrid[-1000:1000:200j,-1000:1000:200j]# plot range 
f=final(x,y) 

plt.figure(figsize=(9,4)) 
plt.subplot(121) 
extent=[np.min(x),np.max(x),np.min(y),np.max(y)] 
cs=plt.contour(f,extent=extent,levels=[0,0.1],colors=["b","r"],linestyles=["solid","dashed"],linewidths=[2,2]) 
plt.show() 

以上是我的代码。我想在一架飞机上绘制final(x,y)= 0。最终(X,Y)是一个小complicated.When我跑我的代码,它会引发提供的函数不返回有效的浮点数

Traceback (most recent call last): 
    File "test.py", line 25, in <module> 
    f=final(x,y) 
    File "test.py", line 22, in final 
    val=dblquad(f3,x_lower, x_upper, lambda x : y_lower, lambda x: y_upper) 
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 433, in dblquad 
    return quad(_infunc,a,b,(func,gfun,hfun,args),epsabs=epsabs,epsrel=epsrel) 
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 252, in quad 
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points) 
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 317, in _quad 
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) 
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 381, in _infunc 
    return quad(func,a,b,args=myargs)[0] 
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 252, in quad 
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points) 
    File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 317, in _quad 
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) 
quadpack.error: Supplied function does not return a valid float. 

功能那么,什么是我的问题吗?谢谢,如果有人能帮助我!

回答

0

不知道你正试图在这里做

dblquad回报2个漂浮,所得到的积分和误差的估计。 通过使一个for循环像下面我们可以导出一个数组,但我不认为这就是你需要

def final(a,b): #final(a,b)=0 to be plotted on a-b plane 
    outcome = [] 
    for i in range(len(a)): 
     #print(i) 
     print(a[i]) 
     m=a[i] 
     n=b[i] 
     def f3(x,y): 
      mass=0 
      for i in range(len(q)): 
       mass+=(L[i]*exp(-(x*x+y*y/(q[i]*q[i]))/(2*sigmak[i]*sigmak[i])))/(2*3.1415926*q[i]*sigmak[i]*sigmak[i]) 
      return mass*(m-x)/((x-m)**2+(y-n)**2)   

     val=dblquad(f3,a[0], a[-1], lambda x : m, lambda x: a[-1]) 

     outcome.append(val) 
    return np.array(outcome) 

y=np.ogrid[-10:10:10j] 
x=np.ogrid[-10:10:10j] 
f=final(x,y) 

是你想达到什么样的更精确。