Bode Plots using Python

 I needed a quick way to plot some Bode plots for a second-order system. I didn't have access to Matlab, instead, I searched for a solution using Python, and I found one. Documentation is a bit sparse so this example might be helpful.

The signals package supports the signal.bode method which turned out to be quite easy to use. Signal is part of the scipy package and is something we bundle with our Tellurium platform.

There is a more comprehensive discussion of Python and Control Theory here.

[code lang="python" fontsize="10"]
#
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

# Coefficients in numerator of transfer function
num = [1]
# Coefficients in denominator of transfer function
# High order to low order, eg 1*s^2 + 0.1*s + 1
den = [1, 0.1, 1]

# Scan over zeta, a parameter for a second-order system
zetaRange = np.arange(0.1,1.1,0.1)

f1 = plt.figure()
for i in range(0,9):
    den = [1, 2*zetaRange[i], 1]
    print den
    s1 = signal.lti(num, den)
    # Specify our own frequency range: np.arange(0.1, 5, 0.01)
    w, mag, phase = signal.bode(s1, np.arange(0.1, 5, 0.01).tolist())
    plt.semilogx (w, mag, color="blue", linewidth="1")
plt.xlabel ("Frequency")
plt.ylabel ("Magnitude")
plt.savefig ("c:\\mag.png", dpi=300, format="png")

plt.figure()

for i in range(0,9):
    den = [1, zetaRange[i], 1]
    s1 = signal.lti(num, den)
    w, mag, phase = signal.bode(s1, np.arange(0.1, 10, 0.02).tolist())
    plt.semilogx (w, phase, color="red", linewidth="1.1")
plt.xlabel ("Frequency")
plt.ylabel ("Amplitude")
plt.savefig ("c:\\phase.png", dpi=300, format="png")

[/code]

Output from Python script:

Magnitude plot:

Phase plot:

Simple RK4 Code

January 4, 2015 10:11 pm 

Thought I'd start off the New Year with a simple freebie, a 4th Order Runge-Kutta class for Delphi, Windows and Mac OS. Should also work on Android and iOS mobile platforms and with Free Pascal. The code below shows how you'd use it. First declare the set of differential equations that need to be solved (TDoubleDynArray type can be found in Types):

 TMySystem = class (TObject)
    class procedure func (time : double; y, p, dydt : TDoubleDynArray);
 end;

class procedure TMySystem.func (time : double; y, p, dydt : TDoubleDynArray);
begin
  dydt[0] := p[0] - p[1]*y[0];
  dydt[1] := p[1]*y[0] - p[2]*y[1];
end;

Next, create a Runge-Kutta object:

// First argument = number of ODEs
// Second argument = number of parameters, if any
// Third argument = ODE Function
rk4 := TRK4.Create (2, 3, TMySystem.func);
rk4.p[0] := vo; rk4.p[1] := k1; rk4.p[2] := k2;

The above code include the option to declare and initialize some parameters that are part of the differential equations. To actually integrate the system by one step use the line:

// Return new time point, assign to startTime ready for next time
startTime := rk4.eval (startTime, y, stepSize);

The y argument will contain the updated solution. Run the line again to do the next step etc. Note that the eval method takes three arguments, that current value for the independent variable (time), an array that contains the current values for the dependent variable, and the step size to use.

Download the code here. RK4.zip The download include a Windows exe that you can try right away.