How to go from SymPy to NumPy

11.1 Success with SymPy

Consider a function:

$f(x)=\sin(x)$

and its derivative:

$f'(x)=\frac{d}{dx}\sin(x)=\cos(x)$

Lets find an extremum for $f(x)$ by solving $f'(x_0)=0$ . Lets also find the crossing between the two curves $f(x)$ and $f'(x)$ by solving $f(x_1)=f'(x_1)$ .

Figure 1

You know how to do that in SymPy:

import sympy as sp

x = sp.symbols('x')
f = sp.sin(x)
fprime = sp.diff(f,x)

x0s = sp.solve(fprime,x)
print(x0s)
x0 = x0s[0]
print('x_0:',x0,'\n')

x1s = sp.solve(f-fprime,x)
print(x1s)
x1 = x1s[1]
print('x_1:',x1)

[pi/2, 3*pi/2]
x_0: pi/2

[-3*pi/4, pi/4]
x_1: pi/4

11.2 Failure with Sympy

Lets now make the problem a little harder:

$f(x)=x^2\sin(x)$

and its derivative:

$f'(x)=\frac{d}{dx}\sin(x)=x^2\cos(x)+2x\sin(x)$

Figure 2

We do as before:

f = sp.sin(x)*x**2
fprime = sp.diff(f,x)
x0s = sp.solve(fprime,x)

But get a long error message ending with:

No algorithms are implemented to solve equation x**2*(1 - tan(x/2)**2)/(tan(x/2)**2 + 1) + 4*x*tan(x/2)/(tan(x/2)**2 + 1)

The same applies when issuing:

x1s = sp.solve(f-fprime,x)

No algorithms are implemented to solve equation -x**2*(1 - tan(x/2)**2)/(tan(x/2)**2 + 1) + 2*x**2*tan(x/2)/(tan(x/2)**2 + 1) - 4*x*tan(x/2)/(tan(x/2)**2 + 1)

We are thus stuck with SymPy and cannot solve this only slightly more difficult problem.

11.3 NumPy to the rescue

The above problems are easily solved numerically in Python using NumPy and SciPy. Here is the hard problem:

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

x = sp.symbols('x')
f = sp.sin(x)*x**2
fprime = sp.diff(f,x)

# convert SymPy functions to NumPy functions
numf = sp.lambdify(x, f, 'numpy')
numfprime = sp.lambdify(x, fprime, 'numpy')

# plot
xs = np.linspace(0,np.pi,30)
plt.rc('font', size=16)
plt.plot(xs,numf(xs))
plt.plot(xs,numfprime(xs))
plt.grid()

x0 = fsolve(numfprime,2)
plt.plot(x0,numf(x0),'o',)

x0 = fsolve(lambda x: numf(x) - numfprime(x),1)
plt.plot(x0,numf(x0),'o',)

plt.legend(['$f(x)$',"$f'(x)$"])
plt.xlabel('$x$')
plt.tight_layout()
plt.savefig('sympy_numpy_2.png')