2013-05-17 42 views
2

我對python非常陌生,作爲一個項目,我決定在Python中編寫我的Mathematica項目,以瞭解它是如何工作的,因此代碼的編寫風格與Mathematica儘可能接近。從插值函數中檢索值

我奮力呼籲從插值函數的值,在一個簡單的工作的例子,我想這樣做:

import numpy as np 
from scipy.interpolate import interp1d 

a = np.linspace(1,10,10) 
b = np.sin(a) 
inter = interp1d(a,b) 

# this is where i get a value from the interpolated function 
s0 = inter(a[0]) 
print(s0) 

這是我MWE不起作用:

import matplotlib 
import numpy as np 
import matplotlib.pyplot as plt 
from sympy import * 
from scipy import integrate 
from scipy.interpolate import interp1d 

E0 = -0.015 
L = 5.5 
Ns = 1000 


# this solves for where the energy E0 intersects the potential y 
def req(E0): 
    L=5.5 
    r = Symbol('r') 
    y = -(2*L**2)/(r**3)+(L**2)/(r**2)-(2)/(r) 
    rr = (E0-y)*(r**4) 
    rreq = Eq(rr, 0) 
    rrt = sorted(solve(rr), key=int) 
    return rrt 

# upper and lower limits on r 
r1 = req(E0)[1] 
r0 = req(E0)[2] 

# initialise the arrays 
a = np.array([1]) 
b = np.array([1]) 

# numerically integrate the function R(r) 
for n in range(2, Ns): 
    # integrate 
    lwlmt = r0 
    uplmt = r0+(n-1)*(r1-r0)/(Ns-1) 

    result, error = integrate.quad(lambda ra: -1/((E0-(-(2*L**2)/(ra**3)+(L**2)/(ra**2)-(2)/(ra)))*(ra**4))**(0.5), r0, uplmt) 

    a = np.append([uplmt],[[a]]) 
    b = np.append([result],[[b]]) 

# chop the 1 from the end 
aa = a[:-1] 
ba = b[:-1] 

# interpolate 
inter = interp1d(aa,ba) 

# this is the problem 
print(inter(110)) 

# this is what i would ideally like to do, 
# get the start and end points however i receive an error 
s0 = inter(aa[0]) 
s1 = inter(aa[len(aa)-1]) 


plt.plot(aa,inter(aa)) 
plt.show() 

奇怪的是,我的MWE只適用於整個數組作爲參數inter(aa)它返回一個插值點列表。 我無法弄清楚爲什麼第一個例子有效,而第二個例子不行。兩個數組看起來都是相同的,但只有第一個例子實際上產生了一個輸出。

編輯:添加返回錯誤

--------------------------------------------------------------------------- 
AttributeError       Traceback (most recent call last) 
/home/nick/Documents/python/<ipython-input-4-b36b3f397c2e> in <module>() 
    45 # this is what i would ideally like to do, 

    46 # get the start and end points 

---> 47 s0 = inter(aa[1]) 
    48 s1 = inter(aa[len(aa)-1]) 
    49 

/usr/lib/python2.7/dist-packages/scipy/interpolate/interpolate.pyc in __call__(self, x_new) 
    364   # from y_new and insert them where self.axis was in the list of axes. 

    365   nx = x_new.ndim 
--> 366   ny = y_new.ndim 
    367 
    368   # 6. Fill any values that were out of bounds with fill_value. 


AttributeError: 'Float' object has no attribute 'ndim' 

這對任何號碼,我把在參數val inter([val])的錯誤,那就是在aa的範圍內。

+1

你能更專注於'不工作':產生錯誤的結果,給出一個錯誤(什麼是追溯)等等? –

回答

1

req()的結果是sympy對象而不是真正的Python浮點數,因此,uplmt也是sympy對象。 numpy的數值例程不知道如何處理這些sympy對象。儘早轉換爲Python浮動對象。

在我的機器上,r0r1的值實際上很複雜,只是帶有微小的虛構分量,它們會比您顯示的錯誤更早造成錯誤。這很容易將它們轉換,雖然:

# upper and lower limits on r 
r1 = complex(req(E0)[1]).real 
r0 = complex(req(E0)[2]).real 

我做出的改變後,你的腳本執行完成對我來說,雖然我不能保證它是給你你想要的數值結果。

+0

謝謝你的回答,這已經完成了我想要的/ –