如何将Sympy与ctypes外部函数结合实现数值计算？

问题核心原因 #

直接将Sympy符号对象传给ctypes函数时，ctypes会尝试把符号转成float类型，而Sympy的Symbol是符号对象而非数值，因此抛出TypeError: Cannot convert expression to float。我们需要让Sympy在数值评估阶段才将符号替换为数值，再调用外部函数，不需要Sympy对外部函数做符号推导。

解决方案步骤 #

1. 先确保ctypes外部函数能正常处理数值参数

以你提供的C++ DLL为例，先完成基础的ctypes调用封装：

C++ DLL代码（simple_math.cpp）

#include <windows.h>

// 导出：返回两数之和
extern "C" __declspec(dllexport) float add(float a, float b) {
    return a + b;
}

// 导出：通过引用返回两数之和（模拟ctRefprop的接口风格）
extern "C" __declspec(dllexport) void add_ref(float a, float b, float* result) {
    *result = a + b;
}

BOOL APIENTRY DllMain(HMODULE hModule, DWORD ul_reason_for_call, LPVOID lpReserved) {
    switch (ul_reason_for_call) {
        case DLL_PROCESS_ATTACH:
        case DLL_THREAD_ATTACH:
        case DLL_THREAD_DETACH:
        case DLL_PROCESS_DETACH:
            break;
    }
    return TRUE;
}

Python端ctypes基础封装

import ctypes
import sympy as sp
from scipy.optimize import root

# 加载DLL
dll = ctypes.CDLL("./simple_math.dll")

# 定义函数参数/返回类型，避免内存错误
dll.add.argtypes = [ctypes.c_float, ctypes.c_float]
dll.add.restype = ctypes.c_float

dll.add_ref.argtypes = [ctypes.c_float, ctypes.c_float, ctypes.POINTER(ctypes.c_float)]
dll.add_ref.restype = None

2. 封装外部函数为Sympy可识别的数值函数

有两种常用方式，根据需求选择：

方式一：自定义Sympy Function类（保留符号表达式能力）

适合需要构建符号方程组，再用Sympy数值求解器求解的场景：

# 包装add函数为Sympy可处理的函数
class AddExternal(sp.Function):
    @classmethod
    def eval(cls, a, b):
        # 若参数是数值，直接调用外部函数返回结果
        if isinstance(a, (int, float)) and isinstance(b, (int, float)):
            return dll.add(float(a), float(b))
        # 若参数是符号，返回符号表达式
        return cls(a, b)

# 包装add_ref函数（处理引用返回）
class AddRefExternal(sp.Function):
    @classmethod
    def eval(cls, a, b):
        if isinstance(a, (int, float)) and isinstance(b, (int, float)):
            result = ctypes.c_float()
            dll.add_ref(float(a), float(b), ctypes.byref(result))
            return result.value
        return cls(a, b)

# 构建方程组并数值求解
x, y = sp.symbols('x y')
eq1 = sp.Eq(x + y, 5)
eq2 = sp.Eq(AddExternal(x, y), 6)

# 使用nsolve数值求解，需传入初始猜测值
solution = sp.nsolve((eq1, eq2), (x, y), (1, 1))
print(f"解：x = {solution[0]}, y = {solution[1]}")

方式二：直接写数值残差函数（更高效，适合纯数值求解）

若不需要保留符号形式，直接将方程组转为数值残差函数，结合scipy求解器：

def residuals(vars):
    x_val, y_val = vars
    # 方程1残差：x + y - 5
    res1 = x_val + y_val - 5
    # 方程2残差：外部add函数结果 - 6
    res2 = dll.add(x_val, y_val) - 6
    return [res1, res2]

# 调用scipy求解器
result = root(residuals, [1, 1])
print(f"解：x = {result.x[0]:.4f}, y = {result.x[1]:.4f}")

针对ctRefprop的适配要点 #

对于ctRefprop的ctypes接口，只需照搬上述封装逻辑：

1. 先确保ctRefprop的ctypes函数能正确接收数值参数、返回/输出数值结果；
2. 将每个需要调用的ctRefprop函数，封装为类似AddExternal的Sympy Function类，或直接写入残差函数中；
3. 用Sympy的nsolve或scipy的数值求解器完成方程组求解。

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