使用SymPy求解含微分边界条件的常微分方程问题

问题分析 #

SymPy的dsolve函数的ics参数仅支持单独指定函数在某点的值或其导数在某点的值，无法直接处理导数与函数的线性组合式边界条件（比如u’(0) - 100000*u(0) = 0）。此时需要手动求解通解后，代入边界条件解出积分常数。

实现步骤与代码 #

1. 定义符号与微分方程
   按需求定义所有符号和微分方程：

   import sympy as sp
   
   # 定义符号变量
   x, EA, l, q0 = sp.symbols('x EA l q0')
   u = sp.Function('u')
   d2u = sp.Derivative(u(x), x, 2)
   
   # 构建微分方程
   eqdiff = sp.Eq(q0 + EA * d2u, 0)
   

2. 求解通解
   先不带初始条件求解微分方程，得到含积分常数的通解：

   # 求通解
   general_sol = sp.dsolve(eqdiff, u(x))
   # 提取u(x)的表达式，方便后续处理
   u_x = general_sol.rhs
   # 计算u(x)的一阶导数u’(x)
   du_x = sp.diff(u_x, x)
   

3. 代入边界条件构建方程组
   把两个边界条件转化为关于积分常数的方程：

   • 第一个边界条件：u(l) = 0
   • 第二个边界条件：u’(0) - 100000*u(0) = 0

   # 提取通解中的积分常数（SymPy默认用C1、C2表示）
   C1, C2 = sp.symbols('C1 C2')
   
   # 代入x=l到u(x)，得到第一个方程
   eq1 = sp.Eq(u_x.subs(x, l), 0)
   # 代入x=0到u’(x)和u(x)，得到第二个方程
   eq2 = sp.Eq(du_x.subs(x, 0) - 100000 * u_x.subs(x, 0), 0)
   

4. 解方程组得到积分常数
   用sp.solve求解关于C1、C2的线性方程组：

   # 解方程组
   constants = sp.solve((eq1, eq2), (C1, C2))
   

5. 代入常数得到特解
   把求解得到的C1、C2代回通解，得到满足边界条件的特解：

   # 代入常数得到特解
   particular_sol = u_x.subs(constants)
   # 打印结果
   sp.pprint(particular_sol)
   

说明 #

• 上述代码会输出包含符号参数EA、l、q0的解析解，你可以根据需要代入具体数值计算。
• 如果SymPy自动生成的积分常数不是C1、C2，可以通过general_sol.free_symbols查看通解中的自由符号，再提取对应的常数变量。

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