gearSoluce.py
# -------------------------------------------------------------------------
#
# PYTHON for DUMMIES
# Probleme 7
#
# Solution detaillée
# Vincent Legat
#
# -------------------------------------------------------------------------
from numpy import *
def stabilityGear(x,y,order):
#
# -1- Coefficients des méthodes de Gear
#
# Source pour les coefficients :
# http://en.wikipedia.org/wiki/Backward_differentiation_formula
#
# C'est bien juste : les calculs savants de quelques tuteurs
# donnent les memes valeurs :-)
# Moralité : ne recalcule pas ce qui a deja été fait :-)
#
if (order > 6) :
return zeros(order+1),zeros(shape(x))
coeffs = {}
coeffs[1] = array([ -1, -1, ])/1
coeffs[2] = array([ -2, -4, 1, ])/3
coeffs[3] = array([ -6, -18, 9, -2, ])/11
coeffs[4] = array([-12, -48, 36, -16, 3, ])/25
coeffs[5] = array([-60, -300, 300, -200, 75, -12, ])/137
coeffs[6] = array([-60, -360, 450, -400, 225, -72, 10 ])/147
coeff = -coeffs[order]
#
# -2- Calcul des racines du polynôme de stabilité
# Le gain est le maximum des modules des racines obtenues
#
z = x + 1j*y
gain = zeros(shape(z))
for i in range(shape(z)[0]) :
for j in range(shape(z)[1]) :
gain[i,j] = max(abs(roots([(coeff[0]*z[i,j] -1),*coeff[1:]])))
return gain,coeff