Python/SymPy primer

0 Setup

The cell below sets up various things. It is included to ensure that this notebook is self-contained. For your own notebooks, it will be better to use from header import * as explained in Section 1 below.

import numpy as np
import sympy as sp
import scipy as sc
import plotly
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import plotly.io as pio
from IPython.display import display, HTML, Math, Latex

def equal_aspect(fig):
    fig.update_yaxes(scaleanchor="x",scaleratio=1)

def use_pi_ticks(fig, k0=-1, k1=1, d = 1, x0=None, x1=None):
    if x0 is not None:
        k0 = sp.Rational(int(np.ceil(x0*d/np.pi)),d)
    if x1 is not None:
        k1 = sp.Rational(int(np.floor(x1*d/np.pi)),d)
    ticks = [sp.Rational(i,d) * sp.pi for i in range(int(k0*d),int(k1*d)+1)]
    tickvals = [float(tick) for tick in ticks]
    ticktext = [f"${sp.latex(tick)}$" for tick in ticks]
    fig.update_xaxes(tickvals=tickvals, ticktext=ticktext)


pio.templates['maths'] = go.layout.Template(
 layout = go.Layout(
  margin=dict(l=0, r=0, t=20, b=20),
  xaxis = dict(zeroline=True, zerolinewidth=1, zerolinecolor='black'),
  yaxis = dict(zeroline=True, zerolinewidth=1, zerolinecolor='black')
 )
)
pio.templates.default = 'plotly+maths'
plotly.offline.init_notebook_mode()
display(HTML(
    '<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_SVG"></script>'
))

1. Getting started

These instructions will assume that you are comfortable with using Python in Jupyter notebooks in VS Code, and that you have installed the packages sympy, numpy, scipy and plotly using pip or conda. You will also need to import these packages and set some options. The most convenient way to do this is as follows. Distributed together with this file is a file called __init__.py.txt. You should rename this to __init__.py and put it in a subfolder of the current folder called header. Then in each of your notebooks you should include the line

 from header import *

at the top.

2. Notation for variables

2.1

By default names like $x$ refer to variables as used in programming languages like Python. They are expected to have a value (which might be a number or a list or something more complicated). If you try to use a variable (for example, by trying to expand $(1+x)^4$) without assigning a value then you will get an error. If you instead want to treat $x$ as an abstract symbol then you need to enter x = sp.symbols('x'). You can do this for several variables at the same time using syntax like x,y,a,b = sp.symbols('x y a b').

The line below will generate an error, because $u$ has not been defined. In VS Code, we get a warning about this: the letter u is shown with a small yellow squiggle underneath. Even if we want to treat $u$ as an abstract symbol with no specific value, we need to tell Python that by entering u = sp.symbols('u').

sp.expand((1 + u) ** 4)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 1
----> 1 sp.expand((1 + u) ** 4)

NameError: name 'u' is not defined

After specifying that $x$ is an abstract symbol, we can expand $(1 + x)^4$. (Note, however, that to raise something to the fourth power we need to enter ** 4, not ^ 4 as in LaTeX or in some other programming languages.)

x = sp.symbols('x')
sp.expand((1 + x) ** 4)

$\displaystyle x^{4} + 4 x^{3} + 6 x^{2} + 4 x + 1$

2.2

If you want variables with indices like $m_0,m_1,m_2$ then you should enter m = sp.IndexedBase('m'), and then enter m[2] for $m_2$, for example.

m = sp.IndexedBase('m')
m[0] + m[1] + m[2]

$\displaystyle {m}_{0} + {m}_{1} + {m}_{2}$

2.3

Most variables will have single-letter names, such as $x$ or $A$. Note that case is significant: Python regards $A$ as different from $a$, and $x$ as different from $X$.

It is perfectly permissible to use long names for variables, such as eqns or sols or SpeedOfLight.

However, you may not use words that already have a special meaning in Python. In particular, you cannot call a variable lambda.

SpeedOfLight = 299792458
print(f"The speed of light is {SpeedOfLight} m/s")

The speed of light is 299792458 m/s

2.4

You can use Greek-letter variables by spelling out their names in full:

alpha = sp.symbols('alpha')
beta = sp.symbols('beta')
display((alpha + beta) / (alpha - beta))

$\displaystyle \frac{\alpha + \beta}{\alpha - \beta}$

Alternatively, there are various ways in which you can set up your computer so that you can enter Greek letters directly (one option is Wincompose). Then you can just use Greek letters like normal English letters.

μ, ν = sp.symbols('μ ν')
sp.expand((μ + ν) * (μ - ν))

$\displaystyle μ^{2} - ν^{2}$

3. Generalities

3.1

The help() command will give you details about any Python or SymPy object that you are interested in. Alternatively, in VS Code you can just enter a name and hover over it with your mouse and some information will appear. However, the output from help() is sometimes very long and overwhelming. It may be better to ask an AI a more specific question about what you are hoping to do.

help(sp.pi)

3.2

See Section 15 for advice on some common problems.

3.3

This performs a calculation and displays the result.

111 * 10000001

1110000111

This performs two calculations, but by default, Python only displays the result of the last one

1111 + 9999
1111 ** 2

1234321

To display the results of both calculations, we need to add explicit print() or display() commands. Here the results are just integers, so there is no difference between print() and display(). For complicated mathematical expressions, print() will give Python notation which you could type if you wanted to enter the expression by hand, whereas display() will give something more similar to traditional mathematical notation.

print(1111 + 9999)
display(1111 ** 2)
x  = (sp.pi ** 2 + 1) / (sp.pi ** 2 - 1)
print(x)
display(x) 
del(x) # Remove the definition of x to avoid problems later

11110

1234321

(1 + pi**2)/(-1 + pi**2)

$\displaystyle \frac{1 + \pi^{2}}{-1 + \pi^{2}}$

3.4

Here the semicolon suppresses the display of the number 55, so no output of any kind is shown.

44
55;

Adding None as the last expression in the cell has the same effect.

44
55
None

3.5

You can use the symbol _ to refer to the last thing that Python calculated and displayed, and __ to refer to the one before that, and so on. However, this mechanism cannot be used to refer to things that were calculated but not displayed.

If you move around a worksheet inserting things in different places, you should remember that _ refers to the most recent result, which may not be the one that you see directly above the line where you are typing. Because this can lead to confusion and error, it is not recommended to use this mechanism in any notebook that you intend to save for future use.

All this is illustrated by the lines below.

1111 * 1 # A simple calculation entered as an example

1111

1111 * 2 # A simple calculation entered as an example

2222

1111 * 3 # A simple calculation entered as an example

3333

1111 * 4 # A simple calculation entered as an example

4444

If we execute all the cells in order, the the symbol __ in the cell below will refer to the result 3333 that appeared two cells above this one, so the result printed will be 33330000.

__ * 10000

33330000

3.6

It is important to understand the difference between = and == and sp.Eq().

The line x = 3 is an assignment statement. It assigns the value 3 to the variable x. Thus, if we then enter print(100 * x) then x gets replaced by 3 so 100 * x becomes 300, and that is what is printed.

x = 3
print(100 * x)

300

The line x == 3 is a comparison statement. It compares the value of x with 3 and returns True if they are equal, and False otherwise. In the previous cell we set x equal to 3, so this comparison will return True. The line x == 33 is another comparison statement, but this one will return False, because x has been set equal to 3, which is not equal to 33.

print(x == 3)  
print(x == 33)  

True
False

Assuming that a and b have been declared to be abstract symbols, the expression sp.Eq(2 * a, b + 1) represents the hypothetical question "Is $2 a$ equal to $b + 1$?", which might be true or false depending on the values of $a$ and $b$. Entering this question does not assign any values to a or b, and does not return a True or False answer. The main use for sp.Eq() is to introduce equations which you later intend to solve; this will be discussed in Section 8. By contrast, the expression 2 * a == b + 1 gives False, because the expression 2 * a is not visibly the same as b + 1.

a, b = sp.symbols('a b')
display(sp.Eq(2 * a , b + 1))
display(2 * a == b + 1)

$\displaystyle 2 a = b + 1$

False

In fact, even (a + 1) * (a - 1) == a ** 2 - 1 will return False, because the two sides of the equation are not visibly the same, even they become the same after expanding out the left hand side. (The right way to test whether $A$ and $B$ are mathematically equivalent is to enter sp.simplify(A-B) and see whether you get zero.)

A = (a + 1) * (a - 1)
B = a**2 - 1
display(A == B) 
display(sp.simplify(A - B))

False

$\displaystyle 0$

3.7

There are various different ways of defining and thinking about functions in Python/SymPy. Here we will just note that the line below gives an error message, and is not the right way to define the function $f(t) = t^2+3$. Correct approaches will be discussed in Section 10.

f(t) = t ** 2 + 3

  Cell In[21], line 1
    f(t) = t ** 2 + 3
    ^
SyntaxError: cannot assign to function call here. Maybe you meant '==' instead of '='?

3.8

Near the top of your VS Code window you should see a menu item marked Restart. If you click this then Python will restart and all values that you might have assigned to variables will be removed. You might want to do this at the beginning of each new exercise, otherwise stray values left over from previous exercises can cause trouble and confusion. Note, however, that most notebooks will have a cell at the top with some import statements. Restarting removes all imports, so you will need to go back to your import cell and execute it again before continuing.

3.9

Here we define values for x, y and z, and then use the del command to remove the values from x and y. Even though we defined z in terms of x and y, and the values of x and y have been removed, the variable z still retains the value that was previously assigned.

x = 1234
y = 5678
z = x + y
del x, y
z

6912

Because x is now undefined, the line below gives an error.

x + 42

4. Numerical approximation

4.1

The expression sp.pi refers to the exact value of the constant π. We can use sp.N() to convert this to an approximate numerical value.

display(sp.pi)
display(sp.N(sp.pi))

$\displaystyle \pi$

$\displaystyle 3.14159265358979$

If A is some kind of algebraic expression, then sp.N(A) will convert any numbers involved in A to floating point form, and leave the symbolic parts unchanged, as illustrated below.

x = sp.symbols('x')
display(x / 3)
display(sp.N(x / 3))

$\displaystyle \frac{x}{3}$

$\displaystyle 0.333333333333333 x$

4.2

By default, sp.N() gives something like 15 significant figures (the precise story is complicated). To get 50 significant figures, we can enter sp.N(...,50).

x = sp.pi ** 2 / 6
display(x)
display(sp.N(x, 50)) 
del(x) # Remove the definition of x to avoid problems later

$\displaystyle \frac{\pi^{2}}{6}$

$\displaystyle 1.6449340668482264364724151666460251892189499012068$

4.3

It is possible to lose accuracy for non-obvious reasons. For example, if you enter sp.N(1/7, 100) then Python will work from the inside out, first calculating $1/7$ to low accuracy and only then calling sp.N() with a request for high accuracy that can no longer be achieved. One way to fix this is to enter sp.N(1,100)/7, which requests high accuracy before attempting to divide by $7$.

x = sp.N(1/7, 100)    # This gives a 100 digit answer, but only the first 16 digits are correct.
y = sp.N(1, 100) / 7  # This gives a 100 digit answer, and all 100 digits are correct.
print(7 * x - 1)      # Because x has some incorrect digits, this gives a non-zero answer, which is not what we expect.
print(7 * y - 1)      # Because y has all correct digits, this gives zero, which is what we expect.

-0.00000000000000005551115123125782702118158340454101562500000000000000000000000000000000000000000000000000000000000000
0

4.4

In some contexts, SymPy will convert fractions to approximate decimals without being asked, which is not what we want when doing symbolic mathematics. To prevent $2/3$ from being converted to 0.66666666667, we need to enter it as sp.Rational(2, 3) instead of 2/3. Sometimes this is not necessary, but the rules for whether it is necessary are complicated, so it is best to err on the safe side.

display(x + 2/3)
display(x + sp.Rational(2, 3))

$\displaystyle x + 0.666666666666667$

$\displaystyle x + \frac{2}{3}$

5. Algebraic operations

For the sake of example, we define $x$, $y$, $n$, $a$, $b$ and $c$ to be abstract symbols.

x, y, n, a, b, c = sp.symbols('x y n a b c')

5.1

Use a * for multiplication; for example, $2xy$ should be entered as 2*x*y, and $(2x+1)\sin(3x)$ should be entered as (2*x+1)*sp.sin(3*x).

display(2*x*y)
display((2*x + 1)*sp.sin(3*x))

$\displaystyle 2 x y$

$\displaystyle \left(2 x + 1\right) \sin{\left(3 x \right)}$

5.2

It is not correct to write xy for the product of x and y. Instead, Python treats xy as a new variable with a two letter name. If xy has not already been assigned some kind of value, this will give an error. If you enter xy in VS Code, you will usually see that it is underlined with a squiggly yellow line to indicate that it is undefined.

print(xy)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[55], line 1
----> 1 print(xy)

NameError: name 'xy' is not defined

5.3

Brackets need to be used for algebraic grouping, just as in ordinary mathematical writing. For the product of $x + 2$ and $x - 2$ we need to enter (x + 2) * (x - 3), not x + 2 * x - 2. Brackets used for algebraic grouping should always be round brackets, not square or curly ones.

display((x + 2) * (x - 2))
display( x + 2  *  x - 2)

$\displaystyle \left(x - 2\right) \left(x + 2\right)$

$\displaystyle 3 x - 2$

5.4

Enter powers using **, for example x ** 6 for $x^6$.

x ** 6

$\displaystyle x^{6}$

The notation x ^ 6 does have a meaning (the exclusive or operator applied to x and 6) but it is almost never what you want; this can unfortunately lead to some confusing error messages.

x ^ 6

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[57], line 1
----> 1 x ^ 6

TypeError: unsupported operand type(s) for ^: 'Symbol' and 'int'

5.5

Bracketing is also important when raising things to powers. To raise $x + 2$ to the power $n+5$, we need to enter (x + 2) ** (n + 5), not x + 2 ** n + 5.

display((x + 2) ** (n + 5)) # This is what we want
display( x + 2  **  n + 5 ) # This is meaningful but different

$\displaystyle \left(x + 2\right)^{n + 5}$

$\displaystyle 2^{n} + x + 5$

The cube root of $x$ is $x^{1/3}$. For this to work properly in SymPy, it is important to enter $1/3$ as sp.Rational(1,3) (which gives $1/3$ exactly) and not $1/3$ (which gives 0.33333).

display(125 ** sp.Rational(1, 3)) # This gives an answer of exactly 5
display(125 ** (1/3))             # This only gives an approximate answer

$\displaystyle 5$

4.999999999999999

5.6

The (nonnegative) square root of $x$ can be entered as sp.sqrt(x) or x ** sp.Rational(1,2). If you enter x ** (1/2) you will only get an approximate answer.

display(sp.sqrt(12))              # exact
display(12 ** sp.Rational(1, 2))  # exact
display(12 ** (1/2))              # approximate

$\displaystyle 2 \sqrt{3}$

$\displaystyle 2 \sqrt{3}$

3.4641016151377544

sp.sqrt(b**2 - 4 * a * c)

$\displaystyle \sqrt{- 4 a c + b^{2}}$

6. Algebraic manipulation

Suppose we have a complicated expression called $A$, which we want to simplify or otherwise manipulate.

x, y, n, a, b, c, d, p, q = sp.symbols('x y n a b c d p q')

6.1

To substitute for some or all of the variables in $A$, use the subs command. For example, if $A$ is $(a+b)^n$, you can change the $b$ to $1$ and the $n$ to $p/q$ by entering A.subs({b:1,n:p/q}), which gives $(a+1)^{p/q}$. (Note that the correct syntax is b:1, not b=1.)

a, b, n, p, q = sp.symbols('a b n p q')
A = (a + b) ** n
display(A)
display(A.subs({b : 1, n : p/q}))

$\displaystyle \left(a + b\right)^{n}$

$\displaystyle \left(a + 1\right)^{\frac{p}{q}}$

6.2

Often the equations to be used in subs(...) will come from the solve(...) command. For example, to find the value of $x^2+y^2$ at the point where $x+y=6$ and $x-y=2$, enter sol=sp.solve([sp.Eq(x+y,6),sp.Eq(x-y,2)]) and then (x**y + y**2).subs(sol).

x, y = sp.symbols('x y')
sol=sp.solve([sp.Eq(x+y,6),sp.Eq(x-y,2)])
display(sol)
display((x**y + y**2).subs(sol))

{x: 4, y: 2}

$\displaystyle 20$

6.3

Enter sp.expand(A) to expand everything out. For example, enter sp.expand((x+y)**3) to expand $(x+y)^3$, giving $x^3+3x^2y+3xy^2+y^3$.

sp.expand((x+y) ** 3)

$\displaystyle x^{3} + 3 x^{2} y + 3 x y^{2} + y^{3}$

x = sp.symbols('x', positive=True)
(x ** 4) ** sp.Rational(1, 2)

$\displaystyle x^{2}$

6.4

To simplify $A$, try entering sp.simplify(A). If this does not do what you want, you can try sp.factor(A). If $A$ involves trigonometric functions, try sp.trigsimp(A) or sp.expand_trig(A).

A = (x + y) ** 2 - (x - y) ** 2
display(A)
display(sp.simplify(A))

$\displaystyle - \left(x - y\right)^{2} + \left(x + y\right)^{2}$

$\displaystyle 4 x y$

A = (x ** 8 - y ** 8) / (x ** 2 - y ** 2)
display(A)
display(sp.factor(A))

$\displaystyle \frac{x^{8} - y^{8}}{x^{2} - y^{2}}$

$\displaystyle \left(x^{2} + y^{2}\right) \left(x^{4} + y^{4}\right)$

A = sp.sin(2 * x) * sp.cos(3 * x)
B = sp.expand(sp.expand_trig(A))
C = sp.trigsimp(B)
display(A)
display(B)
display(C)

$\displaystyle \sin{\left(2 x \right)} \cos{\left(3 x \right)}$

$\displaystyle 8 \sin{\left(x \right)} \cos^{4}{\left(x \right)} - 6 \sin{\left(x \right)} \cos^{2}{\left(x \right)}$

$\displaystyle - \frac{\sin{\left(x \right)}}{2} + \frac{\sin{\left(5 x \right)}}{2}$

6.5

Note that rules like $(a^4)^{1/4}=a$ are not valid when $a$ is negative. The function sp.simplify() will not simplify expressions like this by default. However, if you declare $a$ to be positive using a = sp.symbols('a', positive=True) (instead of just a = sp.symbols('a')) then SymPy will apply rules that are only valid for positive $a$. For example, (a ** 4) ** sp.Rational(1,4) will simplify to $a$. There are also more complex ways to specify that $a$ is an integer or lies in a certain range or has other useful properties that could be used when simplifying; see the SymPy documentation about assumptions.

a = sp.symbols('a')
sp.simplify((a ** 4) ** sp.Rational(1,4))

$\displaystyle \sqrt[4]{a^{4}}$

a = sp.symbols('a', positive=True)
sp.simplify((a ** 4) ** sp.Rational(1,4))

$\displaystyle a$

6.6

If you suspect that $A$ is equal to some other expression (say $B$) and want to check this, the best way is to enter sp.simplify(A-B) and see if you get zero. For example, the cell below checks that the expression $A=\frac{1}{1+\frac{1}{x+1}}$ is the same as $B=1-\frac{1}{x+2}$.

A = 1/(1 + 1/(1 + x))
B = 1 - 1/(x + 2)
display(A)
display(B)
display(sp.simplify(A - B))

$\displaystyle \frac{1}{1 + \frac{1}{x + 1}}$

$\displaystyle 1 - \frac{1}{x + 2}$

$\displaystyle 0$

6.7

To collect together all the terms in $A$ involving the same power of $x$, enter A.collect(x). For example, this converts $1+ax+bx+cx^2+dx^2$ to $1+(a+b)x+(c+d)x^2$.

A = 1+a*x+b*x+c*x**2+d*x**2
display(A)
display(sp.collect(A, x))

$\displaystyle a x + b x + c x^{2} + d x^{2} + 1$

$\displaystyle x^{2} \left(c + d\right) + x \left(a + b\right) + 1$

6.8

If you just want the coefficient of $x^2$ in $A$ then you can instead enter A.coeff(x,2) or A.coeff(x**2). To get the constant term, you must enter A.coeff(x,0) - the syntax A.coeff(x**0) will not work.

A = 1 + 10*x + 111*x**2 + 1000*x**3
print(A.coeff(x, 2))
print(A.coeff(x**3))

111
1000

7. Constants

7.1

If you enter a fraction like 2/9 then it will immediately be converted to an approximate decimal, which is not usually what you want. If you want to treat it as an exact fraction you should enter sp.Rational(2,9) instead.

display(2/9)
display(sp.Rational(2,9))

0.2222222222222222

$\displaystyle \frac{2}{9}$

7.2

The constant $\pi\simeq 3.1415926536$ should be entered as sp.pi. If you just enter pi then that will just be treated as the name of a variable with no special value or properties. Thus, if you have not introduced pi in any way then you will get an error.

display(sp.pi)
display(sp.N(sp.pi, 20))
display(pi) # Gives an error

$\displaystyle \pi$

$\displaystyle 3.1415926535897932385$

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[79], line 3
      1 display(sp.pi)
      2 display(sp.N(sp.pi, 20))
----> 3 display(pi) # Gives an error

NameError: name 'pi' is not defined

7.3

The constant $e\simeq 2.718281828$ must be entered as sp.E (with a capital E) or sp.exp(1).

sp.E, sp.exp(1), sp.N(sp.E)

(E, E, 2.71828182845905)

7.4

The square root of minus one can be entered as sp.I.

display(sp.I)
display(sp.I ** 2)

$\displaystyle i$

$\displaystyle -1$

7.5

In many contexts $\infty$ can be entered as sp.oo. For example, the cell below calculates the limit of $(1+2x)/(3+4x)$ as $x\to\infty$.

x = sp.symbols('x')
sp.limit((1 + 2*x)/(3 + 4*x), x, sp.oo)

$\displaystyle \frac{1}{2}$

8. Solving

8.1

To solve the equation $x^2+6=5x$ (for example), enter sp.solve(sp.Eq(x**2+6,5*x)). The answer is $[2,3]$, indicating that the solutions are $x=2$ and $x=3$.

x = sp.symbols('x')
sp.solve(sp.Eq(x**2+6,5*x))

[2, 3]

If you want to solve for an expression being equal to zero then you do not need sp.Eq(), you can just enter sp.solve(x**2+6-5*x) to solve $x^2+6-5x=0$ for example. Similarly, if you have already defined a function $g$, you can enter solve(g(t),t) to solve the equation $g(t)=0$. It does not work to enter a == b instead of sp.Eq(a, b).

x = sp.symbols('x')
sp.solve(x**2+6-5*x)

[2, 3]

t = sp.symbols('t')

def g(t):
    return (t + 1) ** 2 - (2 * t + 1) ** 2

sp.solve(g(t) , t)

[-2/3, 0]

8.2

The syntax discussed above gives the solution to $x^2+6-5x=0$ in the form [2,3], meaning that $x=2$ or $x=3$. If you prefer to have the answer in the form [{x:2},{x:3}] (which can be more convenient if you want to substitute the result into some other expression), you should instead enter sp.solve({x^2-5*x+6},dict=True).

sp.solve({x**2-5*x+6}, dict=True)

[{x: 2}, {x: 3}]

8.3

The approach discussed above will give complex roots as well as real ones. For example, if we declare x = sp.symbols('x') and then enter sp.solve(x**3+x,x) we get $[0,-i,i]$. If we only want real roots we can either use sp.real_roots(x**3+x,x) or we can declare $x$ using x = sp.symbols('x',real=True) and then just use sp.solve(); this will give $x=0$ as the only real root.

x = sp.symbols('x')
sp.solve(x**3 + x, x) # Includes the complex solutions x = I and x = -I as well as x = 0.
sp.real_roots(x**3 + x, x) # Only includes the real solution x = 0.

[0]

x = sp.symbols('x', real=True)
sp.solve(x**3 + x, x) # Only includes the real solution x = 0.

[0]

8.4

To solve several equations in several variables, put the equations in one set of square brackets, and the variables in another set. For example, to solve $2x + 3y = 5$ and $3x + 4y = 7$ for $x$ and $y$, enter the following:

sp.solve([2*x+3*y-5,3*x+4*y-7],[x,y])

{x: 1, y: 1}

8.5

If there are no solutions, sp.solve() will return the empty list []. For example, this happens if we try to solve the simultaneous equations $x+y=1$ and $2x+2y=10$.

sp.solve([sp.Eq(x+y,1), sp.Eq(2*x+2*y,10)],[x,y])

[]

8.6

If the problem involves trigonometric functions, then there will typically be infinitely many solutions, but sp.solve() will only list a few of them. To get all solutions, you can try using sp.solveset() instead. For example, we can solve the equation $\cos(x)=1$ by entering sp.solveset(sp.cos(x)-1,x). The result will be printed as $\{2n\pi\;\mid\;n\in\mathbb{Z}\}$, indicating that $\cos(x)=1$ if and only if $x$ is an integer multiple of $2\pi$.

display(sp.solve(sp.Eq(sp.cos(x),1), x))     # Gives only the solutions x = 0 and x = 2*pi
display(sp.solveset(sp.Eq(sp.cos(x),1), x))  # Gives all the solutions x = 2*n*pi for any integer n

[0, 2*pi]

$\displaystyle \left\{2 n \pi\; \middle|\; n \in \mathbb{Z}\right\}$

There is quite a long story about how the above result is represented internally in Python; see the SymPy documentation on sets.

Here, as another example, we look for solutions to $\sin(x)=\cos(x)$. The answer says that $x$ either has the form $2n\pi+5\pi/4$ for some integer $n$, or has the form $2n\pi+\pi/4$ for some integer $n$. This is an example where SymPy's answer is correct but not as tidy as it could be: it would be equivalent to say that $x=(m+1/4)\pi$ for some integer $m$.

sp.solveset(sp.Eq(sp.sin(x),sp.cos(x)), x)

$\displaystyle \left\{2 n \pi + \frac{5 \pi}{4}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{\pi}{4}\; \middle|\; n \in \mathbb{Z}\right\}$

8.7

To get approximate numerical solutions, you can use sp.nsolve() instead of sp.solve(). For example, you can enter sp.nsolve(sp.cos(x)-sp.cos(x+1),x,42) to find a number $x$ with $\cos(x)=\cos(x+1)$, starting the search at $x=42$. When using sp.nsolve(), you always need to specify the starting point. If you need to solve a large number of equations numerically then it will usually be best to use the scipy or numpy packages; you should only use sp.nsolve() for occasional numerical solutions when you are mostly working with symbolic equations.

x0 = sp.nsolve(sp.cos(x)-sp.cos(x+1),x,42)
display(x0)
display(sp.cos(x0))
display(sp.cos(x0 + 1))

$\displaystyle 46.6238898038469$

$\displaystyle -0.877582561890374$

$\displaystyle -0.877582561890372$

8.8

Sometimes we need to solve problems like this: find $a$ and $b$ such that $a(x-1)^2+b(x+1)^2=x$ for all $x$. For this we can use sp.solve_undetermined_coeffs() as illustrated below. (The answer is $a=-1/4$ and $b=1/4$, indicating that $\frac{-1}{4}(x-1)^2+\frac{1}{4}(x+1)^2=x$, which you can easily check by expanding it out.)

a, b, x = sp.symbols('a b x')
sp.solve_undetermined_coeffs(a*(x-1)**2+b*(x+1)**2-x,[a,b],x)

{a: -1/4, b: 1/4}

Similarly, to find $u$ such that $\sin(t+u)=\cos(t)$ for all $t$, we could try

 sp.solve_undetermined_coeffs(sp.Eq(sp.sin(t+u),sp.cos(t)),[u],t)

This does not actually work, but we can make it work by giving the solver a hint as follows:

 sp.solve_undetermined_coeffs(sp.Eq(sp.expand_trig(sp.sin(t+u)),sp.cos(t)),[u],t)

t, u = sp.symbols('t u')
display(sp.solve_undetermined_coeffs(sp.Eq(sp.sin(t+u),sp.cos(t)),[u],t))
display(sp.solve_undetermined_coeffs(sp.Eq(sp.expand_trig(sp.sin(t+u)),sp.cos(t)),[u],t))

None

{u: pi/2}

9. Standard functions

9.1

Most functions can be entered using their usual names prefixed with sp., for example sp.sin(x) or sp.ln(y). If you do not want to type the prefix, you can change your initial import statement to include

  from sympy import sin, cos, tan, exp, ln

(for example) then you can just enter sin(x) instead of sp.sin(x).

sp.sin(sp.pi / 6)

$\displaystyle \frac{1}{2}$

9.2

However, you must always use brackets: enter sp.tan(x) instead of tan x, and sp.sin(2*x) instead of sin 2x.

x = sp.symbols('x')
tan x  # Gives an error

  Cell In[99], line 2
    tan x  # Gives an error
        ^
SyntaxError: invalid syntax

9.3

The absolute value of $x$, traditionally written as $|x|$, must be entered as abs(x). (This is an exception to the usual pattern: sp.abs(x) will not work.)

x = sp.symbols('x')
display(abs(-1000))
display(abs( 1000))
display(abs(-x))

1000

1000

$\displaystyle \left|{x}\right|$

9.4

Both sp.ln(x) and sp.log(x) refer to the natural logarithm (to base $e$) of $x$.

display(sp.ln(sp.E ** 99))
display(sp.log(sp.E ** 99))
display(sp.N(sp.log(10)))

$\displaystyle 99$

$\displaystyle 99$

$\displaystyle 2.30258509299405$

9.5

If you want to use the logarithm to base $10$ then you should enter sp.log(x, 10). Similarly, enter sp.log(x, 2) for $\log_2(x)$ (the logarithm to base $2$) and so on.

display(sp.log(1000, 10)) # 1000 = 10 ** 3 so log_10(1000) = 3
display(sp.log(4096,  2)) # 4096 = 2 ** 12 so log_2(4096) = 12

$\displaystyle 3$

$\displaystyle 12$

9.6

The function $e^x$ should be entered as sp.exp(x) or sp.E ** x

display(sp.exp(11))
display(sp.E ** 11)
display(sp.ln(sp.E ** 11))

$\displaystyle e^{11}$

$\displaystyle e^{11}$

$\displaystyle 11$

9.7

The hyperbolic function $\sinh(x)=(e^x-e^{-x})/2$ can be entered as sp.sinh(x), and similarly for the functions $\cosh(x)=(e^x+e^{-x})/2$ and $\tanh(x)=\sinh(x)/\cosh(x)$.

x = sp.symbols('x')
sp.exp(x), sp.sinh(x), sp.cosh(x)

(exp(x), sinh(x), cosh(x))

If we have an expression $A$ which involves hyperbolic functions, we can rewrite it in terms of exponentials using A.rewrite(sp.exp). Below we check the identity $\cosh(x)^2-\sinh(x)^2=1$ by rewriting everything in terms of $e^x$ and expanding it out.

A = sp.cosh(x) ** 2 - sp.sinh(x) ** 2
display(A)
display(A.rewrite(sp.exp))
display(sp.expand(A.rewrite(sp.exp)))

$\displaystyle - \sinh^{2}{\left(x \right)} + \cosh^{2}{\left(x \right)}$

$\displaystyle - \left(\frac{e^{x}}{2} - \frac{e^{- x}}{2}\right)^{2} + \left(\frac{e^{x}}{2} + \frac{e^{- x}}{2}\right)^{2}$

$\displaystyle 1$

9.8

The traditional notation $\sin^2(x)$ refers to the square of $\sin(x)$, which could also be written $(\sin(x))^2$ or $\sin(x)^2$. In Python this must be entered as sp.sin(x)**2, not sp.sin**2(x). Similar comments apply to $\tan^2(x)$ and so on. Note that $\sin(x^2)$ means something different again: it is the sin of the square of $x$, not the square of the sin.

x = sp.symbols('x')
A = sp.sin(x) ** 2 + sp.cos(x) ** 2
display(A)
display(sp.simplify(A))

$\displaystyle \sin^{2}{\left(x \right)} + \cos^{2}{\left(x \right)}$

$\displaystyle 1$

sp.sin ** 2(x) # This is not syntactically correct, and gives an error.

<>:1: SyntaxWarning:

'int' object is not callable; perhaps you missed a comma?

C:\Users\Neil\AppData\Local\Temp\ipykernel_14536\3669833667.py:1: SyntaxWarning:

'int' object is not callable; perhaps you missed a comma?

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[116], line 1
----> 1 sp.sin ** 2(x) # This is not syntactically correct, and gives an error.

TypeError: 'int' object is not callable

9.9

In traditional notation $1/\sin(x)$ is sometimes written as $\text{csc}(x)$ or $\text{cosec}(x)$. In Python sp.csc(x) will work, but sp.cosec(x) will not. Of course, you can also just enter 1/sp.sin(x). Similarly, $1/\cos(x)$ is sp.sec(x).

x = sp.symbols('x')
display(sp.sec(x))
display(sp.sec(x).rewrite(sp.cos))
display(sp.simplify(sp.sec(x) * sp.cos(x)))

$\displaystyle \sec{\left(x \right)}$

$\displaystyle \frac{1}{\cos{\left(x \right)}}$

$\displaystyle 1$

x = sp.symbols('x')
display(sp.csc(x))
display(sp.csc(x).rewrite(sp.sin))
display(sp.simplify(sp.csc(x) * sp.sin(x)))

$\displaystyle \csc{\left(x \right)}$

$\displaystyle \frac{1}{\sin{\left(x \right)}}$

$\displaystyle 1$

9.10

In traditional notation $\tan^{-1}(x)$ or $\arctan(x)$ or $\text{atan}(x)$ refers to the angle $\theta$ such that $\tan(\theta)=x$. This is completely different from the number $\tan(x)^{-1}=1/\tan(x)$. In Python you should enter sp.atan(x) (not tan^(-1)(x) or anything like that). Similar comments apply to $\text{asin}(x)$, $\text{atanh}(x)$ and so on.

x = sp.symbols('x')
θ = sp.atan(x)
sp.simplify(sp.tan(θ)) # As θ = atan(x), we have tan(θ) = x

$\displaystyle x$

sp.atan(sp.sqrt(3))

$\displaystyle \frac{\pi}{3}$

9.11

When $x\geq 0$, the notation sp.sqrt(x) means, by definition, the positive square root of $x$. (There is more to say if $x$ is negative or complex, but we pass over that here.) Thus, for example, sp.sqrt(9) is definitely $3$ and not $-3$.

sp.sqrt(9)

$\displaystyle 3$

10. Defining your own functions

10.1

If you want to define a simple function like $f(x)=x^2$, enter f = lambda x: x^2. It does not work to enter f(x)=x^2 instead.

t = sp.symbols('t')
f = lambda x : x ** 2
display(f(t+1))
display(f(100))

$\displaystyle \left(t + 1\right)^{2}$

10000

10.2

Similarly, to define the sequence $a(n)=(-1)^n/(x + n)$, enter a = lambda n: (-1)**n/(x + n), not a(n)=(-1)**n/(x + n).

x = sp.symbols('x')
k = sp.symbols('k')
a = lambda n : (-1) ** n / (x + n)
display(a(k))
display(list(a(k) for k in range(1, 11)))

$\displaystyle \frac{\left(-1\right)^{k}}{k + x}$

[-1/(x + 1),
 1/(x + 2),
 -1/(x + 3),
 1/(x + 4),
 -1/(x + 5),
 1/(x + 6),
 -1/(x + 7),
 1/(x + 8),
 -1/(x + 9),
 1/(x + 10)]

10.3

It does not work to use notation like f'(x) for the derivative of $f(x)$; see Section 12 instead.

10.4

You can define a function of several variables in a similar way. For example, to define $g(a,u,v)=u^a+v^a$, enter g = lambda a,u,v: u**a+v**a.

x, y = sp.symbols('x y')
g = lambda a, u, v : u**a + v**a
g(2, x, y)

$\displaystyle x^{2} + y^{2}$

10.5

For longer and more complicated functions you may prefer to use the more standard Python function definition syntax, like

def f(x, y, z):
    return (x+y)**8 + (x-y)**8 + (x+z)**8 + (x-z)**8

11. Plotting

In this course we mostly use the sympy library, but that does not, by itself, support plotting. We instead use the numpy and plotly libraries to make plots. (There is also another popular library called matplotlib which could be used instead of plotly, but on balance we have found plotly to be preferable.) More details of everything discussed here can be found on the Plotly website.

11.1

The code below plots the graph $y=x^3-x$ from $x=-2$ to $x=2$.

The line x = np.linspace(-2, 2, 201) sets $x$ to be an array $[x_0,x_1,x_2,\dotsc,x_{200}]=[-2,-1.98,-1.96,\dotsc,2]$ consisting of $201$ equally spaced points between $-2$ and $2$ (inclusive). The line y = x ** 3 - x sets $y$ to be the array $[y_0,\dotsc,y_{200}]=[x_0^3-x_0,\dotsc,x_{200}^3-x_{200}]$ of corresponding $y$ values. The line fig = go.Figure() creates an object called fig representing a blank picture, to which we can add curves, shapes, text and so on. In this case we just add a single curve passing through the points with coordinates $(x_i,y_i)$; this is done by the line fig.add_trace(go.Scatter(x=x, y=y)). Finally, the line fig.show() shows the picture that we have constructed. This is not always necessary, if the last line in a cell is related to plotting then Python will often decide for itself that the resulting plot should be shown, but it does not hurt to make it explicit. The graph has various interactive features which you can discover by hovering and clicking with your mouse.

x = np.linspace(-2, 2, 201)
y = x**3 - x
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y))
fig.show()

11.2

Here is a version with some additional features. We have added options to make the curve into a thick orange dotted line. The option name=r'$y = x^3-x$' specifies the name of the curve, but this is not displayed by default. The option showlegend=True specifies that there should be a box showing the names of all named elements in the plot. In this case we have a single named curve, and the legend is shown at the top right.

Note that the name is written in LaTeX, and we have written name=r'...' instead of just name='...'. The r prefix prevents special treatment of the backslash character, and is generally needed when entering strings containing LaTeX. In this case we have no backslashes so the r is not required, but it is safest to use it consistently.

Note also that the colour must be referred to using the US spelling (line_color rather than line_colour). We have specified the colour as 'orange', and many other normal colour names will also work. However, for more precise control one can use hex codes, such as #40e0d0 for a certain shade of turquoise. Sites such as htmlcolorcodes.com will help you to find appropriate hex codes.

x = np.linspace(-2, 2, 201)
y = x**3 - x
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, name=r'$y = x^3-x$', showlegend=True, line_color='orange', line_dash='dot', line_width=5))
fig.show()

11.3

Various features of the overall picture can be changed using fig.update_layout() (before fig.show()). For example, the code below specifies that the plot should be $800\times 400$ pixels, with the $x$ axis running from $x=-3$ to $x=3$ and the $y$-axis from $y=0$ to $y=10$.

x = np.linspace(-2, 2, 201)
y = x**3 - x
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name=r'$y = x^3-x$', showlegend=True))
fig.update_layout(height=400, width=800, xaxis_range=[-3,3], yaxis_range=[0,10])
fig.show()

11.4

To make the scale on the $x$-axis the same as the scale on the $y$-axis, we can add the options yaxis_scaleanchor="x",yaxis_scaleratio=1 in fig.update_layout(). Alternatively, provided you have done from header import * as suggested at the top of this file, you can enter equal_aspect(fig) to make the scales equal.

x = np.linspace(-2, 2, 201)
y = x**3 - x
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name=r'$y = x^3-x$', showlegend=True))
fig.update_layout(yaxis_scaleanchor="x",yaxis_scaleratio=1)
equal_aspect(fig) # This line has the same effect as the previous one.
fig.show()

11.5

To plot more than one graph in the same picture, we can just call fig.add_trace() several times. For example, the cell below calls fig.add_trace() once to plot the curve $y=x^3-x$, and then a second time to plot the curve $y=x^3+x$ in the same picture. We did not specify any colours, but they were chosen automatically to ensure that the colours of the two curves were different.

x = np.linspace(-2, 2, 201)
y0 = x**3 - x
y1 = x**3 + x
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y0))
fig.add_trace(go.Scatter(x=x, y=y1))
fig.show()

11.6

Here is another example where we have plotted $\cos(t)$ and $\sin(t)$ together. We have also called use_pi_ticks(fig,-3,3,2) to put tick marks on the horizontal axis at positions from $-3\pi$ to $3\pi$ with steps of $\pi/2$. If we did not do that, we would instead get tick marks at positions $-8,-6,-4,\dotsc,8$, which is not very helpful when plotting functions that are periodic with period $2\pi$.

The function use_pi_ticks() will only work if you have done from header import * as suggested at the top of this file.

Note also that we have used he function np.sin() (from the numpy library) rather than sp.sin(). Because t is an array rather than a single number, it would not work to use sp.sin(t). One could instead use [sp.sin(t0) for t0 in t], but that would not be efficient: the numpy library is highly optimised for numerical calculation, and sympy is not.

t = np.linspace(-3*np.pi, 3*np.pi, 201)
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=np.cos(t)))
fig.add_trace(go.Scatter(x=t, y=np.sin(t)))
use_pi_ticks(fig,-3,3,2)
fig.show()

11.7

To make the $x$ and $y$ axes invisible (together with the associated tick marks and labels) add the options xaxis_visible=False, yaxis_visible=False in fig.update_layout().

θ = np.linspace(-2*np.pi, 2*np.pi, 400)
u = np.sin(θ) ** 3
fig = go.Figure()
fig.add_trace(go.Scatter(x=θ, y=u, mode='lines', name=r'$y = \sin^3(\theta)$', showlegend=False))
fig.update_layout(xaxis_visible=False, yaxis_visible=False)

11.8

If we plot a function with discontinuities, like $y=1/x$, then by default we will get spurious vertical lines at the discontinuity points. The cell below shows how to avoid this. The third line finds the values of $x$ where $y$ jumps by more than $0.5$, indicating a discontinuity. (The value $0.5$ might need to be adjusted to give the right outcome in other examples.) The next two lines insert the value np.nan (indicating ``not a number'') at the corresponding places in the arrays $x$ and $y$. The plotting system takes this as an indication that the curve should be broken without spurious vertical lines.

x = np.linspace(-2,2,400)
y = 1/x
jumps = np.where(np.abs(np.diff(y)) >= 0.5)[0]
x = np.insert(x, jumps+1, np.nan)
y = np.insert(y, jumps+1, np.nan)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name=r'$y = 1/x$', showlegend=True))
fig.update_layout(height=400, width=800, yaxis_range=[-8,8], showlegend=True)

11.9

To improve the accuracy of a graph, just increase the third argument in np.linspace. In the two cells below, we plot $y= x\,\sin(\pi/x)$ for $x$ from $-1$ to $1$. This behaves quite wildly when $x$ is close to zero. In the first cell we do the plot with $50$ points, giving a result that is quite jagged. In the second cell, we have $1000$ points, which gives a much smoother picture.

x = np.linspace(-1, 1, 50)
y = x * np.sin(np.pi / x)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name=r'$y = x \sin\left(\frac{\pi}{x}\right)$', showlegend=True))
fig.update_layout(height=400, width=800, showlegend=True)

x = np.linspace(-1, 1, 1000)
y = x * np.sin(np.pi / x)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name=r'$y = x \sin\left(\frac{\pi}{x}\right)$', showlegend=True))
fig.update_layout(height=400, width=800, showlegend=True)

11.10

The example below plots the curve given parametrically by $x=1/(1+t^2)$ and $y=\sin(t)$ from $t=-5$ to $t=5$. We simply define $t$ using np.linspace() then define $x$ and $y$ in terms of $t$ (instead of defining $y$ as a function of $x$, as we did previously).

t = np.linspace(-5, 5, 400)
x = 1/(1 + t ** 2)
y = np.sin(t)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y))
fig.update_layout(height=400, width=400, xaxis_range=[-1.1,1.1], yaxis_range=[-1.1,1.1])
fig.show()

Similarly, the plot below shows the curve given by $x=\sin(2t)$ and $y=\cos(3t)$ for $-\pi\leq t\leq\pi$.

t = np.linspace(-np.pi, np.pi, 400)
x = np.sin(2*t)
y = np.cos(3*t)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode='lines'))
fig.update_layout(height=400, width=400, xaxis_range=[-1.1,1.1], yaxis_range=[-1.1,1.1], xaxis_visible=False, yaxis_visible=False)

11.11

We now explain how to plot a curve given implicitly, by a relation of the form $f(x, y)=c$ instead a formula for $y$ in terms of $x$. We will use the example $f(x, y)=x+y+x^2y^2$.

The first line sets x0 to ba a one-dimensional array of $x$-values, and similarly for the second line. The third line x, y = np.meshgrid(x0, y0) sets $x$ to be a two-dimensional array of $x$-values, and $y$ to be a two-dimensional array of $y$-values. The fourth line uses these to compute a two-dimensional array $z$ of values of the function $x+y+x^2y^2$, then the sixth line draws the curve where $z=1$.

x0 = np.linspace(-3,3,100)
y0 = np.linspace(-4,4,100)
x, y = np.meshgrid(x0, y0)
z = x + y + x**2 * y**2
fig = go.Figure()
fig.add_trace(go.Contour(x=x0, y=y0, z=z, contours=dict(type='constraint', value=1), line_color='red'))
fig.update_layout(height=600, width=600, xaxis_range=[-3,3], yaxis_range=[-4,4])

11.12

Here is an example where we create two separate pictures called fig1 and fig2 (which are not displayed), and then combine them to make a picture called fig, which is displayed.

x = np.linspace(-1,1,100)
y = x * (x ** 2 - 1)
z = 1 - x**4
fig1 = go.Figure()
fig1.add_trace(go.Scatter(x=x, y=y, mode='markers', name=r'$y = x(x^2-1)$', showlegend=True))
fig2 = go.Figure()
fig2.add_trace(go.Scatter(x=x, y=z, mode='lines', name=r'$y = 1 - x^4$', showlegend=True))
fig = go.Figure(data = fig1.data + fig2.data)
fig.show()

11.13

By default, go.Scatter(x=u, y=v) will produce a connected curve passing through the points $(u_i,v_i)$. Sometimes we might instead just want to show the points $(u_i,v_i)$ without connecting lines. For this, we need to add the option mode='markers' to go.Scatter(), as in the example below:

fig = go.Figure()
fig.add_trace(go.Scatter(x=[11,22,33,44,55,66,77,88,99],y=[50,40,30,20,10,20,30,40,50], mode='markers'))

Often we just want to display the $y$-values, and the $x$-values are not important, we just take them to be $0,1,\dotsc,n-1$. In that case we can just omit $x$ when calling go.Scatter(), as in the example below.

fig = go.Figure()
fig.add_trace(go.Scatter(y=[1,2,3,4,5,4,3,2,1,2,3,4,5,4,3,2,1], mode='markers'))

12. Differentiation

x, y, t, θ = sp.symbols('x y t θ')

12.1

To differentiate $x+\sin(x)$ with respect to $x$ (for example), enter diff(x+sin(x),x). This is equivalent to $\frac{d}{dx}(x+\sin(x))$ in traditional notation. Similarly, to differentiate $1/(1-t^2)$ with respect to $t$ we enter diff(1/(1-t^2),t).

x, t = sp.symbols('x t')
display(sp.diff(x + sp.sin(x), x))
display(sp.diff(1/(1 - t**2), t))

$\displaystyle \cos{\left(x \right)} + 1$

$\displaystyle \frac{2 t}{\left(1 - t^{2}\right)^{2}}$

12.2

To differentiate $y$ three times with respect to $x$ (for example), enter diff(y,x,3). This is equivalent to $d^3y/dx^3$ in traditional notation. Similarly, to find the second derivative of $\sin(10\theta)$ with respect to $\theta$ we enter diff(sin(10*theta),theta,2).

x, y = sp.symbols('x y')
y = x ** 101
display(y)
display(sp.diff(y, x, 3))

$\displaystyle x^{101}$

$\displaystyle 999900 x^{98}$

theta = sp.symbols('theta')
sp.diff(sp.sin(10 * theta), theta, 2)

$\displaystyle - 100 \sin{\left(10 \theta \right)}$

12.3

If $x$ and $y$ are related implicitly by an equation (say $x+y+x^2y^2=1$) then you can find $dy/dx$ in terms of $x$ and $y$ by implicit differentiation, like this:

x, y = sp.symbols('x y')
sp.idiff(x + y + x**2*y**2 - 1, y, x)

$\displaystyle - \frac{2 x y^{2} + 1}{2 x^{2} y + 1}$

12.4

Sometimes we just want to mention the derivative without actually working it out. We can do this by writing sp.Derivative() instead of sp.diff(). For example, entering u = sp.Derivative(x**2, x) gives an object $u$ which is displayed as $\frac{d}{dx}x^2$. If you then enter u.doit() then SymPy will evaluate the derivative and return $2x$. Below we display the unevaluated and evaluated derivatives separately, then we illustrate how to make a LaTeX string which displays the equation relating them, namely $\frac{d}{dx}x^2 = 2x$.

x = sp.symbols('x')
u = sp.Derivative(x**2,x)
display(u)
display(u.doit())
display(Latex('$' + sp.latex(u) + ' = ' + sp.latex(u.doit()) + '$'))

$\displaystyle \frac{d}{d x} x^{2}$

$\displaystyle 2 x$

$\frac{d}{d x} x^{2} = 2 x$

Similarly, v = sp.Derivative(t**10,t,2) prints as $\frac{d^2}{dt^2}t^{10}$, whereas v.doit() or sp.diff(t**10,t,2) gives $90 t^8$.

t = sp.symbols('t')
v0 = sp.Derivative(t**10, t, 2)
v = v0.doit()
display(v0)
display(v)
display(Latex('$' + sp.latex(v0) + ' = ' + sp.latex(v) + '$'))

$\displaystyle \frac{d^{2}}{d t^{2}} t^{10}$

$\displaystyle 90 t^{8}$

$\frac{d^{2}}{d t^{2}} t^{10} = 90 t^{8}$

12.5

To find the Taylor series of a function $f(x)$ near a point $x=a$ to order $n$, enter sp.series(f(x),x,a,n). Note that ``order $n$'' means that the highest term allowed is a multiple of $x^{n-1}$, and all terms $x^n$ and higher are discarded. For example, sp.series(1/x,x,2,7) gives the Taylor series for $1/x$ near $x=2$ including terms up to $(x-2)^6$.

x = sp.symbols('x')
sp.series(1/x, x, 2, 7)

$\displaystyle 1 + \frac{\left(x - 2\right)^{2}}{8} - \frac{\left(x - 2\right)^{3}}{16} + \frac{\left(x - 2\right)^{4}}{32} - \frac{\left(x - 2\right)^{5}}{64} + \frac{\left(x - 2\right)^{6}}{128} - \frac{x}{4} + O\left(\left(x - 2\right)^{7}; x\rightarrow 2\right)$

12.6

The command p = sp.series(...) command gives an answer including a term like $O(x^7)$, to indicate that there are infinitely many more terms, starting with a multiple of $x^{7}$. We cannot do very much with this answer until we have converted it to an ordinary polynomial (without the $O(x^7)$ term on the end). We can do this by entering p = p.removeO().

x = sp.symbols('x')
p = sp.series(sp.sin(x), x, 0, 7)
display(p)
display(p.removeO())

$\displaystyle x - \frac{x^{3}}{6} + \frac{x^{5}}{120} + O\left(x^{7}\right)$

$\displaystyle \frac{x^{5}}{120} - \frac{x^{3}}{6} + x$

13. Integration

13.1

To calculate an indefinite integral like $\int\tan(x)\,dx$, enter sp.integrate(sp.tan(x),x).

x = sp.symbols('x')
sp.integrate(sp.tan(x), x)

$\displaystyle - \log{\left(\cos{\left(x \right)} \right)}$

13.2

SymPy's answer will never include an arbitrary constant '$+c$'. In some contexts (particularly, solution of differential equations) the constant makes a difference and must be included. In other contexts it is irrelevant; SymPy cannot help you to decide about this.

13.3

To calculate a definite integral like $\int_1^5t^3\,dt$, enter sp.integrate(t**3,(t,1,5)). Similarly, to calculate $\int_0^2\sqrt{1+x}\,dx$, enter sp.integrate(sp.sqrt(1+x),x=0..2).

t = sp.symbols('t')
sp.integrate(t**3,(t,1,5))

$\displaystyle 156$

x = sp.symbols('x')
sp.integrate(sp.sqrt(1 + x), (x, 0, 2))

$\displaystyle - \frac{2}{3} + 2 \sqrt{3}$

13.4

It is also possible to have infinite limits. For example, $\int_{-\infty}^\infty (1+x^4)^{-1}dx$ can be entered as shown below.

sp.integrate(1/(1+x**4),(x,-sp.oo,sp.oo))

$\displaystyle \frac{\sqrt{2} \pi}{2}$

13.5

As with differentiation, it is sometimes useful to mention an integral without evaluating it. You can do this by writing sp.Integral() instead of sp.integrate(). For example, if you enter u = sp.Integral(x^2,x) then you will get an object that is displayed as $\int x^2\,dx$; if you then enter u.doit(), you will get $x^3/3$.

u0 = sp.Integral(x**2,x)
u = u0.doit()
display(Latex('$' + sp.latex(u0) + ' = ' + sp.latex(u) + '$'))

$\int x^{2}\, dx = \frac{x^{3}}{3}$

14. Sequences and sums

14.1

We might want to talk about the sequence $[1/2, 1/4, 1/6, 1/8, 1/10]$, which is the sequence of terms $1/(2k)$ for $k$ from $1$ to $5$ inclusive. There are two ways to do this. The first is to use an ordinary Python list comprehension, as shown below. For this we need to use the standard Python convention for ranges: the endpoint of the range is never included, so the list $1,\dotsc,5$ is range(1,6) not range(1,5).

Although the cell below mentions a variable k, this is actually a new temporary variable that only has meaning within the list comprehension. It does not matter whether we already had a variable called k, that would not affect anything.

[sp.Rational(1, 2*k) for k in range(1,6)]

[1/2, 1/4, 1/6, 1/8, 1/10]

The second approach is to use sp.sequence(). This is different in various ways:

• Here the variable k needs to be a SymPy symbol with no assigned value.
• The range of the variable (which includes k = 5 but not k = 6) is now written as (k, 1, 5), matching traditional mathematical notation but not the usual Python convention for ranges.
• In this example the upper limit for k is the explicit integer $5$, but we could also do something like sp.sequence(1/(2 * k), (k, 1, n)) where n is a SymPy symbol with no assigned value. This would not be possible with an ordinary list comprehension. We can convert a SymPy sequence to an ordinary list using the list() function.

k = sp.symbols('k')
u = sp.sequence(1/(2 * k), (k, 1, 5))
display(u)
display(list(u))

$\displaystyle \left[\frac{1}{2}, \frac{1}{4}, \frac{1}{6}, \frac{1}{8}, \ldots\right]$

[1/2, 1/4, 1/6, 1/8, 1/10]

14.2

For an ordinary Python list, we can use sum() to find the sum of the terms.

u = list(10 ** k for k in range(5))
display(u)
display(sum(u))
display(sum(10 ** k for k in range(5)))

[1, 10, 100, 1000, 10000]

11111

11111

14.3

For a SymPy sequence in which the upper and/or lower limits are symbols, we need to use sp.Sum() instead of sum(). For example, we can get the geometric progression formula for $\sum_{k=0}^n 3^{-k}$ as follows:

k, n = sp.symbols('k n', integer=True)
u0 = sp.Sum(3**(-k),(k,0,n))
u = sp.simplify(u0.doit())
display(Latex('$' + sp.latex(u0) + ' = ' + sp.latex(u) + '$'))

$\sum_{k=0}^{n} 3^{- k} = \frac{3}{2} - \frac{3^{- n}}{2}$

Note that the naming and behaviour is a bit inconsistent with the naming and behaviour for integrals. The function sp.Sum() always returns an unevaluated sum which is displayed as something like $\sum_{k=0}^n3^{-k}$; it is always necessary to use .doit() to find the actual value of the sum.

14.4

We can also handle some infinite sums. For example, it is a well-known result that $\sum_{k=1}^\infty k^{-2}=\pi^2/6$, and we can see this in SymPy as shown below.

u0 = sp.Sum(1 / n**2, (n, 1, sp.oo))
u = u0.doit()
display(Latex('$' + sp.latex(u0) + ' = ' + sp.latex(u) + '$'))

$\sum_{n=1}^{\infty} \frac{1}{n^{2}} = \frac{\pi^{2}}{6}$

15. Some common problems

15.1

All abstract symbols need to be introduced using something like x = sp.symbols('x'). If you try to use $x$ without having declared it, you will get an error like name 'x' is not defined. It can also happen that you give $x$ a value, perhaps by entering $x = 42$, and later you forget that you did that, and try to treat $x$ as an abstract symbol. This might give you any of a range of different error messages, or it might just give you an answer that is incorrect or meaningless but with no error message. You can check the status of $x$ by entering type(x). If $x$ has been declared as an abstract symbol, then type(x) should give sympy.core.symbol.Symbol. If $x$ has been set equal to $42$, then type(x) should give int. If $x$ has not been introduced in any way, then type(x) will give an error. If $x$ has been given a value that is no longer relevant, you can enter del x to remove it.

x = sp.symbols('x')
display(type(x))
x=42
display(type(x))
del x
display(type(x)) # Gives an error

sympy.core.symbol.Symbol

int

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[207], line 6
      4 display(type(x))
      5 del x
----> 6 display(type(x)) # Gives an error

NameError: name 'x' is not defined

15.2

Remember *'s for multiplication.

• Enter $2(x+y)$ as 2*(x+y). If you instead enter 2(x+y) then you may get a cryptic error message saying 'int' object is not callable.
• Enter $xy$ as $x*y$. If you instead enter xy then you will usually get an error message name 'xy' is not defined.
• For the product of $a$ and $x+y$, enter a*(x+y).
• To evaluate a function $f(t)$ at the point $t=x+y$, enter $f(x+y)$. If $a$ is an object for which a*(x+y) makes sense, then $a(x+y)$ will not make sense and will give an error. If $f$ is an object for which $f(x+y)$ makes sense, then $f*(x+y)$ will not make sense and will give an error.

2(x + y)

<>:1: SyntaxWarning:

'int' object is not callable; perhaps you missed a comma?

C:\Users\Neil\AppData\Local\Temp\ipykernel_14536\1089367336.py:1: SyntaxWarning:

'int' object is not callable; perhaps you missed a comma?

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[14], line 1
----> 1 2(x + y)

TypeError: 'int' object is not callable

xy

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[113], line 1
----> 1 xy

NameError: name 'xy' is not defined

15.3

The fraction $\frac{a+b}{c+d}$ must be entered as (a+b)/(c+d) (not a+b/c+d), and $\frac{a}{(u+v)(x+y)}$ must be entered as a/((u+v)*(x+y)) (not a/(u+v)*(x+y)).

a, b, c, d = sp.symbols('a b c d')
display((a + b)/(c + d))
display( a + b / c + d )

$\displaystyle \frac{a + b}{c + d}$

$\displaystyle a + \frac{b}{c} + d$

a, u, v, x, y = sp.symbols('a u v x y')
display(a / ((u + v) * (x + y)))
display(a /  (u + v) * (x + y) )

$\displaystyle \frac{a}{\left(u + v\right) \left(x + y\right)}$

$\displaystyle \frac{a \left(x + y\right)}{u + v}$

15.4

$a^{1/4}$ must be entered as a ** sp.Rational(1,4), not a ** (1/4) or a ** 1/4.

a = sp.symbols('a')
display(a ** sp.Rational(1, 4))
display(a ** (1 / 4))
display(a **  1 / 4 )

$\displaystyle \sqrt[4]{a}$

$\displaystyle a^{0.25}$

$\displaystyle \frac{a}{4}$

15.5

Always use round brackets (not square or curly ones) for algebraic grouping. Square brackets are used for lists or for subscripts (e.g. $x_5$ is x[5]). Curly brackets are used for sets or Python dictionaries. Using the wrong kind of brackets will usually generate an error.

a, b, c, d = sp.symbols('a b c d')
sp.expand((a + b) * (c + d))
sp.expand({a + b} * {c + d}) # gives an error
sp.expand([a + b] * [c + d]) # gives an error

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[20], line 4
      2 sp.expand((a + b) * (c + d))
      3 # sp.expand({a + b} * {c + d}) # gives an error
----> 4 sp.expand([a + b] * [c + d]) # gives an error

TypeError: can't multiply sequence by non-int of type 'list'

15.6

Remember to remove definitions when you have finished with them, by using del or clicking the Restart button (but remember that you need to re-execute your import cells after restarting). If you enter an expression like $(x+y)^5$ and SymPy converts it to something completely unrelated like $(\sin(\theta)+e^{-u})^5$, the most likely reason is that you previously defined x=sp.sin(theta); and y=sp.exp(-u);, and you forgot to remove these definitions. You can enter x,y = sp.symbols('x y') to reset $x$ and $y$ to be symbols with no assigned value.

u, x, y, θ = sp.symbols('u x y θ')
x = sp.sin(θ)
y = sp.exp(-u)

(x + y) ** 5

$\displaystyle \left(\sin{\left(θ \right)} + e^{- u}\right)^{5}$

x, y = sp.symbols('x y')
(x + y) ** 5

$\displaystyle \left(x + y\right)^{5}$

15.7

Remember the constant $e=2.71828\dots$ must be entered as sp.E or sp.exp(1), and $e^x$ must be entered as exp(x). If you enter e ** x instead of sp.exp(x), you will most often get an error message name 'e' is not defined. However, if you happen to have defined $e$, you might get some other behaviour.

x = sp.symbols('x')
display(sp.exp(1))
display(sp.exp(10))
display(sp.exp(x))
display(sp.E ** x)

$\displaystyle e$

$\displaystyle e^{10}$

$\displaystyle e^{x}$

$\displaystyle e^{x}$

e ** x # Gives an error because e is not defined.

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[27], line 1
----> 1 e ** x # Gives an error because e is not defined.

NameError: name 'e' is not defined

15.8

Remember that the constant $\pi=3.14159\dots$ must be entered as sp.pi, not pi. If you try to use pi instead of sp.pi then you will usually get an error name 'pi' is not defined. If you want to be able to use pi with no prefix, you can enter from sympy import pi at the top of your notebook.

sp.sin(sp.pi / 2)

$\displaystyle 1$

sp.sin(pi / 2) # Gives an error because pi is not defined.

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[31], line 1
----> 1 sp.sin(pi / 2) # Gives an error because pi is not defined.

NameError: name 'pi' is not defined