如何利用jacobian矩阵采用隐式算法GMRES迭代求解固体的非稳态heat equation？ 给出详细的三维算例推导和公式，并给出详细的c++代码范例

在数值计算中，尤其是解决偏微分方程（如非稳态热传导方程）时，Jacobian矩阵常用于优化迭代方法如Generalized Minimal Residual Method (GMRES)。对于三维非稳态热传导方程，该方程通常表示为：

\[
\rho c_p \frac{\partial T}{\partial t} - \nabla \cdot (\kappa \nabla T) = q
\]

其中 \( T \) 是温度，\( \rho \) 是密度，\( c_p \) 是比热容，\( \kappa \) 是热导率，\( t \) 是时间，\( q \) 是热源。

隐式时间积分（如显式-隐式Runge-Kutta法）可以避免直接处理高阶的时间导数，使得非线性问题更容易求解。GMRES适用于大型稀疏系统，通过Jacobian矩阵对每个迭代步进行线性化。

Jacobian矩阵 \( J \) 对于这个方程是关于温度变化的导数，即：

\[
J(T) = \frac{\partial}{\partial T}\left( \rho c_p \frac{\partial T}{\partial t} - \nabla \cdot (\kappa \nabla T) \right)
\]

在隐式算法中，我们有：

\[
J(T^n)(\Delta T^{n+1}) = R^n
\]
其中 \( R^n \) 是残差项，\( \Delta T^{n+1} \) 是下一个时间步的温度增量。

在三维空间上，\( \nabla \cdot (\kappa \nabla T) \) 可以分解为三个方向上的梯度乘积，形成一个三阶偏导数组成的Jacobian矩阵。

C++代码示例过于复杂，这里提供一个简化的思路，真正的代码会涉及到离散网格、矩阵运算库（如Eigen）以及线性代数操作：

// 假设已经有一个求解器类Solver，包含Jacobian矩阵构建和求逆函数
class NonSteadyHeatSolver {
public:
    JacobianMatrix buildJacobian(const Vector3D& T);
    Vector3D solveForTemperatureIncrement(const Vector3D& residual, const Vector3D& oldTemperature);
private:
    // ... 具体的矩阵构建和逆矩阵运算细节
};

// 调用部分
Vector3D temperature;
timeStep++;
residual = computeResidual(timeStep, temperature); // 计算残差

JacobianMatrix jacobian = solver.buildJacobian(temperature);
Vector3D deltaT = solver.solveForTemperatureIncrement(residual, temperature); // 使用GMRES迭代求解

temperature += deltaT; // 更新温度