Skip to content

Instantly share code, notes, and snippets.

@Houdini
Created Feb 26, 2012
Embed
What would you like to do?
cubiс_equation.py
# 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment