使用 Python 求解向量微分方程

如何解决使用 Python 求解向量微分方程

我正在使用 python 和 scipy 来求解具有标量第二个成员的微分方程，当我将其更改为向量时，我无法正确解决。

import numpy as np
import scipy.integrate as si
import matplotlib.pylab as plt

t = np.linspace(0,1,500)
a = 2667613637.976072
phi = [-1,1]
w = [10*np.pi,5*np.pi,10*np.pi]
beta = [1.25e-8,1.12e-8,1.18e-8]

def main():
    plt.plot(t,solver())

def get_g():
    idx = len(w)
    g = 0
    for i in range(0,idx):
        g += w[i]*beta[i]*np.cos(w[i]*(t+phi[i]))    
    return g
    
def diff(f,t):
    # a is a scalar,g is a sum of cos functions
    # the desired equation is dI/dt = a*I(t) + g(t)
    g = get_g()
    return ((f * a)  + g)

def solver():
    sol = si.odeint(diff,t)
    return sol[:,0]

if __name__ == "__main__":
    main()

这是我要解的方程：

解决方法

您要用于评估函数 t 的 get_g 是已传递给 t 函数的 diff，而不是您定义的时间样本数组一开始。实际上不可能（但并非不可能）将时间样本数组的任何值（第一个除外）传递给 diff 函数。

所以只需添加参数

def get_g(t):
    return sum( w[i]*beta[i]*np.cos(w[i]*(t+phi[i])) for i in range(len(w)) )  

然后也在 diff 中使用 g=get_g(t)。

---

对标量参数使用 math 函数以避免不必要的开销。或者将所有数组都设为 numpy，以便在 numpy 数组操作中用隐式 numpy 循环替换显式 python 循环。

phi = np.array([-1,1])
w = np.array([10*np.pi,5*np.pi,10*np.pi])
beta = np.array([1.25e-8,1.12e-8,1.18e-8])

def get_g(t):
    return sum( w*beta*np.cos(w*(t+phi)) )  

def diff(f,t):
    # a is a scalar,g is a sum of cos functions
    # the desired equation is dI/dt = a*I(t) + g(t)
    g = get_g(t)
    return ((f * a)  + g)

def solver():
    sol = si.odeint(diff,t,atol = 1e-18)
    return sol[:,0]

通过这些更改，您应该会得到一个解，但是在到达 t=1e-6 之前它会发散到无穷大。这很可能意味着模型需要做一些工作。