Integral Curves of Vector Fields in SymPy

A week or two ago I implemented some basic functionality for work with integral curves of vector fields. However, I needed to make additional changes in other parts of SymPy in order for the ODE solver to work with systems of equations and with initial conditions. I also wanted to get my plotting module merged so I can show some visualizations if necessary.

Now that all this is ready (even though not everything is merged in SymPy master) I can show you some of the most basic capabilities implemented in the differential geometry module. First, we start with the boilerplate:

A Simple Field

from sympy.diffgeom import *
from sympy.diffgeom.Rn import * # This gives me:
                                    #   - R2_p - the polar coord system
                                    #   - R2_r - the rectangular coord system
                                    #   - x,y,r,theta - the base scalar fields
                                    #   - e_x, ... - the base vector fields
# Define some fields to play with
# (these are the same fields, defined in two different ways):
vector_field_circular_p = R2_p.e_theta
vector_field_circular_r = -R2.y*R2.e_x + R2.x*R2.e_y
# Define the same point in two different ways
point_p = R2_p.point([1,pi/2])
point_r = R2_r.point([0,1])

The r index is for rectangular coordinate systems and the p index is for polar.

Now using intcurve_diffequ we can generate the differential equations for the integral curve. This function also generates the required initial conditions:

#        vector field      free parameter for    starting point
#                  |        the curve     |     /          coord system
#                  v                      v    v      v-- for the equation
intcurve_diffequ(vector_field_circular_p, t, point_p, R2_p)
# output:
#   d            d
#([ ──(f₀(t)),   ──(f₁(t)) - 1  ],
#   dt           dt
#                        π
# [f₀(0) - 1,    f₁(0) - ─      ])
#                        2

intcurve_diffequ(vector_field_circular_p, t, point_p, R2_r)
# output:
#           d                     d
#([ f₁(t) + ──(f₀(t)),   -f₀(t) + ──(f₁(t))  ],
#           dt                    dt
# [ f₀(0),               f₁(0) - 1           ])

Here we have equations for the functions f_0 and f_1 which are by convention the names that intcurve_diffequ gives for the first and second coordinate.

The cool thing is that we can mix the coordinate systems in any way we wish. The code will automatically make the needed coordinate transformation and return the equations in the demanded coordinate system independently of the coordinate systems in which the input objects were defined (at worst you will need to call some simplification routines):

a = intcurve_diffequ(vector_field_circular_p, t, point_p, R2_r)
a == intcurve_diffequ(vector_field_circular_p, t, point_r, R2_r)
# True
a == intcurve_diffequ(vector_field_circular_r, t, point_r, R2_r)
# True
a == intcurve_diffequ(vector_field_circular_r, t, point_p, R2_r)
# True

Solving the equations actually gives (this solver is not yet in SymPy master as of the time of writing):

equ_r, init_r = intcurve_diffequ(vector_field_circular_r, t, point_r, R2_r)
sol_r = dsolve(equ_r+init_r, [Function('f_0')(t), Function('f_1')(t)])
[simplify(s.rewrite(sin)) for s in sol_r] # some simplification
#[f₀(t) = -sin(t), f₁(t) = cos(t)]            # is necessary because
                                              # dsolve returned complex
                                              # exponentials

Even simpler:

equ_p, init_p = intcurve_diffequ(vector_field_circular_p, t, point_p, R2_p)
dsolve(equ_p+init_p, [Function('f_0')(t), Function('f_1')(t)])
# output:
#[f₀(t) = 1,
#             π
# f₁(t) = t + ─
#             2]

This is obviously just a circle (did I mentioned that the vector field that I defined is circular). There is no need to plot it as it is fairly simple. However a slight change will render the field a bit more interesting:

Radial Component

# A circular field that also pushes in radial direction
# towards an equilibrium radius.
v_field = R2.e_theta + (r0 - R2.r)*R2.e_r
# An initial position slightly away from the
# equilibrium one.
start_point = R2_p.point([r0+delta, 0])
equ, init = intcurve_diffequ(v_field, t, start_point)
#                d            d
#[ -r₀ + f₀(t) + ──(f₀(t)),   ──(f₁(t)) - 1 ]
#                dt           dt

#[-δ - r₀ + f₀(0), f₁(0)]

dsolve(equ+init, [Function('f_0')(t), Function('f_1')(t)])
#            -t
#[f₀(t) = δ⋅ℯ   + r₀, f₁(t) = t]

This gives a spiral tending towards the equilibrium radius r_0. Let us extract the coordinates from these equations and plot the resulting curve:

intcurve_coords = [eq.rhs for eq in dsolve(equ+init, [Function('f_0')(t), Function('f_1')(t)])]
#    -t
#[δ⋅ℯ   + r₀, t]

# We need this in Cartesian coordinates for the plot routine.
# We could have solved for Cartesian coordinates since the
# beginning, however our current approach permits us to see
# how to use the `CoordSys` classes to change coordinate systems:
coords_in_cartesian = R2_p.point(intcurve_coords).coords(R2_r)
#⎡⎛   -t     ⎞       ⎤
#⎢⎝δ⋅ℯ   + r₀⎠⋅cos(t)⎥
#⎢                   ⎥
#⎢⎛   -t     ⎞       ⎥
#⎣⎝δ⋅ℯ   + r₀⎠⋅sin(t)⎦

# Substitute numerical values for the plots:
x,y = coords_in_cartesian.subs({delta:0.5, r0:1})
plot(x,y, (t,0,4*pi))
#Plot object containing:
#[0]: parametric cartesian line:
#      ((1 + 0.5*exp(-t))*cos(t), (1 + 0.5*exp(-t))*sin(t))
#      for t over (0.0, 12.566370614359172)

simple integral curve

This is all great, but what happens if one has to work with more complicated fields. For instance the following simple field will not permit analytical solution:

No Analytical Solution

v_field = R2.e_theta + r0*sin(1 - R2.r/r0)*R2.e_r
equ, init = intcurve_diffequ(v_field, t, start_point)
#        ⎛    f₀(t)⎞   d
#- r₀⋅sin⎜1 - ─────⎟ + ──(f₀(t))
#        ⎝      r₀ ⎠   dt       ,
#──(f₁(t)) - 1
#dt           ]

For cases like this one the user can take advantage of one of the numerical ODE solvers from scipy. Or sticking to symbolic work he can use the intcurve_series function that gives the series expansion for the curve:

intcurve_series(v_field, t, start_point, n=1)
#⎡δ + r₀⎤
#⎢      ⎥
#⎣  0   ⎦

intcurve_series(v_field, t, start_point, n=2)
#⎡            ⎛    δ + r₀⎞     ⎤
#⎢δ + r₀⋅t⋅sin⎜1 - ──────⎟ + r₀⎥
#⎢            ⎝      r₀  ⎠     ⎥
#⎢                             ⎥
#⎣              t              ⎦

intcurve_series(v_field, t, start_point, n=4, coeffs=True)
#⎡δ + r₀⎤
#⎢      ⎥
#⎣  0   ⎦,
#⎡        ⎛    δ + r₀⎞⎤
#⎢r₀⋅t⋅sin⎜1 - ──────⎟⎥
#⎢        ⎝      r₀  ⎠⎥
#⎢                    ⎥
#⎣         t          ⎦,
#⎡     2    ⎛    δ + r₀⎞    ⎛    δ + r₀⎞⎤
#⎢-r₀⋅t ⋅sin⎜1 - ──────⎟⋅cos⎜1 - ──────⎟⎥
#⎢          ⎝      r₀  ⎠    ⎝      r₀  ⎠⎥
#⎢                  2                   ⎥
#⎢                                      ⎥
#⎣                  0                   ⎦,
#⎡      ⎛     2                  2            ⎞                ⎤
#⎢    3 ⎜      ⎛    δ + r₀⎞       ⎛    δ + r₀⎞⎟    ⎛    δ + r₀⎞⎥
#⎢r₀⋅t ⋅⎜- sin ⎜1 - ──────⎟ + cos ⎜1 - ──────⎟⎟⋅sin⎜1 - ──────⎟⎥
#⎢      ⎝      ⎝      r₀  ⎠       ⎝      r₀  ⎠⎠    ⎝      r₀  ⎠⎥
#⎢                              6                              ⎥
#⎢                                                             ⎥
#⎣                              0                              ⎦]

However these series do not always converge to the solution, so care should be taken.

There are other amusing possibilities already implemented, however I will write about them another time.

If you want to suggest more interesting examples, please do so in the comments.


Scalar and Vector Fields in SymPy – First Steps

The Differential Geometry module for SymPy already supports some interesting basic operations. However, it would be appropriate to describe its structure before giving any examples.

First of all, there are the Manifold and Patch classes which are just placeholders. They contain all the coordinate charts that are defined on the patch and do not provide, for instance, any topological information. This leads us to the CoordSystem class which contains all the coordinate transformation logic. For example, if I want to define the \mathbb{R}^2 euclidean manifold together with the polar and Cartesian coordinate systems I would do:

R2 = Manifold('R^2', 2)
# Patch and coordinate systems.
R2_origin = Patch('R^2_o', R2)
R2_r = CoordSystem('R^2_r', R2_origin)
R2_p = CoordSystem('R^2_p', R2_origin)

# Connecting the coordinate charts.
x, y, r, theta = [Dummy(s) for s in ['x', 'y', 'r', 'theta']]
R2_r.connect_to(R2_p, [x, y],
                      [sqrt(x**2 + y**2), atan2(y, x)],
                inverse=False, fill_in_gaps=False)
R2_p.connect_to(R2_r, [r, theta],
                      [r*cos(theta), r*sin(theta)],
                inverse=False, fill_in_gaps=False)

All following examples will be about the \mathbb{R}^2 manifold which is already implemented in the code for the module. Also, notice the use of the inverse and fill_in_gaps flags. When they are set to True the CoordSystem classes try to automatically deduce the inverse transformations using SymPy’s solve function.

Now that we have a manifold we would like to create some fields on it and define some points that belong to the manifold. The points are implemented in the Point class. You need to specify some coordinates when you define the point, however after that the object is completely coordinate-system-idependent.

# You need to specify coordinates in some coordinate system
p = Point(R2_p, [r0, theta0])

Then one can define fields. ScalarField takes points to real numbers and VectorField is an operator on ScalarField taking a scalar field to another scalar field by applying a directional derivative. For example, here x and y are the scalar fields taking a point and returning it’s coordinate and d_dx and d_dy are the vector fields \frac{\partial}{\partial x} and \frac{\partial}{\partial y}. R2_r is the Cartesian coordinate system and R2_p is the polar one.

R2_r.x(p) == r0*cos(theta0)
# R2_r.d_dx(R2_r.x) is a also scalar field
R2_r.d_dx(R2_r.x)(p) == 1

Looking at how can these fields be defined:

# For a ScalarField you provide the transformation in some coordinate system
R2_r.x = ScalarField(R2_r, [x0, y0], x0)
#                     /      |        ^-------- the result
#     the coord system     the coordinates

# For a VectorField you provide the components in some coordinate system
R2_r.d_dx = VectorField(R2_r, [x0, y0], [1, 0])
#                        /      |         ^-------- the components
#         the coord system     the coordinates

Obviously one can define much more interesting fields. For instance the potential due to a point charge at the origin is:

potential = ScalarField(R2_p, [r0, thata0], -1/r0)
# And to reiterate, the definition does not limit you
# to use it only in this coordinate system. For instance:
potential(R2_r.point([x0, y0])) == 1/sqrt(x0**2 + y0**2)

However there is another more intuitive way to do it:

# R2_p.r is the scalar field that takes a point and returns the r coordinate
potential2 = 1/R2_p.r
potential2(R2_r.point([x0, y0])) == 1/sqrt(x0**2 + y0**2))

And this new object potential2 is not an instance of ScalarField. It is actually a normal SymPy expression tree that contains a ScalarField somewhere in its leafs (namely in this case it is Pow(R2_p.r, -1)). However, due to the change to one of the base classes of SymPy that I did in this pull request it is now possible for such tree to be a python callable, by recursively applying the argument to each callable leaf in the tree. This change is still debated and it may be reverted.

Vector fields can also be build in this manner. However, they pose a problem. What happens when you multiply a vector field and a scalar field? This operation should give another vector field. And here is a possible problem with the approach of recursively callable expressions trees:

# Naively this operation will call a scalar field on
# another scalar field which is nonsense:
(R2_r.x * R2_r.d_dx)(R2_r.x) == R2_r.x(R2_r.x) * R2_r.d_dx(R2_r.x)
#                         nonsense----^

The current solution is for scalar_field(not_point) to return the callable itself. Thus we have:

(R2_r.x * R2_r.d_dx)(R2_r.x) == R2_r.x * R2_r.d_dx(R2_r.x)
#\________________/ \______/    \_______________________/
#   vector field        ^---scalar fields---^

This way there is no need for complicated logic in __mul__ nor is there need for addition subclasses of Expr in order to accommodate this behavior.

There is not much more to be said about the structure of the module. There are some other nice things already implemented like integral curves, however I will discuss these in a later post.

Among the things that should be done at some point:

  • Should vector fields be callable on points? If yes, what the result should be? An abstract vector, a tuple of coordinates in a certain coordinate system, something else?
  • There are many expressions generated by this code that are not simple enough. I should work on the simplification routines and on the differential geometry module itself in order to get more canonical expressions.
  • The last point is also valid about the solvers: some coordinate transformations are too complicated for the solvers to find the inverse transformation.
  • Manifold and Patch have name attributes. Are these necessary? What is the role of name attributes in SymPy besides printing?
  • Start using Lambda where applicable.
  • Follow better the class structure of SymPy.