use sympy to find gradient and plot vector field

Question

I wrote some code to use sympy to find the gradient of a function f(x,y) = x*y**2, and then to plot the vector field from the gradient. See Below:

%matplotlib inline
import matplotlib.pyplot as plt
import sympy as sp
import numpy as np

sp.init_printing()
x,y = sp.symbols('x y')

def gradient(f):
    return (f.diff(x), f.diff(y))

f = x*y**2
g = gradient(f)
g

X,Y = np.meshgrid(np.linspace(-3,3,15), np.linspace(-3,3,15))

U=[g[0].subs({x:x1, y:y1}) for (x1,y1) in zip(X,Y)]
V=[g[1].subs({x:x1, y:y1}) for (x1,y1) in zip(X,Y)]

plt.quiver(X,Y,U,V, linewidth=1)
plt.title("vector field")
plt.show()

What i'm wondering about is why the sympy "subs" function is not working in this code. Its just returns the expression without inserting a value of X and Y to evaluate to a numerical value, instead its just returning the sympy object without any subsitution.

https://stackoverflow.com/questions/36013063/what-is-the-purpose-of-meshgrid-in-python-numpy — DiscreteMath, Mar 25 '19 at 01:17

score 1 · Answer 1 · answered Aug 16 '22 at 13:51

As an alternative to subs one could use lambdify from sympy.utilities, as in the following reference: https://docs.sympy.org/latest/modules/utilities/lambdify.html

%matplotlib inline
import matplotlib.pyplot as plt
import sympy as sp
import numpy as np
from sympy.utilities.lambdify import lambdify

sp.init_printing()
x,y = sp.symbols('x y')

def gradient(f):
    return (f.diff(x), f.diff(y))

f = x*y**2
g = gradient(f)


X,Y = np.meshgrid(np.linspace(-3,3,15), np.linspace(-3,3,15))

f1 = lambdify([x, y], g[0])
f2 = lambdify([x, y], g[1])

U=[f1(x1, y1) for x1,y1 in zip(X,Y)]
V=[f2(x1,y1) for x1,y1 in zip(X,Y)]


plt.quiver(X,Y,U,V, linewidth=1)
plt.title("vector field")
plt.show()

score 0 · Answer 2 · answered Mar 25 '19 at 01:32

The problem with your code is that you need to access a meshgrid as a 2-dimensional array.

Example: U[i,j] not U[i]

%matplotlib inline
import matplotlib.pyplot as plt
import sympy as sp
import numpy as np

sp.init_printing()
x,y = sp.symbols('x y')

def gradient(f):
    return (f.diff(x), f.diff(y))

f = x*y**2
g = gradient(f)
g

xrange = np.linspace(-3,3,15)
yrange = np.linspace(-3,3,15)
X,Y = np.meshgrid(xrange, yrange)

U=X
V=Y

for i in range(len(xrange)):
    for j in range(len(yrange)):
        x1 = X[i,j]
        y1 = Y[i,j]
        U[i,j] = g[0].subs({x:x1, y:y1})
        V[i,j] = g[1].subs({x:x1, y:y1})

plt.quiver(X,Y,U,V, linewidth=1)
plt.title("vector field")
plt.show()

use sympy to find gradient and plot vector field

2 Answers2