выдержки:
import numpy as np
from scipy.integrate import ode
from scipy.interpolate import interp1d
from pylab import plot, show, gca, savefig, clf
from pprint import pprint as pp
import time
from collections import defaultdict
class XeSys(object):
u"""
"""
def __init__(self,Wr):
self.Wr = Wr
def ystat(self,Wr, y0):
u"""стационарные Значения"""
y0[iXe] = Wr /(1 + Wr)
y0[iId] = Wr
def __call__(self, t, y):
Xe,Id = y
Wr = self.Wr(t)
return [
L_Xe * (-Xe + Id - Wr * Xe) ,
L_Id * (-Id + Wr)
]
sys = XeSys(Wrspline)
r = ode(sys)
r.set_integrator('dopri5')
print "y0=",self.y0
r.set_initial_value(self.y0, self.t0)
for i, t in enumerate(self.time_points):
if t != r.t:
r.integrate(t)
if not r.successful():
print "t", t, "y", r.y, "x", x