python - 二微分方程组的 Runge-Kutta 实现

Question

我有以下微分方程组:

根据他们告诉我的论文，我可以使用 RK 四阶数值求解它。

如您所见，最后两个方程是耦合的，我构造了一个关联 xn 和 yn 的矩阵(概率)，其中 n = 1..(例如，N- 对数，此处 N 等于 4):矢量([x1,x2,...,xn,y1,y2,...,yn])' = Probability.dot(矢量([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

*** diff 的右侧。等式***

# 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)])

***四阶龙格-库塔法的实现***

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'

Answer 1

k1 是一个 python list，其中 multip 表示“重复”，所以

[1,2] * 3 == [1,2,1,2,1,2]

显然这对于​​ float 没有意义

[1,2,3]*2.0
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-118-7d3451361a53> in <module>()
----> 1 [1,2,3]*2.0

TypeError: can't multiply sequence by non-int of type 'float'

您想要 numpy.ndarray 的向量行为，其中

np.array([1,2])*2.0 == np.array([2.0, 4.0])

因此请确保f的返回语句是:

return np.asarray([
    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)
])

<小时/>

根据记录，scipy 有一个 RK4 solver已经不需要自己实现了。