python4oceanographers

Turning ripples into waves

Sympy's function and lambdify example

This post will be a short and simple one. I could not prepare anything more elaborated for this week because I was out on a much deserved holiday.

It is a quick example on how to use sympy, to evaluate and plot functions.

Let's start!

The cell bellow activates Sympy's printing functionality. Now any function object will be displayed as latex equation that gets automagically rendered by the IPython notebook.

In [2]:
from sympy.interactive import printing
printing.init_printing()

The functions we are using contain two variables (x, y) that are represented by sympy's symbols object. We also need sympy's sqrt. (This sqrt is different from either numpy's and math's sqrt functions.) Last but not least we need the Function class to build the function objects.

Note: The S is to simplify the functions. We do not need it unless initializing the equation with quotes. But, by doing so, sympy won't evaluate any number expression.

In [3]:
from sympy import Function, symbols, sqrt, S

x, y = symbols("x y")

Here we create the 5 functions named by color. (The colors will help to identify them later in the plot.)

In [4]:
red = S("(x / 7) ** 2 * sqrt(abs(abs(x) - 3) / (abs(x) - 3)) + (y / 3) ** 2 * sqrt(abs(y + 3 / 7 * sqrt(33)) / (y + 3 / 7 * sqrt(33))) - 1")
red
Out[4]:
$$\frac{x^{2}}{49} \sqrt{\frac{\operatorname{abs}{\left (\operatorname{abs}{\left (x \right )} - 3 \right )}}{\operatorname{abs}{\left (x \right )} - 3}} + \frac{y^{2}}{9} \sqrt{\frac{\operatorname{abs}{\left (y + \frac{3 \sqrt{33}}{7} \right )}}{y + \frac{3 \sqrt{33}}{7}}} - 1$$
In [5]:
orange = S("abs(x / 2) - ((3 * sqrt(33) - 7) / 112) * x ** 2 - 3 + sqrt(1 - (abs(abs(x) - 2) - 1) ** 2) - y")
orange
Out[5]:
$$- x^{2} \left(- \frac{1}{16} + \frac{3 \sqrt{33}}{112}\right) - y + \sqrt{- \left(\operatorname{abs}{\left (\operatorname{abs}{\left (x \right )} - 2 \right )} - 1\right)^{2} + 1} + \operatorname{abs}{\left (\frac{x}{2} \right )} - 3$$
In [6]:
green = S("9 * sqrt(abs((abs(x) - 1) * (abs(x) - 0.75)) / ((1 - abs(x)) * (abs(x) - 0.75))) - 8 * abs(x) - y")
green
Out[6]:
$$- y + 9 \sqrt{\frac{\operatorname{abs}{\left (\left(\operatorname{abs}{\left (x \right )} - 1\right) \left(\operatorname{abs}{\left (x \right )} - 0.75\right) \right )}}{\left(- \operatorname{abs}{\left (x \right )} + 1\right) \left(\operatorname{abs}{\left (x \right )} - 0.75\right)}} - 8 \operatorname{abs}{\left (x \right )}$$
In [7]:
blue = S("3 * abs(x) + 0.75 * sqrt(abs((abs(x) - 0.75) * (abs(x) - 0.5)) / ((0.75 - abs(x)) * (abs(x) - 0.5))) - y")
blue
Out[7]:
$$- y + 0.75 \sqrt{\frac{\operatorname{abs}{\left (\left(\operatorname{abs}{\left (x \right )} - 0.75\right) \left(\operatorname{abs}{\left (x \right )} - 0.5\right) \right )}}{\left(- \operatorname{abs}{\left (x \right )} + 0.75\right) \left(\operatorname{abs}{\left (x \right )} - 0.5\right)}} + 3 \operatorname{abs}{\left (x \right )}$$
In [8]:
pink = S("2.25 * sqrt(abs((x - 0.5) * (x + 0.5)) / ((0.5 - x) * (0.5 + x))) - y")
pink
Out[8]:
$$- y + 2.25 \sqrt{\frac{\operatorname{abs}{\left (\left(x - 0.5\right) \left(x + 0.5\right) \right )}}{\left(- x + 0.5\right) \left(x + 0.5\right)}}$$
In [9]:
brown = S("6 * sqrt(10) / 7 + (1.5 - 0.5 * abs(x)) * sqrt(abs(abs(x) - 1) / (abs(x) - 1)) - (6 * sqrt(10) / 14) * sqrt(4 - (abs(x) - 1) ** 2) - y")
brown
Out[9]:
$$- y + \sqrt{\frac{\operatorname{abs}{\left (\operatorname{abs}{\left (x \right )} - 1 \right )}}{\operatorname{abs}{\left (x \right )} - 1}} \left(- 0.5 \operatorname{abs}{\left (x \right )} + 1.5\right) - \frac{3 \sqrt{10}}{7} \sqrt{- \left(\operatorname{abs}{\left (x \right )} - 1\right)^{2} + 4} + \frac{6 \sqrt{10}}{7}$$

It is fun to play with sympy's functions, just check the docs out and be amazed! There are tons of possibilities.

In the next cell we use lambdify to evaluate the expressions as numpy functions.

In [10]:
from sympy import lambdify

def lamb(func):
    return lambdify((x, y), func, modules='numpy')

red, orange, green, blue, pink, brown = map(lamb, (red, orange, green, blue, pink, brown))

Now it is easy to evaluate the functions at any given grid. (Ignore the warnings, we do not need the imaginary values.)

In [11]:
import numpy as np

spacing = 0.01
x, y = np.meshgrid(np.arange(-7.25, 7.25, spacing),
                   np.arange(-5, 5, spacing))

red = red(x, y)
orange = orange(x, y)
green = green(x, y)
blue = blue(x, y)
pink = pink(x, y)
brown = brown(x, y)
/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in sqrt
  """
/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in sqrt
  """
/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in sqrt
  """
/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in sqrt
  """
/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in sqrt
  """
/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in sqrt
  """

Finally let's plot the figure with matplotlib.

In [12]:
import matplotlib.pyplot as plt

colors = dict(red='#FF0000', orange='#FFA500', green='#008000',
              blue='#003399', pink='#CC0033', brown='#800000')

equations = dict(red=red, orange=orange, green=green,
                 blue=blue, pink=pink, brown=brown)
                 
fig, ax = plt.subplots(figsize=(8, 6))
for key, color in colors.items():
    ax.contour(x, y, equations[key], [0], colors=color, linewidths=3)
    
_ = ax.set_xticks([])
_ = ax.set_yticks([])

You can find these functions at wolfram. But we did not find any bat in the caves we visit during our vacation time... But we did find just this though:

In [13]:
from IPython.display import Image
Image('../../images/poco_encantado.jpg')
Out[13]:
In [14]:
HTML(html)
Out[14]:

This post was written as an IPython notebook. It is available for download or as a static html.

Creative Commons License
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

Comments