## Exploring complex numbers with Python

This notebook is an example of a process for investigating math topics using Python. In this case, I follow the presentation in the Wikipedia page for complex numbers (http://en.wikipedia.org/wiki/Complex_number), translating the examples into Python.

There are a few ways to write complex numbers in Python. One option is the built-in function complex:



In :

c = complex(1, 2)
c




Out:

(1+2j)



Also, the letter j can be appended to any number type to indicate the imaginary unit:



In :

c = 1 + 2j
c




Out:

(1+2j)



Python uses j to denote the complext unit, but you might also want to define a variable:



In :

i = complex(0, 1)
i




Out:

1j



### Overview

The Wikipedia overview gives an example of an equation with no real solution: $(x+1)^2 = -9$. We can write a function to test whether x is a solution to this equation



In :

def func(x):
return (x+1)**2 == -9



We can confirm that 1 is not a solution:



In :

func(1)




Out:

False



And that $-1 + 3i$ is:



In :

x = -1 + 3 * i
func(x)




Out:

True



### Definition

The next section of the Wikipedia article defines the real and imaginary parts of a complex number:



In :

y = -3.5 + 2j



In Python, these parts are accessible as attributes using dot notation:



In :

y.real




Out:

-3.5




In :

y.imag




Out:

2.0



### Complex plane

In the complex plane, complex numbers can be described in polar coordinates.



In :

z = 3 + 4j



Python's built-in function, abs, computes the absolute value, or magnitude, of a complex number, which is the distance from the origin of the complext plane:



In :

abs(z)




Out:

5.0



Sadly, there is no built-in function to compute the angle (also called argument or phase). But the cmath module provides phase and other useful functions for working with complex numbers.



In :

from cmath import *



phase returns the angle of the complex number in radians, relative to the real axis.



In :

phase(z)




Out:

0.9272952180016122



polar expresses a complex number in polar form, returning a tuple of magnitude and angle:



In :

polar(z)




Out:

(5.0, 0.9272952180016122)



### Equality

The Relations section of the Wikipedia page defines equality for complex numbers, which is implemented by the == operator:



In :

y == z




Out:

False



### Ordering

It also indicates that there is no linear ordering on the complex numbers, so the comparator operators raise an exception:



In :

# raise TypeError: no ordering relation is defined for complex numbers
# y > z



### Conjugation

The Elementary Operators section defines the complex conjugate, which is implemented as a method you can invoke on a complex number:



In :

z.conjugate()




Out:

(3-4j)



Conjugation distributes over arithmetic operations like addition. To confirm, I'll create a second complex number:



In :

w = complex(12, 13)



And confirm that $\bar{z+w} = \bar{z} + \bar{w}$:



In :

(z + w).conjugate() == z.conjugate() + w.conjugate()




Out:

True



Python knows how to compute the reciprocal of a complex number:



In :

1/z




Out:

(0.12-0.16j)



But we can also check to confirm that $1/z = \bar{z} / (x^2 + y^2)$:



In :

1 / z == z.conjugate() / (z.real**2 + z.imag**2)




Out:

True



### Multiplication and division

The section on multiplication and division indicates (for the first time, I think) that $i$ is a square root of -1:



In :

i**2




Out:

(-1+0j)



Now we can confirm the behavior of complex multiplication:



In :

a, b, c, d = 1, 2, 3, 4
(a + b*i) * (c + d*i) == (a*c - b*d) + (b*c + a*d)*i




Out:

True



### Square root

Excerise: Use the formulas in the Square root section to compute the square root of $z = 3 - 4i$:



In :

from numpy import sign

# EXERCISE

z = 3 - 4*i
a, b = z.real, z.imag
real = sqrt(a + abs(z))
imag = sign(b) * sqrt(-a + abs(z))
root_z = complex(real, imag) / sqrt(2)
root_z




Out:

(2-1j)



If you check whether your answer is exactly right, it might fail due to floating point approximation:



In :

root_z**2 == z




Out:

True



But you can check it by hand:



In :

root_z**2




Out:

(3-4j)



Or check that the difference is small:



In :

abs(root_z**2 - z)




Out:

0.0



### Polar form

The next section takes us back to polar form. Again, phase computes the angle of a complex number:



In :

phase(z)




Out:

-0.9272952180016122



phase can also be defined in terms of arctan; the best way to compute it is with math.atan2, which correctly figures out which quadrant the number is in:



In :

from math import atan2

atan2(z.imag, z.real)




Out:

-0.9272952180016122



To convert from polar to Cartesian coordinates, you can use trigonometry:



In :

r, phi = polar(z)
r * complex(cos(phi), sin(phi))




Out:

(3.0000000000000004-3.9999999999999996j)



Or use rect, which converts to rectangular coordinates.



In :

rect(r, phi)




Out:

(3.0000000000000004-3.9999999999999996j)



### Multiplication and division

Exercise: use Python to confirm that $z_1 z_2 = r_1 r_2 (\cos(\phi_1 + \phi_2) + i \sin(\phi_1 + \phi_2))$. In other words, when you multiply two complex numbers, you multiply the magnitudes and add the angles.



In [ ]:



### Euler's formula

cmath provides e and pi



In :

e




Out:

2.718281828459045



Which we can use to compute complex exponentials (see the Exponentiation section):



In :

x = 2      # angle in radians
e**(i*x)




Out:

(-0.4161468365471424+0.9092974268256817j)



Using pi, we can test Euler's identity, $e^{i \pi} = -1$, one of the more exquisite results in mathematics:



In :

e**(i*pi)




Out:

(-1+1.2246063538223773e-16j)



Or more generally, we can check Euler's formula:



In :

e**(i*x) == complex(cos(x), sin(x))




Out:

True



### Natural logarithm

We can compute natural logarithms by hand:



In :

log_z = log(r) + phi*i
log_z




Out:

(1.6094379124341003+0.9272952180016122j)



Or using the log function we imported from cmath:



In :

log(z)




Out:

(1.6094379124341003-0.9272952180016122j)



### Integer and fractional exponents

Exercise: test $z^n = r^n(\cos n \phi + i \sin n \phi)$



In [ ]:



Exercise: test $\sqrt[n]{z} = \sqrt[n]{r}(\cos \frac{\phi}{n} + i \sin \frac{\phi}{n})$. Although in this case note that there are multiple possible results, so you and Python might not generate the same one. Hint: for your example, choose a value of $z$ with a small angle.



In [ ]:



## Application

Imaginary numbers turn out to be more useful than we have any reason to expect :) As an example, I'll demonstrate a method for analyzing circuits. Specifically, this circuit: which (as we'll see) is a bandpass filter.

L is the inductance of an inductor in henry, C is the capacitance of a capacitor in farad, and R is the resistance of a resistor in ohms.



In :

R = 1000    # ohm
L = 0.001   # henry



You might already be familiar with Ohm's law, which you can use to compute the current through a resistor as a function of voltage: $I = V/R$.

Ohm's law can be generalized to handle resistors, capacitors, and inductors, by replacing R (which is a real value) with impedance (which is a complex number). For AC signals, impedance depends on the frequency of the signal, $f$, in Hz.



In :

f = 1000    # Hz
Z = R + 2*pi*i*f*L + 1 / (2*pi*i*f*C)



Assuming that the input is a 5V signal with frequency 1000 Hz:



In :

Vin = 5    # Volt



We can use Ohm's law to compute the current:



In :

I = Vin / Z
I




Out:

(0.004999536136154376+4.8157076930153816e-05j)



And then use Ohm's law again to compute the voltage across the resistor, which is the output signal.



In :

Vout = I * R
Vout




Out:

(4.999536136154376+0.048157076930153815j)



The result is a complex number, so it might not be obvious how to interpret it. It turns out to be a phasor, which is a complex number that represents an AC signal: the magnitude of the phasor is the amplitude of the signal in volts; the angle of the phasor is the phase offset between the input and output signals:



In :

abs(Vout)




Out:

4.999768062697697



In this case the amplitude is very close to 5V, so this filter "passes" signals at 1000 Hz with little attenuation.



In :

phase(Vout)




Out:

0.009632011118962535



The phase shift is quite small, so there is only a short delay between the peak of the input signal and the peak of the output signal.

We might be interested to see how this relationship changes for different frequencies. I'll use numpy to make an array of values from $f = 10^0 = 1$ to $f = 10^6$ Hz.



In :

import numpy as np
fs = np.logspace(0, 6, 20, base=10)
fs




Out:

array([  1.00000000e+00,   2.06913808e+00,   4.28133240e+00,
8.85866790e+00,   1.83298071e+01,   3.79269019e+01,
7.84759970e+01,   1.62377674e+02,   3.35981829e+02,
6.95192796e+02,   1.43844989e+03,   2.97635144e+03,
6.15848211e+03,   1.27427499e+04,   2.63665090e+04,
5.45559478e+04,   1.12883789e+05,   2.33572147e+05,
4.83293024e+05,   1.00000000e+06])



Now we can use the array $fs$ anyplace we previously used the scalar value $f$. The result is an array of complex impedances.



In :

zs = R + 2*pi*i*fs*L + 1 / (2*pi*i*fs*C)



And we can use Ohm's law again to compute an array of phasors.



In :

v_outs = Vin / zs * R
v_outs




Out:

array([ 0.01966160+0.31292401j,  0.08310574+0.63923559j,
0.33740507+1.25426599j,  1.18271078+2.1247939j ,
2.85108448+2.47522518j,  4.25208232+1.78331364j,
4.80339240+0.97179393j,  4.95339824+0.48045505j,
4.98977902+0.22583275j,  4.99828460+0.0925963j ,
4.99997947+0.01013141j,  4.99910856-0.06675635j,
4.99348863-0.18031769j,  4.96913188-0.39164747j,
4.86738576-0.80342061j,  4.47506691-1.53268089j,
3.32698186-2.35925857j,  1.58549840-2.32673308j,
0.48919568-1.48548511j,  0.12352323-0.77611735j])



Again, the magnitudes encode the amplitude of the output signal (in volts) for each frequency:



In :

magnitudes = abs(v_outs)
magnitudes




Out:

array([ 0.31354109,  0.64461516,  1.2988554 ,  2.43177999,  3.77563536,
4.61090139,  4.90071036,  4.97664457,  4.9948869 ,  4.99914222,
4.99998974,  4.99955426,  4.99674325,  4.98454205,  4.93324729,
4.73025734,  4.07859158,  2.81558022,  1.56396241,  0.78588557])



And the angles contain the phase offsets (in radians):



In :

offsets = np.angle(v_outs)
offsets




Out:

array([ 1.50804694,  1.44151346,  1.30801119,  1.06288192,  0.7149484 ,
0.39711591,  0.19961965,  0.09669257,  0.0452282 ,  0.0185235 ,
0.00202629, -0.01335286, -0.03609488, -0.07865348, -0.16358701,
-0.32997184, -0.61682642, -0.97264462, -1.2526647 , -1.41296474])



I'll use matplotlib to plot the results.



In :

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(fs, magnitudes)
plt.xscale('log')
plt.xlabel('frequency (Hz)')
plt.ylabel('amplitude (V)')
plt.ylim(0, 5.2)




Out:

(0, 5.2)



Between $10^2 = 100$ Hz and $10^4 = 10$ kHz, the output amplitude is close to 5V. So this circuit "passes" signals in this frequency band (which is roughly the range of audible frequencies). Signals with lower and higher frequencies get attenuated, or "stopped". That's why this circuit is called a "bandpass filter".

The phase offset also depends on the input frequency.



In :

plt.plot(fs, offsets)
plt.xscale('log')
plt.xlabel('frequency (Hz)')
plt.ylim(-pi/2, pi/2)




Out:

(-1.5707963267948966, 1.5707963267948966)



In the pass band, the phase offset is small, so there is not much delay between the peak of the input and the peak of the output. For high frequencies, the offset approaches $\pi/2$, so the output peaks 90 degrees late. For very low frequencies, the output peaks 90 degrees early (or 270 degrees late, if you prefer).

So that's a pretty neat bit of analysis, considering that all we used was complex multiplication and division.

Exercise: try out some different values of R, L, and C, and see what effect they have on the output signal. Can you make a filter with a narrower pass band? Can you shift the pass band to the left or right?



In [ ]: