• API
• FAQ
• Tools
• Archive
daily pastebin goal
28%
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:
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:
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.

Top