使用隐式有限差分法（FDM）求解HJB方程时结果不符的问题排查求助

我正在尝试用隐式有限差分法（FDM）数值求解一个特定的Hamilton-Jacobi-Bellman（HJB）方程，目标是验证已知的解析解，但得到的结果完全不符合预期，实在找不到问题出在哪，希望能得到大家的帮助。

问题背景与方程设定 #

我考虑的是定义在$[0,T] \times \mathbb{R}$上的价值函数$v(t,x)$和增益函数$f(t,x)$，暂时忽略最优控制，先让代码跑通。对应的HJB方程为：
$$-\frac{\partial}{\partial t}v(t,x)-(b_t-x)\frac{\partial}{\partial x}v(t,x)-\frac{\sigma^{2}{2}\frac{\partial}2}{\partial x^2}v(t,x)-f(t,x)=0$$
终端条件为$v(T,x)=0$，其中状态过程是均值回复过程：
$$dX_t= (b_t-X_t)dt + \sigma dW_t, \qquad W_t\sim\mathcal{N}(0,1)$$
这里$b_t$是随时间变化的预测值。

为了验证数值方法的一致性，我选取了已知的解析解$v(t,x)=tx^2$，代入HJB方程后可以得到对应的$f(t,x)$：
$$f(t,x)=-x^{2-(b_{t}-x)2tx-t\sigma}2$$

隐式有限差分的离散化 #

因为是反向时间求解（从终端$T$往初始时刻$0$推进），所以采用隐式格式保证稳定性。对空间导数的离散：

• 一阶导数$\frac{\partial}{\partial x}v$用向后差分
• 二阶导数$\frac{\partial^2}{\partial x^2}v$用中心差分

离散后的方程推导如下：

时间离散：$t_n = n\Delta t$，空间离散：$x_i = x_0 + ih$，其中$\Delta t$是时间步长，$h$是空间步长。

将导数替换为差分后，原方程转化为：
$$-\frac{v_i^{t_n}-v_{i}{t_{n-1}}}{\Delta t}=\frac{\sigma^{2}{2}\frac{v_{i+1}}{t_{n-1}}-2v_{i}^{{t_{n-1}}+v_{i-1}}{t_{n-1}}}{h^{2}+(b_{t_{n-1}}-x_i)\frac{v_{i}}{t_{n-1}}-v_{i-1}^{t_{n-1}}}{h}+f(t_{n-1}, x_i)$$

整理后得到线性方程组的形式：
$$v_i^{t_n} = -\Delta t \big((\frac{\sigma^2}{2h2}-\frac{b_{t_{n-1}}-x_i}{h})v_{i-1}^{{t_{n-1}}-(\frac{\sigma}2}{h^2}+\frac{1}{\Delta t}-\frac{b_{t_{n-1}}-x_i}{h})v_{i}^{{t_{n-1}}+\frac{\sigma}2}{2h^2}v_{i+1}{t_{n-1}}+f(t_{n-1}, x_i)\big)$$

边界条件：因为解析解$v(t,x)=tx^2$已知，所以空间边界$x=a$和$x=b$处的$v(t,a)$、$v(t,b)$可以直接计算得到。

伪代码逻辑 #

我的求解流程大致是：

• 预计算并存储所有时刻的边界真实值到v_a和v_b数组
• 从终端时刻$T$开始反向迭代到$t_1$（避免越界）：
  1. 初始化价值函数矩阵
  2. 构建三对角线性方程组的系数矩阵$A$
  3. 处理边界项，组装右端向量
  4. 求解线性方程组得到当前时间步的内部节点值
  5. 将边界值和内部值拼接成完整的价值函数数组

代码实现 #

以下是我目前的Python代码：

import numpy as np
import scipy.stats as st
import scipy.sparse as spa
import matplotlib.pyplot as plt
import plotly.graph_objs as go

T = 3
d_t = 0.25
t_grid = np.arange(0, T+d_t, d_t)

prodmax = 6
consmax = 6
h = 0.2
x_grid = np.arange(-prodmax, consmax+h, h)

b = np.zeros(len(t_grid))
sig = 0.1
for i, t in enumerate(t_grid):
    if i==0:
        b[0] = 0
    else:
        b[i] = b[i-1]+(6*np.sin(np.pi*t_grid[i-1])-b[i-1])*d_t

# 计算真实价值函数
v_true = np.zeros((len(t_grid), len(x_grid)))
for l, t in enumerate(t_grid):
    for i, x in enumerate(x_grid):
        v_true[l,i] = t*x**2

va_store = v_true[:, 0]
vb_store = v_true[:, -1]
va_store[-1], vb_store[-1] = 0, 0  # 终端条件

v_fdm = np.zeros((len(t_grid), len(x_grid[1:-1])))

for l, t in enumerate(t_grid[1:]):
    # 反向时间迭代的索引
    l_back = len(t_grid)-1-l
    
    # 构建三对角矩阵的三个对角线
    up_diag = sig**2/h**2/2*np.ones(len(x_grid[1:-2]))
    diag = -sig**2/h**2 - 1/d_t + (b[l_back-1]-x_grid[1:-1])/h
    low_diag = sig**2/h**2/2+(b[l_back-1]-x_grid[2:-1])/h
    
    # 计算增益函数f，并处理边界项
    f = np.array([-x**2-(b[l_back-1]-x)*2*t*x-t*sig**2 for x in x_grid[1:-1]])
    f[0] = f[0] + (sig**2/h**2/2-(b[l_back-1]-x_grid[1])/h)*va_store[l_back-1]
    f[-1] = f[-1] + sig**2/h**2/2*vb_store[l_back-1]
    
    # 组装系数矩阵A
    A = np.diag(up_diag, k=1) + np.diag(diag, k=0) + np.diag(low_diag, k=-1)
    
    # 求解线性方程组
    v_fdm[l_back-1, :] = np.linalg.solve(A, -v_fdm[l_back, :]/d_t-f)

# 拼接边界值和内部值
v_fdm = np.column_stack((np.column_stack((va_store, v_fdm)), vb_store))

# 绘制3D结果图
T, X = np.meshgrid(t_grid, x_grid)
fig = go.Figure(data=[go.Surface(x=X, y=T, z=v_fdm.T)])
fig.update_traces(colorscale='Viridis')
fig.update_layout(scene=dict(xaxis_title='x', yaxis_title='t', zaxis_title='v_fdm(t,x)'))
fig.show()

问题现象 #

运行代码后得到的结果完全偏离预期的解析解（$v(t,x)=tx^2$的光滑曲面），得到的是一个混乱的曲面，而预期应该是随时间和x平方增长的光滑结果。

我已经检查了离散化的推导、边界条件的处理、反向迭代的索引，但还是找不到问题所在，希望各位能帮我排查一下哪里出错了。

备注：内容来源于stack exchange，提问作者numbers and me