Carlettos

volterra 3d

Sep 27th, 2020
963
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. from operator import add
  4.  
  5. a = [1, 0.72, 1.53, 1.27]
  6. b = [[1, 1.09, 1.52, 0],
  7.      [0, 1, 0.49, 1.36],
  8.      [2.33, 0, 1, 0.47],
  9.      [1.21, 0.51, 0.35, 1]]
  10. alturas = [0.3013, 0.4586, 0.1307, 0.3557]
  11.  
  12. # cantidad de variables
  13. cv = len(a)
  14. funciones = [lambda n, i=i: a[i]*n[i]*(1 - sum([b[i][j]*n[j] for j in range(cv)])) for i in range(cv)]
  15.  
  16. inicio = 0
  17. final = 1000
  18. N = 10_000
  19.  
  20. tiempos = np.linspace(inicio, final, N)
  21. h = tiempos[1] - tiempos[0]
  22. recorrido = [[] for _ in range(cv)]
  23.  
  24. for t in tiempos:
  25.     for i, altura in enumerate(alturas):
  26.         recorrido[i].append(altura)
  27.     ks = [[] for _ in range(cv)]  # lista de k
  28.     for i, k in enumerate(ks):  # k1
  29.         k.append(h*funciones[i](alturas))
  30.     for j in range(0, 2):  # k2 y k3
  31.         for i, k in enumerate(ks):
  32.             k.append(h*funciones[i](list(map(add, alturas, [ks[w][j] for w in range(cv)]))))
  33.     for i, k in enumerate(ks):  # k4
  34.         k.append(h*funciones[i](list(map(add, alturas, [ks[w][j] for w in range(cv)]))))
  35.     for i, k in enumerate(ks):  # suma final
  36.         alturas[i] += (k[0] + 2*k[1] + 2*k[2] + k[3])/6
  37.     pass
  38.  
  39. for i, r in enumerate(recorrido):
  40.     plt.plot(tiempos, r, label="n{}".format(i))
  41. plt.legend(loc="best")
  42. plt.xlabel("t")
  43. plt.ylabel("n(t)")
  44. plt.show()
  45.  
  46. fig = plt.figure()
  47. ax = fig.add_subplot(111, projection="3d")
  48. ax.plot(recorrido[0], recorrido[1], recorrido[2])
  49. ax.set_xlabel('n1')
  50. ax.set_ylabel('n2')
  51. ax.set_zlabel('n3')
  52. plt.show()
  53.  
RAW Paste Data