algebraic expressions in python

Question

I'm working with some relatively simple algebraic expressions in python and wondered if there's a better way of plotting algebraic expressions than what I'm currently doing:

I have the equation of a notch filter (below image - left side from a book; the right side graph is the fig generated by my code), and so far my code works though is crude.

Is there a better way of plotting the magnitude? As shown the constants are; R = 50.0, C = 470e-12, L = 54e-6, and the desired frequency range is from 0 to 2MHz.

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

R = 50.0
C = 470e-12
L = 54e-6

FREQ = []
DATA = []
for i in range(1, 200):
    f = i*10000.0
    w = 2*np.pi*f

    Ztop = w*L - 1.0/(w*C)
    Zbot = math.sqrt( (math.pow(R,2) + math.pow((w*L) -(1.0/(w*C)),2)) )
    Zout = abs(Ztop / Zbot)
    FREQ.append( f/1e6 )
    DATA.append( Zout )

plt.figure(1)
plt.plot(FREQ,DATA, '-k')
plt.xlabel('Frequency (MHz)')
plt.ylabel('Mag.')
plt.grid()
plt.show()

Answer 1

The main improvement I'd suggest is using Numpy's array broadcasting facility. Basically this will perform a mathematical operation on a whole array at once. Other than that, there's nothing really wrong with your code, though you could calculate the entire magnitude in one expression. That's your call, depending on whether you think it makes the code clearer.

# creates 201 values evenly spaced from 0 to 2e6, inclusive
f = np.linspace(0, 2e6, 201)
w = 2 * np.pi * f
# now calculate the magnitude for all 201 values at once
magnitude = (w * L - 1 / (w * C)) / np.sqrt(R**2 + (w * L - 1 / (w * C))**2)
# now you can plot magnitude vs. frequency
plt.plot(f / 1e6, magnitude, '-k')

Answer 2

You could also use sympy to plot symbolic expressions:

from __future__ import division
from sympy import pi, pprint, sqrt, symbols # $ pip install sympy
from sympy.plotting import plot

f = symbols("f")
L, C, R = symbols("L, C, R", real=True)
w = 2*pi*f
W = w * L - 1 / (w * C)
mag = abs(W / sqrt(R**2 + W**2))
pprint(mag)
plot(mag.subs({R: 50.0, C: 470e-12, L: 54e-6}), (f, 1, 2e6))

Output

│                    1           │
│       2⋅π⋅L⋅f - ───────        │
│                 2⋅π⋅C⋅f        │
│────────────────────────────────│
│     ___________________________│
│    ╱                         2 │
│   ╱   2   ⎛             1   ⎞  │
│  ╱   R  + ⎜2⋅π⋅L⋅f - ───────⎟  │
│╲╱         ⎝          2⋅π⋅C⋅f⎠  │