Ecuaciones Diferenciales
Ā¶

Rafael MartĆ­nez Fonseca (https://rpubs.com/Rafael31416/928321)

Objetivo : Dar soluciĆ³n numerica a sistemas de ecuaciones diferenciales no lineales

Paquetes requeridosĀ¶

InĀ [1]:
from numpy import*
from matplotlib.pyplot import*
from scipy.integrate import*

ODE de primer ordenĀ¶

Problema 1 $\frac{dx}{dt}=ax(t)$Ā¶

InĀ [24]:
a=0.01
T=arange(1,800,0.1)
x0=0.1
InĀ [25]:
def f(x,t):
    return a*x
InĀ [26]:
sol1=odeint(f,x0,T)
InĀ [27]:
plot(T,sol1)
Out[27]:
[<matplotlib.lines.Line2D at 0x282d3ee6280>]
InĀ [28]:
for ci in arange(0.1,2,0.1):
    sol1=odeint(f,ci,T)
    plot(T,sol1)

Problema 2: $\frac{dx}{dt}=ax(1-\frac{x}{k})$Ā¶

InĀ [29]:
k=20
def g(x,t):
    return a*x*(1-(x/k))
InĀ [36]:
sol2=odeint(g,2,T)
InĀ [32]:
odeint(g,2,T)
Out[32]:
array([[ 2.        ],
       [ 2.00180072],
       [ 2.00360289],
       ...,
       [19.93901318],
       [19.93907395],
       [19.93913466]])
InĀ [37]:
plot(sol2)
Out[37]:
[<matplotlib.lines.Line2D at 0x282d58e1880>]
InĀ [38]:
for ci in arange(2,40,2):
    sol2=odeint(g,ci,T)
    plot(sol2)

Sistemas de Ecuaciones DiferencialesĀ¶

Paquetes requeridosĀ¶

InĀ [44]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

Problema 1. EL oscilador ArmĆ³nicoĀ¶

$$\frac{d^2y}{dt^2}=a\frac{dy}{dt}+by+c$$
InĀ [56]:
def f(Y,t):
    y1,y2=Y
    return [y2,-k*y1-b*y2+4*np.cos(2*t)]
InĀ [57]:
k=10.
b=2.
InĀ [58]:
y0=[0.5,0]
InĀ [154]:
T=np.arange(0,20,0.1)
ys=odeint(f,y0,T)
InĀ [70]:
plt.plot(T,ys[:,0],label="PosiciĆ³n")
plt.plot(T,ys[:,1],label="Velocidad")
Out[70]:
[<matplotlib.lines.Line2D at 0x282d7623a90>]

Dibujo del plano fase

InĀ [71]:
plt.plot(ys[:,0],ys[:,1])
Out[71]:
[<matplotlib.lines.Line2D at 0x282d7691070>]

Con distintos valores iniciales

InĀ [94]:
T1=np.arange(0,1000,0.1)
for ci in [0,0.5,1,1.5,2]:
    y0=[ci,0]
    yh=odeint(f,y0,T1)
    plt.plot(yh[:,0],yh[:,1],label="%f"%ci)
plt.legend() #Leyenda de valores iniciales
plt,title("Trayectorias del plano fase")
plt.xlabel("PosiciĆ³n")
plt.ylabel("velocidad")
Out[94]:
Text(0, 0.5, 'velocidad')

Problema 2. Sistema Lotka-VolterraĀ¶

\begin{eqnarray} \frac{dx}{dt}&=&ax-\beta xy \\ \frac{dy}{dt}&=&-\gamma y + \delta xy \end{eqnarray}
InĀ [121]:
def g(Y,t):
    y1,y2=Y
    return [a*y1-b*y1*y2,-c*y1+d*y1*y2]
InĀ [122]:
a=1.8
b=0.9
c=0.81
d=0.58
InĀ [127]:
y0=[4,2]
gs=odeint(g,y0,T)
InĀ [128]:
plt.plot(T,gs[:,0],label="Presas")
plt.plot(T,gs[:,1],label="Depredadores")
plt,legend()
Out[128]:
(<module 'matplotlib.pyplot' from 'C:\\anaconda\\lib\\site-packages\\matplotlib\\pyplot.py'>,
 <matplotlib.legend.Legend at 0x282d9438880>)
InĀ [129]:
plt.plot(gs[:,0],gs[:,1])
Out[129]:
[<matplotlib.lines.Line2D at 0x282d94ba310>]

Problema 3. Sistema de LorenzĀ¶

InĀ [150]:
dt=0.01
num_steps =10000

xs=np.empty(num_steps+1)
ys=np.empty(num_steps+1)
zs=np.empty(num_steps+1)

xs[0],ys[0],zs[0]=(0.,1.,1.05)
InĀ [151]:
for i in range(num_steps):
    x_dot,y_dot,z_dot=lorenz(xs[i],ys[i],zs[i])
    xs[i+1]=xs[i]+(x_dot*dt)
    ys[i+1]=ys[i]+(y_dot*dt)
    zs[i+1]=zs[i]+(z_dot*dt)
InĀ [152]:
plt.plot(xs)
plt.plot(ys)
plt.plot(zs)
Out[152]:
[<matplotlib.lines.Line2D at 0x282da9131c0>]
InĀ [153]:
fig=plt.figure()
ax=fig.gca(projection="3d")
ax.plot(xs,ys,zs,lw=0.5)
plt,title("Atracctor de lorenz")

C:\Users\Rafael Martinez\AppData\Local\Temp\ipykernel_2372\3024403361.py:2: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax=fig.gca(projection="3d")
Out[153]:
(<module 'matplotlib.pyplot' from 'C:\\anaconda\\lib\\site-packages\\matplotlib\\pyplot.py'>,
 Text(0.5, 0.92, 'Atracctor de lorenz'))