sympy

I’ve been using mathematica for all my numerical calculations, but as a web system developer, I’d like to move to python.

First of all, I’d like to start with the basics.

I’ve been using mathematica for all my numerical calculations, but as a web system developer, I’m moving to python.

github

  • The file in jupyter notebook format is here .

google colaboratory

  • If you want to run it in google colaboratory here base/base_nb.ipynb)

Author’s environment

The author’s OS is macOS, and the options are different from Linux and Unix commands.

!sw_vers
ProductName: Mac OS X
ProductVersion: 10.14.6
BuildVersion: 18G6020
python -V
Python 3.7.3

Import the basic libraries and check their versions.

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import matplotlib
import matplotlib.pyplot as plt
import scipy
import numpy as np

print('matplotlib version :', matplotlib.__version__)
print('scipy version :', scipy.__version__)
print('numpy version :', np.__version__)
matplotlib version : 3.0.3
scipy version : 1.4.1
numpy version : 1.16.2

The three numeric types of sympy

sympy is a

  • Real
  • Rational
  • Integer

sympy seems to have three types.

Handling fractions

from sympy import *

a = Rational(1,2)
a
1/2

Handle pi

pi.evalf()
3.14159265358979
pi * 2
2*pi
pi ** 2
pi**2

Handling natural logarithms

exp(1) ** 2
exp(2)
(exp(1) ** 2).evalf()
7.38905609893065

infinity

oo > 999
True
oo + 1
oo

Use evalf to specify the number of digits.

The number of digits is given as an argument.

pi.evalf(1000)
 95505058222317253595942828282848111745252840272070019385210555964642622948954393819644848281964848482828476783165271201909145648566  519415160943357527036575959591953009181617381932161793151515151585480744626379796274912983367336244065664  72414684409012249495343146654958353701050792796892589258923542019951612290219068464043148418159813269774777713099605177072113499999983729780499 51015597317321816096318559502445949455353669088322625222300825334468535353535252619331817070100031378383875282865875332838142061717766914730359825349 0428755468731595628383835378775937151957788185778053271712268606161301927876611195909216420199

logarithm

log(10)
log(10)
log(2,2**4).evalf()
0.250000000000000
E.evalf()
2.71828182845905

Define symbols

Declare variables such as $x and $y to allow algebraic operations.

x = Symbol('x')
y = Symbol('y')
x + y ** 2 + x
2*x + y**2
(x + y) ** 3
(x + y)**3

Algebraic operations

expand and simplify seem to be the same as in mathematica, which is helpful.

expand.

expand((x + y)**3)
x**3 + 3*x**2*y + 3*x*y**2 + y**3

Simplify

simplify((x + x*y) / x)
y + 1

factorization

factor(x**3 + 3*x**2*y + 3*x*y**2 + y**3)
(x + y)**3

Calculus

Extremes

limit(sin(x)/ x, x, 0)
1
limit(x, x, oo)
oo
limit(1 / x, x, oo)
0
limit(x**x, x, 0)
1

Differentiate

diff(sin(x), x)
cos(x)
diff(sin(x) ** 2, x)
2*sin(x)*cos(x)
diff(sin(2 * x), x)
2*cos(2*x)
diff(tan(x), x)
tan(x)**2 + 1

Checking if the derivative is correct

You can take the extremum and check if the derivative is correct. We can take the extremum and check if the derivative is correct, just follow the definition of the derivative.

limit((tan(x+y) - tan(x))/y, y, 0)
tan(x)**2 + 1

Higher-order differentiation is also possible.

diff(sin(x), x, 1)
cos(x)

You can see that the second-order derivative is the original value multiplied by a negative number.

diff(sin(x), x, 2)
-sin(x)
diff(sin(x), x, 3)
-cos(x)
diff(sin(x), x, 4)
sin(x)

Series expansion

Taylor expansion is also possible.

series(exp(x), x)
1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
series(cos(x), x)
1 - x**2/2 + x**4/24 + O(x**6)
series(sin(x), x)
x - x**3/6 + x**5/120 + O(x**6)

I thought it might be possible to specify the order of the expansion as the third argument.

series(exp(x), x, 6)
exp(6) + (x - 6)*exp(6) + (x - 6)**2*exp(6)/2 + (x - 6)**3*exp(6)/6 + (x - 6)**4*exp(6)/24 + (x - 6)**5*exp(6)/120 + O((x - 6)**6, (x, 6))

It looks like the third argument is the number of centers in the numerical calculation, and the fourth argument is the order of the expansion.

series(exp(x), x, 0, 6)
1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)

Integration

If you can do derivatives, integrals are of course supported.

integrate(x**3,x)
x**4/4
integrate(-sin(x),x)
cos(x)
integrate(exp(x),x)
exp(x)
integrate(log(x),x)
x*log(x) - x
integrate(exp(-x**2),x)
sqrt(pi)*erf(x)/2

You can also specify an integration interval

integrate(x**2, (x, -1, 1))
2/3
integrate(sin(x), (x, 0, pi/2))
1

Wide integrals with infinite range are also possible.

integrate(exp(-x), (x, 0, oo))
1

Try to do a Gaussian integral.

integrate(exp(-x**2), (x, -oo, oo))
sqrt(pi)
integrate(exp(-x**2 / 3), (x, -oo, oo))
sqrt(3)*sqrt(pi)

Equation Solving

It is also possible to solve algebraic equations. The interface is almost identical to mathematica.

solve(x**3 - 1, x)
[1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2].
solve(x**4 - 1, x)
[-1, 1, -I, I].

Multivariable linear equations are also supported.

solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
{x: -3, y: 1}

It also seems to calculate Euler’s equation.

solve(exp(x) + 1, x)
[I*pi].

Matrix operations

We’ve been doing matrix operations in numpy for a long time, so we probably won’t be using it for matrices, but let’s run it by hand here.

from sympy import Matrix

a = Matrix([[1,0], [0,1]])
b = Matrix([[1,1], [1,1]])
a
1/2
b
Matrix([
[1, 1],
[1, 1]])
a + b
Matrix([.
[2, 1],
[1, 2]])
a * b
Matrix([.
[1, 1],
[1, 1]])

Differential equations

Ordinary differential equations can also be solved, using dsolve.

f, g = symbols('f g', cls=Function)
f(x)
f(x)
f(x).diff(x,x) + f(x)
f(x) + Derivative(f(x), (x, 2))
dsolve(f(x).diff(x,x) + f(x), f(x))
Eq(f(x), C1*sin(x) + C2*cos(x))

This was a very powerful program that I could see being used for actual numerical calculations.

Of course, there are some parts that are inferior to mathematica, which is a paid service, but I will be actively migrating from mathematica. It’s amazing!

Reference page

I used this page as a reference.