Python-сообщество

 from math import sin, cos, pi
from numpy import matrix, array, arange
from pylab import *
g = 9.78
dt = 0.01
init_conds = [1,2,3,4]
end = 10
def constrain(theta):
    theta = theta % (2*pi)
    if theta > pi:
        theta = -2*pi+theta
        print(theta)
    return theta
def average(x):
    x_i, k1, k2, k3, k4 = x
    return x_i + (k1 + 2.0*(k3 + k4) +  k2) / 6.0
theta = []
class Pendulum(object):
    def __init__(self, dt, init_conds, end):
        self.dt = dt
        self.t = 0.0
        self.x = init_conds[:]
        self.end = end
    def derivative(self, u):
        #x1 = x, x2 = x_dt, x3 = theta, x4 = theta_dt
        x1, x2, x3, x4 = u
        x1_dt, x3_dt = x2, x4
        x2_dt = (sin(x3)*x4**2 - g*cos(x3)*sin(x3)) / (1 - cos(x3)**2)
        x4_dt = (g*sin(x3) - x4**2*cos(x3)*sin(x3) - cos(x3)*x2) / (1 - cos(x3)**2)
        x = [x1_dt, x2_dt, x3_dt, x4_dt]
        return x
    def rk4_step(self, dt):
        dx = self.derivative(self.x)
        k2 = [ dx_i*dt for dx_i in dx ]
        xv = [x_i + delx0_i/2.0 for x_i, delx0_i in zip(self.x, k2)]
        k3 = [ dx_i*dt for dx_i in self.derivative(xv)]
        xv = [x_i + delx1_i/2.0 for x_i,delx1_i in zip(self.x, k3)]
        k4 = [ dx_i*dt for dx_i in self.derivative(xv) ]
        xv = [x_i + delx1_2 for x_i,delx1_2 in zip(self.x, k4)]
        k1 = [self.dt*i for i in self.derivative(xv)]
        self.t += dt
        self.x = list(map(average, zip(self.x, k1, k2, k3, k4)))
        theta.append(constrain(self.x[2]))
        def integrate(self):
            x = []
            while self.t <= self.end:
                self.rk4_step(self.dt)
                x.append([self.t] + self.x)
            return array(x)
p = Pendulum(dt, init_conds, end)
print(p.integrate())