计算科学 - 如何获得有限差分，它是连续的 - 吾爱随笔录

计算科学有限差分 C C++

2021-12-25 10:01:55

我想计算一个有限差分（类似这样的SO Post）。我的数据如下：我的 x 值是 2 的幂（4、8、16、32 和 64）。对应于它们的是 y 值，这样

Y = f (X)

$Y=f(X)$

在哪里 $f$ 是单调函数。我的兴趣是计算有限差分，例如

Diff (x) = (y_{2} - y_{1}) / (x_{2} - x_{1})

$\text{Diff}(x) = (y_2-y_1) /(x_2-x_1)$ 在 X 的不同值，其中 X 在域中（即

4 \leq x \leq 64

$4\le x \le 64$ ）。现在我的问题是，当我采用分段线性近似时，Diff(x) 在 4 到 8 和 8 到 16 之间变化很大。这是不可取的。如果我尝试通过一次取三个点来拟合二次方程，然后计算有限差分，它会在某些点变为负数，这在物理上是不可接受的，因为

f

$f$ 是单调函数。换句话说，我装了

a x^{2} + b x + c

$ax^2+bx+c$ 因此

Diff (X) = 2 a x + b

$\text{Diff}(X)=2ax+b$ . 所以，问题是

2 a x + b

$2ax+b$ 在域中的所有点上都不是正数。

有人可以指出解决此问题的方法。我需要使用 C/C++ 来解决这个问题，我不能使用 matlab。我只有五个 X 值，所以任何复杂的解决方案都是可以接受的。谢谢您的帮助。

3个回答

在这种情况下，一种简单而准确的方法是使用标准差 $\log$ 空间。

\begin{aligned} s (x) & = \log x \\ \partial_{x} f (x) & = \partial_{s} f (e^{s (x)}) \partial_{x} s (x) \\ = \frac{1}{x} \partial_{s} f (e^{s (x)}) \end{aligned}

$\begin{align*} s(x) &= \log x \\ \partial_x f(x) &= \partial_s f(e^{s(x)}) \partial_x s(x) \\ &= \frac 1 x \partial_s f(e^{s(x)}) \\ \end{align*}$

自从 $s(x)$ 等距，很容易使用任何标准差分来近似底线上的导数。

您不想拟合一般多项式，而是考虑使用分段三次样条在每个网格点之间的间隔内插值函数。三次样条的一个好的特性是它们可以最小化拟合曲线的弯曲能量，这可能会导致更少的激进振荡（这可能是您之前尝试过的二次方法遇到的负值的原因）。

您可以通过提出一个单独的三次多项式来获得分段三次样条

$f_i(x)=a_ix^3+b_ix^2+c_ix+d_i$ 在每个区间 $i$

其端点是您选择的网格点。然后，您可以施加以下条件： 1. 所有函数必须与每个节点处的数据一致 2. 函数必须在每个节点点处连续 3. 每个函数的一阶和二阶导数必须在每个节点处连续的节点点

您可以通过施加所谓的“自然样条”条件来获得两个额外的方程，即二阶导数在第一个和最后一个节点处必须为零。这应该为您提供足够的信息来解决每个间隔中每个三次多项式的每个系数。

通过参数化的函数族拟合（如度数 $n$ -多项式，其中有 $n+1$ 自由参数对应于最小化参数空间上的误差项。如果误差项是最小二乘法，那么这会导致二次规划。

如果你希望你的函数是非负的，或者是一个单调函数（不清楚你的意思），那么你会得到一个约束最小化，这本身就是另一个主题。

在您的情况下，低阶多项式就足够了。第一篇维基百科文章应该会带来很多解决问题的方法。

其它你可能感兴趣的问题