Sympy and EMT – Comparison

In Euler Math Toolbox, symbolic expressions are seamlessly integrated by the &-syntax. The expressions and functions are handled externally via the CAS Maxima, but mirrors exists in the numerical side of EMT as simple strings or numerical functions.

Let me give an example.

>function fx(t) &= t*cos(2*pi*t)
 
                             t cos(2 pi t)
 
>function fy(t) &= t*sin(2*pi*t)
 
                             t sin(2 pi t)
 
>plot2d("fx","fy",xmin=0,xmax=1,r=1):

The functions „fx“ and „fy“ are symbolic functions. I.e., Maxima defines the functions in its system, and they can be used in Maxima.

>&fx(z^2+1)
 
                        2                 2
                      (z  + 1) cos(2 pi (z  + 1))

But the functions are also defined in EMT as numerical functions. Thus, they can be used like any other numerical function in EMT.

>fx(linspace(0,1,10))
 [0,  0.0809017,  0.0618034,  -0.0927051,  -0.323607,  -0.5,  -0.48541,
 -0.216312,  0.247214,  0.728115,  1]

Because the functions contain only elementary operations, they can work on vectors too. This is called a „vectorized“ or „mapped“ function. Consequently, they can be called in other functions of EMT. E.g., the „plot2d“ function can plot curves, if two functions (or two expressions) are given.

Continuing our example, we compute the length of the curve using

\(l(f) = \int\limits_0^1 \|f'(t)\| \,dt.\)

>function f(t) &= [fx(t),fy(t)]
 
                     [t cos(2 pi t), t sin(2 pi t)]
 
>function df(t) &= diff(f(t),t)
 
        [cos(2 pi t) - 2 pi t sin(2 pi t), 
                                      sin(2 pi t) + 2 pi t cos(2 pi t)]
 
>integrate(&norm(df(x)),0,1)
 3.3830442855
>&integrate(norm(df(t)),t,0,1), %()
 
                                               2
                   asinh(2 pi) + 2 pi sqrt(4 pi  + 1)
                   ----------------------------------
                                  4 pi
 
 3.3830442855

To compute the length of the curve, we define a vector valued symbolic function „f“. That function does work in Maxima symbolically, or in EMT numerically.

A side remark: If „f“ is applied to a vector instead of one point, it will compute „fx“ and „fy“ and concatenate these vectors, which is not useful most of the time. In future versions, EMT should be able to generate a matrix of values if a vector valued function is applied to a vector of values. In our application, we only need „f“ as a symbolic function.

In the same way, we can then define a function „df(t)“ for the derivative vector simply by symbolic differentiation. „df“ will also work symbolically and numerically. The function „norm“ is a pre-defined symbolic function. You can plug in a symbolic vector and get an expression which represents a function of „t“.

>&norm(v)
 
                         sqrt(v . transpose(v))
 
>&norm(df(t))
 
                                               2
        sqrt((cos(2 pi t) - 2 pi t sin(2 pi t))
                                                                     2
                                 + (sin(2 pi t) + 2 pi t cos(2 pi t)) )

We then use the adaptive integral of EMT to integrate the expression „norm(df(x))“ between 0 and 1.

Finally, Maxima can compute the „exact“ integral in symbolic form, which we evaluate numerically to verify our previous result.

How does this all translate to Python?

The first problem is that there are no symbolic functions in Sympy. There are only expressions, which can be handled by the symbolic features of Sympy. And, of course, there is no automatic parallel definition of a numerical function for Numpy as in EMT. To evaluate a symbolic expression at a point, „subs“ can be used, or it can be made a numerical function by „lambdify“.

Here is an example.

import sympy as sp
import numpy as np

t = sp.symbols('t')
fe = t*sp.cos(2*np.pi*t)
print(fe.subs(t,0.2))
f = sp.lambdify(t,fe)
print(f(0.2))

  -->
  0.0618033988749895
  0.061803398874989326

As you see, „subs“ just substitutes a value for „t“. This is not going to work of „sp.pi“ is in the expression instead of „np.pi“.

t = sp.symbols('t')
fe = t*sp.cos(2*sp.pi*t)
print(fe.subs(t,0.2))
f = sp.lambdify(t,fe)
print(f(0.2))

  -->
  0.2*cos(0.4*pi)
  0.06180339887498949

For Sympy, „pi“ is a mathematical variable with some properties and a value. Because „cos(2pi) „is known, „f(1.0)“ does work as expected. But „f(0.2)“ leaves the expression and is not evaluated numerically.

„lambdify“ translates a symbolic expression so that it can be used by Numpy, effectively translating „sp.pi“ to „np.pi“. Thus, „f“ works in both cases.

Back to our example. We use what we learned to trick Python into doing what we want.

import sympy as sp
import numpy as np

t = sp.symbols('t')
fxe = t*sp.cos(2*sp.pi*t)
fye = t*sp.sin(2*sp.pi*t)
fx = sp.lambdify(t,fxe)
fy = sp.lambdify(t,fye)

vt = np.linspace(0,1,400)
x = fx(vt)
y = fy(vt)
fig,ax = plt.subplots(figsize=(8,8))
ax.set_xlim(-1,1)
ax.set_ylim(-1,1)
ax.grid(True)
ax.plot(x,y);

It might be a good idea to define a function „plot_curve“ to produce the plot, like the following.

def plot_curve (fx,fy,xmin,xmax,n=100,figsize=(8,8),r=1):
    vt = np.linspace(0,1,400)
    x=fx(vt); y=fy(vt);
    fig,ax = plt.subplots(figsize=figsize)
    ax.set_xlim(-r,r); ax.set_ylim(-r,r)
    ax.grid(True)
    return ax.plot(x,y)

plot_curve(fx,fy,0,1,400);

Such functions can go into a Python module. I already started a module „renestools.py“ with functions I need to make my life in Python easier.

Now, let us compute the length. The integration is done using an adaptive Gauss quadrature in my package. I explained in a previous posting, how this is working. Note that we need to pass a vectorized function to the integration. „fnorm“ is not vectorized (in contrast to the expression „norm(f(t))“ which we used in EMT). „np.vectorize“ transforms any function for single real values into vectorized functions which will work on Numpy matrices or vectors.

fe = sp.Matrix([fxe,fye])
fde = sp.diff(fe,t)

def norm(v):
    return np.sqrt(np.sum(v*v))

fd = sp.lambdify(t,fde)
def fnorm(t):
    return norm(fd(t))

import renestools as rt
rt.integrate(np.vectorize(fnorm),0,1)

  --> 3.3830442855028062

To compute the integral in symbolic form, we need to help Sympy with „trigsimp“. We evaluate the result to 20 digits. (We could have done that in Maxima too, by the way.)

res = sp.integrate(sp.trigsimp(sp.sqrt(fde[0]**2+fde[1]**2)),(t,0,1))
print(res)
print(res.evalf(20))

  -->
  asinh(2*pi)/(4*pi) + sqrt(1 + 4*pi**2)/2
  3.3830442855028069597

As you see, Maxima and EMT are much more intuitive than Sympy and Numpy. I say this in the hope that many of you find their way to Euler Math Toolbox and help the development. Currently, Python has a lot more potentially useful toolboxes than EMT.

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht.

Diese Website verwendet Akismet, um Spam zu reduzieren. Erfahre mehr darüber, wie deine Kommentardaten verarbeitet werden.