У меня есть следующая система дифференциальных уравнений:
И согласно бумаге, они сказали, что я могу решить это численно, используя РК 4-го порядка.
Как видите, два последних уравнения связаны, и я построил матрицу (Вероятность), которая связывает xn и yn, где n = 1 .. (N- количество пар, например, здесь N равно 4): vector ([ x1, x2, ..., xn, y1, y2, ..., yn]) '= Probability.dot (vector ([x1, x2, ..., xn, y1, y2, ..., yn] )), где штрих - дифференцирование по времени. Но, с другой стороны, для каждого шага у меня есть дополнительный член суммы (un * xn и то же самое для yn), это была первая проблема, с которой я столкнулся и не знал, как с ней справиться.
Я написал код и возникло много ошибок, с которыми я не справился.
Пока я пытаюсь справиться с этим самостоятельно, буду очень благодарен за любую помощь в этом.
Выше показан мой код:
Импорт библиотек
import numpy as np
import math
import scipy.constants as sc
from scipy.sparse import diags
from scipy.integrate import ode
import matplotlib.pyplot as plt
from matplotlib import mlab
Исходные данные и константы
dimen_paramT0 = [0,0,0,0]
step = 0.00001
Mn = 1e-21 # mass of pairs
thau = 1e-15 # character time
wn_sqr = 1e-2 # to 10e-6
wn_prime = 3e-2 # to 10e-5
n = 4 #count of repetition; forinstance 4
gamma_n = 3e-9
Dn = 4e-2
an = 4.45
alphaPrime_n = 0.13
Volt = 0.4
hn = 0
hn_nInc = 1.276 #hn,n+1
hn_nDec = 1.276 #hn,n-1
Un = thau * alphaPrime_n / Mn
ksi_n = an * Un
Omega_n = 2 * Dn * an * (thau ** 2) / (Mn * Un)
* Построение матрицы для ассоциаций с diff / eq / для вероятности *
k = np.array([hn_nInc*np.ones(n-1),hn*np.ones(n),hn_nDec*np.ones(n-1)])
offset = [-1,0,1]
Probability = diags(k,offset).toarray() # bn(tk)=xn(tk)+iyn(tk)
xt0_list = [0] * n
yt0_list = [0] * n
* Правая сторона разн. экв. *
# dimen_param = [un,vn,zn,vzn] [tn]
# x_list = [x1,...,xn] [tn]
# y_list = [y1,...,yn] [tn]
def fun(dimen_param, x_list, y_list):
return dimen_param[1]
def fvn(dimen_param, x_list, y_list):
return -(x_list[len(x_list)-1]**2 + y_list[len(y_list)-1]**2)- wn_prime*dimen_param[1] + Omega_n * (1-np.e ** (-ksi_n * dimen_param[0]))*np.e ** (-ksi_n * dimen_param[0])
def fzn(dimen_param, x_list, y_list):
return dimen_param[3]
def fvzn(dimen_param, x_list, y_list):
return -wn_prime * dimen_param[3]-(wn_sqr ** 2) * dimen_param[2] - 1
def fxn(dimen_param, x_list, y_list):
return Probability.dot(y_list)
def fyn(dimen_param, x_list, y_list):
return -Probability.dot(x_list)
#xv = [dimen_param, x_list, y_list]
def f(xv):
k_d = xv[0:4]
k_x = xv[4:4+len(xt0_list)]
k_y = xv[4+len(xt0_list):4+len(xt0_list)+len(yt0_list)]
return ([fun(k_d, k_x, k_y),fvn(k_d, k_x, k_y),fzn(k_d, k_x, k_y),fvzn(k_d, k_x, k_y),fxn(k_d, k_x, k_y),fyn(k_d, k_x, k_y)])
* Реализация метода Рунге - Кутта 4-го порядка *
def RK4(f, dimen_paramT0, xt0_list, yt0_list):
T = np.linspace(0, 1. / step, 1. / step +1)
xvinit = np.concatenate([dimen_paramT0, xt0_list, yt0_list])
xv = np.zeros( (len(T), len(xvinit)) )
xv[0] = xvinit
for i in range(int(1. / step)):
k1 = f(xv[i])
k2 = f(xv[i] + step/2.0*k1)
k3 = f(xv[i] + step/2.0*k2)
k4 = f(xv[i] + step*k3)
xv[i+1] = xv[i] + step/6.0 *( k1 + 2*k2 + 2*k3 + k4)
return T, xv
* Бег *
print RK4(f, dimen_paramT0, xt0_list, yt0_list)
Проблема на данный момент:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-104-be8ed2e50d37> in <module>()
----> 1 RK4(f, dimen_paramT0, xt0_list, yt0_list)
<ipython-input-103-8c48cf5efe73> in RK4(f, dimen_paramT0, xt0_list, yt0_list)
7 for i in range(int(1. / step)):
8 k1 = f(xv[i])
----> 9 k2 = f(xv[i] + step/2.0*k1)
10 k3 = f(xv[i] + step/2.0*k2)
11 k4 = f(xv[i] + step*k3)
TypeError: can't multiply sequence by non-int of type 'float'
k1
в указанной строке должен быть массивом numpy? Обычные списки и кортежи Python не поддерживают поэлементные математические операции. - person Aaron   schedule 04.04.2018numba
и своевременной компиляции. - person Aaron   schedule 04.04.2018