Kernel: Python 3 (system-wide)
In [2]:
import numpy as np import matplotlib.pyplot as plt def rungekutta2_fg(f,g,t0,x0,y0,h,muestras): tamaño = muestras tabla = np.zeros(shape= (tamaño, 3), dtype=float) tabla[0] = [t0,x0,y0] ti = t0 xi = x0 yi = y0 for i in range(1, tamaño, 1): K1x = h*f(ti,xi,yi) K1y = h*g(ti,xi,yi) K2x = h*f(ti+h, xi + K1x, yi+K1y) K2y = h*g(ti+h, xi + K1x, yi+K1y) xi = xi + (1/2)*(K1x + K2x) yi = yi + (1/2)*(K1y + K2y) ti = ti + h tabla[i] = [ti,xi,yi] tabla = np.array(tabla) return(tabla) a = 0.5471 b = 0.0481 c = 0.001 d = 0.0466 e = 0.8439 s = 0.001 f = lambda t,x,y : a*x - b*x*y - c*x**2 g = lambda t,x,y : d*x*y - e*y - s*y**2 t0 = 0 x0 = 0 y0 = 0 h = 0.0001 muestras = 1000 tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras) ti = tabla[: , :] xi = tabla[: , :] yi = tabla[: , :] np.set_printoptions(precision=6) print(' [ ti, xi, yi] ') print(tabla)
Out[2]:
[ ti, xi, yi]
[[0. 0. 0. ]
[0.0001 0. 0. ]
[0.0002 0. 0. ]
...
[0.0997 0. 0. ]
[0.0998 0. 0. ]
[0.0999 0. 0. ]]
In [3]:
#grafica tiempo vs población plt.plot(ti,xi, label='xi presa') plt.plot(ti,yi, label='yi depredador') plt.title('Presa-Depredador') plt.xlabel('t tiempo') plt.ylabel('población') plt.legend() plt.grid() plt.show()
Out[3]:
In [4]:
#grafica xi vs yi plt.plot(xi,yi) plt.title('Presa-Depredador [xi,yi]') plt.xlabel('x presa') plt.ylabel('y depredador') plt.grid() plt.show()
Out[4]:
In [0]:
In [40]:
def rungekutta4_fg(f,g,t0,x0,y0,h,muestras): tamaño = muestras tabla = np.zeros(shape= (tamaño, 3), dtype=float) tabla[0] = [t0,x0,y0] ti = t0 xi = x0 yi = y0 for i in range(1, tamaño, 1): K1x = h*f(ti,xi,yi) K1y = h*g(ti,xi,yi) K2x = h*f(ti+0.5*h, xi + h*0.5*K1x, yi + h*0.5*K1y) K2y = h*g(ti+0.5*h, xi + h*0.5*K1x, yi + h*0.5*K1y) K3x = h*f(ti+0.5*h, xi + h*0.5*K2x, yi + h*0.5*K2y) K3y = h*g(ti+0.5*h, xi + h*0.5*K2x, yi + h*0.5*K2y) K4x = h*f(ti + 1, xi + h*K3x, yi + h*K3y) K4y = h*g(ti + 1, xi + h*K3x, yi + h*K3y) xi = xi + (h*(K1x + 2*K2x + 2*K3x + K4x))/6.0 yi = yi + (h*(K1y + 2*K2y + 2*K3y + K4y))/6.0 ti = ti + h tabla[i] = [ti,xi,yi] tabla = np.array(tabla) return(tabla) a = 0.5471 b = 0.0481 c = 0.001 d = 0.0466 e = 0.8439 s = 0.001 f = lambda t,x,y : a*x - b*x*y - c*x**2 g = lambda t,x,y : d*x*y - e*y - s*y**2 t0 = 0 x0 = 1 y0 = 2 h = 0.0001 muestras = 1000 tabla = rungekutta4_fg(f,g,t0,x0,y0,h,muestras) ti = tabla[: , :] xi = tabla[: , :] yi = tabla[: , :] np.set_printoptions(precision=6) print(' [ ti, xi, yi] ') print(tabla)
Out[40]:
[ ti, xi, yi]
[[0.000000e+00 1.000000e+00 2.000000e+00]
[1.000000e-04 1.000000e+00 2.000000e+00]
[2.000000e-04 1.000000e+00 2.000000e+00]
...
[9.970000e-02 1.000004e+00 1.999984e+00]
[9.980000e-02 1.000004e+00 1.999984e+00]
[9.990000e-02 1.000004e+00 1.999984e+00]]
In [41]:
plt.plot(ti,xi, label='xi presa') plt.plot(ti,yi, label='yi depredador') plt.title('Presa-Depredador') plt.xlabel('t tiempo') plt.ylabel('población') plt.legend() plt.grid() plt.show()
Out[41]:
In [39]:
plt.plot(xi,yi) plt.title('Presa-Depredador [xi,yi]') plt.xlabel('x presa') plt.ylabel('y depredador') plt.grid() plt.show()
Out[39]:
In [0]:
In [30]:
import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy import integrate import ipywidgets as ipw
In [31]:
alpha = 0.5471 #mortality rate due to predators beta = 0.0481 gamma = 0.001 delta = 0.0466 epsilon = 0.8439 zeta = 0.001 x0 = 10 y0 = 5 def derivative(X, t, alpha, beta, gamma, delta, epsilon, zeta): x, y = X dotx = alpha*x - beta*x*y - gamma*x**2 doty = delta*x*y - epsilon*y - zeta*y**2 return np.array([dotx, doty])
In [42]:
Nt = 1000 tmax = 30 t = np.linspace(0.,tmax, Nt) X0 = [x0, y0] res = integrate.odeint(derivative, X0, t, args = (alpha, beta, gamma, delta, epsilon, zeta)) x, y = res.T
In [43]:
plt.figure() plt.grid() plt.title("odeint method") plt.plot(t, x, 'xb', label = 'Depredadores') plt.plot(t, y, '+r', label = "Presas") plt.xlabel('Time t, [days]') plt.ylabel('Population') plt.legend() plt.show()
Out[43]:
In [0]:
In [10]:
import matplotlib.pyplot as plt from numpy import arange
In [33]:
def f(x,t): return alpha*x - beta*x*t - gamma*x**2 a, b = 0.0, 50.0 N = 1000 h = 0.0001 x = 1.0 tpoints = arange(a, b, h) rk4_xpoints = [] for t in tpoints: rk4_xpoints.append(x) k1 = h*f(x, t) k2 = h*f(x+0.5*k1, t+0.5*h) k3 = h*f(x+0.5*k2, t+0.5*h) k4 = h*f(x+k3, t+h) x += (k1+2*k2+2*k3+k4)/6.0 plt.plot(tpoints, rk4_xpoints, label="RK4") plt.xlabel("$t$") plt.ylabel("$x(t)$") plt.legend() plt.show()
Out[33]:
In [0]:
for t in tpoints: rk2_xpoints.append(x) k1 = h*f(x, t) k2 = h*f(x+0.5*k1, t+0.5*h) x += k2 x = 10.0 plt.plot(tpoints, rk2_xpoints, label="RK2")
In [0]:
for t in tpoints: euler_xpoints.append(x) x += h*f(x, t) x = 5.0 plt.plot(tpoints, euler_xpoints, label="Euler")
In [0]:
In [34]:
delta = 0.0466 epsilon = 0.8439 zeta = 0.001 def f(y,t): return delta*t*y - epsilon*y - zeta*y**2 a, b = 0.0, 60.0 N = 1000 h = 0.0001 x = 1.0 tpoints = arange(a, b, h) rk4_xpoints = [] for t in tpoints: rk4_xpoints.append(x) k1 = h*f(x, t) k2 = h*f(x+0.5*k1, t+0.5*h) k3 = h*f(x+0.5*k2, t+0.5*h) k4 = h*f(x+k3, t+h) x += (k1+2*k2+2*k3+k4)/6.0 plt.plot(tpoints, rk4_xpoints, label="RK4") plt.xlabel("$t$") plt.ylabel("$x(t)$") plt.legend() plt.show()
Out[34]:
In [0]: