This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.
Suppose \( x^2=2 \). Then most of us are able to find out that \( x=\sqrt{2} \) is a solution to the equation. The more mathematically interested reader will also remark that \( x=-\sqrt{2} \) is another solution. But faced with the equation \( x^2=-2 \), very few are able to find a proper solution without any previous knowledge of complex numbers. Such numbers have many applications in science, and it is therefore important to be able to use such numbers in our programs.
On the following pages we extend the previous material on computing with real numbers to complex numbers. The text is optional, and readers without knowledge of complex numbers can safely drop this section and jump to the section Summary.
A complex number is a pair of real numbers \( a \) and \( b \), most often written as \( a+bi \), or \( a+ib \), where \( i \) is called the imaginary unit and acts as a label for the second term. Mathematically, \( i=\sqrt{-1} \). An important feature of complex numbers is definitely the ability to compute square roots of negative numbers. For example, \( \sqrt{-2}=\sqrt{2}i \) (i.e., \( \sqrt{2}\sqrt{-1} \)). The solutions of \( x^2=-2 \) are thus \( x_1=+\sqrt{2}i \) and \( x_2=-\sqrt{2}i \).
There are rules for addition, subtraction, multiplication, and division between two complex numbers. There are also rules for raising a complex number to a real power, as well as rules for computing \( \sin z \), \( \cos z \), \( \tan z \), \( e^z \), \( \ln z \), \( \sinh z \), \( \cosh z \), \( \tanh z \), etc. for a complex number \( z = a+ib \). We assume in the following that you are familiar with the mathematics of complex numbers, at least to the degree encountered in the program examples. $$ \begin{equation*} \hbox{let } u=a+bi\hbox{ and }v=c+di\end{equation*} $$
The following rules reflect complex arithmetics: $$ \begin{align*} u &= v\quad\Rightarrow\quad a=c, \ b=d\\ -u &= -a -bi\\ u^* &\equiv a - bi\quad\hbox{(complex conjugate)}\\ u+ v &= (a+c) + (b+d)i\\ u - v &= (a-c) + (b-d)i\\ uv &= (ac - bd) + (bc+ad)i\\ u/v &= \frac{ac +bd}{c^2 + d^2} + \frac{bc-ad}{c^2+d^2}i\\ |u| &= \sqrt{a^2 + b^2}\\ e^{iq} &= \cos q + i\sin q \end{align*} $$
Python supports computation with complex numbers. The imaginary unit
is written as j in Python, instead of \( i \) as in mathematics. A
complex number \( 2-3i \) is therefore expressed as (2-3j) in Python.
We remark that the number \( i \) is written as 1j, not just j. Below
is a sample session involving definition of complex numbers and some
simple arithmetics:
>>> u = 2.5 + 3j # create a complex number
>>> v = 2 # this is an int
>>> w = u + v # complex + int
>>> w
(4.5+3j)
>>> a = -2
>>> b = 0.5
>>> s = a + b*1j # create a complex number from two floats
>>> s = complex(a, b) # alternative creation
>>> s
(-2+0.5j)
>>> s*w # complex*complex
(-10.5-3.75j)
>>> s/w # complex/complex
(-0.25641025641025639+0.28205128205128205j)
A complex object s has functionality for extracting
the real and imaginary
parts as well as computing the complex conjugate:
>>> s.real
-2.0
>>> s.imag
0.5
>>> s.conjugate()
(-2-0.5j)
Taking the sine of a complex number does not work:
>>> from math import sin
>>> r = sin(w)
Traceback (most recent call last):
File "<input>", line 1, in ?
TypeError: can't convert complex to float; use abs(z)
The reason is that the sin function from the math module only
works with real (float) arguments, not complex. A similar module,
cmath, defines functions that take a complex number as argument and
return a complex number as result. As an example of using the cmath
module, we can demonstrate that the relation \( \sin(ai)=i\sinh a \)
holds:
>>> from cmath import sin, sinh
>>> r1 = sin(8j)
>>> r1
1490.4788257895502j
>>> r2 = 1j*sinh(8)
>>> r2
1490.4788257895502j
Another relation, \( e^{iq} = \cos q + i\sin q \), is exemplified next:
>>> q = 8 # some arbitrary number
>>> exp(1j*q)
(-0.14550003380861354+0.98935824662338179j)
>>> cos(q) + 1j*sin(q)
(-0.14550003380861354+0.98935824662338179j)
The cmath functions always return complex numbers. It would be nice
to have functions that return a float object if the result is a real
number and a complex object if the result is a complex number. The
Numerical Python package has such versions of the basic mathematical
functions known from math and cmath. By taking a
from numpy.lib.scimath import *
one obtains access to these flexible versions of mathematical functions. The functions also get imported by any of the statements
from scipy import *
from scitools.std import *
A session will illustrate
what we obtain.
Let us first use the sqrt function in the
math module:
>>> from math import sqrt
>>> sqrt(4) # float
2.0
>>> sqrt(-1) # illegal
Traceback (most recent call last):
File "<input>", line 1, in ?
ValueError: math domain error
If we now import sqrt from cmath,
>>> from cmath import sqrt
the previous sqrt function is overwritten by the new one.
More precisely, the name sqrt was previously bound to a function
sqrt from the math module, but is now bound to another
function sqrt from the cmath module. In this case, any
square root results in a complex object:
>>> sqrt(4) # complex
(2+0j)
>>> sqrt(-1) # complex
1j
If we now take
>>> from numpy.lib.scimath import *
we import (among other things) a new sqrt function. This function
is slower than the versions from math and cmath, but it
has more flexibility since the returned object is float if that is
mathematically possible, otherwise a complex is returned:
>>> sqrt(4) # float
2.0
>>> sqrt(-1) # complex
1j
As a further illustration of the need for flexible treatment of both complex and real numbers, we may code the formulas for the roots of a quadratic function \( f(x)=ax^2 + bx + c \):
>>> a = 1; b = 2; c = 100 # polynomial coefficients
>>> from numpy.lib.scimath import sqrt
>>> r1 = (-b + sqrt(b**2 - 4*a*c))/(2*a)
>>> r2 = (-b - sqrt(b**2 - 4*a*c))/(2*a)
>>> r1
(-1+9.94987437107j)
>>> r2
(-1-9.94987437107j)
Using the up arrow, we may go back to the definitions of the coefficients and change them so the roots become real numbers:
>>> a = 1; b = 4; c = 1 # polynomial coefficients
Going back to the computations of r1 and r2 and performing
them again, we get
>>> r1
-0.267949192431
>>> r2
-3.73205080757
That is, the two results are float objects. Had we applied
sqrt from cmath, r1 and r2 would always be
complex objects, while sqrt from the math module
would not handle the first (complex) case.