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)
equ
# d d
#[ -r₀ + f₀(t) + ──(f₀(t)), ──(f₁(t)) - 1 ]
# dt dt
init
#[-δ - 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)])]
intcurve_coords
# -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)
coords_in_cartesian
#⎡⎛ -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)
[
][]
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)
equ
#[
# ⎛ f₀(t)⎞ d
#- r₀⋅sin⎜1 - ─────⎟ + ──(f₀(t))
# ⎝ r₀ ⎠ dt ,
#
#d
#──(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.
[]: https://krastanov.files.wordpress.com/2012/06/integral_curve.png