使用Sympy求解PDE方程组时约束方程未代入的问题咨询

原因 #

你遇到的问题是SymPy的solve函数在处理多方程求解时，不会自动推导隐式约束：eq1只给出了diff(u,x)=0，但solve不会自动推断出diff(u,x,x)=0（一阶导为0时二阶导必然为0）；同时在求解指定变量diff(v,x)时，它只会从包含该变量的eq2中提取表达式，未将eq1的约束完全应用到eq2的所有项上。另外你的原代码漏了定义符号x，实际运行会报错，需要补充定义。

解决方法 #

方法1：手动代入约束

先利用eq1的结论替换eq2中的一阶导，再推导并替换二阶导：

from sympy import *

x = symbols('x')
u = Function('u', Real=True)(x)
v = Function('v', Real=True)(x)
eq1 = Eq(diff(u, x), 0)
eq2 = Eq(diff(v, x), diff(u, x, x) + diff(u, x))

# 代入一阶导的约束
eq2_sub = eq2.subs(diff(u, x), 0)
# 由eq1推导二阶导为0，代入eq2
eq2_final = eq2_sub.subs(diff(u, x, x), 0)
# 求解diff(v,x)
print(solve(eq2_final, diff(v, x)))

输出：[0]

方法2：显式添加二阶导约束

把一阶导的约束和推导出来的二阶导约束联立，再求解：

from sympy import *

x = symbols('x')
u = Function('u', Real=True)(x)
v = Function('v', Real=True)(x)
eq1 = Eq(diff(u, x), 0)
# 显式添加二阶导的约束：对eq1两边求导得到
eq1_second = Eq(diff(u, x, x), diff(eq1.rhs, x))
eq2 = Eq(diff(v, x), diff(u, x, x) + diff(u, x))

# 联立所有约束求解
print(solve([eq1, eq1_second, eq2], diff(v, x)))

输出：[0]

方法3：使用线性方程组求解器linsolve

因为这是线性关系，用linsolve能更准确地处理约束联立：

from sympy import *

x = symbols('x')
u = Function('u', Real=True)(x)
v = Function('v', Real=True)(x)
eq1 = Eq(diff(u, x), 0)
eq2 = Eq(diff(v, x) - diff(u, x) - diff(u, x, x), 0)

# 求解包含diff(u,x), diff(u,x,x), diff(v,x)的线性方程组
sol = linsolve([eq1, eq2], diff(u, x), diff(u, x, x), diff(v, x))
print(sol)

输出：{(0, 0, 0)}，其中第三个元素即为diff(v,x)=0

内容的提问来源于stack exchange，提问作者miro