Lorenz Equations with Python

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117

import matplotlib.pyplot as plt import numpy as np   vSigma = 10 vBeta = 8/3 vRho = 28     def f1(t,x,y,z):     f1 = vSigma*y - vSigma*x     return f1   def f2(t,x,y,z):     f2 = vRho*x - y - x*z     return f2   def f3(t,x,y,z):     f3 = -vBeta*z + x*y     return f3     print("Sample input: lorenzEquation(0,20,200,5,5,5)")   def lorenzEquation(first,second,N,alpha1,alpha2, alpha3):     h = (second-first)/N     t = first       w1 = alpha1     w2 = alpha2     w3 = alpha3       A = [w1]     B = [w2]     C = [w3]     tempTime = alpha1     time = [alpha1]           A1 = [w1]     B1 = [w2]     C1 = [w3]         for num in range(1,N):         k1x = h*f1(t,w1,w2,w3)         k1y = h*f2(t,w1,w2,w3)         k1z = h*f3(t,w1,w2,w3)                   k2x = h*f1(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))         k2y = h*f2(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))         k2z = h*f3(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))                   k3x = h*f1(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))         k3y = h*f2(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))         k3z = h*f3(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))                             k4x = h*f1(t + h, w1 + k3x, w2 + k3y, w3 + k3z)         k4y = h*f2(t + h, w1 + k3x, w2 + k3y, w3 + k3z)         k4z = h*f3(t + h, w1 + k3x, w2 + k3y, w3 + k3z)                   w1 = w1 + (1/6)*(k1x + 2*k2x + 2*k3x + k4x)         w2 = w2 + (1/6)*(k1y + 2*k2y + 2*k3y + k4y)         w3 = w3 + (1/6)*(k1z + 2*k2z + 2*k3z + k4z)           A.append(w1)         B.append(w2)         C.append(w3)         tempTime = alpha1 + num*h         time.append(tempTime)       w1 = alpha1 + 0.001     w2 = alpha2     w3 = alpha3         for num in range(1,N):         k1x = h*f1(t,w1,w2,w3)         k1y = h*f2(t,w1,w2,w3)         k1z = h*f3(t,w1,w2,w3)                   k2x = h*f1(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))         k2y = h*f2(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))         k2z = h*f3(t + (h/2), w1 + (k1x/2), w2 + (k1y/2), w3 + (k1z/2))                   k3x = h*f1(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))         k3y = h*f2(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))         k3z = h*f3(t + (h/2), w1 + (k2x/2), w2 + (k2y/2), w3 + (k2z/2))                             k4x = h*f1(t + h, w1 + k3x, w2 + k3y, w3 + k3z)         k4y = h*f2(t + h, w1 + k3x, w2 + k3y, w3 + k3z)         k4z = h*f3(t + h, w1 + k3x, w2 + k3y, w3 + k3z)                   w1 = w1 + (1/6)*(k1x + 2*k2x + 2*k3x + k4x)         w2 = w2 + (1/6)*(k1y + 2*k2y + 2*k3y + k4y)         w3 = w3 + (1/6)*(k1z + 2*k2z + 2*k3z + k4z)           A1.append(w1)         B1.append(w2)         C1.append(w3)     #     plt.plot(time, A)     plt.plot(time, A1)     plt.legend(['Initial condition [5,5,5]', 'Initial condition [5.001,5,5]'], loc='upper left')     plt.show()       plt.plot(A, B)     plt.show()           plt.plot(A, C)     plt.show()       return     <h2> Lotka Voltera Equations with Python</h2>

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

import matplotlib.pyplot as plt import numpy as np   a = 1.2 b = 0.6 c = 0.8 d = 0.3   def f1(t,x,y):     f1 = a*x - b*x*y     return f1   def f2(t,x,y):     f2 = -c*y + d*x*y     return f2     print("Sample input:  predatorPrey(0, 30, 300, 2, 1)")     def predatorPrey(first,second,N,alpha1,alpha2):     h = (second-first)/N     t = first           w1 = alpha1     w2 = alpha2           A = [w1]     B = [w2]     tempTime = alpha1     time = [alpha1]           for num in range(1,N):         k1x = h*f1(t,w1,w2)         k1y = h*f2(t,w1,w2)                   k2x = h*f1(t + (h/2), w1 + (k1x/2), w2 + (k1y/2))         k2y = h*f2(t + (h/2), w1 + (k1x/2), w2 + (k1y/2))                   k3x = h*f1(t + (h/2), w1 + (k2x/2), w2 + (k2y/2))         k3y = h*f2(t + (h/2), w1 + (k2x/2), w2 + (k2y/2))                             k4x = h*f1(t + h, w1 + k3x, w2 + k3y)         k4y = h*f2(t + h, w1 + k3x, w2 + k3y)                   w1 = w1 + (1/6)*(k1x + 2*k2x + 2*k3x + k4x)         w2 = w2 + (1/6)*(k1y + 2*k2y + 2*k3y + k4y)                   A.append(w1)         B.append(w2)         tempTime = alpha1 + num*h         time.append(tempTime)                     plt.plot(time, A)     plt.plot(time, B)     plt.legend(['x, prey', 'y, predator'], loc='upper left')     ax.grid()     ax.set_xlabel("Time (h)")     plt.show()           plt.plot(A, B)     plt.show()           return