Python数值分析--------模拟Henon-Heiles系统

代码段

#Henon-Heiles系统
#V(x,y)=1/2(x^2+y^2)+x^2y-1/3y^3
import matplotlib.pyplot as plt
from math import *
import decimal

decimal.getcontext().prec = 10000000

time = 20
h = 0.001
step = int(time/h)
x0,y0 = 0,0
px0,py0 = -0.1,0.1

x = [x0] + [0. for i in range(step)]
y = [y0] + [0. for i in range(step)]
px = [px0] + [0. for i in range(step)]
py = [py0] + [0. for i in range(step)]


order_4x = [x0] + [0 for i in range(step)]
order_4y = [y0] + [0 for i in range(step)]
order_4px = [px0] + [0 for i in range(step)]
order_4px = [py0] + [0 for i in range(step)]

t = [i*h for i in range(1+step)]
for i in range(step):
    try:
        x[i+1] = x[i] + h * px[i]
        y[i+1] = y[i] + h * py[i]
        px[i+1] = px[i] - h * (x[i] + 2 * x[i] * y[i])
        py[i+1] = py[i] - h * (y[i] + x[i]**2 - y[i]**2)
    except:
        print(i)
        print(x[i])
        break

    
plt.plot(x,y,label = "euler method")
plt.legend()
plt.show()