daily pastebin goal
70%
SHARE
TWEET

Untitled

a guest Dec 9th, 2018 58 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. from math import sqrt
  2. from collections import OrderedDict
  3.  
  4. import numpy as np
  5. import matplotlib.pyplot as plt
  6. from prettytable import PrettyTable
  7.  
  8. if __name__ == '__main__':
  9.     solutions_table = PrettyTable('x1 x2 x3 x4 p1 p2 p3 p4'.split())
  10.     solutions_table.align = 'l'
  11.  
  12.     verification_table = PrettyTable('x1 dx1/dt dx2/dt dx3/dt dx4/dt Jacobian'.split())
  13.     verification_table.align = 'l'
  14.  
  15.     plt.title('Bifurcation Diagram')
  16.     plt.xlabel('p3')
  17.     plt.ylabel('p2')
  18.     plt.grid(True)
  19.  
  20.     points = OrderedDict()
  21.     p1 = 2
  22.  
  23.     for x1 in (val * 0.1 for val in range(1, 40)):
  24.         x3 = 2 * p1 - x1
  25.  
  26.         if x1 == p1:
  27.             solutions_table.add_row('-' * 8)
  28.             continue
  29.  
  30.         alpha1 = 20 * (x1 ** 2 + x3 ** 2)
  31.         beta1 = 2 * x1 ** 2 * x3 ** 2 + 10 * (x1 ** 2 + x3 ** 2)
  32.         gamma1 = -20 * x1 * x3
  33.         delta1 = x1 ** 2 * x3 ** 2
  34.  
  35.         alpha2 = alpha1 - (80 * p1 * x3 * (x1 - p1)) / x1
  36.         beta2 = beta1 - 4 * x3 * (x1 - p1) * (x1 ** 2 + (10 * p1) / x1)
  37.         gamma2 = 10 * (x1 ** 2 + x3 ** 2 - 4 * p1 * x3)
  38.         delta2 = delta1 - 2 * x1 ** 2 * x3 * (x1 - p1)
  39.  
  40.         alpha = alpha1 * gamma2 - alpha2 * gamma1
  41.         beta = beta1 * gamma2 - beta2 * gamma1
  42.         delta = delta1 * gamma2 - delta2 * gamma1
  43.  
  44.         p3 = (-beta + sqrt(beta ** 2 - 4 * alpha * delta)) / (2 * alpha)
  45.  
  46.         if not p3:
  47.             solutions_table.add_row('-' * 8)
  48.             continue
  49.  
  50.         p4 = 10 * p3
  51.         p2 = -(alpha2 * p3 + beta2 + delta2 / p3) / gamma2
  52.         x2 = (p2 * x1 + (2 * p3 + 1) * (x1 - p1)) / (x1 ** 2)
  53.         x4 = x2 + (2 * p3 + 1) * (x1 - p1) / p4
  54.  
  55.         points[p3] = p2
  56.         solutions_table.add_row((x1, x2, x3, x4, p1, p2, p3, p4))
  57.  
  58.         jacobian_matrix = np.array((2 * x1 * x2 - p2 - p3 - 1, x1 ** 2, p3, 0,
  59.                                     -2 * x1 * x2 + p2, -x1 ** 2 - p4, 0, p4, p3, 0,
  60.                                     2 * x3 * x4 - p2 - p3 - 1, x3 ** 2,
  61.                                     0, p4, -2 * x3 * x4 + p2, -x3 ** 2 - p4))
  62.         jacobian_matrix = jacobian_matrix.reshape(4, 4)
  63.         fmt = '{:.0e}'
  64.         dx1 = fmt.format(p1 - (p2 + 1) * x1 + x1 ** 2 * x2 + p3 * (x3 - x1))
  65.         dx2 = fmt.format(p2 * x1 - x1 ** 2 * x2 + p4 * (x4 - x2))
  66.         dx3 = fmt.format(p1 - (p2 + 1) * x3 + x3 ** 2 * x4 - p3 * (x3 - x1))
  67.         dx4 = fmt.format(p2 * x3 - x3 ** 2 * x4 - p4 * (x4 - x2))
  68.         det = fmt.format(np.linalg.det(jacobian_matrix))
  69.         verification_table.add_row((x1, dx1, dx2, dx3, dx4, det))
  70.  
  71.     print 'Solution:\n', solutions_table
  72.     print '\n Verification:\n', verification_table
  73.  
  74.     plt.plot(points.keys(), points.values())[0].axes.set_xscale('log')
  75.     points.clear()
  76.  
  77.     for p3 in (val * 0.1 for val in range(-10, 5) if val != 0):
  78.         points[p3] = (2 * p3 + 1) * (p1 ** 2 + 20 * p3) / (20 * p3)
  79.  
  80.     plt.plot(points.keys(), points.values())
  81.     plt.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top