Copyright Allen Downey 2015, MIT license
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 [56]:
c = complex(1, 2)
c
Out[56]:
Also, the letter j can be appended to any number type to indicate the imaginary unit:
In [57]:
c = 1 + 2j
c
Out[57]:
Python uses j to denote the complext unit, but you might also want to define a variable:
In [145]:
i = complex(0, 1)
i
Out[145]:
In [146]:
def func(x):
return (x+1)**2 == -9
We can confirm that 1 is not a solution:
In [148]:
func(1)
Out[148]:
And that $-1 + 3i$ is:
In [147]:
x = -1 + 3 * i
func(x)
Out[147]:
In [63]:
y = -3.5 + 2j
In Python, these parts are accessible as attributes using dot notation:
In [64]:
y.real
Out[64]:
In [65]:
y.imag
Out[65]:
In [150]:
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 [151]:
abs(z)
Out[151]:
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 [152]:
from cmath import *
phase returns the angle of the complex number in radians, relative to the real axis.
In [69]:
phase(z)
Out[69]:
polar expresses a complex number in polar form, returning a tuple of magnitude and angle:
In [70]:
polar(z)
Out[70]:
In [71]:
y == z
Out[71]:
In [72]:
# raise TypeError: no ordering relation is defined for complex numbers
# y > z
In [73]:
z.conjugate()
Out[73]:
Conjugation distributes over arithmetic operations like addition. To confirm, I'll create a second complex number:
In [74]:
w = complex(12, 13)
And confirm that $\bar{z+w} = \bar{z} + \bar{w}$:
In [75]:
(z + w).conjugate() == z.conjugate() + w.conjugate()
Out[75]:
Python knows how to compute the reciprocal of a complex number:
In [153]:
1/z
Out[153]:
But we can also check to confirm that $1/z = \bar{z} / (x^2 + y^2)$:
In [77]:
1 / z == z.conjugate() / (z.real**2 + z.imag**2)
Out[77]:
In [154]:
i**2
Out[154]:
Now we can confirm the behavior of complex multiplication:
In [155]:
a, b, c, d = 1, 2, 3, 4
(a + b*i) * (c + d*i) == (a*c - b*d) + (b*c + a*d)*i
Out[155]:
In [184]:
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[184]:
If you check whether your answer is exactly right, it might fail due to floating point approximation:
In [185]:
root_z**2 == z
Out[185]:
But you can check it by hand:
In [186]:
root_z**2
Out[186]:
Or check that the difference is small:
In [187]:
abs(root_z**2 - z)
Out[187]:
In [167]:
phase(z)
Out[167]:
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 [168]:
from math import atan2
atan2(z.imag, z.real)
Out[168]:
To convert from polar to Cartesian coordinates, you can use trigonometry:
In [169]:
r, phi = polar(z)
r * complex(cos(phi), sin(phi))
Out[169]:
Or use rect, which converts to rectangular coordinates.
In [171]:
rect(r, phi)
Out[171]:
In [ ]:
In [188]:
e
Out[188]:
Which we can use to compute complex exponentials (see the Exponentiation section):
In [191]:
x = 2 # angle in radians
e**(i*x)
Out[191]:
Using pi, we can test Euler's identity, $e^{i \pi} = -1$, one of the more exquisite results in mathematics:
In [194]:
e**(i*pi)
Out[194]:
Or more generally, we can check Euler's formula:
In [193]:
e**(i*x) == complex(cos(x), sin(x))
Out[193]:
In [90]:
log_z = log(r) + phi*i
log_z
Out[90]:
Or using the log function we imported from cmath:
In [174]:
log(z)
Out[174]:
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 [ ]:
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:
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 [108]:
R = 1000 # ohm
L = 0.001 # henry
C = 0.00001 # farad
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 [109]:
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 [110]:
Vin = 5 # Volt
We can use Ohm's law to compute the current:
In [111]:
I = Vin / Z
I
Out[111]:
And then use Ohm's law again to compute the voltage across the resistor, which is the output signal.
In [112]:
Vout = I * R
Vout
Out[112]:
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 [118]:
abs(Vout)
Out[118]:
In this case the amplitude is very close to 5V, so this filter "passes" signals at 1000 Hz with little attenuation.
In [119]:
phase(Vout)
Out[119]:
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 [176]:
import numpy as np
fs = np.logspace(0, 6, 20, base=10)
fs
Out[176]:
Now we can use the array $fs$ anyplace we previously used the scalar value $f$. The result is an array of complex impedances.
In [132]:
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 [133]:
v_outs = Vin / zs * R
v_outs
Out[133]:
Again, the magnitudes encode the amplitude of the output signal (in volts) for each frequency:
In [139]:
magnitudes = abs(v_outs)
magnitudes
Out[139]:
And the angles contain the phase offsets (in radians):
In [140]:
offsets = np.angle(v_outs)
offsets
Out[140]:
I'll use matplotlib to plot the results.
In [180]:
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[180]:
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 [181]:
plt.plot(fs, offsets)
plt.xscale('log')
plt.xlabel('frequency (Hz)')
plt.ylabel('phase offset (radian)')
plt.ylim(-pi/2, pi/2)
Out[181]:
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 [ ]: