Skip to content
GitLab
Commit 9cbab1b3 authored by AMKalinin's avatar AMKalinin
Browse files

create

parents
Hide whitespace changes
Inline Side-by-side
Integ.py 0 → 100644
View file @ 9cbab1b3
import math
import Lin
import model
class TIntegrator():
def __init__(self, h, t0, tk, model:model.CentralGravity):
self.h = h
self.t0 = t0
self.tk = tk
self.model = model
def integrate(self):
pass
class TRunge(TIntegrator):
def __init__(self, h, t0, tk, model):
super().__init__(h, t0, tk, model)
def integrate(self):
C0 = 0.0
C1 = 1.0 / 5.0
C2 = 3.0 / 10.0
C3 = 4.0 / 5.0
C4 = 8.0 / 9.0
C5 = 1.0
C6 = 1.0
O = 1
b1 = O * (1 + O * (-1337 / 480 + O * (1039 / 360 + O * (-1163 / 1152))))
b2 = 0
b3 = 100 * (O ** 2) * (1054 / 9275 + O * (-4682 / 27825 + O * (379 / 5565))) / 3
b4 = -5 * (O ** 2) * (27 / 40 + O * (-9 / 5 + O * (83 / 96))) / 2
b5 = 18225 * (O ** 2) * (-3 / 250 + O * (22 / 375 + O * (-37 / 600))) / 848
b6 = -22 * (O ** 2) * (-3 / 10 + O * (29 / 30 + O * (-17 / 24))) / 7
A00 = 0.0
A10 = 1.0 / 5.0
A20 = 3.0 / 40.0
A21 = 9.0 / 40.0
A30 = 44.0 / 45.0
A31 = -56.0 / 15.0
A32 = 32.0 / 9.0
A40 = 19372.0 / 6561.0
A41 = -25360.0 / 2187.0
A42 = 64448.0 / 6561.0
A43 = -212.0 / 729.0
A50 = 9017.0 / 3168.0
A51 = -355.0 / 33.0
A52 = 46732.0 / 5247.0
A53 = 49.0 / 176.0
A54 = -5103.0 / 18656.0
A60 = 35.0 / 384.0
A61 = 0.0
A62 = 500.0 / 1113.0
A63 = 125.0 / 192.0
A64 = -2187.0 / 6784.0
A65 = 11.0 / 84.0
t1 = self.t0
t = self.t0
x0 = self.model.v.Clone()
self.model.v[0].append(0.0)
while (t <= self.tk):
self.InitH = self.h + 0
tmp = x0.Clone()
tmp = self.model.calc(t, tmp)
k1 = tmp.Clone() * self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i] * (0.2)
tmp = self.model.calc(t + (0.2) * self.h, tmp)
k2 = tmp.Clone() * self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i] * (A20) + k2[0][i] * (A21)
tmp = self.model.calc(t + (0.3) * self.h, tmp)
k3 = tmp.Clone() * self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i] * (A30) + k2[0][i] * (A31) + k3[0][i] * (A32)
tmp = self.model.calc(t + (0.8) * self.h, tmp)
k4 = tmp.Clone() * self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i] * (A40) + k2[0][i] * (A41) + k3[0][i] * (A42) + k4[0][i] * (A43)
tmp = self.model.calc(t + (C4) * self.h, tmp)
k5 = tmp.Clone() * self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i] * (A50) + k2[0][i] * (A51) + k3[0][i] * (A52) + k4[0][i] * (A53) + \
k5[0][i] * (A54)
tmp = self.model.calc(t + 1 * self.h, tmp)
k6 = tmp.Clone() * self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i] * (A60) + k2[0][i] * (A61) + k3[0][i] * (A62) + k4[0][i] * (A63) + \
k5[0][i] * (A64) + k6[0][i] * (A65)
tmp = self.model.calc(t + 1 * self.h, tmp)
k7 = tmp.Clone() * self.h
xt1 = x0 + k1 * (b1) + k2 * (b2) + k3 * (b3) + k4 * (b4) + k5 * (b5) + k6 * (b6) + k7 * (0)
x0 = xt1.Clone()
xt1[0].append(t)
self.model.v.add(xt1)
t += self.h
return self.model.v
class TDormanPrinceIntegrate(TIntegrator):
def __init__(self, h, t0, tk, model,Emax):
super().__init__( h, t0, tk, model)
self.Emax = Emax
def E(self, x0,x4,x5):
E1 = 0
for i in range(self.model.p):
E1+= (self.InitH*(x4[0][i]-x5[0][i])/max( max( math.pow(10,-5), math.fabs(x4[0][i])), max(math.fabs(x0[0][i]), 2*self.U()/self.Emax)))**2
return math.sqrt(E1/self.model.p)
def StepNew(self, x0,x4,x5):
return (self.h/max(0.1, min(5,math.pow((self.E(x0,x4,x5)/self.Emax),1/5)/0.9)))
def U (self):
v = 1
while (1+v>1):
u = v
v = v/2
return u
def integrate(self):
C0 = 0.0
C1 = 1.0 / 5.0
C2 = 3.0 / 10.0
C3 = 4.0 / 5.0
C4 = 8.0 / 9.0
C5 = 1.0
C6 = 1.0
B10 = 35.0 / 384.0
B11 = 0.0
B12 = 500.0 / 1113.0
B13 = 125.0 / 192.0
B14 = -2187.0 / 6784.0
B15 = 11.0 / 84.0
B16 = 0.0
B20 = 5179.0 / 57600.0
B21 = 0.0
B22 = 7571.0 / 16695.0
B23= 393.0 / 640.0
B24 = -92097.0 / 339200.0
B25 = 187.0 / 2100.0
B26 = 1.0 / 40.0
A00 = 0.0
A10 = 1.0 / 5.0
A20 = 3.0 / 40.0
A21 = 9.0 / 40.0
A30 = 44.0 / 45.0
A31 = -56.0 / 15.0
A32 = 32.0 / 9.0
A40= 19372.0 / 6561.0
A41 = -25360.0 / 2187.0
A42 = 64448.0 / 6561.0
A43 = -212.0 / 729.0
A50 = 9017.0 / 3168.0
A51 = -355.0 / 33.0
A52 = 46732.0 / 5247.0
A53 = 49.0 / 176.0
A54 = -5103.0 / 18656.0
A60 = 35.0 / 384.0
A61 = 0.0
A62 = 500.0 / 1113.0
A63 = 125.0 / 192.0
A64 = -2187.0 / 6784.0
A65 = 11.0 / 84.0
x =[]
y = []
t1 = self.t0
t = self.t0
CH = self.h+0
x0 = self.model.v.Clone()
while (t <= self.tk):
while True:
self.InitH = self.h+0
tmp = x0.Clone()
tmp = self.model.calc(t,tmp)
k1 = tmp.Clone()*self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i]*(0.2)
tmp = self.model.calc(t+(0.2)*self.h, tmp)
k2 = tmp.Clone()*self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i]*(A20) + k2[0][i]*(A21)
tmp = self.model.calc(t+(0.3)*self.h, tmp)
k3 = tmp.Clone()*self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i]*(A30) + k2[0][i]*(A31)+k3[0][i]*(A32)
tmp = self.model.calc(t+(0.8)*self.h, tmp)
k4 = tmp.Clone()*self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i]*(A40) + k2[0][i]*(A41) + k3[0][i]*(A42) + k4[0][i]*(A43)
tmp = self.model.calc(t+(C4)*self.h ,tmp)
k5 = tmp.Clone()*self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i]*(A50) + k2[0][i]*(A51) + k3[0][i]*(A52) + k4[0][i]*(A53) + k5[0][i]*(A54)
tmp = self.model.calc(t+1*self.h ,tmp)
k6 = tmp.Clone()*self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i]*(A60) + k2[0][i]*(A61) + k3[0][i]*(A62) + k4[0][i]*(A63) + k5[0][i]*(A64) + k6[0][i]*(A65)
tmp = self.model.calc(t+1*self.h ,tmp)
k7 = tmp.Clone()*self.h
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i]*(B10) + k2[0][i]*(B11) + k3[0][i]*(B12) + k4[0][i]*(B13) + k5[0][i]*(B14) + k6[0][i]*(B15) #+ k7[0][i]*(B16) b16 = 0
x4 = tmp.Clone()
for i in range(self.model.p):
tmp[0][i] = x0[0][i] + k1[0][i]*(B20) + k2[0][i]*(B21) + k3[0][i]*(B22) + k4[0][i]*(B23) + k5[0][i]*(B24) + k6[0][i]*(B25) +k7[0][i]*(B26)
x5 = tmp.Clone()
self.h = self.StepNew(x0, x4, x5)
if (self.E(x0,x4,x5) < self.Emax):
break
if (t1 < self.InitH +t):
while(t1 < self.InitH +t):
O = (t1 - t)/self.InitH
b1 = O*(1+O*(-1337/480+O*(1039/360+O*(-1163/1152))))
b2 = 0
b3 = 100*(O**2)*(1054/9275+O*(-4682/27825+O*(379/5565)))/3
b4 = -5*(O**2)*(27/40+O*(-9/5+O*(83/96)))/2
b5 = 18225*(O**2)*(-3/250+O*(22/375+O*(-37/600)))/848
b6 = -22*(O**2)*(-3/10+O*(29/30+O*(-17/24)))/7
xt1 = x0 + k1*(b1) + k2*(b2) + k3*(b3) + k4*(b4) + k5*(b5) + k6*(b6) +k7*(0)
xt1[0].append(t1)
self.model.v.add(xt1)
t1 += CH
x0 = x4.Clone()
t += self.h
self.model.v[0].append(0.0)
return self.model.v
Lin.py 0 → 100644
View file @ 9cbab1b3
import math
class Vector:
def __str__(self):
return '(' + ', '.join(map(str, self.v)) + ')'
def __init__(self, *par):
self.v= list(par)
# Функция получения значения элемента
def __getitem__(self, key):
return self.v[key]
# установка значения элемента
def __setitem__(self,key,value):
self.v[key] = value
# Получение размера вектора
def __len__(self):
return len(self.v)
# Получение индекса последнего елемента
def GetHigh(self):
return len(self) - 1
# Создание копии объекта
def Clone(self):
return Vector(*self[:])
# Изменение размера вектора
def Resize(self, size):
self
# Сложение векторов
def __add__(self, other):
tmp = []
for i in range(len(self)):
tmp.append(self[i] + other[i])
return Vector(*tmp)
# Вычетание векторов
def __sub__(self, other):
tmp = []
for i in range(len(self)):
tmp.append(self[i] - other[i])
return Vector(*tmp)
# Умножение вектора на число / матрицу / скалярное
def __mul__(self, other):
tmp = []
t=0
if( type(other) == Vector ):
for i in range(len(self)):
t += other[i] * self[i]
return t
elif (type(other) == Matrix):
s=0 #сумма
t=[] #временная матрица
m3=[] # конечная матрица
r1 = len(self) #количество строк в первой матрице
c2 = other.ColCount # количество столбцов во 2ой матрице
if (r1 == c2):
c1 = 1 #Количество столбцов в 1
for z in range(0,r1):
for j in range(0,c2):
for i in range(0,c1):
s=s+self[z]*other.m[i][j]
t.append(s)
s=0
m3.append(t)
t=[]
return Matrix(*m3)
else :
for i in range(len(self)):
tmp.append(other * self[i])
return Vector(*tmp)
#Векторное произведение векторов
def Cross(self,other):
tmp = [[1,1,1]]
n = len(self)
for i in range(n-1):
t = []
for j in range(n):
if i==0:
t.append(self[j])
else:
t.append(other[j])
tmp.append(t)
z = []
m = Matrix(*tmp)
for i in range(n):
z.append(((-1)**i)*(m.Minor(0,i).Det()))
return Vector(*z)
# Получение модуля вектора
def Length(self):
tmp = 0
for i in range(len(self)):
tmp += self.v[i]**2
tmp = tmp**(1/2)
return tmp
# Нормирование ввектора
def Norm(self):
tmp = []
for i in range(len(self)):
tmp.append( self[i] / self.Length())
self.v = tmp
def rotateByRodrig(self, vec, phi):
phi1 = phi
vec.Norm()
tmp = vec * (vec * self) * (1-math.cos(phi1)) + vec.Cross(self) * math.sin(phi1) + self * math.cos(phi1)
return tmp
def rotateByQuaternion(self, quat):
tmp =(quat * self * quat.conj())
return Vector(tmp[1], tmp[2], tmp[3],)
# Свойство для чтения размера вектора
Size = property(__len__)
# Для чтения последнего элемента
High = property(GetHigh)
class Matrix:
def __init__(self, *par):
self.m = list(par)
def __str__(self):
return '\n'.join(map(str, self.m))
# Обращение матриц
def Inv(self):
tmp = self.Stac()
for k in range(tmp.RowCount):
swap_row = tmp.pick_nonzero_row(k)
if swap_row != k:
tmp[k][:] , tmp[swap_row][:] = tmp[swap_row][:] , tmp.Clone()[k][:]
if tmp[k][k] != 1:
a = tmp[k][k]
for i in range(tmp.ColCount):
tmp[k][i] = tmp[k][i] / a
for row in range(k + 1, tmp.RowCount):
a = tmp[row][k]
for i in range(tmp.ColCount):
tmp[row][i] -= tmp[k][i] * a
for k in range(tmp.RowCount-1, 0, -1):
for row in range(k-1, -1, -1):
if tmp[row][k]:
a = tmp[row][k]
for i in range(tmp.ColCount):
tmp[row][i] = tmp[row][i] - tmp[k][i] * a
tmp.Split()
return tmp
def Stac(self):
tmp = self.Clone()
for i in range(self.RowCount):
for j in range(self.ColCount):
if (i == j):
tmp[i].append(1)
else:
tmp[i].append(0)
return tmp
def Split(self):
n = self.RowCount
for i in range(n):
j = n -1
while (j>=0):
self[i].pop(j)
j -= 1
def pick_nonzero_row(self, k):
while k < self.RowCount and not self[k][k]:
k += 1
return k
# Функция получения значения элемента
def __getitem__(self, key):
return self.m[key]
# Процедура установки значения элемента
def __setitem__(self,key,value):
self.m[key] = value
# Функция получения количества строк
def GetRowCount(self):
return len(self.m)
# Функция получения количества столбцов
def GetColCount(self):
return len(self.m[0])
# Получение индекса последней строки
def GetRowHigh(self):
return len(self.m) - 1
# Получение индекса последнего столбца
def GetColHigh(self):
return len(self.m[0]) - 1
# Создание копии объекта
def Clone(self):
t = []
for i in range(self.RowCount):
tmp = []
for j in range (self.ColCount):
tmp.append(self[i][j])
t.append(tmp)
return Matrix(*t[:])
# Изменение размерности
def Resize(self, n, m):
return
# Добавить строку
def add(self,a):
self.m.append(*a)
# Сложение матриц
def __add__(self, other):
tmp = self.Clone()
for i in range(self.RowCount):
for j in range (self.ColCount):
tmp[i][j] += other[i][j]
return tmp
# Вычитание матриц
def __sub__(self, other):
tmp = self.Clone()
for i in range(self.RowCount):
for j in range (self.ColCount):
tmp[i][j] -= other[i][j]
return tmp
# Умножение на матрицу / число / вектор
def __mul__(self, other):
if (type(other) == Vector):
tmp = []
t = self.Clone()
for i in range(self.RowCount):
c = 0
for j in range(self.ColCount):
t[i][j] *= other[j]
c += t[i][j]
tmp.append(c)
return Vector(*tmp)
elif(type(other) == Matrix):
s=0 #сумма
t=[] #временная матрица
m3=[] # конечная матрица
if other.RowCount == self.ColCount:
r1 = self.RowCount #количество строк в первой матрице
c1 = self.ColCount #Количество столбцов в 1
c2 = other.ColCount # количество столбцов во 2ой матрице
for z in range(0,r1):
for j in range(0,c2):
for i in range(0,c1):
s=s+self[z][i]*other[i][j]
t.append(s)
s=0
m3.append(t)
t=[]
return Matrix(*m3)
else:
tmp = self.Clone()
for i in range(self.RowCount):
for j in range(self.ColCount):
tmp[i][j] *= other
return tmp
# Вычисление определителя
def Det(self):
c = self.ColCount
r = self.RowCount
if (c != r):
return 'error'
if c == 2:
return self[0][0] * self[1][1] - self[0][1] * self[1][0]
else:
return sum((-1) ** j * self[0][j] * self.Minor(0, j).Det()
for j in range(c))
# Минор матрицы
def Minor(self, i, j):
mas = []
c = self.ColCount
r = self.RowCount
for k in range (r):
tmp = []
if k != i :
for l in range(c):
if l != j:
tmp.append(self[k][l])
mas.append(tmp)
return Matrix(*mas)
# Транспонирование матрицы
def T(self):
mas = []
for i in range (self.ColCount):
tmp = []
for j in range (self.RowCount):
tmp.append(self[j][i])
mas.append(tmp)
return Matrix(*mas)
# Переcтановка строк
def SwapRows(self,i,j):
tmp = []
for k in range (self.ColCount):
tmp.append( self[i][k] )
self[i][k] = self[j][k]
self[j][k] = tmp[k]
return self
# Св-во для чтения кол-ва строк
RowCount = property(GetRowCount)
# Св-во для чтения кол-ва столбцов