多元数值微分的张量积差分公式及其在非均匀网格上的构造
题目描述
考虑一个定义在矩形区域 \([a,b] \times [c,d]\) 上的二元函数 \(f(x,y)\),我们已知它在非均匀网格点 \((x_i, y_j)\) 上的函数值 \(f_{ij}\),其中 \(i=0,1,\dots,m\),\(j=0,1,\dots,n\)。目标是构造一个高精度的张量积型差分公式,用以计算该函数在区域内部某个点 \((x^*, y^*)\)(不一定是网格点)处的混合偏导数 \(\frac{\partial^{p+q} f}{\partial x^p \partial y^q}\) 的数值近似。特别地,公式需要适应非均匀网格分布,并给出截断误差的表达式。
解题过程
第一步:理解问题核心与挑战
- 多元数值微分:与一元函数不同,多元函数的偏导数需要沿多个方向分别近似,并且混合偏导数涉及交叉差分。
- 张量积思想:对于矩形区域上的规则网格(即使是非均匀的),我们可以将二维问题分解为两个一维问题的组合。即先沿 \(x\) 方向构造一维差分公式,再沿 \(y\) 方向构造一维差分公式,然后通过张量积(外积)的方式组合。
- 非均匀网格:网格节点在 \(x\) 和 \(y\) 方向上可以不均匀分布,这要求我们使用基于非均匀节点的一维插值多项式来构造差分公式,而不是简单的等距差分系数。
- 目标导数:对于混合偏导数 \(\frac{\partial^{p+q} f}{\partial x^p \partial y^q}\),我们可以先对 \(x\) 求 \(p\) 阶偏导,再对 \(y\) 求 \(q\) 阶偏导(顺序可交换)。
第二步:一维非均匀节点的差分公式构造(以 \(x\) 方向为例)
假设我们在 \(x\) 方向上有 \(m+1\) 个非均匀节点 \(x_0 < x_1 < \dots < x_m\),对应函数值 \(f(x_i, y)\)(注意此时 \(y\) 暂时视为固定)。我们想要在点 \(x^*\)(可能不在网格点上)近似 \(p\) 阶导数。
- 选择局部插值节点:为平衡精度与稳定性,通常选取 \(x^*\) 附近的 \(k+1\) 个节点(\(k \ge p\)),记为 \(x_{i_0}, x_{i_1}, \dots, x_{i_k}\),满足 \(x^*\) 位于这些节点的最小闭区间内。
- 构造 Lagrange 插值多项式:
\[ L_x(x) = \sum_{r=0}^{k} f(x_{i_r}, y) \cdot \ell_r(x), \quad \ell_r(x) = \prod_{\substack{s=0 \\ s \ne r}}^{k} \frac{x - x_{i_s}}{x_{i_r} - x_{i_s}}. \]
- 对插值多项式求导:将 \(L_x(x)\) 对 \(x\) 求 \(p\) 阶导数,然后计算在 \(x=x^*\) 的值:
\[ \frac{\partial^p f}{\partial x^p}(x^*, y) \approx \sum_{r=0}^{k} f(x_{i_r}, y) \cdot \ell_r^{(p)}(x^*), \]
其中 \(\ell_r^{(p)}(x^*)\) 是 Lagrange 基函数 \(\ell_r(x)\) 在 \(x^*\) 处的 \(p\) 阶导数值。这些系数可以通过递推公式或直接对基函数求导得到。
- 截断误差:根据多项式插值误差的导数,一维近似误差为:
\[ E_x = \frac{f^{(k+1)}(\xi_x, y)}{(k+1)!} \cdot \prod_{s=0}^{k} (x^* - x_{i_s}), \]
其中 \(\xi_x\) 位于节点所张区间内。注意,对于导数近似,实际截断误差还涉及对误差项求导,通常与节点间距的 \(k+1-p\) 次幂相关。
第三步:推广到二维(张量积构造)
现在考虑 \(y\) 方向。假设在 \(y\) 方向选取 \(l+1\) 个非均匀节点 \(y_{j_0}, y_{j_1}, \dots, y_{j_l}\)(\(l \ge q\))用于逼近 \(y\) 方向的导数。
- 先固定 \(y\) 方向节点:对于每个选定的 \(y_{j_t}\)(\(t=0,\dots,l\)),我们已经有了 \(x\) 方向的近似:
\[ \frac{\partial^p f}{\partial x^p}(x^*, y_{j_t}) \approx \sum_{r=0}^{k} f(x_{i_r}, y_{j_t}) \cdot a_r, \quad a_r = \ell_r^{(p)}(x^*). \]
这里 \(a_r\) 是仅依赖于 \(x\) 方向节点和 \(x^*\) 的常数系数。
- 再沿 \(y\) 方向进行导数逼近:将上述近似值视为关于 \(y\) 的函数(在节点 \(y_{j_t}\) 处已知),构造 \(y\) 方向的 Lagrange 插值多项式并求 \(q\) 阶导数:
\[ \frac{\partial^{p+q} f}{\partial x^p \partial y^q}(x^*, y^*) \approx \sum_{t=0}^{l} \left[ \sum_{r=0}^{k} f(x_{i_r}, y_{j_t}) \cdot a_r \right] \cdot b_t, \]
其中 \(b_t = m_t^{(q)}(y^*)\),\(m_t(y)\) 是关于 \(y\) 方向节点的 Lagrange 基函数。
- 合并为张量积形式:
\[ \frac{\partial^{p+q} f}{\partial x^p \partial y^q}(x^*, y^*) \approx \sum_{r=0}^{k} \sum_{t=0}^{l} f(x_{i_r}, y_{j_t}) \cdot a_r b_t. \]
这就是最终的数值微分公式,它是一个双求和形式,系数 \(a_r b_t\) 是 \(x\) 和 \(y\) 两个方向一维差分系数的外积。
第四步:公式特点与误差分析
- 节点选择:通常取 \(k = p + d_x\),\(l = q + d_y\),其中 \(d_x, d_y\) 为额外节点数以提高精度(如 \(d_x = d_y = 2\) 可得到较高阶近似)。
- 误差来源:
- 插值误差的导数:总误差由两个方向的误差共同贡献。对于混合偏导数,截断误差的主要项与两个方向节点间距的幂次有关。若 \(x\) 方向节点最大间距为 \(h_x\),\(y\) 方向为 \(h_y\),则误差通常为 \(O(h_x^{k+1-p} + h_y^{l+1-q})\),交叉项可能更高阶。
- 数值稳定性:非均匀节点可能导致 Lagrange 基函数的导数系数很大(特别是高阶导数时),从而放大函数值的噪声。应确保节点分布不太极端,或考虑使用重心 Lagrange 形式提高稳定性。
- 在非网格点 \((x^*, y^*)\) 处求导:公式允许在任意点求导,只需在该点附近选取局部节点即可。如果 \((x^*, y^*)\) 恰好是网格点,可以简化计算。
第五步:具体计算步骤示例(以计算 \(\frac{\partial^2 f}{\partial x \partial y}\) 为例)
假设 \(p=1, q=1\),选取 \(x^*\) 附近3个节点(\(k=2\)),\(y^*\) 附近3个节点(\(l=2\))。
- 确定局部节点:设 \(x\) 方向节点为 \(x_{i-1}, x_i, x_{i+1}\),\(y\) 方向节点为 \(y_{j-1}, y_j, y_{j+1}\),且 \((x^*, y^*)\) 位于这些节点围成的矩形内。
- 计算一维差分系数:
- 对 \(x\) 方向,计算基函数在 \(x^*\) 处的一阶导数 \(\ell_{i-1}^{(1)}(x^*), \ell_i^{(1)}(x^*), \ell_{i+1}^{(1)}(x^*)\)(公式略)。
- 对 \(y\) 方向类似,得到 \(m_{j-1}^{(1)}(y^*), m_j^{(1)}(y^*), m_{j+1}^{(1)}(y^*)\)。
- 构造张量积系数表:系数 \(c_{r,t} = \ell_{x_r}^{(1)}(x^*) \cdot m_{y_t}^{(1)}(y^*)\),其中 \(r,t \in \{-1,0,1\}\) 对应节点偏移索引。
- 计算近似值:
\[ \frac{\partial^2 f}{\partial x \partial y}(x^*, y^*) \approx \sum_{r=-1}^{1} \sum_{t=-1}^{1} f(x_{i+r}, y_{j+t}) \cdot c_{r,t}. \]
第六步:应用注意事项
- 边界点处理:若 \((x^*, y^*)\) 靠近区域边界,可能无法取到对称的局部节点,此时应使用非对称的节点集合,精度可能会下降。
- 高阶混合偏导数:方法可直接推广到更高阶混合偏导数,只需增加每个方向的节点数(\(k > p, l > q\))以提高精度。
- 与均匀网格对比:在均匀网格下,系数 \(a_r, b_t\) 可预先算出(成为标准中心/前向/后向差分系数),计算更简单。非均匀网格需要根据节点位置实时计算系数。
- 多维推广:对于 \(d\) 元函数,可类似通过 \(d\) 个一维差分公式的张量积构造,但计算量随维度指数增长(维度灾难)。
总结
本题目展示了一种基于张量积思想的多元数值微分方法,适用于非均匀网格。其核心是将高维问题分解为一维问题的组合,利用 Lagrange 插值多项式求导得到差分系数,再通过外积方式组合。该方法具有系统性和可推广性,但需要注意非均匀节点带来的稳定性问题以及维度增长时的计算成本。