import sys
#The following are my extensions to sympy
from printer import *
from asycode import *
from piecewise import *
from dop import *
from sympy import *
Format(True,False) # Format is called to suppress arguments in printing general functions (F instead of F(x,y,z) for example)
(x,xp,w,a,y) = symbols(r'x xp \omega a y',real=True)
half = S(1)/S(2)
This notebook gives examples of several classes I have written as extensions of sympy. The printer
class includes the functions Format
, gprint
, and xpdf
. Format
enables the gprint
function which is an analog of the print
function but for $\LaTeX$ output. it's main use in a Jupyter Notebook is to allow one to annotate lines of output as demonstated below and to print more than one expression on a single line. Also it lets one print multiple annotated lines of output in a single cell. In a python script these functions along with xpdf
allow one to generate and view $\LaTeX$ output in a pdf file.
The piecewise
class implements piecewise functions on a fixed grid. I implemented this class when trying to calculate the inverse Fourier Transforms of powers of the $\f{\sinc}{\w/2}$ function. The interest in this was because I knew the the Fourier Transform of a gate function was a $\f{\sinc}{\w/2}$ function and the Fourier Transform of a triangular pulse was a $\f{\sinc^2}{\w/2}$ function. I wished to compute exactly the inverse Fourier Transform of $\f{\sinc^n}{\w/2}$ where $n>2$. Knowing that products in Fourier Transform space are given by convolution products in the time domain I wrote a piecewise function where I could calculate convolutions of piecewise functions. The piecewise function class included with sympy
could usually do the first convolution but not higher order convolutions. I think it suffered from being too general. Show below are the convolution powers of the gate function which correspond to powers of the Fourier Transforms of the gate function. One possible use of piecewise
is calculating either Fourier Transforms of inverse Fourier Transforms of functions which can be factored into products of functions that you can directly calculate the Fourier Transform or the inverse Fourier Transform. Note that $h^{*n}$ indicates the $n^{th}$ convolution product of $h$.
h1 = PieceWiseFunction([S(0),S(1),S(0)],[-half,half])
h2 = h1 & h1
h3 = h1 & h2
h4 = h1 & h3
h5 = h1 & h4
h6 = h1 & h5
def F(g):
return r'\mathcal{F}\{'+str(g)+r'\} = '
gprint(r'#Test of piecewise function convolution using repeated convolutions of the gate function $h$ and Fourier Transforms')
gprint(r'#of the convolutions which are powers of the sinc function.')
gprint('h = ',h1)
gprint(r'\FT{h} =',h1.FT(),r' = \f{\sinc}{\omega/2}')
gprint('h^{*2} = ',h2)
gprint(F('h^{*2}'),h2.FT().subs(w,2*y).expand(trig=True).subs(cos(y)**2,1-sin(y)**2).expand().subs(y,(w/2)),r' = \f{\sinc^2}{\omega/2}')
gprint('h^{*3} = ',h3)
gprint(F('h^{*3}'),h3.FT().subs(w,2*y).expand(trig=True).subs(cos(y)**2,1-sin(y)**2).expand().subs(y,(w/2)),r' = \f{\sinc^3}{\omega/2}')
gprint('h^{*4} =',h4)
gprint(F('h^{*4}'),h4.FT().subs(w,2*y).expand(trig=True).subs(cos(y)**2,1-sin(y)**2).expand().subs(y,(w/2)),r' = \f{\sinc^4}{\omega/2}')
gprint('h^{*5} =',h5)
gprint(F('h^{*5}'),h5.FT().subs(w,2*y).expand(trig=True).subs(cos(y)**2,1-sin(y)**2).expand().subs(y,(w/2)),r' = \f{\sinc^5}{\omega/2}')
gprint('h^{*6} =',h6)
gprint(F('h^{*6}'),h6.FT().subs(w,2*y).expand(trig=True).subs(cos(y)**2,1-sin(y)**2).expand().subs(y,(w/2)),r' = \f{\sinc^6}{\omega/2}')
The asycode
package converts python/sympy
code to Asymptote code. The code below plots the first six convolution powers of the gate function and a simple function to see how Asymptote can be used within a notebook or a python script. In this exercise the piecewise class has it's own Asymptote code generator. While it is also shown how to plot the sympy expression given by fofx
. The Asy.ExprAsy('F',fofx,'F,x,real')
code along with the sympy expression fofx
defines the Asymptote function to be plotted. Since Asymptote is a strongly typed langauge the ExprAsy
code is needed to define the variable types in the function. Note that the Asymptote plotting commands are in the asy_str
variable. Since one has to learn the commands anyway to plot something I did not see any utility in wrapping them in python functions. The main effort is in converting the python/sympy code for functions you wish to plot to Asymptote code. Since I have not used mathplotlib
I would be interested in a comparison how my Asymptote code compares to what would be required to generate a similar plot using mathplotlib
.
Asy.init() # Initialize Asymptote code file including LaTeX Macros and predefined function definitions
H = h1.asy('h1')+h2.asy('h2')+h3.asy('h3')+h4.asy('h4')+h5.asy('h5')+h6.asy('h6')
Asy.add(H) # Add piece wise functions to Asymptote code
fofx = sin(1/(x**2+0.05))**2
Asy.add(Asy.ExprAsy('F',fofx,'F,x,real')) # Add fofx function to Asymptote code with name 'F'
#Asymptote drawing commands
asy_str = r"""
picture Pic;
size(Pic,20cm,15cm,IgnoreAspect);
draw(Pic,graph(Pic,h1.f,-3,3,n=1000),black,legend="$h$");
draw(Pic,graph(Pic,h2.f,-3,3,n=500),red,legend="$h^{*2}$");
draw(Pic,graph(Pic,h3.f,-3,3,n=250),green,legend="$h^{*3}$");
draw(Pic,graph(Pic,h4.f,-3,3,n=250),blue,legend="$h^{*4}$");
draw(Pic,graph(Pic,h5.f,-3,3,n=250),orange,legend="$h^{*5}$");
draw(Pic,graph(Pic,h6.f,-3,3,n=250),cyan,legend="$h^{*6}$");
draw(Pic,graph(Pic,F,-3,3,n=5000),purple,legend="$\f{\sin^2}{(x^2+0.05)^{-1}}$");
xaxis(Pic,"$x$",Bottom,Ticks,xmin=-3,xmax=3);
yaxis(Pic,Ticks,ymin=0,ymax=1.5);
attach(Pic,legend(Pic,1),(-2.75,1.45),UnFill);
shipout(Pic);
"""
Asy.add(asy_str)
Asy.display('test_pwf_plots') # Show Asymptote plots in notebook
The dop
class creates a additional API for differentiation. It creates a partial derivative operator class, Pdop
, and a scalar differential operator class, Sdop
. The Pdop
class allows one to define a partial differntial operator say Dx = Pdop(x)
which is $\pdiff{}{x}$. Then Dx**3
is $\npdiff{}{x}{3}$ and if f
is a sympy function then (Dx**3)*f
is $\npdiff{f}{x}{3}$. A Sdop
(scalar differential operator) is a linear combination of sympy expressions (on the left) with Pdop
's on the right. For example x**2*Dx+x*(Dx**2)
. The following code demostrates the properties of Pdop
and Sdop
.
(x,y,z) = symbols('x y z',real=True)
r2 = x**2+y**2+z**2
Dx = Pdop(x)
Dy = Pdop(y)
Dz = Pdop(z)
nabla2 = Dx**2+Dy**2+Dz**2
gprint('r^2 = ',r2)
gprint(r'\nabla^2 = ',nabla2)
gprint(r'\nabla^2 r^2 = ',nabla2*r2)
gprint(r'r^2 \nabla^2 = ',r2*nabla2)
gprint(r'r^2 \nabla^2 r^2 = ',r2*nabla2*r2)
gprint(r'r^2+ \nabla^2 = ',r2+nabla2)
rdotgrad = x*Dx+y*Dy+z*Dz
gprint(r'\bs{r}\cdot\bs{\nabla} = ',rdotgrad)
gprint(r'(\bs{r}\cdot\bs{\nabla}) r^2= ',rdotgrad*r2)
gprint(r'(\bs{r}\cdot\bs{\nabla})(xyz) = ',rdotgrad*(x*y*z))
F = Function('F',real=True)(x,y,z)
G = Function('G',real=True)(x,y,z)
gprint(r'\nabla^2 F = ',nabla2*F)
gprint(r'\nabla^2 (FG) = ',nabla2*(F*G))
gprint(r'\nabla^2 (F/G) = ',nabla2*(F/G))
gprint(r'\nabla^2 (\sqrt{F}) = ',nabla2*(sqrt(F)))
With little effort the API could be extened to vector operators such as $\nabla$ in rectangular coordinate or other types of coordinate systems.