Python Planeta.avanza_euler方法代码示例

# 需要导入模块: from planeta import Planeta [as 别名]
# 或者: from planeta.Planeta import avanza_euler [as 别名]
N_steps = 100000
dt=4000./N_steps
t=np.linspace(0,4000,N_steps)

#creamos los arrays en que se guardara la informacion
x= np.zeros(N_steps)
y= np.zeros(N_steps)
energia = np.zeros(N_steps)

x[0]= 10
y[0]= 0
energia[0]=p.energia_actual

for i in range(1, N_steps):
    p.avanza_euler(dt)
    x[i]=p.y_actual[0]
    y[i]=p.y_actual[1]
    p.energia_total()
    energia[i]=p.energia_actual

fig = plt.figure(1)
fig.clf()

ax1 = fig.add_subplot(211)
plt.suptitle('Trayectoria y energia vs tiempo con $v_{y}(t=0)=0.4$ y  ' r'$\alpha=0$ (Euler explicito)')
fig.subplots_adjust(hspace=.3)
ax1.plot(x,y)
ax1.grid(True)
ax1.set_xlabel('x')
ax1.set_ylabel('y')

# 需要导入模块: from planeta import Planeta [as 别名]
# 或者: from planeta.Planeta import avanza_euler [as 别名]
t_final =  1025.
numero_pasos = 10000
dt= t_final / (float)(numero_pasos)
t_values = np.linspace(0,t_final,numero_pasos)
x = np.zeros(numero_pasos)
y = np.zeros(numero_pasos)
vx = np.zeros(numero_pasos)
vy = np.zeros(numero_pasos)

energia = np.zeros(numero_pasos)

[x[0],y[0],vx[0],vy[0]] = condicion_inicial
energia[0] = Jupiter.energia_total()
for i in range (1,numero_pasos):
    #pdb.set_trace()
    Jupiter.avanza_euler(dt)
    resultados = Jupiter.y_actual
    x[i] = resultados[0]
    y[i] = resultados[1]
    vx[i] = resultados[2]
    vy[i] = resultados[3]
    energia[i] = Jupiter.energia_total()

fig=plt.figure(1)
plt.subplot(2, 1, 1)
fig.subplots_adjust(hspace=.5)
plt.plot(x , y)
plt.axis('equal')
plt.title(' 'u'Ó''rbita descrita por planeta, m'u'é''todo de Euler')
plt.xlabel("x [m]")
plt.ylabel("y [m]")

# 需要导入模块: from planeta import Planeta [as 别名]
# 或者: from planeta.Planeta import avanza_euler [as 别名]
def orbitar(CI,solucion, prt=True):
    '''
    Esta función  es cool, porque le entregas las condiciones iniciales (CI)
    en forma de x0,y0,vx0,vy0, además de un string con el tipo de solución
    que quieres, como 'euler', 'rk4' o 'verlet' y nos devuelve las coordenadas
    x,y y la energia total del planeta. Además de un ploteo de la órbita, apropi-
    adamente labelado y un ploteo de la energía vs el tiempo.
    La ultima parte es para que plotee, el valor default es que
    plotee (true or false).
    Ejemplo de uso:
        orbitar([10,0,0,0.1],'verlet')
        x=[0,0.1,...]
        y=[0,0.3,...]
        energia=[7,8,...]
    '''
    if solucion=='euler':
        n= 2500 #grande y dsps se arregla
    if solucion=='rk4':
        n= 5000 #grande y dsps se arregla
    if solucion=='verlet':
        n=800
    Aiur= Planeta(CI) #se crea el planeta
    if solucion=='verlet_reloaded': #caso especial para la ultima parte
        n=6000
        Aiur= Planeta(CI,alpha=alpha2) #si es el verlet bkn se parcha
    dt= 1 #error si es muy chico
    tmax=n*dt
    x=[] #listas vacias a llenar
    y=[]
    pylon=[] #energia

    x= np.append(x,CI[0]) #llenamos con el primer termino
    y= np.append(y,CI[1])
    pylon= np.append(pylon,Aiur.energia_total())

    for i in range(tmax):   #aca se hace la iteracion del paso
        if solucion=='euler':
            Aiur.avanza_euler(dt)
        if solucion=='rk4':
            Aiur.avanza_rk4(dt)
        if solucion=='verlet' or solucion=='verlet_reloaded':
            Aiur.avanza_verlet(dt)
        xn,yn,vxn,vyn=Aiur.y_now
        x= np.append(x,xn)
        y= np.append(y,yn)
        pylon= np.append(pylon,Aiur.energia_total())
        times=np.linspace(0,tmax,tmax+1)
    if prt==True:
        fig1=plt.figure(1)  #aca se plotea, muy standar
        fig1.clf()
        if solucion=='euler':
            plt.plot(x,y,'r',label= 'Trayectoria metodo Euler')
            plt.title('Trayectoria metodo Euler')
        if solucion=='rk4':
            plt.plot(x,y,'r',label= 'Trayectoria metodo RK4')
            plt.title('Trayectoria metodo RK4')
        if solucion=='verlet' or solucion=='verlet_reloaded':
            plt.plot(x,y,'r',label= 'Trayectoria metodo Verlet')
            plt.title('Trayectoria metodo Verlet')

        plt.xlabel('X')
        plt.ylabel('Y')
        plt.savefig('verlet_rld.png')
        plt.draw()
        plt.show()

        fig2=plt.figure(2)
        fig2.clf()
        plt.plot(times, pylon,'g')
        plt.xlabel('Tiempo')
        plt.ylabel('Energia')
        plt.title('Tiempo vs Energia')
        plt.savefig('verlet_rld_energia.png')
        plt.draw()
        plt.show()
    return x,y,pylon

# 需要导入模块: from planeta import Planeta [as 别名]
# 或者: from planeta.Planeta import avanza_euler [as 别名]
import numpy as np
import matplotlib.pyplot as plt

#Solucion a trayectoria
condicion_inicial = [10.0, 0.0, 0.0, 0.3]

Mars = Planeta(condicion_inicial)
x_pos = []
y_pos = []
energia = []
dt=0.1
for i in range(10000):
    x_pos.append(Mars.y_actual[0])
    y_pos.append(Mars.y_actual[1])
    energia.append(Mars.energia_total())
    Mars.avanza_euler(dt)

x_pos = np.array(x_pos)
y_pos = np.array(y_pos)
energia = np.array(energia)

#Resultados
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.plot(x_pos,y_pos)
ax.plot(0,0,'o',color='r',linewidth=2.0)
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_xlim([-10,12])
ax.set_ylim([-12,12])
ax.set_aspect('equal')

# 需要导入模块: from planeta import Planeta [as 别名]
# 或者: from planeta.Planeta import avanza_euler [as 别名]
N = 1000
t_paso= t_total/N

x = np.zeros(N); y = np.zeros(N)
vx = np.zeros(N); vy = np.zeros(N)

E = np.zeros(N)

[x[0],y[0],vx[0],vy[0]] = CI

E[0] = p.energia_total()
t= np.linspace(1,t_total,N)

for i in range (1,N):

    p.avanza_euler(t_paso)
    x[i] = p.y_actual[0]
    y[i] = p.y_actual[1]
    vx[i] = p.y_actual[2]
    vy[i] = p.y_actual[3]
    E[i] = p.energia_total()

fig=plt.figure(1)
fig.subplots_adjust(hspace=0.4)

plt.subplot(2, 1, 1)
plt.plot(x, y, 'c-')
plt.title("Trayectoria potencial central (Metodo Euler)")
plt.xlabel("x")
plt.ylabel("y")

# 需要导入模块: from planeta import Planeta [as 别名]
# 或者: from planeta.Planeta import avanza_euler [as 别名]
condicion_inicial = [10, 0, 0, 0.2]

p = Planeta(condicion_inicial,alpha)
print "Energia Total", p.energia_total()
contador = 0

x=np.zeros(100001)
y=np.zeros(100001)
t=np.linspace(0,2000,100001)
E_total = np.zeros(100001)

while  contador<=100000:
    E_total[contador] = p.energia_total()
    x[contador] = p.y_actual[0]
    y[contador] = p.y_actual[1]
    p.avanza_euler(0.02)
    contador +=1


plt.xlabel("X")
plt.ylabel("Y")
plt.title("Orbitas con Metodo de Euler")
plt.plot(x,y,'-')
plt.savefig("Orbitas con el Metodo de Euler")
plt.show()

plt.xlabel("Tiempo")
plt.ylabel("Energia Total")
plt.title("Energia total del sistema como funcion del Tiempo")
plt.plot(t,E_total,'-')
plt.savefig("Energia con el Metodo de Euler")