使用PETSc4py求解ODE遇异常：输出全0及报错求助

你的代码存在几个关键错误，导致求解结果异常，以下是问题分析和修复方案：

核心问题分析 #

• RHS函数签名不匹配：PETSc4py的TS求解器要求右端函数（RHS）必须遵循f(t, y, ydot)的签名，你额外添加的k参数会导致求解器无法正确调用函数，进而无法计算导数项。
• 向量赋值错误：直接执行ydot = -k * y只是替换了Python变量的引用，并没有修改PETSc向量ydot的实际内存数据，求解器无法获取正确的导数信息。
• 问题类型设置错误：y'=-ky是线性ODE，应设置为TS.ProblemType.LINEAR，而非NONLINEAR，错误的类型会触发不合适的求解逻辑。
• 错误信息被屏蔽：PETSc.Sys.popErrorHandler()会关闭PETSc的错误提示，导致你无法看到底层报错信息，增加排查难度。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
import petsc4py
from petsc4py import PETSc
from functools import partial

# 初始化PETSc
petsc4py.init()

# ODE右端函数：修正签名，正确修改向量内容
def f(t, y, ydot, k):
    ydot[:] = -k * y  # 用切片赋值修改PETSc向量的实际数据

k = 0.5
# 用partial绑定参数k，让函数符合TS要求的签名
f_bound = partial(f, k=k)

# 创建时间步进器
ts = PETSc.TS().create()
ts.setType('bdf')
# 设置为线性问题类型
ts.setProblemType(PETSc.TS.ProblemType.LINEAR)
# 绑定处理后的RHS函数
ts.setRHSFunction(f_bound)

# 设置时间参数
t_max = 100
num_steps = 30
ts.setMaxTime(t_max)
ts.setTimeStep(t_max / num_steps)
ts.setFromOptions()

# 设置初始条件
y0 = np.array([1.0])
y = PETSc.Vec().createWithArray(y0)
ts.setTime(0.0)
ts.setSolution(y)

# 求解并记录结果
solution_times = []
solution_values = []

while ts.getTime() < ts.getMaxTime():
    ts.step()
    solution_times.append(ts.getTime())
    solution_values.append(y[0])

# 绘制结果
plt.plot(solution_times, solution_values, label='PETSc4py Solution')
# 绘制解析解作为对比
t_arr = np.array(solution_times)
plt.plot(t_arr, np.exp(-k*t_arr), '--', label='Analytical Solution')
plt.xlabel('Time')
plt.ylabel('y(t)')
plt.legend()
plt.show()

额外说明 #

• 使用functools.partial是PETSc4py中传递额外参数到RHS函数的标准方式，确保函数签名符合求解器要求。
• 移除PETSc.Sys.popErrorHandler()后，如果后续遇到问题，PETSc会输出详细的错误日志，帮助快速定位问题。
• 加入了解析解的对比曲线，可以直观验证求解结果的正确性。

内容的提问来源于stack exchange，提问作者A.Be.