符号看得有点花…… 如果你想做的是「在 x_1 和 x_2 分别趋近 x_0 後所得的极限 = Taylor 多项式」， 那你需要的就是 MVT of divided differences。 https://en.wikipedia.org/wiki/Mean_value_theorem_(divided_differences) 直接套上去就马上做完了。 也不必去算新多项式的导数。 但是如果要一步一步来就没那麽好算了。 (x_n - x_0)lim_{x_1→x_0} f[x_0,...,x_n] = f[x_0,x_2,...,x_n] - lim_{x_1→x_0} f[x_0,...,x_{n-1}] 上面这条递回式是用来算极限的。 本来的插值多项式是 f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)^2。 在 x_1→x_0 之後，变成 f(x_0) + f'(x_0)(x-x_0) + {f[x_0,x_2]-f'(x_0)}/(x_2-x_0) * (x-x_0)^2。 然後 {f[x_0,x_2] - f'(x_0)}/(x_2 - x_0) 在 x_2→x_0 下的极限 = f"(x_0)/2， 所以多项式的极限就变成 f(x_0) + f'(x_0)(x-x_0) + f"(x_0)/2 * (x-x_0)^2。 不过我本来在想的是用 Lagrange 观点。 e_0(x) := Π_{i=1}^n (x-x_i)/(x_0-x_i)，其他 e_j 类推。 还是先用 n = 2 来观察， 插值多项式 = f(x_0)e_0(x) + f(x_1)e_1(x) + f(x_2)e_2(x) 然後也先让 x_1→x_0，多项式变成 f(x_0){e_0(x)+e_1(x)} + {f(x_1)-f(x_0)}e_1(x) + f(x_2)e_2(x)。 所以我们分成三项来看： 1. e_0(x)+e_1(x) = (x-x_1)(x-x_2)/(x_0-x_1)(x_0-x_2) + (x-x_0)(x-x_2)/(x_1-x_0)(x_1-x_2) 公因式 = (x-x_2)/(x_1-x_0) * (x_1-x_0)(x_0-x_2+x_1-x)/(x_0-x_2)(x_1-x_2) = (x-x_2)(x_0-x_2+x_1-x)/(x_0-x_2)(x_1-x_2) → (x-x_2)(2x_0-x_2-x)/(x_0-x_2)^2 = 1 - (x-x_0)^2/(x_0-x_2)^2 最後这个多项式，他代 x_0 得 1、导数得 0，而代 x_2 得 0。 2. {f(x_1)-f(x_0)}e_1(x) = {f(x_1)-f(x_0)}(x-x_0)(x-x_2)/(x_1-x_0)(x_1-x_2) → f'(x_0)(x-x_0)(x-x_2)/(x_0-x_2) (x-x_0)(x-x_2)/(x_0-x_2) 代 x_0 得 0、导数得 1，而代 x_2 得 0。 3. e_2(x) = (x-x_0)(x-x_1)/(x_2-x_0)(x_2-x_1) → (x-x_0)^2/(x_2-x_0)^2 最後这个多项式也是代 x_0 得 0、导数得 0，而代 x_2 得 1。 我们得到三个可以各自突显 f(x_0), f'(x_0), f(x_2) 的多项式， 刚好跟 Lagrange 观点有谋而合。 最後再让 x_2→x_0， f(x_0){1-(x-x_0)^2/(x_0-x_2)^2} + f'(x_0)(x-x_0)(x-x_2)/(x_0-x_2) + f(x_2)(x-x_0)^2/(x_2-x_0)^2 = f(x_0)+f'(x_0)(x-x_0) + { f(x_2) - f(x_0) - f'(x_0)(x_2-x_0) }(x-x_0)^2/(x_2-x_0)^2 → f(x_0) + f'(x_0)(x-x_0) + f"(x_0)(x-x_0)^2/2 其实仔细看，1, x-x_0, (x-x_0)^2/2 也是 在函数值、一阶导数、二阶导数之中各自突显一项，而消灭其他两项的多项式函数， 同样符合 Lagrange 观点的插值概念。 真正麻烦的还是 general case： 有资料的点是 x_0, ..., x_n，每个点的高阶导数已知阶数不尽相同。 像是已知 f(-1), f(0), f'(0), f"(0), f(5), f(100), f'(100) 这样。 然後先用 -1, 0, a, b, 5, 100, c 插值，再让 a,b→0 和 c→100， 之後要检查在 x = 0 的一阶二阶导数和在 x = 100 的一阶导数。 不过我想，应该也是这样一步步算极限就好。 但是那个 general form 就真的很难看，所以平常都是给 algorithm。 --

※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 1.162.224.247 (台湾) ※ 文章网址: https://webptt.com/cn.aspx?n=bbs/Math/M.1691321340.A.CC3.html 先修一点 cap 的错漏。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/07/2023 04:56:47 把 divided differences 算一算，感觉很有趣。 在 f is smooth enough 的前提下，divided differences 其实都可以用极限延拓。 说的是类似 f[a,a] 这种东西可以自然定义成 f'(a)。 不过这样一来，如果 f 是 C^1，就保证 f[x,y] 有 C^0。 而要让 f[x,y,z] 能处处连续，f 至少也得是 C^2。 总之，f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)^2 的极限， 就是 f[x_0] + f[x_0,x_0](x-x_0) + f[x_0,x_0,x_0](x-x_0)^2。 f[x_0,x_0] 就是 f'(x_0) 没问题，而 f[x_0,x_0,x_0] 则是 f"(x_0)/2。 有公式为证： d/dx f[z_1,z_2,...,z_n,x] = f[z_1,z_2,...,z_n,x,x] 用数学归纳法甚至可以得到： (D_x^3/3!)(D_y^2/2!)(D_z^2/2!) f[x,y,z] = f[x,x,x,x,y,y,y,z,z,z]。 f[x_0,x_0,x_0] 的情况比较简单，就是 f[x_0] 在 x_0 的二阶导数再除以 2!。 如果是已知 f, f', f" 在 -1, 0, 1 上的值， x f(x) f'(x) f"(x) -1 2 -8 56 0 1 0 0 1 2 8 56 那这种 p(x) 的 Newton form 就是 f[-1] + f[-1,-1](x+1) + f[-1,-1,-1](x+1)^2 + f[-1,-1,-1,0](x+1)^3 + f[-1,-1,-1,0,0]x(x+1)^3 + f[-1,-1,-1,0,0,0]x^2(x+1)^3 + f[-1,-1,-1,0,0,0,1]x^3(x+1)^3 + f[-1,-1,-1,0,0,0,1,1]x^3(x+1)^3(x-1) + f[-1,-1,-1,0,0,0,1,1,1]x^3(x+1)^3(x-1)^2 然後 wiki 上那张表，就只是很理所当然地计算这些系数而已。 所以你应该是卡在： 1. f[1,2,3,4,5] 我会，但 f[-1,-1,-1,0,0] 到底是什麽？ 2. 为什麽 f[-1,b,c,0,e] 在 b,c→-1 且 e→0 的时候会收敛到 f[-1,-1,-1,0,0]？ 根据前面的脉络，两个问题的回答是一起的： f[a,b,c,d,e] 在 Ｒ^5 上无法直接用 divided difference 写下的位置， 以其极限取代之。则 f[a,b,c,d,e] 在 Ｒ^5 上连续。 总之就是要确认 a divided difference with repeated arguments is well defined。 f[a,a] = lim_{x→a} { f(x)-f(a) }/{ x-a } = f'(a) 而一般一点的情况，建议从 expanded form 下手。 以 f[0,0,0,1,1] 为例， f[0,b,c,1,e] = f(0)/(-b)(-c)(-1)(-e) + f(b)/b(b-c)(b-1)(b-e) + f(c)/c(c-b)(c-1)(c-e) + f(1)/1(1-b)(1-c)(1-e) + f(e)/e(e-b)(e-c)(e-1) 後二项在 b,c→0 的时候，会趋近於 -f(1)/(e-1) + f(e)/e^3(e-1)。 然後在 e→1 的时候会再趋近於 f'(1)-3f(1)。 前三项在 e→1 的时候，会→ f(0)/bc - f(b)/b(c-b)(b-1)^2 + f(c)/c(c-b)(c-1)^2。 最後在 b,c→0 下会趋近於 3f(0) + 2f'(0) + f"(0)/2。 关於连续性，适合的参考资料应该是 https://ftp.cs.wisc.edu/Approx/deboor2.pdf。 f[x,y] 在 Ｒ^2 上连续，这直接做就好，没有很难做。 更多变数的情况下就要一些技巧了。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/08/2023 13:45:45 要 f[x,y] 连续至少要 C^1，因为沿着 y=x 靠近 (a,a) 的时候要 f' 连续。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/08/2023 17:31:27 https://i.imgur.com/kDaTMXM.png 对，我知道因为 MVT 的要求是 "f conti. on [a,b], f is diff. on (a,b)."， 所以如果只是要 f[0,b,c,1,e] 会收敛到 f[0,0,0,1,1] 上，可能可以放宽条件， 甚至上面这个应该不用到四阶导数存在，只要二阶可导就可以了。 可是大家应该也都知道…… (partially) derivable, differentiable, continuously differentiable 很烦人。 但我的确是想先建构函数再考虑，毕竟，f[x,y,z,u,v] 那种函数很美嘛。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/09/2023 17:34:42 找到比较一劳永逸的方式： 首先，这次既不是用 Newton form，也不是用 Lagrange form， 毕竟 p(x) = Σ_{i=1}^n (c_i x^i) 这个 standard form 还是最容易求导数的长相。 已知 f(x_i) for i = 0, 1, ... , n，而且 x_i 各不相同。 那麽，https://i.imgur.com/sPmLIUo.png。 因为 Vandermonde determinant 不是 0，所以这组 c_i 有唯一解。 下一步是让 x_1, x_2, ... , x_{m_0} 都趋近於 x_0， 不过在那之前，要先处理一下前 1 + m_0 列。 把上图的等式拿来做以下列运算： for(int i=1; i<=1+m_0; i++) for(int j=1+m_0; j>=i; j--) (第 j 列 -= 第 j-1 列) /= (x_j - x_{j-i}); 总之，经过这串列运算以後，等式会变成 https://i.imgur.com/7UsrbCz.png。 其中的「１」代表 1 函数，各个「ｘ^k」则各自代表 k 次方函数。 其实１[x_0,x_1] = 0，上部矩阵的下三角都是 0。 然後ｘ[x_0,x_1] = 1，上部矩阵的主对角线上都是 1。 下部矩阵没有动到，照抄。 接下来要确认一下他的行列式值。 虽然长相很凶恶，但是因为我们之前有纪录列运算的过程， 所以实际上是 Vand. det. / sub-Vand. det. of {x_0, ... , x_{m_0}}， 所以这个行列式 = Π_{j>i>m_0} (x_j - x_i) * Π_{j>m_0≧i} (x_j - x_i)， 在 x_1, x_2, ... , x_{m_0} 都趋近於 x_0 的时候， 收敛到 Π_{j>i>m_0} (x_j - x_i) * Π_{j>m_0} (x_j - x_0)^{1+m_0} ≠ 0。 这表示如果那个方阵的极限存在的话，行列式值非零。 终於要算极限了，前面那个方程式左侧的方阵和右侧的行矩阵各自都是收敛的， 而且方阵极限的行列式值非零，所以 c 那一个行矩阵也收敛。 整个方程式的极限是 https://i.imgur.com/AoSFHbG.png。 跟刚刚一样，其实上部矩阵的下三角都是 0，而且主对角线上都是 1。 只是为了能有个通式的长相就拉他们下水。 下部矩阵没有动到，照抄。 然後反覆把想拿来简并的 x_k 并在一起，我们就得到了退化多项式 p(x) 的系数 c_i。 前 1 + m_0 列已经不会再被动到了。 p(x_0) = f(x_0) 可直接参照矩阵方程式的第一列。 观察第二列可得 p'(x_0) = 1 + 2x_0 + 3x_0^2 + ... + nx_0^{n-1} = f'(x_0)， 同理可得其他高阶导数的等式。 这个作法就是从最初的插值多项式直接退化成 p(x)， 并且证明了直到第「简并数」阶之前，p 在简并点上的导数 = f 在简并点上的导数。 ------------------------- 至於 f[a,b,c] 收敛到 f"(a)/2!，似乎真的可以用 f" 存在来证。 总之先对 b,c 做 MVT： 如果 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 收敛到 f"(a)/2， 那 f[a,b,c] 就收敛到 f"(a)/2!。 可是罗下去会卡在 f" 的连续性上，所以虽然我爱罗但是不能罗， 因为罗後不收敛不代表罗前也不收敛。 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 = [ f'(ξ)-f'(a) ]/(ξ-a) - [ f(ξ)-f(a)-f'(a)(ξ-a) ]/(ξ-a)^2 →　f"(a) - f"(a)/2 = f"(a)/2 前项直接算极限，後项则是先罗一次。 前面做 MVT 也很辛苦，f[a,x] 的连续性是显然的， 而 f[a,x] 的可微性则要考虑是否在 a 微分。 不在 a 的时候，很简单，商法则可搞定。 在 a 的话，[ f[a,x]-f'(a) ]/(x-a) 把分母处理好以後就是先罗一次， 其实跟上面那个先罗一次的後项是一样的东西，所以也是收敛到 f"(a)/2。 f[a,b,c] = d/dx f[a,x] |x=ξ for some ξ between b and c， 这个 ξ 有可能是 a， 所以 f[a,b,c] 可能是 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 或 f"(a)/2。 最後当 b,c→a 的时候，也让 ξ→a，然後就回到前头的计算了。 可是再多一个 d 的时候…… 感觉上应该是可以做 MIT 的，但我觉得累。 f[a,b,c,d] = (f[a,b,c]-f[a,b,d])/(c-d) 且 c≠d，进行 MVT： 1. 检查 f[a,b,x] 的连续性。 虽然 b≠a，但是 x 可能是两者之一。 所以对於 x 的连续性是依赖於 f \in C^1。 f 三阶可微，所以 f 当然 \in C^1。 2. 检查 f[a,b,x] 的可微性。 当 x 不是 a 也不是 b 的时候，直接用公式算导函数。 然後在 a 上计算 (f[a,b,x]-f[a,b,a])/(x-a) 的极限， 也就是 f[a,a,x,b] 的极限，这个应该可以由 C^2 保证。 （f[b,b,x,a] 在 b 的情形是类似的，因为此时 a,b 有对称性。） 到此，应用 MVT 得 f[a,b,c,d] = d/dx f[a,b,x] |x=ξ = lim_{x→ξ} (f[a,b,x]-f[a,b,ξ])/(x-ξ) = lim_{x→ξ} f[a,b,ξ,x] 而这个 ξ 只是介於 c,d 之间而已，可能是两数之间的任何数。 接下来分成 ξ 刚好是 a：此时 f[a,b,c,d] = f[a,b,a,a] （∵C^2） ξ 刚好是 b：此时 f[a,b,c,d] = f[a,b,b,b] （∵C^2） ξ 两者皆非：此时 f[a,b,c,d] = f[a,b,ξ,ξ] （∵C^1） 所以综合来说，f[a,b,c,d] = f[a,b,ξ,ξ] 都是对的。 然後计算 b,ξ→a 的时候 f[a,b,c,d] 的极限， 基於 b,ξ 的异同情况，f[a,b,c,d]=f[a,b,b,b] 或 f[a,η,η,η]， 统一写作 f[a,b,c,d] = f[a,η,η,η]，其中 η 介於 b,ξ 之间或者就是 b， 最後利用三阶可微能算出他收敛到 f[a,a,a,a] = f"'(a)/3!。 感觉上应该还是要善用 divided difference 的符号， 不过扯上简并的连续性跟可微性都得验证。 或者可以针对 b,c,d 直接做某种推广版的 MVT： f[a,b,c,d] = f[a,η,η,η], where min(b,c,d)<η<max(b,c,d). ※ 编辑: Vulpix (1.160.12.97 台湾), 08/13/2023 07:26:34 f[x_0,x_1] = [ f(x_1)-f(x_0) ]/(x_1-x_0) f 代入 1 函数，得 (1-1)/(x_1-x_0) = 0。 f 代入 k 次方函数，得 (x_1^k-x_0^k)/(x_1-x_0) = Σ_{i=0}^{k-1} x_0^i x_1^{k-1-i}。 例如 k = 2 的 2 次方函数， ｘ^2[x_0,x_1] = (x_1^2-x_0^2)/(x_1-x_0) = x_1+x_0， ｘ^2[x_1,x_2] = (x_2^2-x_1^2)/(x_2-x_1) = x_2+x_1， ｘ^2[x_0,x_1,x_2] = [ (x_2+x_1) - (x_1+x_0) ]/(x_2-x_0) = 1 这样。 所以其实（上部矩阵的）下三角都是 0、（上部矩阵的）主对角线都是 1。 那就把次方都改成 max(0, 原本的次方) 好了…… 或者乾脆改成原次方的绝对值。 不过反正负次方项的系数都是 0 啦，当成 notation 看就好，没有什麽代值的功能。 首先，不管 Newton、Lagrange 插出来的多项式都跟 https://i.imgur.com/sPmLIUo.png 解出来的多项式是同一个。 所以虽然我算的是这个多项式的退化，但可以直接当成 Newton form 的极限，没问题。 你把这个矩阵方程式展开，其实每一条都只是 p(x_i) = f(x_i) 而已。 因为我想要都用 divided difference 写，所以会比 CV 多除以一些阶乘。 也是因为这样，我的系数才是组合数 C 而不是 CV 的排列数 P。 对，但我不要。这样的行列式里面会有一大堆阶乘。 我的证明过程须要算行列式，不想让他变丑。 数学归纳法。 不过我还在想到底有哪些东西要定义跟证明XD 以这个流程来说，至少就是 f[x_1, ..., x_n] 收敛到 f^{(n)}(x_0)/n! 吧。 是啊。 ※ 编辑: Vulpix (1.160.12.97 台湾), 08/13/2023 22:58:39 既然反矩阵麻烦，那为什麽不用克拉玛公式呢？ 作者: Vulpix (Sebastian) 看板: Math 标题: Re: [分析] Hermite内插演算法的证明 时间: Sun Aug 6 19:28:57 2023 符号看得有点花…… 如果你想做的是「在 x_1 和 x_2 分别趋近 x_0 後所得的极限 = Taylor 多项式」， 那你需要的就是 MVT of divided differences。 https://en.wikipedia.org/wiki/Mean_value_theorem_(divided_differences) 直接套上去就马上做完了。 也不必去算新多项式的导数。 但是如果要一步一步来就没那麽好算了。 (x_n - x_0)lim_{x_1→x_0} f[x_0,...,x_n] = f[x_0,x_2,...,x_n] - lim_{x_1→x_0} f[x_0,...,x_{n-1}] 上面这条递回式是用来算极限的。 本来的插值多项式是 f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)^2。 在 x_1→x_0 之後，变成 f(x_0) + f'(x_0)(x-x_0) + {f[x_0,x_2]-f'(x_0)}/(x_2-x_0) * (x-x_0)^2。 然後 {f[x_0,x_2] - f'(x_0)}/(x_2 - x_0) 在 x_2→x_0 下的极限 = f"(x_0)/2， 所以多项式的极限就变成 f(x_0) + f'(x_0)(x-x_0) + f"(x_0)/2 * (x-x_0)^2。 不过我本来在想的是用 Lagrange 观点。 e_0(x) := Π_{i=1}^n (x-x_i)/(x_0-x_i)，其他 e_j 类推。 还是先用 n = 2 来观察， 插值多项式 = f(x_0)e_0(x) + f(x_1)e_1(x) + f(x_2)e_2(x) 然後也先让 x_1→x_0，多项式变成 f(x_0){e_0(x)+e_1(x)} + {f(x_1)-f(x_0)}e_1(x) + f(x_2)e_2(x)。 所以我们分成三项来看： 1. e_0(x)+e_1(x) = (x-x_1)(x-x_2)/(x_0-x_1)(x_0-x_2) + (x-x_0)(x-x_2)/(x_1-x_0)(x_1-x_2) 公因式 = (x-x_2)/(x_1-x_0) * (x_1-x_0)(x_0-x_2+x_1-x)/(x_0-x_2)(x_1-x_2) = (x-x_2)(x_0-x_2+x_1-x)/(x_0-x_2)(x_1-x_2) → (x-x_2)(2x_0-x_2-x)/(x_0-x_2)^2 = 1 - (x-x_0)^2/(x_0-x_2)^2 最後这个多项式，他代 x_0 得 1、导数得 0，而代 x_2 得 0。 2. {f(x_1)-f(x_0)}e_1(x) = {f(x_1)-f(x_0)}(x-x_0)(x-x_2)/(x_1-x_0)(x_1-x_2) → f'(x_0)(x-x_0)(x-x_2)/(x_0-x_2) (x-x_0)(x-x_2)/(x_0-x_2) 代 x_0 得 0、导数得 1，而代 x_2 得 0。 3. e_2(x) = (x-x_0)(x-x_1)/(x_2-x_0)(x_2-x_1) → (x-x_0)^2/(x_2-x_0)^2 最後这个多项式也是代 x_0 得 0、导数得 0，而代 x_2 得 1。 我们得到三个可以各自突显 f(x_0), f'(x_0), f(x_2) 的多项式， 刚好跟 Lagrange 观点有谋而合。 最後再让 x_2→x_0， f(x_0){1-(x-x_0)^2/(x_0-x_2)^2} + f'(x_0)(x-x_0)(x-x_2)/(x_0-x_2) + f(x_2)(x-x_0)^2/(x_2-x_0)^2 = f(x_0)+f'(x_0)(x-x_0) + { f(x_2) - f(x_0) - f'(x_0)(x_2-x_0) }(x-x_0)^2/(x_2-x_0)^2 → f(x_0) + f'(x_0)(x-x_0) + f"(x_0)(x-x_0)^2/2 其实仔细看，1, x-x_0, (x-x_0)^2/2 也是 在函数值、一阶导数、二阶导数之中各自突显一项，而消灭其他两项的多项式函数， 同样符合 Lagrange 观点的插值概念。 真正麻烦的还是 general case： 有资料的点是 x_0, ..., x_n，每个点的高阶导数已知阶数不尽相同。 像是已知 f(-1), f(0), f'(0), f"(0), f(5), f(100), f'(100) 这样。 然後先用 -1, 0, a, b, 5, 100, c 插值，再让 a,b→0 和 c→100， 之後要检查在 x = 0 的一阶二阶导数和在 x = 100 的一阶导数。 不过我想，应该也是这样一步步算极限就好。 但是那个 general form 就真的很难看，所以平常都是给 algorithm。 --

※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 1.162.224.247 (台湾) ※ 文章网址: https://webptt.com/cn.aspx?n=bbs/Math/M.1691321340.A.CC3.html 先修一点 cap 的错漏。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/07/2023 04:56:47 把 divided differences 算一算，感觉很有趣。 在 f is smooth enough 的前提下，divided differences 其实都可以用极限延拓。 说的是类似 f[a,a] 这种东西可以自然定义成 f'(a)。 不过这样一来，如果 f 是 C^1，就保证 f[x,y] 有 C^0。 而要让 f[x,y,z] 能处处连续，f 至少也得是 C^2。 总之，f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)^2 的极限， 就是 f[x_0] + f[x_0,x_0](x-x_0) + f[x_0,x_0,x_0](x-x_0)^2。 f[x_0,x_0] 就是 f'(x_0) 没问题，而 f[x_0,x_0,x_0] 则是 f"(x_0)/2。 有公式为证： d/dx f[z_1,z_2,...,z_n,x] = f[z_1,z_2,...,z_n,x,x] 用数学归纳法甚至可以得到： (D_x^3/3!)(D_y^2/2!)(D_z^2/2!) f[x,y,z] = f[x,x,x,x,y,y,y,z,z,z]。 f[x_0,x_0,x_0] 的情况比较简单，就是 f[x_0] 在 x_0 的二阶导数再除以 2!。 如果是已知 f, f', f" 在 -1, 0, 1 上的值， x f(x) f'(x) f"(x) -1 2 -8 56 0 1 0 0 1 2 8 56 那这种 p(x) 的 Newton form 就是 f[-1] + f[-1,-1](x+1) + f[-1,-1,-1](x+1)^2 + f[-1,-1,-1,0](x+1)^3 + f[-1,-1,-1,0,0]x(x+1)^3 + f[-1,-1,-1,0,0,0]x^2(x+1)^3 + f[-1,-1,-1,0,0,0,1]x^3(x+1)^3 + f[-1,-1,-1,0,0,0,1,1]x^3(x+1)^3(x-1) + f[-1,-1,-1,0,0,0,1,1,1]x^3(x+1)^3(x-1)^2 然後 wiki 上那张表，就只是很理所当然地计算这些系数而已。 所以你应该是卡在： 1. f[1,2,3,4,5] 我会，但 f[-1,-1,-1,0,0] 到底是什麽？ 2. 为什麽 f[-1,b,c,0,e] 在 b,c→-1 且 e→0 的时候会收敛到 f[-1,-1,-1,0,0]？ 根据前面的脉络，两个问题的回答是一起的： f[a,b,c,d,e] 在 Ｒ^5 上无法直接用 divided difference 写下的位置， 以其极限取代之。则 f[a,b,c,d,e] 在 Ｒ^5 上连续。 总之就是要确认 a divided difference with repeated arguments is well defined。 f[a,a] = lim_{x→a} { f(x)-f(a) }/{ x-a } = f'(a) 而一般一点的情况，建议从 expanded form 下手。 以 f[0,0,0,1,1] 为例， f[0,b,c,1,e] = f(0)/(-b)(-c)(-1)(-e) + f(b)/b(b-c)(b-1)(b-e) + f(c)/c(c-b)(c-1)(c-e) + f(1)/1(1-b)(1-c)(1-e) + f(e)/e(e-b)(e-c)(e-1) 後二项在 b,c→0 的时候，会趋近於 -f(1)/(e-1) + f(e)/e^3(e-1)。 然後在 e→1 的时候会再趋近於 f'(1)-3f(1)。 前三项在 e→1 的时候，会→ f(0)/bc - f(b)/b(c-b)(b-1)^2 + f(c)/c(c-b)(c-1)^2。 最後在 b,c→0 下会趋近於 3f(0) + 2f'(0) + f"(0)/2。 关於连续性，适合的参考资料应该是 https://ftp.cs.wisc.edu/Approx/deboor2.pdf。 f[x,y] 在 Ｒ^2 上连续，这直接做就好，没有很难做。 更多变数的情况下就要一些技巧了。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/08/2023 13:45:45 要 f[x,y] 连续至少要 C^1，因为沿着 y=x 靠近 (a,a) 的时候要 f' 连续。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/08/2023 17:31:27 https://i.imgur.com/kDaTMXM.png 对，我知道因为 MVT 的要求是 "f conti. on [a,b], f is diff. on (a,b)."， 所以如果只是要 f[0,b,c,1,e] 会收敛到 f[0,0,0,1,1] 上，可能可以放宽条件， 甚至上面这个应该不用到四阶导数存在，只要二阶可导就可以了。 可是大家应该也都知道…… (partially) derivable, differentiable, continuously differentiable 很烦人。 但我的确是想先建构函数再考虑，毕竟，f[x,y,z,u,v] 那种函数很美嘛。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/09/2023 17:34:42 找到比较一劳永逸的方式： 首先，这次既不是用 Newton form，也不是用 Lagrange form， 毕竟 p(x) = Σ_{i=1}^n (c_i x^i) 这个 standard form 还是最容易求导数的长相。 已知 f(x_i) for i = 0, 1, ... , n，而且 x_i 各不相同。 那麽，https://i.imgur.com/sPmLIUo.png。 因为 Vandermonde determinant 不是 0，所以这组 c_i 有唯一解。 下一步是让 x_1, x_2, ... , x_{m_0} 都趋近於 x_0， 不过在那之前，要先处理一下前 1 + m_0 列。 把上图的等式拿来做以下列运算： for(int i=1; i<=1+m_0; i++) for(int j=1+m_0; j>=i; j--) (第 j 列 -= 第 j-1 列) /= (x_j - x_{j-i}); 总之，经过这串列运算以後，等式会变成 https://i.imgur.com/7UsrbCz.png。 其中的「１」代表 1 函数，各个「ｘ^k」则各自代表 k 次方函数。 其实１[x_0,x_1] = 0，上部矩阵的下三角都是 0。 然後ｘ[x_0,x_1] = 1，上部矩阵的主对角线上都是 1。 下部矩阵没有动到，照抄。 接下来要确认一下他的行列式值。 虽然长相很凶恶，但是因为我们之前有纪录列运算的过程， 所以实际上是 Vand. det. / sub-Vand. det. of {x_0, ... , x_{m_0}}， 所以这个行列式 = Π_{j>i>m_0} (x_j - x_i) * Π_{j>m_0≧i} (x_j - x_i)， 在 x_1, x_2, ... , x_{m_0} 都趋近於 x_0 的时候， 收敛到 Π_{j>i>m_0} (x_j - x_i) * Π_{j>m_0} (x_j - x_0)^{1+m_0} ≠ 0。 这表示如果那个方阵的极限存在的话，行列式值非零。 终於要算极限了，前面那个方程式左侧的方阵和右侧的行矩阵各自都是收敛的， 而且方阵极限的行列式值非零，所以 c 那一个行矩阵也收敛。 整个方程式的极限是 https://i.imgur.com/AoSFHbG.png。 跟刚刚一样，其实上部矩阵的下三角都是 0，而且主对角线上都是 1。 只是为了能有个通式的长相就拉他们下水。 下部矩阵没有动到，照抄。 然後反覆把想拿来简并的 x_k 并在一起，我们就得到了退化多项式 p(x) 的系数 c_i。 前 1 + m_0 列已经不会再被动到了。 p(x_0) = f(x_0) 可直接参照矩阵方程式的第一列。 观察第二列可得 p'(x_0) = 1 + 2x_0 + 3x_0^2 + ... + nx_0^{n-1} = f'(x_0)， 同理可得其他高阶导数的等式。 这个作法就是从最初的插值多项式直接退化成 p(x)， 并且证明了直到第「简并数」阶之前，p 在简并点上的导数 = f 在简并点上的导数。 ------------------------- 至於 f[a,b,c] 收敛到 f"(a)/2!，似乎真的可以用 f" 存在来证。 总之先对 b,c 做 MVT： 如果 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 收敛到 f"(a)/2， 那 f[a,b,c] 就收敛到 f"(a)/2!。 可是罗下去会卡在 f" 的连续性上，所以虽然我爱罗但是不能罗， 因为罗後不收敛不代表罗前也不收敛。 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 = [ f'(ξ)-f'(a) ]/(ξ-a) - [ f(ξ)-f(a)-f'(a)(ξ-a) ]/(ξ-a)^2 →　f"(a) - f"(a)/2 = f"(a)/2 前项直接算极限，後项则是先罗一次。 前面做 MVT 也很辛苦，f[a,x] 的连续性是显然的， 而 f[a,x] 的可微性则要考虑是否在 a 微分。 不在 a 的时候，很简单，商法则可搞定。 在 a 的话，[ f[a,x]-f'(a) ]/(x-a) 把分母处理好以後就是先罗一次， 其实跟上面那个先罗一次的後项是一样的东西，所以也是收敛到 f"(a)/2。 f[a,b,c] = d/dx f[a,x] |x=ξ for some ξ between b and c， 这个 ξ 有可能是 a， 所以 f[a,b,c] 可能是 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 或 f"(a)/2。 最後当 b,c→a 的时候，也让 ξ→a，然後就回到前头的计算了。 可是再多一个 d 的时候…… 感觉上应该是可以做 MIT 的，但我觉得累。 f[a,b,c,d] = (f[a,b,c]-f[a,b,d])/(c-d) 且 c≠d，进行 MVT： 1. 检查 f[a,b,x] 的连续性。 虽然 b≠a，但是 x 可能是两者之一。 所以对於 x 的连续性是依赖於 f \in C^1。 f 三阶可微，所以 f 当然 \in C^1。 2. 检查 f[a,b,x] 的可微性。 当 x 不是 a 也不是 b 的时候，直接用公式算导函数。 然後在 a 上计算 (f[a,b,x]-f[a,b,a])/(x-a) 的极限， 也就是 f[a,a,x,b] 的极限，这个应该可以由 C^2 保证。 （f[b,b,x,a] 在 b 的情形是类似的，因为此时 a,b 有对称性。） 到此，应用 MVT 得 f[a,b,c,d] = d/dx f[a,b,x] |x=ξ = lim_{x→ξ} (f[a,b,x]-f[a,b,ξ])/(x-ξ) = lim_{x→ξ} f[a,b,ξ,x] 而这个 ξ 只是介於 c,d 之间而已，可能是两数之间的任何数。 接下来分成 ξ 刚好是 a：此时 f[a,b,c,d] = f[a,b,a,a] （∵C^2） ξ 刚好是 b：此时 f[a,b,c,d] = f[a,b,b,b] （∵C^2） ξ 两者皆非：此时 f[a,b,c,d] = f[a,b,ξ,ξ] （∵C^1） 所以综合来说，f[a,b,c,d] = f[a,b,ξ,ξ] 都是对的。 然後计算 b,ξ→a 的时候 f[a,b,c,d] 的极限， 基於 b,ξ 的异同情况，f[a,b,c,d]=f[a,b,b,b] 或 f[a,η,η,η]， 统一写作 f[a,b,c,d] = f[a,η,η,η]，其中 η 介於 b,ξ 之间或者就是 b， 最後利用三阶可微能算出他收敛到 f[a,a,a,a] = f"'(a)/3!。 感觉上应该还是要善用 divided difference 的符号， 不过扯上简并的连续性跟可微性都得验证。 或者可以针对 b,c,d 直接做某种推广版的 MVT： f[a,b,c,d] = f[a,η,η,η], where min(b,c,d)<η<max(b,c,d). ※ 编辑: Vulpix (1.160.12.97 台湾), 08/13/2023 07:26:34 f[x_0,x_1] = [ f(x_1)-f(x_0) ]/(x_1-x_0) f 代入 1 函数，得 (1-1)/(x_1-x_0) = 0。 f 代入 k 次方函数，得 (x_1^k-x_0^k)/(x_1-x_0) = Σ_{i=0}^{k-1} x_0^i x_1^{k-1-i}。 例如 k = 2 的 2 次方函数， ｘ^2[x_0,x_1] = (x_1^2-x_0^2)/(x_1-x_0) = x_1+x_0， ｘ^2[x_1,x_2] = (x_2^2-x_1^2)/(x_2-x_1) = x_2+x_1， ｘ^2[x_0,x_1,x_2] = [ (x_2+x_1) - (x_1+x_0) ]/(x_2-x_0) = 1 这样。 所以其实（上部矩阵的）下三角都是 0、（上部矩阵的）主对角线都是 1。 那就把次方都改成 max(0, 原本的次方) 好了…… 或者乾脆改成原次方的绝对值。 不过反正负次方项的系数都是 0 啦，当成 notation 看就好，没有什麽代值的功能。 首先，不管 Newton、Lagrange 插出来的多项式都跟 https://i.imgur.com/sPmLIUo.png 解出来的多项式是同一个。 所以虽然我算的是这个多项式的退化，但可以直接当成 Newton form 的极限，没问题。 你把这个矩阵方程式展开，其实每一条都只是 p(x_i) = f(x_i) 而已。 因为我想要都用 divided difference 写，所以会比 CV 多除以一些阶乘。 也是因为这样，我的系数才是组合数 C 而不是 CV 的排列数 P。 对，但我不要。这样的行列式里面会有一大堆阶乘。 我的证明过程须要算行列式，不想让他变丑。 数学归纳法。 不过我还在想到底有哪些东西要定义跟证明XD 以这个流程来说，至少就是 f[x_1, ..., x_n] 收敛到 f^{(n)}(x_0)/n! 吧。 是啊。 ※ 编辑: Vulpix (1.160.12.97 台湾), 08/13/2023 22:58:39 既然反矩阵麻烦，那为什麽不用克拉玛公式呢？ 你误会了，x_0^{-N} 之类的只是为了让通式漂亮硬塞的而已。 「C(n,m), n<m」才是 0。 组合解释：从少数物品中，任取多数物品，方法数 = 0。 或者也可以把 (n-m)! 想像成 Γ(n-m+1) 这些发散到无限大的东西。 所以才麻烦，因为一直有长相不合群的在捣蛋。 不过毕竟要假设 3 阶可微了，那 f[x,y] 就是可微分的函数应该没有错， 至少 f[x,y] 在所有位置上的两个偏导数都存在。 那要从 f[a,a] 做出 f[a,a,c] 或者从 f[a,b] 做出 f[a,b,a] 都不是难事。 Ex. 1) f[a,a,c] f[a] 〉f[a,a] f[a] 〉f[a,a,c] 〉f[a,c] f[c] 一开始我们的表只有白字，其中 f[a,a]=f'(a)。 为了算 f[a,a,c] 加了一条红字，因为 a 与 c 不同， 所以用普通的均差递回关系就好。 Ex. 2) f[a,b,a] f[a] 〉f[a,b] f[b] 〉f[a,b,a] 〉f[b,a] f[a] 这次因为 a 与 a 相同，所以要算 f[b,x] 的导数了。 f[a,b,a] = lim_{x→a} (f[b,x]-f[a,b])/(x-a) = lim_{x→a} f[a,b,x]。 虽然看起来很美好，但因为我们的假设很弱，所以这个极限突然下不了手。 幸好均差是对称函数，所以这个极限还可以不用爆开来算。 f[a,b,a] = lim_{x→a} f[a,b,x] = lim_{x→a} f[a,x,b] = (f[a,b]-f[a,a])/(b-a) = f[a,a,b] 回到了跟前一个例子相同的情况。 我是这样想的：f[x,a_1,...,a_n,y] 都由那个均差表来定义。 如果 x 与 y 不同，那就是 (f[a_1,...,a_n,y]-f[x,a_1,...,a_n])/(y-x)。 但要是他们一样，那就改成上式在 y→x 下的极限。 当然综合来说都是 lim_{z→y} (f[a_1,...,a_n,z]-f[x,a_1,...,a_n])/(z-x)。 也因为是这样定义的，所以至少 f[x,a_1,...,a_n,z] 对 z 是连续的。 然後靠对称性（我怀疑这对於有简并情形的均差是要证明一下的）移动 z 的位置， 就能知道均差对所有自变数都连续，只是要注意一点：都是视为单变数函数。 所以如果不要太贪心，一口气让 x_i 都趋近於 x_0， 而是一个一个来，那证明应该会轻松许多。 照上面的定义来看第三个例子。 Ex. 3) f[a,a,a] f[a,a,a] = lim_{x→a} (f[a,x]-f[a,a])/(x-a) = lim_{x→a} [ (f(x)-f(a))/(x-a) - f'(a) ]/(x-a) = lim_{x→a} (f(x)-f(a)-f'(a)(x-a))/(x-a)^2 = lim_{x→a} (f'(x)-f'(a))/2(x-a) = f"(a)/2 最後一步不能罗，要直接用导数定义。 反正怎麽搞这一步都跑不掉，索性就不要把 f[a,a] 定义成 f'(a) 了。 看着会觉得有趣的地方在分子那边，其实会是 f(x)-泰勒多项式， 所以罗起来自然非常流畅。 m_1 个点退化到另一个点上、m_2 个点退…… 完成这个 modified CV 就可以了。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/14/2023 10:06:09 这个是直接插多项式出来，所以会在比较早的位置看到各阶导数。 不过，就算先为了方便起见而把 x_i 依序排好， 他也跟均差一样，遇到高阶项就要背公式。 均差表的公式是 f[a,a,a] = f"(a)/2! 这些， 而 Neville's algorithm 的是 P_{0,2}(x) = f(a)+f'(a)(x-a)+f"(a)/2! (x-a)^2。 这都不是单纯的「把前一层的结果代入递回公式并让 x_1,x_2→a」。 要算这类极限，在算到这一层之前都不可以先算极限。 不过的确只要把这点先处理好，就可以比较轻松地检查各点上的导数。 毕竟 Newton form 最後有「选」一条从零阶均差直到 n 阶均差的路径：top row， 所以只有 x_0 那一点的导数好算，其他点的导数要找「其他路径」去看。 所以使用 Newton form 的时候，要随时记得「终点一样的路就代表同一个多项式」。 也就是虽然我们只写出一个 Newton form， 但是在把均差表完成的当下就已经有了同一个多项式的 2^n 个 Newton form。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/14/2023 16:49:24 直接上实例：f(x) = x^3，选相异四数 a,b,c,d 算插值多项式。 P_aa 〉P_ab P_bb 〉P_ac 〉P_bc 〉P_ad = P P_cc 〉P_bd 〉P_cd P_dd a^3 〉(a^2+ab+b^2)x-ab(a+b) b^3 〉(a+b+c)x^2-(ab+bc+ca)x+abc 〉(b^2+bc+c^2)x-bc(b+c) 〉x^3 c^3 〉(b+c+d)x^2-(bc+cd+db)x+bcd 〉(c^2+cd+d^2)x-cd(c+d) d^3 然後让 b,c→a： a^3 〉 {3a^2}x-2a^3 a^3 〉 {3a}x^2-{3a^2}x+a^3 〉 {3a^2}x-2a^3 〉x^3 a^3 〉(2a+d)x^2-(a^2+2ad)x+{a^2}d 〉(a^2+ad+d^2)x-ad(a+d) d^3 其中 {3a}x^2-{3a^2}x+a^3 ≠ lim_{t→a} [(x-a)({3t^2}x-2t^3)-(x-t)({3a^2}x-2a^3)]/(t-a) ＝ {6a}x^2 - {9a^2}x + 4a^3 三项系数分别是正确值的 2、3、4 倍。 不像 {3a^2}x-2a^3 = lim_{t→a} [(x-a)(t^3)-(x-t)(a^3)]/(t-a) 那麽单纯。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/14/2023 21:27:43 不是，就是乾脆不要先定义成这种形式。 因为 https://i.imgur.com/lylZcfd.png 会更方便。 至少在 f 有一定平滑性的时候这样定义反而比较方便。 （对於不连续的 f，应该是有办法把他也包容进来才对。） 1. 存在即连续，所以不用检查太多次连续性。 2. 要是 x_n = x_0，还能利用 n 个变数的均差的对称性快速得到 https://i.imgur.com/88E2qLe.png 这条我们想要的微分公式。 这样一来，每次均差的变数变多，首先就要证明 x_i 都可以对调。 pf): 1) x_n 可与 x_0 对调。 如果 x_n 与 x_0 相同，那两者对调确实不改变均差。 即便两者不同，定义中的极限变回普通的差分商， https://i.imgur.com/Mn5WvDa.png，两者一样可以对调。 2) x_0 可与 x_1 对调。 如果 x_0 与 x_1 相同，那两者对调确实不改变均差。 如果两者不同，https://i.imgur.com/qt1duqL.png 说明还是可对调。 其中有先假设变数不到 n+1 个的时候可以任意排列变数。 以下补述少变数的情况，剩下的部份由数学归纳法递回可得。 i) 双变数的情况 如果 a≠b：f[a,b] = (f(b)-f(a))/(b-a) = (f(a)-f(b))/(a-b) = f[b,a] 至於 a=b：trivial case。 ii) 三变数的情况 如果 a,b,c 相异，可以直接写出对称式： f[a,b,c] = f(a)/(a-b)(a-c) + f(b)/(b-a)(b-c) + f(c)/(c-a)(c-b) 如果 a＝b≠c： f[a,a,c] = f[c,a,a] = [(f(c)-f(a))/(c-a) -f'(a)]/(c-a) = ( f(c)-f(a)-f'(a)(c-a) )/(c-a)^2 f[a,c,a] = lim_{b→a} (f[c,b]-f[a,c])/(b-a) = d/da f[a,c] = ( f(c)-f(a)-f'(a)(c-a) )/(c-a)^2 至於 a = b = c：trivial case。 3) x_1,...,x_{n-1} 可任意排列。 任取一个 1,...,n-1 之间的排列 σ，https://i.imgur.com/7SJV8cT.png。 这些排列可以生成任意一个在 0,1,...,n 之间的排列，done。 没事，那个你没空可以暂时不去管，因为我觉得那个也有那个的麻烦。 不过解决 Newton form 的 top row 过於偏心的方案可能要了解一下。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/15/2023 02:55:44 是指在後面要计算 f[a,...,a] 的时候，该有的麻烦还是一样在。 top row 的问题是你另一篇回文一开始的问题。 因为 top row 就真的不能让你直接看到 p(0)、p'(1) 那些东西。 这个只要看那篇的推文就好。 对 b,c 做，那当然要先换成 f[b,a,c] 再做。 令 g(x) = f[a,x] = f[x,a] 因为 f'(a) 存在，所以 g 在 a 连续。 又 g[a,x] = f[a,a,x] = (f(x)-f(a)-f'(a)(x-a))/(x-a)^2 → f"(a)/2 其他 x 就不用说了，g(x) 肯定连续，用商法则就能算 g'(x) = f[a,x,x]。 总之 f 只要二阶可微，那 g 就符合 MVT 的使用资格。 所以 f[b,a,c] = (f[a,c]-f[b,a])/(c-b) = g[b,c] = g'(ξ) g'(ξ) 有两种可能长相： 1. ξ = a，此时 g'(ξ) = g'(a) = f"(a)/2。 2. ξ≠a，此时 g'(ξ) = f[a,ξ,ξ] = (f'(ξ)(ξ-a)-f(ξ)+f(a))/(ξ-a)^2 在 b,c→a 的过程中，这两个情况也可能交错发生， 但只要 f[a,ξ,ξ] 也收敛到 f[a,a,a] = f"(a)/2 就好。 证明方式就是罗，上面两串紫字都要罗一次。 不过 f[a,a,x] 可以直接罗，但 f[a,ξ,ξ] 要先换成 f'[a,ξ] - f[a,a,ξ] 才有罗。 这些该罗的，一个都跑不掉。 再开题外：f'[a,ξ] = f[a,ξ,ξ] + f[a,a,ξ] 这种形状的公式， 好像真的应该先做出来放着备用。 ※ 编辑: Vulpix (163.13.112.58 台湾), 08/15/2023 15:15:14 这应该是 f'[x_0,...,x_n] = Σ_{r=0}^n f[x_0,...,x_n, x_r]。 => f"[x_0,...,x_n] = 2! Σ_{0≦r≦s≦n} f[x_0,...,x_n, x_r,x_s] => f"'[x_0,...,x_n] = 3! Σ_{0≦r≦s≦t≦n} f[x_0,...,x_n, x_r,x_s,x_t] and so on. ※ 编辑: Vulpix (163.13.112.58 台湾), 08/15/2023 16:26:02 对，大概本来是要打 f[a,x,x] 吧？ 直接 C^n 就没有那麽多事了。只要求 Diff^n 是会厘清一些事情啦…… Def): The maximal degeneracy of (x_0,...,x_n) is the frequency of the mode of the sequence {x_i}_{i=0}^n. 总之就是，最大简并数 := 众数的出现次数。 f \in Diff^k => f[x_0,...,x_n] is well-defined at the points where the maximal degeneracy ≦ k+1. In particular, f[x_0,...,x_k] is well-defined everywhere. f \in C^k => f[x_0,...,x_n] is continuous at the points where the maximal degeneracy ≦ k+1. In particular, f[x_0,...,x_k] is continuous everywhere. ※ 编辑: Vulpix (163.13.112.58 台湾), 08/16/2023 04:38:19 其实放弃 Lagrange form 也是挺可惜的。 令 u_{ik}(x) 为 Π_{j≠i} {(x-x_j)/(x_i-x_j)}^{-1-m_j} 在 x = x_i 展开至 m_i-k 次项的多项式。 那麽 u_{ik}(x) * (x-x_i)^k/k! * Π_{j≠i} {(x-x_j)/(x_i-x_j)}^{1+m_j} 们 就是我们 Lagrange form 的基底……没有很漂亮，但也没有很不漂亮啦。 具体来说，https://i.imgur.com/ItFFMG9.png。 整个计算现在变成要算 https://i.imgur.com/LXZTetm.png， 而且因不看超过 m_i 次的项，最後一式也可变成 https://i.imgur.com/q85ceOt.png。 其实 u_{ik}(x) 可以用辗转相除、长除法、泰勒展开等方式去做，方便就好。 ※ 编辑: Vulpix (1.160.14.56 台湾), 08/19/2023 08:45:44