11 Using symbolic math within Julia

This page briefly describes the use of symbolic math for the Calculus II topics of MTH 232 within Julia utilizing the SymPy library of Python. Either the SymPy or SymPyPythonCall packages provide access, the difference being the package used for interop between Julia and Python.

The SymPy package is loaded when the MTH229 package is. We also add the Plots package to handle our plotting needs.¹

using MTH229 # loads SymPy (v2.0 or greater)
using Plots
plotly()

┌ Warning: Failed to load integration with PlotlyBase & PlotlyKaleido.
│   exception =
│    ArgumentError: Package PlotlyBase not found in current path.
│    - Run `import Pkg; Pkg.add("PlotlyBase")` to install the PlotlyBase package.
│    Stacktrace:
│      [1] macro expansion
│        @ ./loading.jl:1772 [inlined]
│      [2] macro expansion
│        @ ./lock.jl:267 [inlined]
│      [3] __require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1753
│      [4] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│      [5] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│      [6] require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1746
│      [7] top-level scope
│        @ ~/.julia/packages/Plots/ju9dp/src/backends.jl:565
│      [8] eval
│        @ ./boot.jl:385 [inlined]
│      [9] _initialize_backend(pkg::Plots.PlotlyBackend)
│        @ Plots ~/.julia/packages/Plots/ju9dp/src/backends.jl:564
│     [10] backend
│        @ ~/.julia/packages/Plots/ju9dp/src/backends.jl:244 [inlined]
│     [11] plotly()
│        @ Plots ~/.julia/packages/Plots/ju9dp/src/backends.jl:87
│     [12] top-level scope
│        @ In[2]:3
│     [13] eval
│        @ ./boot.jl:385 [inlined]
│     [14] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
│        @ Base ./loading.jl:2076
│     [15] softscope_include_string(m::Module, code::String, filename::String)
│        @ SoftGlobalScope ~/.julia/packages/SoftGlobalScope/u4UzH/src/SoftGlobalScope.jl:65
│     [16] execute_request(socket::ZMQ.Socket, msg::IJulia.Msg)
│        @ IJulia ~/.julia/packages/IJulia/Vo51o/src/execute_request.jl:67
│     [17] #invokelatest#2
│        @ ./essentials.jl:892 [inlined]
│     [18] invokelatest
│        @ ./essentials.jl:889 [inlined]
│     [19] eventloop(socket::ZMQ.Socket)
│        @ IJulia ~/.julia/packages/IJulia/Vo51o/src/eventloop.jl:8
│     [20] (::IJulia.var"#15#18")()
│        @ IJulia ~/.julia/packages/IJulia/Vo51o/src/eventloop.jl:38
└ @ Plots ~/.julia/packages/Plots/ju9dp/src/backends.jl:574

Plots.PlotlyBackend()

There are a few basic ideas behind the SymPy and SymPyPythonCall packages

Symbolic variables and values can be easily created
Generic operations, such as math operations, have specialized methods for symbolic values and expressions that utilize the underlying Python library to manipulate or create new symbolic expressions.
Functions specific to the Python library are accessed through the sympy object, such as sympy.expand_trig.
Methods of objects in the Python library are called with the same dot call syntax, e.g. obj.subs(...) to call the subs method of the object

11.1 Symbolic variables

There are numerous ways to create symbolic objects in SymPy (symbols, sympify, …) but we use the @syms macro to create variables with a given name:

@syms x

(x,)

One or more variables can be created. In this example, the last one creates a range of values:

@syms x y zs[1:5]

(x, y, Sym{PyCall.PyObject}[zs₁, zs₂, zs₃, zs₄, zs₅])

Variables can have basic assumptions placed on them (positive, real), e.g.

@syms a::positive, b::real

(a, b)

Symbolic functions can be naturally specified

@syms u()
u(x)

$u{\left(x \right)}$

11.2 Symbolic numbers

If x is symbolic, the expression 2x will also be symbolic

2x

$2 x$

Behind the scenes, before multiplying, 2 and x are promoted to a common type, which will be symbolic. So most regular Julia expressions involving one or more symbolic values will simply yield symbolic answers without modification. Integers and rational numbers promote to exact counterparts; floating point numbers convert to floating point values.

Still, exactness can be lost when the floating point values are created before becoming symbolic:

1/10 * x + 2/10 * x - 3/10 * x

$5.55111512312578 \cdot 10^{-17} x$

The above “should” be $0$, but it isn’t, due to the conversion to floating point prior to become symbolic. For example, the expression 1/10 * x first divides 10 into 1, producing a floating point approximation for 1/10 and this is promoted to symbolic when multiplied by x.

Using rational numbers or re-expressing for earlier promotion to symbolic will keep the exactness:

1//10 * x + 2//10 * x - 3//10*x, 1x/10 + 2x/10 - 3x/10

(0, 0)

Sometimes symbolic values must be created before passing to a function, so dispatch to the symbolic function occurs. Compare these two calls, where the Sym constructor is used to create symbolic numbers:

log(Sym(2)), log(2)

(log(2), 0.6931471805599453)

When in doubt, use Sym

But most of the time, just being mindful to avoid non-symbolic divisions and function calls with Julia numbers is enough.

11.3 Substitution

Substituting one value in for another, is done many ways, for example there is subs, replace, and xreplace. We discuss subs.

The subs function (or object method) replaces an old value for an new value after sympifying² the arguments. The old/new pairs can be specified with a tuple (e.g. (old, new)), a dictionary, or the more Julian syntax, a pairs expression old => new, with multiple substitutions possibly combined in a tuple.

@syms a b c x y
ex = a*x^2 + b* x + c - y
subs(ex, x => y-1)

$a \left(y - 1\right)^{2} + b \left(y - 1\right) + c - y$

subs(ex, a => 2)

$b x + c + 2 x^{2} - y$

subs(ex, a=>3, b=>2, c=>1)

$3 x^{2} + 2 x - y + 1$

The latter is akin to function evaluation. As such, the call notation for symbolic expressions is defined to use subs:

ex(a=>3, b=>2, c=>1)

$3 x^{2} + 2 x - y + 1$

11.4 Numerical evaluation

An exact symbolic value may be the result of a computation, but it may be of interest to convert that value to a number (e.g. 1.414... instead of $\sqrt{2}$).

SymPy objects have an evalf method that finds a floating-point value.

u = exp(Sym(2)) - exp(-Sym(2))
u

$- \frac{1}{e^{2}} + e^{2}$

The value of u is exact. In this next cell, we call the evalf method and compare to a computation within Julia:

u.evalf(), exp(2) - exp(-2)

(7.25372081569404, 7.253720815694038)

The default for .evalf() is 15 decimal points; similar but not the same as for Julia. The method allows a specification of the number of digits in its first argument. (The verbose=true argument illustrates some of the algorithm employed.)

u.evalf(50)

$7.2537208156940375353364279656025234097726840246423$

The result of evalf looks numeric, but internally the values are symbolic values within SymPy.³ In SymPy.jl the N function converts the exact value to a Julia number type with an attempt to match the type in Python with a corresponding type in Julia. In this example, that is Rational.

u = Sym(1) / 1000
u.evalf(), N(u)

(0.00100000000000000, 1//1000)

For repeated conversions of expressions to numeric values, the lambdify method creates a Julia function from an expression, which can then be called as any other function.

@syms x
ex = cos(x)
l =  lambdify(ex)
ex(PI/3), ex(pi/3), l(pi/3)

(1/2, 0.500000000000000, 0.5000000000000001)

11.5 Solving an equation

SymPy can be used to solve many algebraically solvable equations. The solve and solveset functions are used.

By default, an expression is interpreted as an equation set equal to 0. However, we will always model the use of rhs ~ lhs to separate two sides of an equation.⁴

Here the equation $x^2 = 2$ is solved for:

@syms x y::real z::positive
solve(x^2 ~ 2)

2-element Vector{Sym{PyCall.PyObject}}:
 -√2
  √2

Simple assumptions on the variables are respected:

solve(y^2 + 1 ~ 0), solve(z + 1 ~ 0)

(Any[], Any[])

Equations can have symbolic variables included. In which case, passing in the second argument the variable(s) to solve for is done:

@syms a b c
solve(a*x^2 + b*x + c ~ 0, x)

2-element Vector{Sym{PyCall.PyObject}}:
 (-b - sqrt(-4*a*c + b^2))/(2*a)
 (-b + sqrt(-4*a*c + b^2))/(2*a)

The solve function has many different types of return values. The solveset function always returns a set. For finite sets, a Julia Set object is returned. For infinite sets, this is not the case.

solveset(a*x^2 + b*x + c ~ 0, x)

Set{Sym} with 2 elements:
  -b/(2*a) - sqrt(-4*a*c + b^2)/(2*a)
  -b/(2*a) + sqrt(-4*a*c + b^2)/(2*a)

This may look the same as the last call using solve, but a finite set is unordered, so indexing doesn’t work unless the set is collected (using collect) into a vector.

For infinite sets, the answer is queried. In this command we check for inclusion:⁵

u = solveset(sin(x) ~ 0, x)
pi in u

true

Intersecting with a finite interval will create a finite set:⁶

S = sympy.Interval(-5, 5)
intersect(S, u)

Set{Sym} with 3 elements:
  0
  -pi
  pi

11.5.1 Numerically solving an equation

Not all equations can be solved symbolically. For example, finding the solution(s) to $\cos(x) = x$ will result in an error:

try solve(cos(x) ~ x) catch err "No algorithms are implemented to solve equation -x + cos(x)" end

"No algorithms are implemented to solve equation -x + cos(x)"

Within SymPy, there is nsolve:

sympy.nsolve(cos(x) ~ x, pi/4)

$0.739085133215161$

Newer version of Roots

The find_zero function of the Roots package (if new enough) accepts an equation for input:

find_zero(cos(x) ~ x, (0, pi/2))  # using [0,π/2] as a bracketing interval

Roots is re-exported by the MTH229 package.

11.5.2 Systems of equations

The solve function can solve systems of equations. These can be specified by wrapping one or more equations into a container (using parentheses to form a tuple is suggested) and specifying the symbols similarly:

@syms x::real y::real
solve((x^2 + y^2 ~ 1, y ~ x^2), (x,y))

2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
 (-sqrt(-1/2 + sqrt(5)/2), -1/2 + sqrt(5)/2)
 (sqrt(-1/2 + sqrt(5)/2), -1/2 + sqrt(5)/2)

11.6 Calculus specific functions

SymPy provides function for specific topics of Calculus II.

11.6.1 Limits

Limits may be taken symbolically using limit. The specification $x \rightarrow c$ can be specified with a tuple ((x,c)) or pair (x => c). Directional limits are specified with dir="-", dir="+" (the default), or dir="+-".

@syms a::positive b::positive x::real
ex = ((a^x - x*log(a)) / (b^x - x*log(b)))^(1/x^2)
L = limit(ex, x=>0)

$e^{\frac{\log{\left(a \right)}^{2}}{2} - \frac{\log{\left(b \right)}^{2}}{2}}$

For a given set of values of a and b, we can see different answers:

N(L(a=>2, b=>3))

0.6954139679176876

Watch out for early conversion

The above example will fail if instead of making a and b symbolic and then substituting values of 2 and 3, the expression were substituted in at the outset:

@syms x::real
a, b = 2, 3
ex = ((a^x - x*log(a)) / (b^x - x*log(b)))^(1/x^2)
L = limit(ex, x=>0)

$\infty$

This is because both log(a) and log(b) are then inexact values for $\log(2)$ and $\log(3)$ respectively.

Simplification can be necessary

In this example we wish to take the limit as $x \rightarrow \infty$ of

\[ f(x) = \frac{\cos{\left(2 x \right)} + 1}{\left(x + \sin{\left(x \right)} \cos{\left(x \right)}\right) e^{\sin{\left(x \right)}} \cos{\left(x \right)} + \left(- \sin^{2}{\left(x \right)} + \cos^{2}{\left(x \right)} + 1\right) e^{\sin{\left(x \right)}}} \]

A straightforward approach doesn’t work:

den(x) = 1//2*sin(2x) + x
num(x) = exp(sin(x))*(cos(x)*sin(x) + x)
up, vp = diff(den(x),x), diff(num(x),x)
limit(up/vp, x=> oo)

$\lim_{x \to \infty}\left(\frac{\cos{\left(2 x \right)} + 1}{\left(x + \sin{\left(x \right)} \cos{\left(x \right)}\right) e^{\sin{\left(x \right)}} \cos{\left(x \right)} + \left(- \sin^{2}{\left(x \right)} + \cos^{2}{\left(x \right)} + 1\right) e^{\sin{\left(x \right)}}}\right)$

As seen, SymPy returns an unevaluated limit expression, as no answer can be found. However, if the expression is simplified, SymPy will return the correct answer:

limit(simplify(up/vp), x => oo)

$0$

11.6.2 Differentiation

Derivatives are taken through diff(ex, var, ...) with variants for multiple or mixed derivatives. The basic usage is pretty straightforward. The derivative with respect to $x$ of an expression is found as follows:

@syms x a b
ex = sin(a*x - b)
diff(ex, x)

$a \cos{\left(a x - b \right)}$

Second derivatives can be more succinctly expressed by adding more variables:

diff(ex, x, x) # 2nd derivative

$- a^{2} \sin{\left(a x - b \right)}$

11.6.3 Integration

In Calculus II there are techniques of integration to be learned:

integration by parts
trigonometric integrals
partial fractions

The definition of integration can be extended to incorporate infinities:

improper integrals

Further, there are several formulas where an integral takes on a geometric meaning beyond the area under a curve:

Area between 2 curves: $A = \int_a^b (f(x) - g(x)) dx$.
volume of revolution: $V = \int_a^b \pi r(x)^2 dx$.
cylindrical shell: $V = \int_a^b 2\pi x f(x) dx$.
arc length: $L = \int_a^b \sqrt{1 + f'(x)^2} dx$,
surface area: $SA = \int_a^b 2\pi f(x) \sqrt{1+ f'(x)^2} dx$.

In SymPy, the integrate function computes integrals symbolically by finding an anti-derivative. Bear in mind, not all integrands have an antiderivative that can be algebraically expressed! SymPy uses the tricks of integration above, and in addition implements an algorithm (in part) due to Risch.

The indefinite integral, $\int f(x) dx$, is computed with integrate(f(x), x) where f(x) is some expression depending on x and may have symbolic parameters.

Definite integrals, $\int_a^b f(x)dx$, are computed with integrate(f(x), (x, a, b)).

For example:

@syms x c a b
integrate(x * exp(x^2), x) # integration by parts

$\frac{e^{x^{2}}}{2}$

integrate(sin(x)^2, x)  # integration by parts

$\frac{x}{2} - \frac{\sin{\left(x \right)} \cos{\left(x \right)}}{2}$

integrate(x^5/(36x^2 + 1)^(3//2), (x, 0, 1//6)) # trig subs with x = tan(theta)/6

$\frac{1}{17496} - \frac{11 \sqrt{2}}{279936}$

ex = (x^4 + 1) / ( (x^2 + 1)^2 * (x^2 - 4)^2)
integrate(ex, x) # partial fraction

$\frac{- 9 x^{3} - 49 x}{200 x^{4} - 600 x^{2} - 800} - \frac{37 \log{\left(x - 2 \right)}}{4000} + \frac{37 \log{\left(x + 2 \right)}}{4000} - \frac{\operatorname{atan}{\left(x \right)}}{125}$

We can see how this is done by taking the partial fraction decomposition provided by apart:

us = apart(ex)

$- \frac{6}{125 \left(x^{2} + 1\right)} + \frac{2}{25 \left(x^{2} + 1\right)^{2}} + \frac{37}{4000 \left(x + 2\right)} + \frac{17}{400 \left(x + 2\right)^{2}} - \frac{37}{4000 \left(x - 2\right)} + \frac{17}{400 \left(x - 2\right)^{2}}$

and then integrate term-by-term:

[integrate(a, x) for a ∈ Introspection.arguments(us)] # newer SymPy

6-element Vector{Sym{PyCall.PyObject}}:
    -37*log(4000*x - 8000)/4000
                 -6*atan(x)/125
 2*x/(50*x^2 + 50) + atan(x)/25
              -17/(400*x - 800)
              -17/(400*x + 800)
     37*log(4000*x + 8000)/4000

This finds the area bounded by two parabola:

f(x) = 1 - 2x^2
g(x) = x^2
as = solve(f(x) ~ g(x), x)
a,b = sort(as)
integrate(f(x) - g(x), (x, a, b))

$\frac{4 \sqrt{3}}{9}$

Following an example from AP, we look at the solid formed by rotating the region bounded by $y = \sqrt{x+2}$ and $y=e^x$ about the line $y = -2$.

We first find the intersection points, solutions to $f(x) = g(x)$. We will solve for these, and will have success, as SymPy uses some special functions

@syms x
f(x) = exp(x)
g(x) = sqrt(x + 2)
a,b = sort(solve(f(x) ~ g(x), x))

2-element Vector{Sym{PyCall.PyObject}}:
     -2 - LambertW(-2*exp(-4))/2
 -2 - LambertW(-2*exp(-4), -1)/2

The integral requires two radii from the line $y=-2$. We have on $(a,b)$ that $g(x) > f(x) > 0$, so the distance of $g(x)$ from $y=-2$ is greater than that of $f(x)$, hence the ordering below:

r1, r2 = (g(x) - (-2)), (f(x) - (-2))
v = integrate(PI * (r1^2 - r2^2) , (x, a, b))

$- \pi \left(-4 + \frac{4 \sqrt{2} \sqrt{- W\left(- \frac{2}{e^{4}}\right)} \left(-2 - \frac{W\left(- \frac{2}{e^{4}}\right)}{2}\right)}{3} - \frac{4}{e^{\frac{W\left(- \frac{2}{e^{4}}\right)}{2} + 2}} - \frac{1}{2 e^{W\left(- \frac{2}{e^{4}}\right) + 4}} - W\left(- \frac{2}{e^{4}}\right) + \frac{8 \sqrt{2} \sqrt{- W\left(- \frac{2}{e^{4}}\right)}}{3} + \frac{\left(-2 - \frac{W\left(- \frac{2}{e^{4}}\right)}{2}\right)^{2}}{2}\right) + \pi \left(- 4 e^{-2 - \frac{W_{-1}\left(- \frac{2}{e^{4}}\right)}{2}} - 4 - \frac{e^{-4 - W_{-1}\left(- \frac{2}{e^{4}}\right)}}{2} + \frac{\left(-2 - \frac{W_{-1}\left(- \frac{2}{e^{4}}\right)}{2}\right)^{2}}{2} + \frac{4 \sqrt{2} \sqrt{- W_{-1}\left(- \frac{2}{e^{4}}\right)} \left(-2 - \frac{W_{-1}\left(- \frac{2}{e^{4}}\right)}{2}\right)}{3} - W_{-1}\left(- \frac{2}{e^{4}}\right) + \frac{8 \sqrt{2} \sqrt{- W_{-1}\left(- \frac{2}{e^{4}}\right)}}{3}\right)$

This isn’t a satisfying answer if you want to know the scale. We call evalf to give the value:

v.evalf()

$19.7247527844568$

This finds the length of the graph of $x^2$ between $0$ and $1$:

f(x) = x^2
dL = sqrt(1 + diff(f(x),x)^2)
p = integrate(dL, (x, 0, 1))

$\frac{\operatorname{asinh}{\left(2 \right)}}{4} + \frac{\sqrt{5}}{2}$

With numeric value

p.evalf()

$1.4789428575446$

To find the surface area of the volume formed by rotating the graph of $f(x) = \sqrt{9 - x^2}$ between $-1 \leq x \leq 2$ we have:

f(x) = sqrt(9 - x^2)
dSA = 2 * PI * f(x) * sqrt(1 + diff(f(x), x)^2)

$2 \pi \sqrt{9 - x^{2}} \sqrt{\frac{x^{2}}{9 - x^{2}} + 1}$

SymPy can’t see how to do this, we need to help. By squaring, we can see terms will cancel:

dSA = sqrt(cancel(dSA^2))

$6 \pi$

And

integrate(dSA, (x, -1, 2))

$18 \pi$

11.6.4 Algebraic manipulation

The integration of rational functions can always be done, as is known since the early days of calculus. The ability to compute a partial fraction decomposition of any rational function (of a single variable in this particular case), allows this integration as each possible term in this decomposition permits an antiderivative.

SymPy provides some tools to algebraically manipulate expressions, with some addressing rational expressions. Common algebraic manipulations include simpliy/expand along with various specialized functions primarily involving polynomial expressions: factor, cancel, apart, and together.

Symbolic expressions are routinely simplified. However, as simplification can be expensive, only light simplification is done by default. For example, x^3 * x^2 is reduced to x^5; x/3 is reduced to (1/3) * x. Further manipulations are possible. The simplify function is an interface that iterates through dozens of specific simplification steps to completion.

In this example, the gamma function is $\Gamma(n) = (n-1)!$ for integer $n > 0$. We see by simplifying, the cancellation occurs:

@syms n::integer
u = gamma(n) / gamma(n-1)
u, simplify(u)

(gamma(n)/gamma(n - 1), n - 1)

The expand function takes an expression and basically multiplies it out. If new default cancellations occur, expand can actually result in shorter expressions:

expand((x-1)^3)

$x^{3} - 3 x^{2} + 3 x - 1$

ex = (x + 1) * (x - 2) - (x - 1) * x
ex, expand(ex)

(-x*(x - 1) + (x - 2)*(x + 1), -2)

The factor function applied to polynomials is the opposite of expand. This function uses a multivariate factorization algorithm over the rational numbers. This point means not everything factorable is factored:

factor(x^2 - 4), factor(x^2 - 3)

((x - 2)*(x + 2), x^2 - 3)

Whereas, the latter factoring could be achieved through

prod(x - u for u in solve(x^2 - 3))

$\left(x - \sqrt{3}\right) \left(x + \sqrt{3}\right)$

The collect function of sympy collects common powers of a term in an expression. The collect function of Julia takes an iterable and returns a vector or array of values. As these are different concepts, the collect function must be called through sympy.collect.

ex = x*y + x - 3 + 2*x^2 - z*x^2 + x^3
sympy.collect(ex, x)

$x^{3} + x^{2} \cdot \left(2 - z\right) + x \left(y + 1\right) - 3$

For rational expressions, the cancel function re-expresses in a canonical form, no common factors and expanded:

ex = (x - 3)*(x^2 + 2*x + 1)/(x^2 + x)  # common factor of x + 1
ex, cancel(ex)

((x - 3)*(x^2 + 2*x + 1)/(x^2 + x), (x^2 - 2*x - 3)/x)

The apart function finds a partial fraction decomposition for a rational expression. Such decompositions allow the integration of rational polynomials, as the result denominators are of a specific type (powers of linear or quadratic terms).

apart(ex)

$x - 2 - \frac{3}{x}$

The together function re-expresses such a decomposition as a rational expression.

11.6.5 Sequences and series

A finite sequence can be generated in a Julian manner using a comprehension. For example,

[1/Sym(i)^2 for i in 1:10]

10-element Vector{Sym{PyCall.PyObject}}:
     1
   1/4
   1/9
  1/16
  1/25
  1/36
  1/49
  1/64
  1/81
 1/100

The sympy.sequence function can also be used to generate sequences, in this case possibly infinite ones. The specification of the values to iterate over follow the (var, start, end) pattern of integration. In the following oo is used for

@syms i::(integer, nonnegative)
sympy.sequence(1/i^2, (i, 1, oo))

$\left[1, \frac{1}{4}, \frac{1}{9}, \frac{1}{16}, \ldots\right]$

There are some methods available for working with sequences which are not pursued here.

The Sum function of SymPy represents symbolic sums, both finite and infinite. It is called through sympy.Sum as sympy.Sum(ex, (i, a, b)) where ex is an expression of symbolic function in the variable i. The sum is from i=a to i=b where there is convention when a is not less than b and both are finite.

@syms i::(integer, nonnegative) r::real u()
s = sympy.Sum(u(i), (i, 0, 5))

$\sum_{i=0}^{5} u{\left(i \right)}$

As seen, s is an unevaluated sum. It has a doit method to carry out the evaluation:

s.doit()

$u{\left(0 \right)} + u{\left(1 \right)} + u{\left(2 \right)} + u{\left(3 \right)} + u{\left(4 \right)} + u{\left(5 \right)}$

The Sum function is aware of some common sums:

@syms n
s1 = sympy.Sum(1/i, (i, 1, n))
s2 = sympy.Sum(1/i^2, (i, 1, n))
s3 = sympy.Sum((-1)^i /i, (i, 1, n))
s1.doit(), s2.doit(), s3.doit()

(harmonic(n), harmonic(n, 2), -(-1)^(n + 1)*lerchphi(exp_polar(I*pi), 1, n + 1) - log(2))

Infinite sums are specified with oo for $\infty$. (Inf as well, but oo is more fun to type.)

s4 = sympy.Sum(1/i^4, (i, 1, oo))
s5 = sympy.Sum(1/r^i, (i, 0, oo))
s4.doit(), s5.doit()(r => 3)

(pi^4/90, 3/2)

The methods is_convergent and is_absolutely_convergent are useful:

s3.is_convergent(), s3.is_absolutely_convergent()

(True, False)

Suppose $s_n = (n^2 + 4n)\cdot e^{-2n}$ is $\Sigma s_n$ convergent?

sₙ(n) = (n^2 + 4n) * exp(-2n)
@syms n::(integer, nonnegative)
sympy.Sum(sₙ(n), (n, 1, oo)).is_convergent()

$\text{True}$

This is an example on the integral test where $s_n = \ln(n)/n$.

sₙ(n) = log(n)/n
sympy.Sum(sₙ(n), (n, 1, oo)).is_convergent()

$\text{False}$

However,

sₙ(n) = log(n)/n^2
sympy.Sum(sₙ(n), (n, 1, oo)).is_convergent()

$\text{True}$

11.7 Parametric description of functions

Consider a person on a ferris wheel with position coordinates $x(t)$ and $y(t)$. Both of these are cyclical, one side-by-side (e.g., x(t) = r(t)$, one up-and-down $(y(t) = R - r\cos(\omega t))$. A parametric description of the motion, just combines the two into a time-dependent vector

\[ r(t) = (r\sin(\omega\cdot t), R - r\cos(\omega \cdot t)). \]

The function $r$ is vector-valued, not scalar valued. Within Julia, vectors are created using [] with commas separating components. The above might be:

@syms R r t ω
rho(t) = [r*sin(ω*t), R - r*cos(ω*t)]

rho (generic function with 1 method)

We can plot for specific values:

d = (R => 35, r=>30, ω => 2PI/120)
rt = subs.(rho(t), d...) # dots to broadcast and "splat"
plot(rt..., 0, 100)      # dots to splat to two variable

11.7.1 Polar coordinates

Polar coordinates are a special case of a representation. The are an alternate to Cartesian descriptions and use $(r,\theta)$ as an alternate to $(x,y)$. The translation is straightforward: $(x,y) = (r\cdot \cos(\theta), r\cdot\cos(\theta))$.

There isn’t much special support for polar coordinates within base SymPy or Julia. The plot function has a projection=:polar argument that plots polar plots where $r=r(\theta)$.

For example,

rho(theta) = 1 + cos(theta) * sin(theta)^2
@syms θ
plot(rho(θ), 0, 2pi; projection=:polar, legend=false)

11.7.2 Differential equations

A differential equation is an equation involving a variable, a function, and its derivatives. In the following we define an operator, using Differential to simplify the specification of the equation $y'(x) = y(x)$:

@syms x, u()
∂ = Differential(x)
ex = ∂(u(x)) ~ u(x)

$\frac{d}{d x} u{\left(x \right)} = u{\left(x \right)}$

(The use of Differential just cleans up visually what can be achieved through diff(u(x), x) above.)

The dsolve function is then used to solve differential equations:

dsolve(ex)

$u{\left(x \right)} = C_{1} e^{x}$

The constants come from integration. Initial conditions reduce the number of constants. These are specified to the ics keyword argument using a dictionary:

dsolve(ex, ics=Dict(u(0) => 3))

$u{\left(x \right)} = 3 e^{x}$

The initial conditions can also involve constraints on derivatives, etc.

We follow with a much more involved example. The following equations model projectile motion using Newton’s laws and adding a drag proportional to velocity. This proportion is given by a constant $\gamma$. We use u and v to model the $x$ and $y$ coordinates of motion over time, $t$. Inverting $x$ to get $t$, allows the solution of $y$ in terms of $t$.

@syms x0::real y0::real v0::real γ::nonnegative 𝑔::real
@syms t::positive x u() v()
Dₜ = Differential(t)
D2ₜ = Dₜ ∘ Dₜ
eq₁ = Dₜ(Dₜ(u))(t) ~    - γ * Dₜ(u)(t)
eq₂ = Dₜ(Dₜ(v))(t) ~ -𝑔 - γ * Dₜ(v)(t)

$\frac{d^{2}}{d t^{2}} v{\left(t \right)} = - γ \frac{d}{d t} v{\left(t \right)} - 𝑔$

If we set $\gamma = 0$ and solve, we get the $y$ values as a function of $t$:

a1, a2 = dsolve.((eq₁(γ=>0), eq₂(γ=>0)), (u(t), v(t)))
tₓ = only(solve(x - rhs(a1), t)) # invert; only one solution
yₓ = rhs(a2)(t => tₓ)

$x - \frac{𝑔 \left(- C_{1} + x\right)^{2}}{2 C_{2}^{2}}$

We can see the form of a parabola, with two parameters adjusting its shape related to the initial conditions. This is anticipated from physics.

Here, we specify initial values so that we can ultimately plot:

@syms α::real v₀::real
icx = Dict(u(0) => Sym(0), diff(u(t),t)(0) => v₀*cos(α))
icy = Dict(v(0) => Sym(0), diff(v(t),t)(0) => v₀*sin(α) )

a1 = dsolve(eq₁(γ => 0), u(t), ics=icx)
a2 = dsolve(eq₂(γ => 0), v(t), ics=icy)

tₓ = only(solve(x - rhs(a1), t)) # only one solution
yₓ = rhs(a2)(t => tₓ)

$\frac{x \sin{\left(α \right)}}{\cos{\left(α \right)}} - \frac{x^{2} 𝑔}{2 v₀^{2} \cos^{2}{\left(α \right)}}$

If there is a positive $\gamma$, there will be drag and the formulas will change to reflect that.

𝑎1 = dsolve(eq₁, u(t), ics=icx)
𝑎2 = dsolve(eq₂, v(t), ics=icy)

𝑡ₓ = only(solve(x - rhs(𝑎1), t)) # only one solution
𝑦ₓ = rhs(𝑎2)(t => 𝑡ₓ)

$- \frac{𝑔 \log{\left(\frac{v₀ \cos{\left(α \right)}}{v₀ \cos{\left(α \right)} - x γ} \right)}}{γ^{2}} + \frac{v₀ γ \sin{\left(α \right)} + 𝑔}{γ^{2}} + \frac{\left(v₀ \cos{\left(α \right)} - x γ\right) \left(- v₀ γ \sin{\left(α \right)} - 𝑔\right)}{v₀ γ^{2} \cos{\left(α \right)}}$

The difference can be visualized through plotting. The proper way to get the values for the constants is to specify numeric values in the initial conditions:

a = yₓ(𝑔 => 32, v₀ => 200, α => PI/4)
M = maximum(solve(a ~ 0))
plot(a, 0, N(M))

𝑎 = 𝑦ₓ(𝑔 => 32, v₀ => 200, α => PI/4, γ => 1//2)
M = nsolve(𝑎 ~ 0, 2*100*sqrt(2)-10) # avoid starting at asymptote
plot!(𝑎, 0, N(M))

Other plotting packages can be used, of course, but there are a few recipes to easily plot SymPy expressions↩︎
The sympify function takes an expression and makes it symbolic, which can be used to convert numbers to symbolic numbers or strings into variables.↩︎
The alert reader may notice a different type face in the output below; one used for simpler symbolic expressions.↩︎
We can’t use =, the mathematical indication of an equations, as that is used for assignment within Julia, as it is with most programming languages.↩︎
There is a subtlety, pi is promoted to the symbolic value of PI before checking, as the floating point value associated to $\pi$, 3.141592653589793, is not in the set.↩︎
To run this example requires a newer version of SymPy.↩︎