Gabriel Vasconcelos
Gabriel Vasconcelos

Reputation: 606

Trouble with SymPy diff() function

I'm working on implementation of various numerical methods on python, including backward Euler. Thus, I decided to implement a quick Newton-Raphson method that would give me the desired root for variable equations my calculator will have to deal with. Here it is:

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

The code excerpt that involves the call of Newton() function:

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

Where f is a 'sympified' string. In my failing test case it was: f = sympify('(y**2 + t)/(y - t)')

As you can see, I'm printing aux content to trace exactly where the diff function fails. It was in the first iteration for y(0) = 1 with h = 0.1, where aux had: -y + 1 + 0.1*(y**2 + 0.1)/(y - 0.1). When passing aux as f in Newton(), and then running diff(f, y), it gives me:

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

Surprisingly, when I manually 'simpify' that same equation, and 'diff' it, it works. Newton's method also works as a plume. Still, I'm clueless on what can be wrong.

Upvotes: 4

Views: 603

Answers (1)

smichr
smichr

Reputation: 19029

It looks like you are using the mpmath diff instead of the SymPy one. Try "from sympy import diff" before you call your routine. You can check which one you are using by asking for help:

>>> 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:

When I run your program with SymPy's diff I get

>>> 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]

Upvotes: 4

Related Questions