2014-11-06 118 views
3

我正在使用python實現各種數值方法,包括後向Euler。因此,我決定實施一種快速的Newton-Raphson方法,該方法將爲我提供我的計算器必須處理的變量方程所需的根。那就是:SymPy diff()函數的問題

def Newton(y0, f, tol=1e-10): 
    y = Symbol('y') 
    df = diff(f, y) 
    while f.subs(y, y0) > tol: 
     y0 -= f.subs(y, y0)/df.subs(y, y0) 
    return y0 

的代碼摘錄涉及Newton()函數的調用:

def InverseEuler(f, y0, t0, tf, h): 
    coords = [[t0, y0]] 
    t, y = symbols('t y') 

    while round(t0, 7) < round(tf, 7): 
     aux = sympify(y0 - y + h*f.subs(t, t0 + h)) 
     print aux 
     y0 = utils.Newton(y0, aux) 
     t0 += h 
     coords.append([t0, y0]) 
    return coords 

f是 'sympified' 字符串。在我的失敗測試案例中,它是:f = sympify('(y**2 + t)/(y - t)')

正如您所看到的,我正在打印aux內容以準確跟蹤diff函數失敗的位置。 y(0)= 1的第一次迭代中h = 0.1,其中aux爲:-y + 1 + 0.1*(y**2 + 0.1)/(y - 0.1)。 當Newton()傳遞auxf,然後運行diff(f, y),它給了我:

Traceback (most recent call last): 
    File "C:\Users\Gabriel Vasconcelos\Documents\python\Metodos\metodos.py", line 97, in <module> 
    GPHandler.replot(InverseEuler(f, 1, 0, 4, 0.1)) 
    File "C:\Users\Gabriel Vasconcelos\Documents\python\Metodos\metodos.py", line 31, in InverseEuler 
    y0 = utils.Newton(y0, aux) 
    File "C:\Users\Gabriel Vasconcelos\Documents\python\Metodos\utils.py", line 6, in Newton 
    df = diff(f, y) 
    File "C:\Python27\lib\site-packages\sympy\mpmath\calculus\differentiation.py", line 188, in diff 
    values, norm, workprec = hsteps(ctx, f, x, n, prec, **options) 
    File "C:\Python27\lib\site-packages\sympy\mpmath\calculus\differentiation.py", line 61, in hsteps 
    values = [f(x+k*h) for k in steps] 
TypeError: 'Add' object is not callable 

出人意料的是,當我手動「simpify」同一方程,「差異」它,它的工作原理。牛頓的方法也可以用作羽流。儘管如此,我仍然對可能出錯的東西毫無頭緒。

回答

4

它看起來像你正在使用mpmath diff而不是SymPy。在你打電話給你的例程之前,試試「從sympy import diff」。您可以通過尋求幫助檢查您所使用的一個:

>>> help(diff) 
Help on method diff in module sympy.mpmath.calculus.differentiation: 

VS

>>> from sympy import diff 
>>> help(diff) 
Help on function diff in module sympy.core.function: 

當我與SymPy的DIFF運行您的程序,我得到

>>> InverseEuler(f,1, 0, 1, 0.1)[-1] 
-y + 1 + 0.1*(y**2 + 0.1)/(y - 0.1) 
-y + 1.13404207225556 + 0.1*(y**2 + 0.2)/(y - 0.2) 
-y + 1.30637177730098 + 0.1*(y**2 + 0.3)/(y - 0.3) 
-y + 1.52036600337727 + 0.1*(y**2 + 0.4)/(y - 0.4) 
-y + 1.77886566311632 + 0.1*(y**2 + 0.5)/(y - 0.5) 
-y + 2.08466045990054 + 0.1*(y**2 + 0.6)/(y - 0.6) 
-y + 2.44089877798719 + 0.1*(y**2 + 0.7)/(y - 0.7) 
-y + 2.85134771320082 + 0.1*(y**2 + 0.8)/(y - 0.8) 
-y + 3.32053168504081 + 0.1*(y**2 + 0.9)/(y - 0.9) 
-y + 3.85380349564978 + 0.1*(y**2 + 1.0)/(y - 1.0) 
[0.99999999999999989, 4.45738956363267] 
+0

謝謝!我需要習慣這些衝突......自然,切換位置我進口sympy和我導入mpmath的地方使它工作得很好! – 2014-11-06 23:59:42