2017-09-07 6 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を平面上にプロットしたいと思います。 final(x、y)はちょっと複雑な関数です。コードを実行すると、それは発生します提供された関数が有効な浮動小数点を返さない

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台の山車、結果としての積分と誤差の推定値をやろうとしているか確認してください。以下我々は、配列をエクスポートすることができますが、私はそれはあなたが達成しようとしているもので、より正確に

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) 

を必要とするものだとは思わないようにループのために作ることで 。

関連する問題