我正在致力于在 python 上实现各种数值方法,包括向后 Euler。因此,我决定实现一种快速的 Newton-Raphson 方法,该方法将为我的计算器必须处理的变量方程提供所需的根。这里是:
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
是一个“符号化”字符串。在我失败的测试用例中,它是: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)
.路过时
aux
如 f
在 Newton()
,然后运行 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
令人惊讶的是,当我手动“简化”相同的等式并“区分”它时,它起作用了。牛顿法也可用作羽流。尽管如此,我对什么可能是错的一无所知。
最佳答案
看起来您正在使用 mpmath 差异而不是 SymPy 差异。在调用例程之前尝试“from sympy import diff”。您可以通过寻求帮助来检查您使用的是哪一种:
>>> help(diff)
Help on method diff in module sympy.mpmath.calculus.differentiation:
对比
>>> 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]
关于python - SymPy diff() 函数的问题,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/26790774/