如何用Sympy提取二阶齐次线性ODE的两个基解函数

问题场景 #

我有一个二阶齐次线性常微分方程（ODE），其通解为ur(r) = C1*u1(r) + C2*u2(r)，希望通过设置(C1=1, C2=0)和(C1=0, C2=1)分别获取u1(r)和u2(r)这两个基解函数。

尝试了如下代码，但运行后subs操作未生效，输出仍为C1/r + C2*r：

import IPython.display as disp
import sympy as sym
sym.init_printing(latex_mode='equation')

r = sym.symbols('r', real=True, positive=True)
ur = sym.Function('ur', real=True)
def Naxp(ur):
    return r**2 * sym.diff(ur(r), r, 2) + r * sym.diff(ur(r), r) - ur(r)
Naxp_ur = Naxp(ur)
EqNaxpur = sym.Eq(Naxp_ur)
#print('Radial ODE:')
#display(Naxp_ur)
ur_sol = sym.dsolve(Naxp_ur, ur(r))
ur_sol1 = ur_sol.rhs
#display(sym.Eq(ur(r), ur_sol1))

# I still have to find how to extract the two basis functions
#u1 = sym.Function('u1', real=True)
C1, C2 = sym.symbols('C1, C2', real=True)
u1 = ur_sol1.subs([(C1, 1), (C2, 0)])
print(u1)

问题原因 #

核心问题是符号定义顺序错误：先调用dsolve得到通解，之后才手动定义C1、C2符号。此时通解中的C1、C2是SymPy自动生成的内部符号，和后续手动定义的C1、C2并非同一个对象，导致subs无法匹配替换。

解决方法 #

方法1：提前定义符号 #

在调用dsolve前先定义C1、C2，让SymPy求解时直接使用这些预定义符号，subs即可正常工作：

import sympy as sym
sym.init_printing(latex_mode='equation')

# 提前定义所有符号，包括C1、C2
r = sym.symbols('r', real=True, positive=True)
C1, C2 = sym.symbols('C1, C2', real=True)
ur = sym.Function('ur', real=True)

def Naxp(ur):
    return r**2 * sym.diff(ur(r), r, 2) + r * sym.diff(ur(r), r) - ur(r)

Naxp_ur = Naxp(ur)
ur_sol = sym.dsolve(Naxp_ur, ur(r))
ur_sol1 = ur_sol.rhs

# 替换符号获取基解
u1 = ur_sol1.subs([(C1, 1), (C2, 0)])
u2 = ur_sol1.subs([(C1, 0), (C2, 1)])

print("u1(r):", u1)
print("u2(r):", u2)

方法2：直接提取基函数 #

无需手动定义符号，直接从通解的线性组合项中分离出基函数：

import sympy as sym
sym.init_printing(latex_mode='equation')

r = sym.symbols('r', real=True, positive=True)
ur = sym.Function('ur', real=True)

def Naxp(ur):
    return r**2 * sym.diff(ur(r), r, 2) + r * sym.diff(ur(r), r) - ur(r)

Naxp_ur = Naxp(ur)
ur_sol = sym.dsolve(Naxp_ur, ur(r))
ur_sol1 = ur_sol.rhs

# 提取通解中的线性项，分离基函数
u1 = ur_sol1.args[0].coeff(ur_sol1.args[0].free_symbols.pop())
u2 = ur_sol1.args[1].coeff(ur_sol1.args[1].free_symbols.pop())

print("u1(r):", u1)
print("u2(r):", u2)

方法3：直接构造基函数（针对本例） #

由于本例通解形式简单，可直接写出基函数：

u1 = sym.Rational(1)/r
u2 = r

内容的提问来源于stack exchange，提问作者sancho.s ReinstateMonicaCellio