如何将给定Mathematica微分方程组代码转换为Python代码

原Mathematica代码 #

NDSolve[{1/(1 - f)^2.5*g''''[y] + 
       G*Sin[\[Alpha]]*((1 - f) + f*(p1*e1)/(p2*e2))*
        t'[y] - (((M^2)*A)/(1 + (A*m)^2) + 1/(k*(1 - f)^2.5))*
        g''[y] == 0, 
     b1*t''[y] + Br/(1 - f)^2.5*(g''[y])^2 + 
       Br*(((M^2)*A)/(1 + (A*m)^2) + 1/(k*(1 - f)^2.5)) (g'[y] + 
          1)^2 + \[Epsilon] == 0, g[h1] == F/2, g'[h1] == -1, 
     g[h2] == -F/2, g'[h2] == -1, t[h1] == -1/2, t[h2] == 1/2}, {g, 
     t}, {y, h1, h2}]

转换思路与Python实现 #

这是一个耦合的四阶+二阶边值问题（BVP），Python中常用scipy.integrate.solve_bvp求解这类问题，核心步骤是将高阶微分方程降阶为一阶方程组，再定义残差函数与边界条件。

1. 前置准备 #

首先给所有符号参数赋值（Python求解数值解需要具体数值，必须与Mathematica中使用的参数完全一致），示例参数如下：

import numpy as np
from scipy.integrate import solve_bvp

# 替换为你Mathematica中使用的实际数值
G = 1.0
alpha = np.pi/6
sin_alpha = np.sin(alpha)
p1, e1, p2, e2 = 2.0, 0.5, 1.0, 0.3
M = 0.8
A = 1.2
m = 0.5
k = 0.7
b1 = 0.4
Br = 0.9
epsilon = 0.1
F = 1.0
h1 = -1.0
h2 = 1.0
f = 0.2

2. 降阶为一阶方程组 #

将高阶方程拆解为一阶方程组，是数值求解的必要步骤：

• 对于四阶函数g(y)：
  令 g0 = g(y)、g1 = g'(y)、g2 = g''(y)、g3 = g'''(y)，则原四阶方程转化为：

  g0' = g1
  g1' = g2
  g2' = g3
  g3' = [ ((M²A)/(1+(Am)²) + 1/(k(1-f)^2.5))*g2 - G*sin_alpha*((1-f)+f*(p1e1)/(p2e2))*t1 ] * (1-f)^2.5
  

• 对于二阶函数t(y)：
  令 t0 = t(y)、t1 = t'(y)，则原二阶方程转化为：

  t0' = t1
  t1' = [ -Br/(1-f)^2.5 * g2² - Br*((M²A)/(1+(Am)²) + 1/(k(1-f)^2.5))*(g1+1)² - epsilon ] / b1
  

3. 定义残差函数与边界条件 #

# 定义微分方程组残差
def fun(y, z):
    # z = [g0, g1, g2, g3, t0, t1]
    g0, g1, g2, g3, t0, t1 = z
    
    # 计算g3'（即g''''）
    term_g = ((M**2 * A)/(1 + (A*m)**2) + 1/(k*(1 - f)**2.5)) * g2
    term_t = G * sin_alpha * ((1 - f) + f*(p1*e1)/(p2*e2)) * t1
    g3_prime = (term_g - term_t) * (1 - f)**2.5
    
    # 计算t1'（即t''）
    term1 = Br/(1 - f)**2.5 * (g2)**2
    term2 = Br * ((M**2 * A)/(1 + (A*m)**2) + 1/(k*(1 - f)**2.5)) * (g1 + 1)**2
    t1_prime = (-term1 - term2 - epsilon) / b1
    
    return np.vstack([g1, g2, g3, g3_prime, t1, t1_prime])

# 定义边界条件残差
def bc(za, zb):
    # za是y=h1处的z值，zb是y=h2处的z值
    g0_a, g1_a, _, _, t0_a, _ = za
    g0_b, g1_b, _, _, t0_b, _ = zb
    
    return np.array([
        g0_a - F/2,    # g(h1) = F/2
        g1_a + 1,      # g'(h1) = -1
        g0_b + F/2,    # g(h2) = -F/2
        g1_b + 1,      # g'(h2) = -1
        t0_a + 0.5,    # t(h1) = -1/2
        t0_b - 0.5     # t(h2) = 1/2
    ])

4. 求解并验证 #

# 生成y的采样点
y = np.linspace(h1, h2, 100)

# 初始猜测（需贴近真实解，否则可能不收敛）
z0 = np.zeros((6, y.size))
# 给g0设线性分布匹配边界条件
z0[0] = F/2 - F*(y - h1)/(h2 - h1)
# g1初始值设为边界条件的-1
z0[1] = -1
# t0设线性分布匹配边界条件
z0[4] = -0.5 + (y - h1)/(h2 - h1)

# 求解BVP
sol = solve_bvp(fun, bc, y, z0, tol=1e-6)

# 检查求解状态
if sol.success:
    print("求解成功")
    # 获取结果
    g_sol = sol.sol(y)[0]
    g_prime_sol = sol.sol(y)[1]
    t_sol = sol.sol(y)[4]
else:
    print(f"求解失败：{sol.message}")

输出不一致的常见原因及解决 #

• 参数不匹配：确保Python中所有参数的数值、精度与Mathematica完全一致，比如用2.5还是5/2。
• 初始猜测不合理：solve_bvp对初始猜测敏感，若偏离真实解过远，可能收敛到不同解或不收敛，可参考Mathematica的输出调整初始猜测。
• 求解器精度差异：调整solve_bvp的tol参数（比如设为1e-8），或增加max_nodes提升计算精度，匹配Mathematica的求解精度。
• 降阶错误：仔细核对降阶后的一阶方程组，确保符号、系数与原方程完全等价。

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