I have the following system of differential equations:
And according to the paper they told I can solve it numerically by using RK 4th order.
As you can see, two last equations are coupled and I constructed a matrix (Probability) that associates xn and yn, where n = 1..(N- count of pairs, for instance, here N is equal 4) : vector([x1,x2,...,xn,y1,y2,...,yn])' = Probability.dot(vector([x1,x2,...,xn,y1,y2,...,yn])), where prime is time differentiation. But on the ither hand for each step I have additional term of sum (un * xn and the same for yn), it was the first problem I faced and didn't know how to deal with it.
I wrote a code and a lot of errors which I couldn't manage arose.
While I am trying to cope it by myself, I'll be very thankful for any help in it.
Above show you my code:
Importing Libraries
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
Initial data and constants
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)
*** Constructing matrix for associations with dif/ eq/ for probability ***
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
*** Right side of dif. eq. ***
# 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)])
*** Realisation of Runge - Kutta method of 4th order ***
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
*** Running ***
print RK4(f, dimen_paramT0, xt0_list, yt0_list)
The problem at this moment is:
---------------------------------------------------------------------------
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'
