计算科学 - 求解具有多个自变量的耦合微分方程 - 吾爱随笔录

我目前正在使用scipy.integrate.odeint 的包装器odeintw 来求解我的方程，因为它们是复值的。

目前，我有 3 个具有 2 个自变量的耦合一阶微分方程

\begin{aligned} \frac{\partial C}{\partial t} & = \frac{i}{a} (A B - c d) - \frac{C}{t_{1}} \\ \frac{\partial A}{\partial t} & = \frac{2 i b^{2}}{a} c C - \frac{A}{t_{2}} \\ \frac{\partial B}{\partial t} & = i e d + f \end{aligned}

$\begin{align} \frac{\partial C}{\partial t} &= \frac{i}{a}(AB - cd) - \frac{C}{t_1}\\ \frac{\partial A}{\partial t} &= \frac{2ib^2}{a}cC - \frac{A}{t_2}\\ \frac{\partial B}{\partial t} &= ied + f \end{align}$

其中大写字母 A、B 和 C 是我的函数；小写字母 a、b、c、t1、t2 和 d 只是参数；z 和 t 是我的自变量。我已经能够通过简单地使 B 等于一个常数（即忽略第三个方程）并输入某些初始条件/参数值来解决前两个方程。代码如下所示，分别提供解的实部和虚部（注意：有些变量名称不同，但方程是等价的）：

from odeintw import odeintw
import numpy as np
import matplotlib.pyplot as plt

def Wfunc(W, t, hbar, Pm, Em, Ep, T1, T2, d):
    N, Pp = W
    return [(1j/hbar)*(Pp*Ep - Em*Pm) - N/T1,
            (2j*(d**2)/hbar)*Em*N - Pp/T2]

W0 = np.array([1+2j, 3+4j])
t = np.linspace(0, 5, 1001)
hbar = 1.
Pm = 4 - 2j
Em = 2.5
Ep = 10.
T1 = 100.
T2 = 10
d = 0.1
W, infodict = odeintw(Wfunc, W0, t, args=(hbar, Pm, Em, Ep, T1, T2, d),
                  full_output=True)

plt.figure(1)
plt.clf()
color1 = (0.5, 0.4, 0.3)
color2 = (0.2, 0.2, 1.0)
plt.plot(t, W[:, 0].real, color=color1, label='N.real', linewidth=1.5)
plt.plot(t, W[:, 0].imag, '--', color=color1, label='N.imag', linewidth=2)
plt.plot(t, W[:, 1].real, color=color2, label='Pp.real', linewidth=1.5)
plt.plot(t, W[:, 1].imag, '--', color=color2, label='Pp.imag', linewidth=2)
plt.xlabel('t')
plt.grid(True)
plt.legend(loc='best')

plt.show()

这适用于具有相同自变量的 2 个方程，但现在我想介绍我的另一个自变量 z，并且遇到了一些困难。是否有任何有效的方法可以同时解决这 3 个方程？扩展涉及自变量和更复杂微分方程的大区间的情况的代码将是有意义的。