Using a computer algebra system (CAS), you can quicklycarry out mathematical calculations like: computing mathematical expressions, solving equations, performing calculus procedures, and simulating physics systems. Computer algebra systems are great when you're learning math, because they will help you simplify complicated expressions and solve hairy equations.

All computer algebras systems offer essentially the same functionality so it doesn't matter which one you'll use: the free ones like sympy, Magma, and Octave, or the commercial ones Maple, MATLAB, and Mathematica. This tutorial is an introduction to sympy, which is a symbolic computer algebra system written in python. In a symbolic CAS numbers and operations are represented precisely so the answers you obtain are exact. For example, in sympy the number $\sqrt{2}$ is represented as the object Pow(2,1/2), whereas in numerical computer algebra systems like octave the number $\sqrt{2}$ is represented as the approximation $1.41421356237310$ (a float). For most purposes the approximation is okay, but sometimes approximations can lead to problems: float(sqrt(2))*float(sqrt(2)) = 2.00000000000000044 $\neq 2$. Using sympy, you won't run into such problems: Pow(2,1/2)*Pow(2,1/2)=2 $= 2$.

Do not use sympy to cheat on your homework! The whole point of homework is for you to suffer. Mathematical skill is developed through mathematical suffering—only when trying to solve a problem that you haven't solved before will you be forced to think and practice your skills. You can use sympy to check your answers though. In fact, one of the best ways to learn math is to invent practice problems for yourself, solve them by hand, and then check your answers using sympy.

This tutorial is organized as follows. We'll begin by introducing the basics of using sympy and the bread-and-butter functions used for calculating, manipulating expressions, and solving equations. Afterward, we'll discuss the sympy functions that implement calculus operations like differentiation and integrations. We'll also introduce the sympy functions used to deal with vectors and complex numbers. Later we'll see how to use vectors and calculus operations to understand Newtonian mechanics. In the last section, we'll introduce the linear algebra functions available in sympy.

The easiest way to use sympy, provided you're connected to the Internet, is to visit ''live.sympy.org''. You'll be presented with an interactive prompt into which you can enter your commands—right from your browser.

If you want to use sympy on your own computer, you must install python and the sympy python package. You can then open a command prompt and start a sympy session using:

you@host$ python
Python X.Y.Z 
[GCC a.b.c (Build Info)] on platform
Type "help", "copyright", or "license" for more information.
>>> from sympy import *
>>> 

The »> prompt indicates you are in the python shell which accepts python commands. The python command from sympy import * imports all the sympy functions into the current namespace. All the sympy functions are now available to you.

To exit the python shell press CTRL+D:

>>>  ( press CTRL + D )
you@host$  

I highly recommend that you also install ipython, which is an improved interactive python shell. If you have ipython and sympy installed, you can start an ipython shell with sympy pre-imported using the command isympy. For an even better experience, you can try ipython notebook, which is a web frontend for the ipython shell.

Each section in this tutorial begins with a specific python import statement, which indicates the functions used in that section. If you use the statement from sympy import * in the beginning of your code, you don't need to run these individual import statements, but it's a good idea to look at them so you'll know what part of your sympy vocabulary will be developed in that section.

In this section we'll learn about the basic sympy objects and the math operations we can carry out on them. We'll learn the sympy equivalents of many of the math verbs we used throughout the book: “to solve” (an equation), “to simplify” (an expression), “to expand” (an expression), etc.

This section serves a dual purpose. It serves to introduce readers familiar with the fundamentals of math to the power of computers, so they will never have to suffer through another tedious arithmetic calculation done by hand again. This section also serves as a review of the fundamental concepts of mathematics for computer people, who may have forgotten their high school math.

The basic number objects in sympy are the python numbers: ints and floats:

>>> 3
3            # an int
>>> 3.0
3            # a float 

Special care is required when specifying rational numbers, because the int division doesn't return the right answer.

>>> 1/7
0                  # an int 

If you force float division by using 1.0 instead of 1 you'll obtain:

1.0/7
0.14285714285714285    # a float

This is better, we only have a float approximation of the number $\frac{1}{7} \in \mathbb{Q}$. To obtain an exact representation of $\frac{1}{7}$ you need to create a sympy expression. You can sympify an expression using the shortcut function S():

S('1/7')
1/7         #  Rational(1,7)
  

We could have achieved the same result using S('1')/7 since, the sympy object divided by an int is a sympy object.

Apart from the tricky division operator, the other math operators: +, -, and * work as you would expect. The syntax *EXP* is used in python to denoted exponents:

>>> 2**10
1024

When solving math problems, it is always preferable to work with sympy objects, and only compute the numeric answer in the end. To obtain a numeric approximation of a sympy object as a float, call its .evalf() method:

>>> pi
pi
>>> pi.evalf()
3.14159265358979

The method .n() and global function N() are equivalent to .evalf(). You can also specify the number of digits of precision you want as an argument to the above functions: pi.n(400) computes an approximation of $\pi$ to 300 decimals.

  from sympy import Symbol, symbols

Python is a civilized language so you don't need to define variables before you use them. When you write

>>> a = 3

you define a new name a and set it the value 3. From now on you can use a in your calculations.

Most interesting sympy calculations require us to define some symbols, which are the basic math objects in sympy. For your convenience, when [http://live.sympy.org|live.sympy.org] starts, it runs the following commands for you:

>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)

THe first statements instructs python to convert 1/7 to 1.0/7 when dividing, thus potentially saving you from the int division confusion. The second statement imports all the sympy functions. The remaining statements define some generic symbols x, y, z, and t, and some other symbols with special properties (to be discussed later).

Note the difference between the following two statements:

>>> x + 2
x + 2       # an Expr 
>>> p + 2 
NameError: name 'p' is not defined

The name x was defined to be a symbol so sympy knows that x + 2 is an expression, but the variable p was not defined so sympy doesn't know what to make of p + 2. If you want to use p in expressions, you must first define it as a symbol:

>>> p = Symbol('p')     # the same as p = symbols('p')
>>> p+2 
p + 2

You can define a sequence of variables using the following notation:

>>> a0, a1, a2, a3, a4 = symbols('a0:4')

You can use any name you want for a variable, but it's best if you avoid the letters in the sting QCOSINE because they have a special use in sympy: I is the unit imaginary number $i\equiv\sqrt{-1}$, E is the base of the natural logarithm, S() is the sympify function, N() is used to obtain numeric approximations, O is used for big-O notations, etc.

The underscore symbol _ is a special variable that contains the result of the last return value. The variable _ is analogous to the ans button on certain calculators, and is useful in multi-step calculations:

>>> 3+3
6
>>> _*2
12

>>> from sympy import simplify, factor, expand, collect

You define sympy expressions by combining symbols, the basic math operations, and other functions:

>>> expr = 2*x + 3*x - sin(x) - 3*x+42
>>> simplify(expr)
2*x - sin(x) + 42

The function simplify can be used with any expression to simplify it. The examples below illustrate some other useful sympy functions that correspond to common mathematical operations on expressions:

>>> factor( x**2-2*x-8 )
(x - 4)*(x + 2)
>>> expand( (x-4)*(x+2) )
x**2 - 2*x - 8
>>> collect(x**2 + x*b + a*x + a*b, x)
x**2 + (a+b)*x + a*b

To substitute a given value into an expression, we call the .subs() method, passing in a python dictionary object with the variable–value substitutions you want to make:

>>> expr  = sin(x) + cos(y)
>>> expr
sin(x) + cos(y)
>>> expr.subs({x:1, y:2})
sin(1) + cos(2)
>>> expr.subs({x:1, y:2}).n()
0.425324148260754

Note we used the .n() method to obtain a the numeric value of the expression.

>>> from sympy import solve

The function solve is the workhorse of sympy. This incredibly powerful function knows how to solve all kinds of equations. You can throw pretty much any equation at solve to find the solutions. The function solve is so powerful that it can be somewhat disheartening for high school students to find out about it—they'll wonder why they wasted 5 years of their life learning about various functions and ways to solve equations when solve can do all the work for them. Don't worry,

The function solve takes two arguments. The first argument specifies the equation you want to solve expressed in the form expr==0. The second argument is the variable you want to solve for. For example, to solve the quadratic equation $x^2+2x-8=0$, we use

>>> solve( x**2 + 2*x - 8, x)
[2, -4]

In this case the equation has two solutions so solve returns a list. Check that $x=2$ and $x=-4$ satisfy the equation $x^2+2x-8=0$.

The best part about solve and sympy is that you can obtain symbolic answers when solving equations. Instead of solving one specific quadratic equation, we can solve all possible equations of the form $ax^2 + bx+c=0$ using

>>> a, b, c = symbols('a b c')
>>> solve( a*x**2 + b*x + c, x)
[(-b + sqrt(b**2 - 4*a*c))/(2*a), (-b-sqrt(b**2-4*a*c))/(2*a)]

In this case solve calculated the solution in terms of the symbols a, b, and c. You should be able to recognize the solution in this case corresponds to the quadratic formula.

You can verify that substiuting the coefficients $a=1$, $b=2$, and $c=-8$ into general solutions results in the same answer as we obtained above:

>>> gen_sol = solve( a*x**2 + b*x + c, x)
>>> [ gen_sol[0].subs({'a':1,'b':2,'c':-8}), 
      gen_sol[1].subs({'a':1,'b':2,'c':-8})  ]
[2, -4]

To solve a system of equations, you can feed solve with the list of equations as the first argument and specify a list of unknowns to solve for as the second argument. For example to solve for $x$ and $y$ in the system of equations $x+y =3$ and $3x-2y=0$, use:

>>> solve([x + y - 3, 3*x - 2*y], [x, y])
{x: 6/5, y: 9/5}

The function solve is like a swiss army knife you can use to solve all kind of problems. For example, suppose you want to completing the square in the expression $x^2-4x+7$, that is, you want to find constants $h$ and $k$ such that $x^2-4x+7 = (x-h)^2 +k$. There is no special purpose “complete the square” function in sympy, but you call solve on the equation $(x-h)^2 +k \ - \ (x^2-4x+7) = 0$ to find the unknowns $h$ and $k$:

>>> h, k = symbols('h k')
>>> solve( (x-h)**2 + k  - (x**2-4*x+7), [h,k] )
[(2, 3)]                               # so h = 2 and k = 3
>>> ((x-2)**2+3).expand()              # verify...
x**2 - 4*x + 7

>>> from sympy import together, apart

By default, sympy will not put together or split apart rational expressions. You have to use together to calculate symbolically the addition of fractions:

>>> togeter( 1 + 1/(x + 2) )
(x + 3)/(x + 2)

Alternately, if you have a rational expression an you want to divide the numerator by the denominator, use the apart function:

>>> apart( (x**2+x+4)/(x+2)  )
x - 1  +  6/(x + 2)

Euler's constant $e=2.71828\ldots$ in defined one of several ways: \[ e \equiv \lim_{n\to \infty} \left( 1 + \frac{1}{n}\right)^{n} \equiv \lim_{\epsilon \to 0} \left( 1 + \epsilon\right)^{1/\epsilon} \equiv \sum_{n=0}^\infty \frac{1}{n!}, \] is denoted E in sympy. The funtion exp(x) is equivalent to E*EXP*x.

The function log computes the natural logarithm (logarithm base $e$):

>>> log(E**3)    # same as ln(E**3)
3

By default, the sympy assumes the inputs to functions like exp and log are complex numbers, so it will not expand certain logarithmic expressions. However, indicating to sympy that the inputs are real numbers will make the epxations work:

>>> x, y = symbols('x y')
>>> log(x*y).expand()
log(x*y)
>>> a, b = symbols('a b', positive=True)
>>> log(a*b).expand()
log(a) + log(b)

Let's define a polynomial $P$ that has roots at $x=1$, $x=2$ and $x=3$.

>>> P = (x-1)*(x-2)*(x-3)
>>> P
(x - 1)*(x - 2)*(x - 3)

To see the expaned version of the polynomial, call its expand method:

>>> P.expand()
x**3 - 6*x**2 + 11*x - 6

When you see the polynomial in it's expanded form $P(x)=x^3-6^2 + 11x - 6$, you can't immediately tell what it's roots are. This is why the factored form $P(x)=(x-1)(x-2)(x-3)$ is preferable. To factor a polynomial simply call its factor method:

>>> P.factor()
(x - 1)*(x - 2)*(x - 3)

Or, you can call the more generic simplify which is used to “simplify” any expression:

>>> P.simplify()
(x - 1)*(x - 2)*(x - 3)

Recall that the roots of the polynomial are defined as the values of $x$ where $P(x)=0$. We can use the solve function to find the roots of the polynomial:

>>> roots = solve(P,x)
>>> roots
[1, 2, 3]
# let's check if equal 
>>> simplify( P  -  (x-roots[0])*(x-roots[1])*(x-roots[2]) )   
0  

Note we used the simplify function to check whether two expressions are equal in the previous section. This way of checking equality works because $P=Q$ if and only if $P-Q=0$. This is the best for if two expeession are equal in sympy because it takes care of all possible simplifications while comparing the expressions. Below is a list of other ways for checking if two quantities are equal, and some examples where they work and where they fail:

>>> p = (x-5)*(x+5)
>>> q = x**2 - 25
>>> p == q                        # fail
False
>>> p - q == 0                    # fail 
False
>>> expand(p - q) == 0         
True
>>> factor(p - q) == 0          
True
>>> sin(x)**2 + cos(x)**2  == 1   # fail
False
>>> simplify( sin(x)**2 + cos(x)**2  - 1) == 0
True

from sympy import sin, cos, tan, trigsimp, expand_trig

The trigonometric functions sin and cos take inputs in radians:

>>> sin(pi/6)
1/2
>>> cos(pi/6)
sqrt(3)/2

For angles in degrees, you need the conversiton factor of $\frac{\pi}{180}$[rad/$\circ$] like this:

>>> sin(30*pi/180)
1/2

The inverse trigonometric functions are $\sin^{-1}(x)\equiv \arcsin(x)\equiv$asin(x) and $\cos^{-1}(x)\equiv \arccos(x)\equiv$acos(x):

>>> asin(1/2)
pi/6
>>> acos(sqrt(3)/2)
pi/6

Recall that $\tan(x) \equiv \frac{\sin(x)}{\cos(x)}$. The inverse function of $\tan(x)$ is $\tan^{-1}(x)\equiv \arctan(x)\equiv$atan(x)

>>> tan(pi/6)
1/sqrt(3)             # == ( 1/2 )/( sqrt(3)/2 )
>>> atan( 1/sqrt(3) )
pi/6

Here are some trigonometric identities that sympy knows about:

>>> sin(x) == cos(x - pi/2)      
True

>>> simplify( sin(x)*cos(y)+cos(x)*sin(y) )
sin(x + y)
>>> trigsimp(sin(x)*cos(y)+cos(x)*sin(y))
sin(x + y)

>>> e = 2*sin(x)**2 + 2*cos(x)**2
>>> trigsimp(e)
2
>>> trigsimp(log(e))
log(2*sin(x)**2 + 2*cos(x)**2)
>>> trigsimp(log(e), deep=True)
log(2)

>>> simplify(sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4)
cos(4*x)/2 + 1/2

The function trigsimp does essentially the same job as simplify.

If instead of simplifying you want to expand a trig expression, you should use expand_trig because the default expand won't do trig functions:

>>> expand( sin(2*x) )         # = (sin(2*x)).expand()
sin(2*x)
>>> expand_trig( sin(2*x) )    # = (sin(2*x)).expand(trig=True)
2*sin(x)*cos(x)

The hyperbolic sine and cosine in sympy are denoted sinh and cosh respectively and sympy is smart enough to recognize them when simplifying expressions:

>>> simplify( (exp(x)+exp(-x))/2 )
cosh(x)
>>> simplify( (exp(x)-exp(-x))/2 )
sinh(x)

Recall that $x=\cosh(\mu)$ and $y=\sinh(\mu)$ are defined as $x$ and $y$ coordinates of a point on the the hyperbola with equation $x^2 - y^2 = 1$ and therefore satisfy the identity $\cosh^2 x - \sinh^2 x =1$:

>>> simplify( cosh(x)**2 - sinh(x)**2 )
1

>>> from sympy import I, re, im, Abs, arg, conjugate

Ever since Newton, the work “number” has referred to all of the following sets: the naturals $\mathbb{N}$, the integers $\mathbb{Z}$, the rationals $\mathbb{Q}$, and the real numbers $\mathbb{R}$. Each set of numbers is associated with a different class of equations. The natural numbers $\mathbb{N}$ appear as solutions of the equation $m+n=x$, where $m$ and $n$ are natural numbers (denoted $m,n \in \mathbb{N}$). The integers $\mathbb{Z}$ are the solutions to equations of the form $x+m=~n$, where $m,n \in \mathbb{N}$. The rational numbers $\mathbb{Q}$ are necessary to solve for $x$ in $mx=n$, with $m,n \in \mathbb{Z}$. To find the solutions of $x^2=2$, we need the real numbers $\mathbb{R}$. New types of numbers are requited for solving more complicated types of equations.

There are no real solutions to the quadratic equaiton $x^2=-1$, but if we define an imaginary number $i =\sqrt{-1}$ (denoted I in sympy) that satisfies this equation:

 >>> solve( x**2 + 1 , x)
 [-I, I]

The solutions are $x=-i$ and $x=i$, and indeed we can verify that $i^2+1=0$ and $(-i)^2+1=0$.

The set of complex numbers $\mathbb{C}$ is defined as follows $\{ a+bi \vert a,b \in \mathbb{R} \}$. Complex numbers have a real part and an imaginary part:

>>> z = 4 + 3*I
>>> z 
4 + 3*I
>>> re(z1)
4
>>> im(z1)
3 

Complex numbers can also be expressed in polar form: $z=|z|e^{i\theta}$. For a complex number $z=a+bi$, the quantity $|z|=\sqrt{a^2+b^2}$ is known as the absolute value of $z$, and $\theta$ is called the phase or the argument of $z$:

>>> Abs(z)
5
>>> arg(z)
atan(3/4)

The complex conjugate of the number $z=a+bi$ is the number $\bar{z} = a-bi$:

>>> conjugate( z )
4 - 3*I

Complex conjugation is important for computing the absolute value of $z$: $|z|=\sqrt{ z\bar{z} }$ and for division: $\frac{1}{z} \equiv \frac{\bar{z}}{|z|^2}$.

>>> from sympy import expand, rewrite

Euler's formula shows a important relation between the exponential function $e^x$ and the trigonometric functions $\sin(x)$ and $\cos(x)$: \[ e^{ix} = \cos x + i \sin x. \] To obtain this result in sympy, you must specify that the number $x$ is real and also tell expand that you're interested in complex expansions:

>>> x = symbols('x', real=True)
>>> exp(I*x).expand(complex=True)
cos(x) + I*sin(x) 
>>> re( exp(I*x) ) 
cos(x)
>>> im( exp(I*x) ) 
sin(x)

So basically, $\cos(x)$ is the real part of $e^{ix}$, and $\sin(x)$ is the imaginay part of $e^{ix}$. Whaaat? I know it's weird, but weird things are to be expected when one inputs imaginary values into a funciton!

Euler's identity is Euler's formula with $x=\pi$: $e^{i\pi} = \cos \pi + i\sin \pi = - 1$, which we can rewrite as $e^{i\pi} + 1=0$. Let's check that sympy knows about Euler's identity:

>>> (exp(I*pi) + 1).expand(complex=True)
0

Euler's formula is often used to rewrite the functions $\sin$ and $\cos$ in terms of complex exponentials. For example,

>>> (cos(x)).rewrite(exp)
exp(I*x)/2 + exp(-I*x)/2

The oprations of calculus are used to descibe the limit behaviour of functions, calculate the rates of chage, calculate areas and volumes, etc. In this section we'll learn about the sympy functions that correspond to the operations of calculus: limits, derivatives, integrals, sequences, and series.

from sympy import oo

The infinity symbol in sympy is denoted oo (two lowercase o's) in sympy. Infiniy is not a number but a process: the process of counting forever. Thus, $\infty + 1 = \infty$, $\infty$ is greater than any finite number, and $1/\infty$ is an infinitely small number. Sympy knows how to correctly treat infiniy in expressions:

>>> oo+1
oo
>>> oo > 5000
True
>>> 1/oo
0

from sympy import limit

We use limits to describe with mathematical precision infinitely large quantities, infinitely small quantities, and precedues with infinitely many of steps.

For example, the number $e=2.71828\ldots$ in defined as the limit $e \equiv \lim_{n\to \infty} \left( 1 + \frac{1}{n}\right)^{n}$:

>>> limit( (1+1/n)**n, n, oo)
E

This limit expression desribes the annual growth rate of a loan with interest rate $100\%$, and compoinding performed infinitely frequenctly.

Limits are also useful to describe the behaviour of functions. Consider for example the function $f(x)=\frac{1}{x}$. We can see what happens to $f(x)$ as $x$ goes to infitnity, and near $x=0$ using the limit command:

>>> limit( 1/x, x, oo)
0
>>> limit( 1/x, x, 0, dir="+")
oo
>>> limit( 1/x, x, 0, dir="-")
-oo

Indeed, as $x$ becomes larger and larger, the fraction $\frac{1}{x}$ becomes smaller ans smaller, and in the limit $\lim_{x\to \infty} \frac{1}{x} = 0$. On the other hand, when $x$ taks on smaller and smaller positive values, the expression $\frac{1}{x}$ becmes infinite: $\lim_{x \to 0^+} \frac{1}{x} = \infty$. When $x$ approaches $0$ from the left, we have $\lim_{x \to 0^-} \frac{1}{x} = -\infty$. If these calculations are not clear to you, you should study the the graph of $f(x)=\frac{1}{x}$.

Here are some other examples of limits

>>> limit(sin(x)/x, x, 0)
1

Limits are also used to define the derivative and the integral operations.

The derivative function, denoted $f'(x)$, $\frac{d}{dx}f(x)$, $\frac{df}{dx}$, or $\frac{dy}{dx}$, describes the rate of change of the function $f(x)$. The sympy function diff computes the derivative of any expression:

>>> diff(x**3, x)
3*x**2

The differentiation operation knows about the product rule $[f(x)g(x)]^\prime=f^\prime(x)g(x)+f(x)g^\prime(x)$, the chain rule $f(g(x))' = f'(g(x))g'(x)$, and the quatient rule $\left[\frac{f(x)}{g(x)}\right]^\prime = \frac{f'(x)g(x) - f(x)g'(x)}{g(x)^2}$:

>>> diff( x**2*sin(x), x )
2*x*sin(x) + x**2*cos(x)
>>> diff( sin(x**2),   x )
cos(x**2)*2*x
>>> diff( x**2/sin(x), x )
(2*x*sin(x) - x**2*cos(x))/sin(x)**2

The second derivative of a function can be obtained using the diff(f, x, 2):

>>> diff(x**3, x, 2)       # same as diff(diff(x**3, x), x)
6*x

Recall the exponential function $f(x)=e^x$ is special because it is equal to its derivative:

>>> diff( exp(x), x)     # same as diff( E**x, x)
exp(x)                   # same as E**x

A differential equation is an equation that relates some unknown function $f(x)$ to its derivatives. An example of a differential equation is $f'(x)=f(x)$. What is the function $f(x)$ which is equal to its derivative? You can either try to guess what $f(x)$ is or use the dsolve function:

>>> x = symbols('x')
>>> f = symbols('f', cls=Function)
>>> dsolve( f(x) - diff(f(x),x), f(x) )
f(x) == C1*exp(x)

We'll discuss dsolve in more detail in the section on mechanics.

The tangent line to the function $f(x)$ at $x=x_0$ corresponds to the line that passes through the point $(x_0, f(x_0))$ and has the same slope as the function at that point. The tangent line to the function $f(x)$ at the point $x=x_0$ is described by the equation \[ T_1(x) = f(x_0) \ + \ f'(x_0)(x-x_0). \] What is the equation of a the tangent line to $f(x)=\frac{1}{2}x^2$ at $x_0=1$?

>>> f = S('1/2')*x**2
>>> f
x**2/2
>>> df = diff(f,x)
>>> df
x
>>> T_1 = f.subs({x:1}) + df.subs({x:1})*(x - 1)
>>> T_1
x - 1/2
>>> T_1.subs({x:1}) == f.subs({x:1})
True
>>> diff(T_1,x).subs({x:1}) == diff(f,x).subs({x:1})
True

The tangent line $T_1(x)$ has the same value as $f(x)$ and the same slope as $f(x)$ at $x=1$.

Optimization is about choosing an input for a function $f(x)$ that results in the best value for $f(x)$. The best value usually means the maximum value (if the function measures something desirable like profits) or the minimum value (if the function describes something undesirable like costs).

The derivative $f'(x)$ encodes the information about the slope of $f(x)$. The second derivative $f^{\prime\prime}(x)$ encodes the information about the curvature of $f(x)$. Positive curvature means the function looks like $x^2$, negative curvature means the function looks like $-x^2$. The critical points of a function $f(x)$ are the solutions to the equation $f'(x)=0$. Each critical point is a candidate to be either a maximum or a minimum of the function.

As an example, let's find the where the maximum of the function $f(x)=x^3-2x^2+x$ occurs on the interval $x \in [0,1]$. You should plot to see what it looks like.

>>> x = Symbol('x')
>>> f = x**3-2*x**2+x
>>> diff(f, x)
3*x**2 - 4*x + 1
>>> sols = solve( diff(f,x),  x)
>>> sols
[1/3, 1]
>>> diff(diff(f,x), x).subs( {x:sols[0]} )
-2
>>> diff(diff(f,x), x).subs( {x:sols[1]} )
2

The point $x=\frac{1}{3}$ is a local maximum because it is critical point of $f(x)$ where the curvature is negative curvature, meaning $f(x)$ looks like the peak of a mountain. The maximum value of $f(x)$ is $f\!\left(\frac{1}{3}\right)=\frac{4}{27}$. The point $x=1$ is a local minimum because it is a critical point with positive curvature, meaning $f(x)$ looks like the bottom of a valley there.

The \emph{integral} of $f(x)$ corresponds to the computation of the area under the graph of $f(x)$. The area under $f(x)$ between the points $x=a$ and $x=b$ is denoted as follows: \[ A(a,b) = \int_a^b f(x) \: dx. \] The integral function $F(c)$ corresponds to the area calculation as a function of the upper limit of integration: \[ F(c) \equiv \int_0^c \! f(x)\:dx\,. \] The area under $f(x)$ between $x=a$ and $x=b$ is obtained by calculating the change in the integral function: \[ A(a,b) = \int_a^b \! f(x)\:dx = F(b)-F(a). \]

In sympy we use the integrate(f, x) to obtain the integral function $F(x) = \int_0^x f(u)\,du$:

>>> integrate(x**3, x)
x**4/4
>>> integrate(sin(x), x)
-cos(x)
>>> integrate(ln(x), x)
x*log(x) - x

To compute the area under $f(x)$ between $x=a$ and $x=b$ use integrate(f, (x,a,b)) and sympy will compute the integral $\int_a^b f(x) \, x$ for you:

>>> integrate(x**3, (x,0,1))    
1/4

Note we can obtain the same result by fist finding the integral function $F(x)$, and then use the formula $A(a,b) = F(b) - F(a)$:

>>> F = integrate(x**3, x)
>>> F.subs({x:1}) - F.subs({x:0})   
1/4

Integrals correspond to signed area calculations:

>>> integrate(sin(x), (x, 0, pi))
2
>>> integrate(sin(x), (x, pi, 2*pi))
-2
>>> integrate(sin(x), (x, 0, 2*pi))
0

During the first half of its $2\pi$ cycle, the function $\sin(x)$ is above the $x$-axis and has positive area under the curve. During the second half of its cycle (from $x=\pi$ to $x=2\pi$), the $\sin(x)$ is below the $x$-axis so contributes negative area. Draw the graph of $\sin(x)$ to see what is going on.

The integral is the “inverse operation” of the derivative. If you perform the integral operation followed by the derivative operation on some function, you'll obtain the same function: \[ \left(\frac{d}{dx} \circ \int dx \right) f(x) \equiv \frac{d}{dx} \int_c^x f(u)\:du = f(x). \]

>>> f = x**2
>>> F = integrate(f, x)
>>> F
x**3/3           # + C
>>> diff(F,x)
x**2

Alternately, if you compute the derivative of a function followed by the integral, you will obtain the original function $f(x)$ (up to a constant): \[ \left( \int dx \circ \frac{d}{dx}\right) f(x) \equiv \int_c^x f'(u)\;du = f(x) + C. \]

>>> f = x**2
>>> df = diff(f,x)
>>> df
2*x
>>> integrate(df, x)
x**2    # + C

The fundamental theorem of calculus is important because it tells us how to solve differential equations. If we have to solve for $f(x)$ in the differential equation $\frac{d}{dx}f(x) = g(x)$, we can take the integral on both sides of the equation to obtain $f(x) = \int g(x)\,dx + C$.

Sequences are functions that take integers as inputs. Instead of continuous inputs $x\in \mathbb{R}$, a sequences take natural numbers $n\in\mathbb{N}$ as inputs. We denote the sequences as $a_n$ instead of the usual function notation $a(n)$.

We define a sequence by specifying an expression for it's $n$th term:

>>> a_n = 1/n
>>> b_n = 1/factorial(n)

We can see the value of the $n$th term by subsituting the desrired value of $n$:

>>> a_n.subs({n:5})
1/5

We can use the python list comprehention syntaxt [item for item in list] to print the sequence values for some range of indices:

>>> [ a_n.subs({n:i}) for i in range(0,8) ]
[oo, 1, 1/2, 1/3, 1/4,  1/5,   1/6,   1/7]  
>>> [ b_n.subs({n:i}) for i in range(0,8) ]
[1,  1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040]

Observe that $a_n$ is not properly defined for $n=0$ since $\frac{1}{0}$ is a division-by-zero error, so to be precise we should say $a_n$s input space is only the positive naturals $a_n:\mathbb{N}^+ \to \mathbb{R}$. Observe also how quickly the factorial function grows: $n!=1\cdot2\cdot3\cdots(n-1)\cdot n$.

We're often interested in calculating the limits of sequences as $n\to \infty$. What happens to the terms in the sequence when $n$ takes on very large values?

>>> limit(a_n, n, oo)
0
>>> limit(b_n, n, oo)
0

Both $a_n$ and $b_n$ converge to $0$ as $n\to\infty$.

Many important math quantities are defined as the limit of some expression. An interesting example to consider is the number $\pi$, which is defined as the area of a circle of radius $1$. We can approximate the area of the unit circle by drawing a many-sided regular polygon around the circle. Splitting the $n$-sided regular polygon into identical triangular splices, we can obtain a formula for its area $A_n$ (see problem~REFref{mathprob:octagon-approx-to-circle}). In the limit as $n\to \infty$, the $n$-sided-polygon approximation to the area of the unit-circle becomes exact:

>>> A_n = n*tan(2*pi/(2*n))
>>> limit(A_n, n, oo)
pi

Suppose we're given a sequence $a_n$ and we want to compute the sum of all the values in this sequence. Series are sums of sequences. Summing the values of a sequence $a_n:\mathbb{N}\to \mathbb{R}$ is analogous to taking the integral of a funciton $f:\mathbb{R}\to \mathbb{R}$.

In sympy we can use the summation notation whose syntax is analogous to the intagrate command:

>>> a_n = 1/n
>>> b_n = 1/factorial(n)
>>> summation(a_n, [n, 1, oo])
oo
>>> summation(b_n, [n, 0, oo])
E

We say the series $\sum a_n$ diverges (or is divergent) while the series $\sum b_n$ converges (or is convergent). As we sum together more and more terms of the sequence $b_n$, the total becomes closer and closer to some finite number. In this case, the infinite sum $\sum_{n=0}^\infty \frac{1}{n!}$ converges to the number $e=2.71828\ldots$.

Using the summation command is useful because is allows us to compute infinite summations, but for most practical applications we don't need to take an infinite number of terms in a series to obtain a good approximation. Indeed, this is what is so neat about series: they represent a great way to obtain approxmations.

Using the standard python tools, we can obtain an approxmation to $e$ that is accurate to six decimals as follows:

>>> import math
>>> def b_nf(n): 
    return 1.0/math.factorial(n)
>>> sum([b_nf(n) for n in range(0,10)] )
2.718281 52557319
>>> E.evalf()
2.718281 82845905       # true value

But wait there is more! Not only can series be used to approximate numbers, series can also be used to approximate functions.

A power series is a series whose terms contain different powers of the variable $x$. The $n$th term in a a power series is a function of both the sequence index $n$ and the input variable $x$.

For example the power series of the function $f(x)=\exp(x)=e^x$ is \[ \exp(x) \equiv 1 + x + \frac{x^2}{2} + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots = \sum_{n=0}^\infty \frac{x^n}{n!}. \] This is, IMHO, one of the most important ideas in calculus: you can compute the value of $\exp(5)$ by taking the infinite sum of the terms in the power series with $x=5$:

>>> exp_xn = x**n/factorial(n)
>>> summation( exp_xn.subs({x:5}), [n, 0, oo] ).evalf()
148.413159102577
>>> exp(5).evalf()
148.413159102577        # the true value

Note that sympy is actually smart enough to recognize the infinite series you're computing corresponds to $\exp(5)$:

>>> summation( exp_xn.subs({x:5}), [n, 0, oo])
exp(5)

That taking as few as 35 terms in the series is sufficient to obtain a good approximation, so series are not some abstract thing for mathematicians but a practical trick you can when you code:

>>> summation( exp_xn.subs({x:5}), [n, 0, 35]).evalf()
148.413159102577

>>> import math                  # redo using sympy-less python 
>>> def exp_xnf(x,n): 
    return x**n/math.factorial(n)
>>> sum([exp_xnf(5,i) for i in range(0,35)] ).evalf()
148.413159102577

The coefficients in the power series of a function (also known as the Taylor series or Maclaurin series) depend on the value of the higher derivatives of the function. The formula for the $n$th term of the Taylor series of $f(x)$ expanded at $x=0$ is $a_n(x) = \frac{f^{(n)}(0)}{n!}x^n$, where $f^{(n)}(0)$ is the value of the $n$th derivative of $f(x)$ evaluated at $x=0$.

The sympy function series is a convenient way to obtain the series of any function. Calling series(expr,var,at,nterms) will show you the series of expr expanded around var=at up to nterms:

>>> series( sin(x), x, 0, 8)
x - x**3/6 + x**5/120 - x**7/5040 + O(x**8)
>>> series( cos(x), x, 0, 8)
1 - x**2/2 + x**4/24 - x**6/720 + O(x**8)
>>> series( sinh(x), x, 0, 8)
x + x**3/6 + x**5/120 + x**7/5040 + O(x**8)
>>> series( cosh(x), x, 0, 8)
1 + x**2/2 + x**4/24 + x**6/720 + O(x**8)

Note the power series of $\sin$ and $\sinh$ contain only odd powers of $x$ while the power series of $\cos$ and $\cosh$ contain only even powers.

Some functions are not defined at $x=0$, so we expand them at a different value of $x$. For example the power series of $\ln(x)$ expanded at $x=1$ is

>>> series(ln(x), x, 1, 6)     # Taylor series of ln(x) at x=1
x - x**2/2 + x**3/3 - x**4/4 + x**5/5  + O(x**6)

In this case the result retuned by sympy is a bit misleading, the Taylor series of $\ln(x)$ expanded at $x=1$ actually contains terms of the form $(x-1)^n$: \[ \ln(x) = (x-1) - \frac{(x-1)^2}{2} + \frac{(x-1)^3}{3} - \frac{(x-1)^4}{4} + \frac{(x-1)^5}{5} + \cdots. \] You can verify this is the correct formula by substituting $x=1$. The sympy developers have simply chosen not to show this fact, because when dealing with series in general we're mostly interested in the coefficients.

Instead of expanding $\ln(x)$ around $x=1$, we can obtain an equivalent expression if we expand $\ln(x+1)$ around $x=0$:

>>> series(ln(x+1), x, 0, 6)   # Maclaurin series of ln(x+1)
x - x**2/2 + x**3/3 - x**4/4 + x**5/5 + O(x**6)

The term Taylor series applies to all series expansions of functions. The term Maclaurin series applies specifically to Taylor series expansions at $x=0$, which are generally easier to work with.

A vector $\vec{v} \in \mathbb{R}^n$ is an $n$-tuple of real numbers. For example, consider a vector that has three components: \[ \vec{v} = (v_1,v_2,v_3) \ \in \ (\mathbb{R},\mathbb{R},\mathbb{R}) \equiv \mathbb{R}^3. \] To specify the vector $\vec{v}$, we specify the values for its three components $v_1$, $v_2$, and $v_3$.

A matrix $A \in \mathbb{R}^{m\times n}$ is a rectangular array of real numbers with $m$ rows and $n$ columns. A vector is a special type of matrix; we can think of a vector $\vec{v}\in \mathbb{R}^n$ either as a row vector (a $1\times n$ matrix) or a column vector (an $n \times 1$ matrix). Because of this equivalence between vectors and matices, there is no need for a special vector object in sympy, and Matrix objects are used for vectors as well.

Here is how to define a vector in sympy and a demo of some basic vectors operations:

>>> u = Matrix([4,5,6])  # defines new col. vector u
>>> u    
Matrix([                 # col vec = 3x1 matrix 
[4],
[5],
[6]])
>>> u.T                  # use the transpose operation to 
Matrix([[4, 5, 6]])      # convert a col vec to a row vec

>>> u[0]                 # 0-based indexing for vector entries 
4
>>> u.norm()             # length of u 
sqrt(77)
>>> uhat = u/u.norm()    # unit length vec in dir of u
>>> uhat.T
[4/sqrt(77), 5/sqrt(77), 6/sqrt(77)]
>>> uhat.norm()
1

The dot product of two $3$-vectors $\vec{u}$ and $\vec{v}$ is \[ \vec{u}\cdot\vec{v}\equiv u_xv_x+u_yv_y+u_zv_z \equiv \|\vec{u}\|\|\vec{v}\|\cos(\varphi) \quad \in \mathbb{R}, \] where $\varphi$ is the angle between the two vectors.

>>> u = Matrix([4,5,6])
>>> v = Matrix([-1,1,2])
>>> u.dot(v)
13

We can combine the algebraic and the geometric formulas for the dot product to obtain the cosine of the angle between the vectors \[ \cos(\varphi) = \frac{ \vec{u}\cdot\vec{v} }{ \|\vec{u}\|\|\vec{v}\| } = \frac{ u_xv_x+u_yv_y+u_zv_z }{ \|\vec{u}\|\|\vec{v}\| }, \] and use the acos function to find the measure of the angle:

>>> acos(u.dot(v)/(u.norm()*v.norm())).evalf()
0.921263115666387       # in radians 

Just looking at the coordinates of the vectors $\vec{u}$ and $\vec{v}$, it would be difficult to get an intuition of their relative direction, but thanks to the above calculation we know the angle between the vectors is $52.76^\circ$, which means they kind of point in the same direction. Vectors that are at an angle $\varphi=90^\circ$ are called orthogonal, meaning at right angles with each other. Vectors for which $\varphi > 90^\circ$ have a negative dot product.

The notion of “angle between vectors” generalizes to vectors with any number of dimensions. The dot product for $n$-dimensional vectors is: $\vec{u}\cdot\vec{v}=\sum_{i=1}^n u_iv_i$. This means you can talk about “the angle between” 1000-dimensional vectors. That's pretty crazy if you think about it—there is no way we could possibly “visualize” 1000-dimensional vectors, yet given two such vectors we can tell if they point mostly in the same direction, in perpendicular directions, or mostly in opposite directions.

The dot product is commutative $\vec{u}\cdot\vec{v} = \vec{v}\cdot\vec{u}$:

>>> u.dot(v) == v.dot(u)
True

Dot products are very useful for computing projections. Assume you have two $\vec{u}$ and $\vec{n}$ and you want to find the component of $\vec{u}$ that points in the $\vec{n}$ direction. The following formula based on the dot product will give you the answer: \[ \Pi_{\vec{n}}( \vec{u} ) = \frac{ \vec{u} \cdot \vec{n} }{ \| \vec{n} \|^2 } \vec{n}. \]

This is how to implement this formula in sympy:

>>> u = Matrix([4,5,6])
>>> n = Matrix([1,1,1])
>>> (u.dot(n) / n.norm()**2)*n
Matrix([5, 5, 5])

In the case where the direction vector $\hat{n}$ is of unit length $\|\hat{n}\| = 1$, the projection formula simplifies to $(\vec{u}\cdot\hat{n})\hat{n}$.

Consider now the plane $P$ defined by the equation $(1,1,1)\cdot[(x,y,z)-(0,0,0)]=0$. A plane is a two dimensional subspace of $\mathbb{R}^3$. We can decompose any vector $\vec{u} \in \mathbb{R}^3$ into two parts $\vec{u}=\vec{v} + \vec{w}$ such that $\vec{v}$ lies inside the plane and $\vec{w}$ is perpendicular to the plane (parallel to $\vec{n}=(1,1,1)$).

To obtain the perpendicular-to-$P$ component of $\vec{u}$, compute the projection of $\vec{u}$ in the direction $\vec{n}$:

>>> w = (u.dot(n) / n.norm()**2)*n
Matrix([5, 5, 5])  

To obtain the in-the-plane-$P$ component of $\vec{u}$, start with $\vec{u}$ and subtract the perpendicular-to-$P$ part:

>>> v = u - (u.dot(n)/n.norm()**2)*n      # same as u - w
Matrix([ -1, 0, 1])

You should check on your own that $\vec{v}+\vec{w}=\vec{u}$ as claimed.

The cross product, denoted $\times$, takes two vectors as inputs and produces another vector as the output. The cross products of individual basis elements are defined as follows: \[ \hat{\imath}\times\hat{\jmath} =\hat{k}, \ \ \ \hat{\jmath}\times\hat{k} =\hat{\imath}, \ \ \ \hat{k}\times \hat{\imath}= \hat{\jmath}. \]

Here is how to compute the cross product of two vectors in sympy:

>>> u = Matrix([4,5,6])
>>> v = Matrix([-1,1,2])
>>> u.cross(v)
Matrix([[4, -14, 9]])

The vector $\vec{u}\times \vec{v}$ is orthogonal to both $\vec{u}$ and $\vec{v}$. The norm of the cross product $|\vec{u}\times \vec{v}|$ is proportional to the lengths of the vectors and the sine of the angle between them:

(u.cross(v).norm()/(u.norm()*v.norm())).n()
0.796366206088088    # = sin(0.921..) 

The name “cross product” is well suited for this operation, because the entry for each dimension in the cross product is computed by “cross-multiplying” the entries in the other dimensions: \[ \vec{u}\times\vec{v}= \left( u_yv_z-u_zv_y, \ u_zv_x-u_xv_z, \ u_xv_y-u_yv_x \right). \]

By individual symbols for the entries of two vectors, we can ask sympy to show us the cross-product formula:

>>> u1,u2,u3 = symbols('u1:4')
>>> v1,v2,v3 = symbols('v1:4')
>>> Matrix([u1,u2,u3]).cross(Matrix([v1,v2,v3]))
[ u2*v3 - u3*v2]
[-u1*v3 + u3*v1]
[ u1*v2 - u2*v1]

The dot product is anti-commutative $\vec{u}\times\vec{v} = -\vec{v}\times\vec{u}$:

>>> u.cross(v)
Matrix([[4, -14, 9]])
>>> v.cross(u)
Matrix([[-4, 14, -9]])

Watch out for this: the product of two numbers and the dot product of two vectors are commutative operations. The cross product isn't: $\vec{u}\times\vec{v} \neq \vec{v}\times\vec{u}$.

The module called sympy.physics.mechanics contains elaborate tools for describing mechanical systems, reference frames, forces, torques, etc. These specialized functions are not necessary for a first-year mechanics course though. The basic sympy functions like solve (to solve equations), and the vector functions are powerful enough to solve most Newtonian mechanics problems.

The net force acting on an object is the sum of all the external forces acting on it.

Let $x(t)$ denote the position of an object, $v(t)\equiv x'(t)$ denote its velocity, and $a(t)\equiv v'(t)$ denote its acceleration. Together $x(t)$, $v(t)$, and $a(t)$ are knowns as the \emph{equations of motion} of the object. The equations of motion are related by the derivative operation: \[ a(t) \overset{\frac{d}{dt} }{\longleftarrow} v(t) \overset{\frac{d}{dt} }{\longleftarrow} x(t). \] Assuming we know the initial position $x_i\equiv x(0)$ and the initial velocity of the object $v_i\equiv v(0)$, we want to find $x(t)$ for all later times. We can do this starting from the dynamics of the problem—the forces acting on the object.

Newton's second law $\vec{F}_{\textrm{net}} = m\vec{a}$ states that a net force $\vec{F}_{\textrm{net}}$ applied on an object of mass $m$ produces acceleration $\vec{a}$. Thus, we can obtain an objects acceleration calculating if we know the net force acting on it. Starting from the knowledge of $a(t)$, we can obtain $v(t)$ by integrating then find $x(t)$ by integrating $v(t)$: \[ a(t) \ \ \ \overset{v_i+ \int\!dt }{\longrightarrow} \ \ \ v(t) \ \ \ \overset{x_i+ \int\!dt }{\longrightarrow} \ \ \ x(t). \] The reasoning follows from from the fundamental theorem of calculus: if $a(t)$ represents the change in $v(t)$, then the total of $a(t)$ accumulated between $t=t_1$ and $t=t_2$ is equal to the total change in $v(t)$ between these times: $\Delta v = v(t_2) - v(t_1)$. Similarly, the integral of $v(t)$ from $t=0$ until $t=\tau$ is equal to $x(\tau) - x(0)$.

Let's analyze the case where the net force on object is constant. A constant force causes a constant acceleation $a = \frac{F}{m} = \textrm{constant}$. The acceleration function is constant over time $a(t)=a$. We find $v(t)$ and $x(t)$ as follows:

»> t, v_i, x_i = symbols('t v_i x_i')

>>> a = Symbol("a")
>>> v = v_i + integrate(a, (t, 0,t) )
>>> v
a*t + v_i
>>> x = x_i + integrate(v, (t, 0,t) )
>>> x
a*t**2/2 + v_i*t + x_i

You may remember these equations from your high school mechanics class. These are known as the uniform accelerated motion (UAVM) equations: \begin{align} a(t) = a,
v(t) = v_i + at,
x(t) = x_i + v_it + \frac{1}{2}a^2. \end{align} In high school, you probably had to memorize these equations. Now you know how to derive them yourself starting from first principles.

For the sake of completeness, we'll now derive the fourth UAM equation, which relates the findal velocity to the initial velocity, the displacement, and the acceleration without reference to time:

>>> (v*v).expand()
a**2*t**2 + 2*a*t*v_i + v_i**2
>>> ((v*v).expand() - 2*a*x).simplify()
-2*a*x_i + v_i**2

The above calculation shows $v_f^2 - 2ax_f = -2ax_i + v_i^2$. After moving the term $2ax_f$ to the other side of the equation, we obtain \begin{align} (x(t))^2 \ = \ x_f^2 = v_i^2 + 2a\Delta x \ = \ v_i^2 + 2a(x_f-x_i). \end{align} The fourth equation is important for practical purposes because allows you to solve physics problems in a time-less manner. It's also important for theoretical purposes since it resembles an energy calculation. Indeed, if you multiply the fourth equation by $\frac{1}{2}m$, you'll obtain conservation of energy equation: $K_f = K_i + mg\Delta h$.

The procedure $a(t) \ \overset{v_i+ \int\!dt }{\longrightarrow} \ v(t) \ \overset{x_i+ \int\!dt }{\longrightarrow} \ x(t)$ can be used to obtain the position function $x(t)$ even when the acceleration is not constant in time. For example, suppose we have $a(t)=\sqrt{k t}$, then the position as a function of time will be

>>> k = Symbol("k")
>>> a = sqrt(k*t)
>>> v = v_i + integrate(a, (t, 0,t) )
>>> v
v_i + 2*(k*t)**(3/2)/(3*k)
>>> x = x_i + integrate(v, (t, 0,t) )
>>> x
x_i + v_i*t + (4/15)*(k*t)**(5/2)/k**2 

Find the position function of an object at time $t =3$[s], if it starts from $x_i=20$[m], with $v_i=10$[m/s] and undergoes a constant acceleration of $a=5$[m/s$^2$].

>>> x_i = 20   # initial position
>>> v_i = 10   # initial velocity
>>> a   = 5    # acceleration (constant during motion)
>>> x = x_i + integrate( v_i+integrate(a,(t,0,t)),   (t,0,t) )   
>>> x
5*t**2/2 + 10*t + 20
>>> x.subs({t:3})
145/2   # [m]

Do you see what I mean when math and physics knowledge combined with computer skills is like superpowers?

Instead of working with the kinematic equations of motion $x(t)$, $v(t)$, and $a(t)$, we can solve physics problems using energy calculations. A key connection between the world of forces and the world of energy is the concept of potential energy. If you move an object against a conservative force (think raising a ball in the air against the force of gravity), you can think of the work you did as stored as potential energy.

For each force $\vec{F}(x)$ there is a corresponding potential energy function $U_F(x)$. The change in potential energy associated with the force $\vec{F}$ is defined as the negative of the work done by the force during the displacement: $U_F(x) = - W = - \int_{\vec{d}} \vec{F}\cdot d\vec{x}$.

The potential energies associated with the gravitational force $\vec{F}_g = -mg\hat{\jmath}$ and the force of a spring $\vec{F}_s = -k\vec{x}$ are

>>> x, y, z, t = symbols('x y z t')
>>> m, g, k, h = symbols('m g k h')
>>> F_g = -m*g
>>> U_g = - integrate( F_g, (y,0,h) )
>>> U_g
g*h*m
>>> F_s = -k*x
>>> U_s = - integrate( F_s, (x,0,x) )
>>> U_s
k*x**2/2

Note the negative sign in the definition of potential, which gets canceled because the dot product $\vec{F}\cdot d\vec{x}$ is negative in both calculations. When the force acts in the opposite direction of the displacement, the work done by the force is negative.

from sympy import Function, dsolve

The force exerted by a spring is given by the formula $F=-kx$. If the only force acting on a mass $m$ is the force of a spring, we can use Newton's second law to obtain the following equation \[ F=ma \quad \Rightarrow \quad -kx = ma \quad \Rightarrow \quad -kx(t) = m\frac{d^2}{dt^2}\Big[x(t)\Big]. \] The motion of a mass-spring system is described by the differential equation $\frac{d^2}{dt^2}x(t) + \omega^2 x(t)=0$, where the constant $\omega = \sqrt{\frac{k}{m}}$ is the angular frequency. We can find the position function $x(t)$ using the dsolve method:

>>> t = Symbol('t')                   # time            t
>>> x = Function('x')                 # position function x(t)
>>> w = Symbol('w', positive=True)    # angular frequency w

>>> sol = dsolve( diff(x(t),t,t) + w**2*x(t), x(t) )
>>> sol
x(t) == C1*sin(w*t) + C2*cos(w*t)  
>>> x = sol.rhs           
>>> x 
C1*sin(w*t) + C2*cos(w*t)

Note the solution $x(t)=C_1\sin(\omega t)+C_2 \cos(\omega t)$ is equivalent to $x(t) = A\cos(\omega t + \phi)$, which is more commonly used to describe simple harmonic motion. We can use the expand function with the argument trig=True to convince ourselves:

>>> A, phi = symbols("A phi")
>>> (A*cos(w*t - phi)).expand(trig=True)
A*sin(phi)*sin(w*t) + A*cos(phi)*cos(w*t)

If we define $C_1=A\sin(\phi)$ and $C_2=A\cos(\phi)$, we obtain the form $x(t)=C_1\sin(\omega t)+C_2 \cos(\omega t)$ that sympy found.

We can also verify that the total energy of the mass-spring system is conserved: $E_T(t) = U_s(t) + K(t) = \textrm{constant}$:

>>> x = sol.rhs.subs({"C1":0,"C2":A}) 
>>> x
A*cos(t*w)
>>> v = diff(x, t)
-A*w*sin(t*w)
>>> E_T = (0.5*k*x**2 + 0.5*m*v**2).simplify()
>>> E_T
0.5*A**2*(k*cos(t*w)**2 + m*w**2*sin(t*w)**2)
>>> E_T.subs({k:m*w**2}).simplify()
0.5*m*(w*A)**2                     # = K_max
>>> E_T.subs({w:sqrt(k/m)}).simplify()
0.5*k*A**2                         # = U_max 

from sympy import Matrix

A matrix $A \in \mathbb{R}^{m\times n}$ is a rectangular array of real numbers with $m$ rows and $n$ columns. To specify a matrix $A$, we specify the values for its $mn$ components $a_{11}, a_{12}, \ldots, a_{mn}$.

>>> A = Matrix( [[2,-3,-8, 7],
                 [-2,-1,2,-7],
                 [1 ,0,-3, 6]])

Use the square brackets to access the matrix elements or to obtain a submatrix:

>>> A[0,1]
-3
>>> A[0:2, 0:3] 
[ 2, -3, -8]
[-2, -1,  2]

Some commonly used matrices can be created with shortcut methods:

>>> eye(2)
[1, 0]
[0, 1]
>>> zeros((2, 3))
[0, 0, 0]
[0, 0, 0]

The transpose operation flips the matrix through its diagonal:

>>> A.transpose()       # the same as A.T
[ 2, -2,  1]
[-3, -1,  0]
[-8,  2, -3]
[ 7, -7,  6]

Recall that the transpose is used to convert row vectors into column vectors and vice versa.

The Gauss–Jordan elimination procedure is a sequence of row operations you can perform on any matrix to bring it to its reduced row echelon form (RREF). In sympy, matrices have a rref method that computes the RREF:

>>> A = Matrix( [[2,-3,-8, 7],
                 [-2,-1,2,-7],
                 [1 ,0,-3, 6]])
>>> A.rref()
([1, 0, 0,  0]  # RREF of A
 [0, 1, 0,  3]                    # locations of pivots
 [0, 0, 1, -2],                   [0, 1, 2]              )

Note the rref() method returns a tuple of values, the first value is the RREF of $A$, while the second tells you the indices of the leading ones (also known as pivots) in the RREF of $A$. To get just the RREF of $A$, select the 0th entry form the tuple:

>>> Arref = A.rref()[0]
>>> Arref
[1, 0, 0,  0]
[0, 1, 0,  3]
[0, 0, 1, -2]

The fundamental spaces of a matrix are its column space $\mathcal{C}(A)$, its null space $\mathcal{N}(A)$, and its row space $\mathcal{R}(A)$. These vector spaces are important when you consider the matrix product $A\vec{x}=\vec{y}$ as “applying” the linear transformation $T_A:\mathbb{R}^4 \to \mathbb{R}^3$ to an input vector $\vec{x} \in \mathbb{R}^4$ to produce the output vector $\vec{y} \in \mathbb{R}^3$.

The set of linear transformations $T_A:\mathbb{R}^n \to \mathbb{R}^m$ (vector functions) is equivalent to the set of $m\times n$ matrices. This is one of the fundamental ideas in linear algebra. You can think of $T_A$ as the abstract description of the transformation and $A \in \mathbb{R}^{m\times n}$ as a concrete implementation of $T_A$. By this equivalence, the fundamental spaces of a matrix $A$ tell us facts about the domain and image of the linear transformation $T_A$. The columns space $\mathcal{C}(A)$ is the same as the image space space $\textrm{Im}(T_A)$ (the set of all possible outputs). The null space $\mathcal{N}(A)$ is the same as the null space $\textrm{Null}(T_A)$ (the set of inputs that $T_A$ maps to the zero vector). The row space $\mathcal{R}(A)$ is the orthogonal complement of the null space. Input vectors in the row space of $A$ are in one-to-one correspondence with the output vectors in the column space of $A$.

Okay, enough theory! Let's see how to compute the fundamental spaces of the matrix $A$ defined above. The non-zero rows in the reduced row echelon form of $A$ form a basis for its row space:

>>> [ A.rref()[0][r,:] for r in A.rref()[1] ]   # R(A)
[  [1, 0, 0, 0],  [0, 1, 0, 3],  [0, 0, 1, -2]   ]

The column space of $A$ is the span of the columns of $A$ that contain the pivots in the reduced row echelon form of $A$:

>>> [ A[:,c] for c in  A.rref()[1] ]   # C(A)
[  [ 2]    [-3]   [-8]
   [-2],   [-1],  [ 2]
   [ 1]    [ 0]   [-3]  ]
   

Note we took the columns form the original matrix $A$ and not its RREF.

To find the null space of $A$ call its nullspace method:

>>> A.nullspace()          # N(A)
[  [0, -3, 2, 1]  ]

The determinant of a matrix, denoted $\det(A)$ or $|A|$, is a particular way to multiply the entries of the matrix to produce a single number.

>>> M = Matrix( [[1, 2, 3], 
                 [2,-2, 4],
                 [2, 2, 5]] )
>>> M.det()   
2

We use determinants for all kinds of tasks: to compute areas and volumes, to solve systems of equations, to check whether a matrix is invertible or not, etc.

For every invertible matrix $A$, there exists an inverse matrix $A^{-1}$ which undoes the effect of $A$. The cumulative effect of the product of $A$ and $A^{-1}$ in any order is the identity matrix: $AA^{-1}= A^{-1}A=\mathbbm{1}$.

>>> A = Matrix( [[1,2], 
                 [3,9]] ) 
>>> A.inv()                        
[ 3, -2/3]
[-1,  1/3]
>>> A.inv()*A
[1, 0]
[0, 1]
>>> A*A.inv()
[1, 0]
[0, 1]

The fundamental eigenvalue equation is \[ A\vec{e}_\lambda =\lambda\vec{e}_\lambda, \] where $\lambda$ is an eigenvalue of $A$ and $\vec{e}_\lambda$ is the corresponding eigenvector. When we multiply $A$ on of its eigenvectors $\vec{e}_\lambda$, the result is the same vector scaled by the constant $\lambda$.

The matrix $A$ can be written as the product of three matrices: $A=Q\Lambda Q^{-1}$. This is called the eigendecomposition of $A$. The matrix $Q$ contains the eigenvectors of $A$ as columns. The matrix $\Lambda$ contains the eigenvalues of $A$ on its diagonal.

To diagonalize a matrix $A$ is to find its $Q$ and $\Lambda$ matrices.

>>> A = Matrix([ [ 9, -2],
                 [-2,  6]])
>>> A.eigenvals()
{5: 1, 10: 1}   # eigenvalues 5 and 10 with multiplicity 1
>>> A.eigenvects()   
[(5, 1, [ 1]  
        [ 2]  ),      (10, 1, [-2]
                              [ 1]  )]
>>> Q, Λ = A.diagonalize()
>>> Q               # the matrix of eigenvectors 
[1,  -2]            # as columns 
[2,   1]
>>> Q.inv()
[ 1/5, 2/5]
[-2/5, 1/5]
>>> Λ               # the matrix of eigenvalues
[5,  0]
[0, 10]
>>> Q*Λ*Q.inv()     # eigendecomposition of A
[ 9, -2]
[-2,  6]
>>> Q.inv()*A*Q     # obtain Λ from A and Q 
[5,  0]
[0, 10]

Not all matrices are diagonalizable. You can check if a matrix is diagonalizable by calling its is_diagonalizable method:

>>> A.is_diagonalizable()
True
>>> B = Matrix([[ 1, 3],
  	          [ 0, 1]])
>>> B.is_diagonalizable()
False
>>> B.eigenvals()
{1: 2}  
>>> B.eigenvects()
[(1, 2,  [1]
         [0]  )]

The matrix $B$ is not diagonalizable because it doesn't have a full set of eigenvectors. To diagonalize a $2\times 2$ matrix, we need two orthogonal eigenvectors but $B$ has only a single eigenvector. Therefore, we can't construct the matrix of eigenvectors $Q$ (we're missing a column!) and so $B$ is not diagonalizable.

Non-square matrices do not have eigenvectors and therefore don't have an eigendecomposition. Instead, we can use the singular value decomposition to brake up a non-square matrix $A$ into left singular vectors, right singular vectors, and a diagonal matrix of of singular values. Use the singular_values method on any matrix to find its singular values.

It is possible to write a matrix $A$ as the product of an orthogonal matrix $Q$ and an upper triangular matrix $R$. This is known as the QR-decomposition.

>>> A=Matrix([[12,-51,4], 
              [6,167,-68], 
              [-4,24,-41] ])
>>> Q,R = A.QRdecomposition()
>>> Q
[ 6/7, -69/175, -58/175]
[ 3/7, 158/175,   6/175]
[-2/7,    6/35,  -33/35]
>>> Q*Q.T
[1, 0, 0]         
[0, 1, 0]         # Q is orthogonal 
[0, 0, 1]  
>>> R             
[14,  21, -14]    
[ 0, 175, -70]    # R is upper triangular
[ 0,   0,  35]
>>> Q*R
[12, -51,   4]
[ 6, 167, -68]
[-4,  24, -41]

Each sumpy matrix is also equipped with ''LUdecomposition'' and ''cholesky'' decomposition methods.

I would like to conclude with some words of caution about the overuse of computers. Computer technology is very powerful and is everywhere around us, but let's not forget that computers are actually very dumb: computers are mere calculators and they depend on your knowledge to direct them. It's important that you know how to do complicated math by hand in order to be able to instruct computers to do math for you and to check results of your computer calculations. I don't want you to use the tricks you leaned in this tutorial to avoid math problems from now on and simply rely blindly on sympy for all your math needs. I want both you and the computer to become math powerhouses! The computer will help you with tedious calculations (they're good at that) and you'll helping the computer by fixing when it gets stuck (humans are good at that).

The coverage of the math and physics topics in this tutorial were not sufficient to do justice to the subjects. The goal here is to quickly introduce you to the useful sympy commands.

you should consider some of my other tutorials on mechanics If you're interested in learning more about calculus and mechanics, you should consider the No bullshit guide to math and physics—a short textbook like no other. The book starts with an introductory chapter that reviews high school math concepts in order to make math and physics ideas accessible to everyone.

[ The official sympy tutorial ]
http://docs.sympy.org/0.7.2/tutorial.html

[ A list of python gotchas ]
http://docs.sympy.org/dev/gotchas.html