cubiс_equation.py

gistfile1.txt

# encoding: utf8
import pylab, math

# Задание переменных
delta_omega_1 = 0.1
delta_omega_2 = 0.2
omega_2_omega_1 = 10.0
gamma = 2.0

# левая и правая границы графика
x_range = [-15, 15]

# Что варьировать?
# var = ['delta_omega_1', 'delta_omega_2', 'omega_2_omega_1', 'gamma']
var = ['gamma']

# Задание пределов варьирования
delta_omega_1_range = [0.01, 0.1]
delta_omega_2_range = [0.01, 0.1]
omega_2_omega_1_range = [0.1, 10]
gamma_range = [0.1, 2]

# -------------- Реализация --------------
def make_range(array):
  n = 20
  diff = array[-1] - array[0]
  step = diff / float(n)
  res = [array[0] + step*i for i in range(n+1)]
  return res

def calculate_coef(_delta_omega_1, _delta_omega_2, _omega_2_omega_1, _gamma):
  c = 1.0 + _gamma * _omega_2_omega_1**2
  b1 = 1.0/_gamma**(1.0/3.0) * (_delta_omega_1)**(4.0/3.0)
  b2 = _gamma**(2.0/3.0) / (_delta_omega_1)**(2.0/3.0)
  b3 = (_omega_2_omega_1)**(4.0/3.0) / (_delta_omega_2)**(2.0/3)
  b = b1 + b2 + b3
  a = (1 + _gamma)*(_delta_omega_1 / _gamma)**(2.0/3.0)
  return [a, b, c]

# x^3 - ax^2 - bx + c = 0
def solve(_a, _b, c):
  a = -_a
  b = -_b

  print 'x^3 - ax^2 - bx + c = 0'
  print 'a = ', a
  print 'b = ', b
  print 'c = ', c

  q = (a**2 - 3*b)/9.0
  r = (2*a**3 - 9*a*b + 27*c)/54.0

  s = q**3 - r**2

  if s > 0:
    sqrt_q = q**0.5
    fi = 1.0/3.0*math.acos(r/q**(3.0/2.0))

    x1 = -2 * sqrt_q * math.cos(fi) - a/3.0
    x2 = -2 * sqrt_q * math.cos(fi + 2.0/3.0 * math.pi) - a/3.0
    x3 = -2 * sqrt_q * math.cos(fi - 2.0/3.0 * math.pi) - a/3.0
    return [x1, x2, x3]
  elif s < 0:
    fi = 1.0/3.0 * math.acosh( abs(r)/q**(3.0/2.0) )
    sign_r = math.copysign(1, r)
    x1 = -2 * sign_r * q**0.5 * math.cosh(fi) - a/3.0
    return [x1]
  else:
    return []

def fx(x, a, b, c):
  return x**3 - a*x**2 - b*x + c

x = pylab.arange(x_range[0], x_range[1], 0.01)
pylab.plot(x, fx(x, *calculate_coef(delta_omega_1, delta_omega_2, omega_2_omega_1, gamma)), 'k--', linewidth = 5, color='r')

print 'Решение:'
solution = solve(*calculate_coef(delta_omega_1, delta_omega_2, omega_2_omega_1, gamma))
print solution

# Построение графиков
delta_omega_1_range = make_range(delta_omega_1_range)
delta_omega_2_range = make_range(delta_omega_2_range)
omega_2_omega_1_range = make_range(omega_2_omega_1_range)
gamma_range = make_range(gamma_range)

if 'delta_omega_1' in var:
  for i in range(len(delta_omega_1_range) - 1):
    pylab.plot(x, fx(x, *calculate_coef(delta_omega_1_range[i], delta_omega_2, omega_2_omega_1, gamma)), color='g')
  pylab.plot(x, fx(x, *calculate_coef(delta_omega_1_range[-1], delta_omega_2, omega_2_omega_1, gamma)), color='g', label='do1')
  pylab.legend()

if 'delta_omega_2' in var:
  for i in range(len(delta_omega_2_range) - 1):
    pylab.plot(x, fx(x, *calculate_coef(delta_omega_1, delta_omega_2_range[i], omega_2_omega_1, gamma)), color='b')
  pylab.plot(x, fx(x, *calculate_coef(delta_omega_1, delta_omega_2_range[-1], omega_2_omega_1, gamma)), color='b', label='do2')
  pylab.legend()

if 'omega_2_omega_1' in var:
  for i in range(len(omega_2_omega_1_range) - 1):
    pylab.plot(x, fx(x, *calculate_coef(delta_omega_1, delta_omega_2, omega_2_omega_1_range[i], gamma)), color='pink')
  pylab.plot(x, fx(x, *calculate_coef(delta_omega_1, delta_omega_2, omega_2_omega_1_range[-1], gamma)), color='pink', label='o2_o1')
  pylab.legend()

if 'gamma' in var:
  for i in range(len(gamma_range) - 1):
    pylab.plot(x, fx(x, *calculate_coef(delta_omega_1, delta_omega_2, omega_2_omega_1, gamma_range[i])), color='brown')
  pylab.plot(x, fx(x, *calculate_coef(delta_omega_1, delta_omega_2, omega_2_omega_1, gamma_range[-1])), color='brown', label='g')
  pylab.legend()

pylab.grid(True)
pylab.show()