★ FEATURED ARTICLE 简介面向数值分析课程学习与上机实践的作业文档围绕 MATLAB 编程实现三类经典数值算法适合正在做课程上机题、需要参照完整求解流程与结果验证思路的本科生与研究生使用。文档以三次样条插值函数为主线之一给出第一型样条通用程序的编制方法、车门曲线型值点算例与端点条件处理并打印节点间的插值结果重积分部分采用复化梯形公式配合逐次二分与外推思想给出通用程序结构及精度 ε 控制下的计算实例常微分方程初值问题则覆盖 RK4、AB4 及 AB4-AM4 预测校正方法的程序实现并将计算结果与精确解对照列出四种方法的误差范数与精度排序。压缩包共 1 个 docx 文件约 566KB内容集中于算法公式推导、程序思路与运行结果整理。已有 904 人学习可作为上机报告撰写与算法复现的参考。1. 三次样条为什么“分段低次”反而比“整体高次”更接近真曲线数值分析上机题里最容易被低估的一道是 3 次样条插值。直觉告诉人给 11 个型值点用一次高阶多项式拟合岂不比分段拼凑更“整体”、更精确真跑一遍就翻车了——10 次多项式在区间两端剧烈摆动中间倒是贴着数据走这就是 Runge 现象的经典现场。同样这 11 个点改用分段三次、只要求一阶和二阶导数在节点处连续曲线在整段区间上几乎没有折角两端导数还能精确落在给定的 0.8 和 0.2 上。这个反差揭示了一条准则逼近精度不来自多项式次数而来自约束的结构。想通“约束怎么选、方程组怎么闭合”这两件事样条就不再是背公式而是一道可以自己推出来的题。2. 三弯矩方程把“处处二阶连续”压成一组三对角方程三次样条的定义听起来很松每个小区间上是一个三次多项式拼起来要二阶连续。松到什么程度如果每段三次多项式有 4 个待定系数n 段就是 4n 个未知量而插值条件只有 2n 个每段两个端点值再加 n-1 个一阶导连续、n-1 个二阶导连续一共 4n-2 个。差 2 个——这 2 个自由度就是端点条件要补的位置。数一遍自由度比背“第一型、第二型、周期型”有用得多因为它告诉你端点条件不是可选项是方程组的闭合条件少一个就解不出来多一个就过约束。2.1 换一组未知量从“每段四个系数”到“每个节点的二阶导”直接拿 4n 个系数列方程是灾难矩阵又大又不稀疏。真正让样条变得可算的一步是把未知量从“每段的系数”换成“每个节点上的二阶导数值” $M_j S(x_j)$。这一步的妙处在于在每个小区间上$S(x)$ 是线性的两端值一给整段的二阶导就被拉格朗日插值唯一确定再积分两次一阶导和函数值也只剩两个积分常数而这两个常数恰好被两端函数值 $y_{j-1}$、$y_j$ 定死。也就是说$M$ 一旦确定整条曲线就唯一确定一个多余的自由度都没有。未知量从 4n 个塌缩成 n1 个方程数也正好是 n1 个n-1 个内点的一阶导连续条件加上 2 个端点条件。矩阵从稠密变成三对角求解复杂度从 $O(n^3)$ 掉到 $O(n)$。这是整个推导里最关键的一次“换元”值得单独记一笔。2.2 μ、α、d 三个系数是怎么从差商里长出来的把一阶导连续条件代进三次 Hermite 形式整理后得到每个内点 $j1,\dots,n-1$ 的一条方程$$\mu_j M_{j-1} 2M_j \alpha_j M_{j1} d_j$$其中$$\mu_j \frac{h_{j-1}}{h_{j-1}h_j},\quad \alpha_j \frac{h_j}{h_{j-1}h_j} 1-\mu_j,\quad h_j x_j - x_{j-1}$$$$d_j 6,f[x_{j-1},x_j,x_{j1}] \frac{6}{h_{j-1}h_j}\left(\frac{y_{j1}-y_j}{h_j} - \frac{y_j-y_{j-1}}{h_{j-1}}\right)$$这三个量里$\mu_j$ 和 $\alpha_j$ 只跟网格有关$d_j$ 才是真正携带数据的那个。注意 $d_j$ 的量纲它是二阶差商乘 6也就是“斜率的变化率”。$d_j$ 为正说明曲线在此处向上弯$d_j$ 为负说明向下弯。这解释了后面那张 $d$ 向量表为什么在 $i6$ 的位置突然变负——$y$ 从 5.78 掉到 5.40斜率由正转负曲率必然出现一个大负数。选用汽车门轮廓这组数据$x_i 0,1,\dots,10$$y$ 从 2.51 到 5.80作为验算靶子是因为它的网格是均匀的$h_j \equiv 1$所有 $\mu_j \alpha_j 0.5$$d_j$ 退化成 $3(\Delta_j - \Delta_{j-1})$手算两三个就能对上。这组数我建议每个做这道题的人都手推一遍作为调代码的锚点$j$012345$d_j$-0.06-0.15-0.24-0.42-0.60-0.24$j$678910$d_j$-1.861.65-0.12-0.090.60其中 $d_0$ 和 $d_{10}$ 用的是端点导数$d_0 \frac{6}{h_0}\left(\frac{y_1-y_0}{h_0} - y_0\right) 6(0.79-0.8) -0.06$$d_{10} \frac{6}{h_9}\left(y_{10} - \frac{y_{10}-y_9}{h_9}\right) 6(0.2-0.10)0.60$。如果你算出来的 $d$ 跟这张表整体差一个 2 倍因子几乎可以断定是忘了除 $(h_{j-1}h_j)$或者把 6 抄成了 3——这是这道题最高频的翻车点后面单独说。2.3 第一型边界条件把开区间方程组闭合的两把锁内点方程只覆盖 $j1$ 到 $jn-1$首尾两个 $M_0$ 和 $M_n$ 悬空。第一型边界条件给的是端点一阶导 $S(x_0)y_0$、$S(x_n)y_n$把它们代入端点处的导数表达式得到两条只涉及两个未知量的方程$$2M_0 M_1 \frac{6}{h_0}\left(\frac{y_1-y_0}{h_0} - y_0\right),\qquad M_{n-1} 2M_n \frac{6}{h_{n-1}}\left(y_n - \frac{y_n-y_{n-1}}{h_{n-1}}\right)$$形状和内点方程完全一致只是 $\mu_00,\alpha_01$ 和 $\mu_n1,\alpha_n0$。于是整个方程组可以写成一个统一的 $(n1)\times(n1)$ 三对角系统主对角全是 2副对角是 $\mu_j$ 和 $\alpha_j$。这个矩阵有一个非常重要的性质主对角元 2副对角元之和 $\mu_j \alpha_j 1 2$严格对角占优。严格对角占优保证追赶法Thomas 算法消元过程中不会出现主元趋近于零的情况也就意味着不需要选主元、不会数值失稳。样条插值之所以在实践中几乎“从不失败”根源就在这条对角占优上——它跟数据长什么样无关只跟 $h_j0$ 有关。2.4 组装、追赶法与中点闭式解下面这段是把上面三节直接翻译成代码。注意 MATLAB 是 1 基索引而公式里 $j$ 从 0 到 $n$我在注释里标了对应关系免得读代码时把 $h_{j-1}$ 和h(j-1)搞混function M spline3_M(x, y, dy0, dyn) % x, y : 型值点长度 n1 的行/列向量需同向 % dy0,dyn: 第一型边界给定的两端一阶导 % M : 各节点的二阶导值长度 n1 x x(:); y y(:); n numel(x) - 1; h diff(x); % h(j) x(j1)-x(j)1 基对应公式的 h_{j} mu zeros(n1,1); la zeros(n1,1); d zeros(n1,1); % 内点公式 j1..n-1 对应代码下标 j1 2..n for j 2:n mu(j) h(j-1) / (h(j-1) h(j)); la(j) 1 - mu(j); d(j) 6 * ( (y(j1)-y(j))/h(j) - (y(j)-y(j-1))/h(j-1) ) / (h(j-1)h(j)); end % 左端点2M0 M1 d0即 mu0, la1 mu(1) 0; la(1) 1; d(1) 6/h(1) * ((y(2)-y(1))/h(1) - dy0); % 右端点M_{n-1} 2M_n d_n即 mu1, la0 mu(n1) 1; la(n1) 0; d(n1) 6/h(n) * (dyn - (y(n1)-y(n))/h(n)); % 追赶法前向消元 b 2*ones(n1,1); for k 2:n1 w mu(k) / b(k-1); b(k) b(k) - w * la(k-1); d(k) d(k) - w * d(k-1); end % 回代 M zeros(n1,1); M(n1) d(n1) / b(n1); for k n:-1:1 M(k) (d(k) - la(k)*M(k1)) / b(k); end end两个循环加起来恰好 $2n$ 次乘除没有任何矩阵求逆。b是消元后被修改的主对角w是消元因子因为严格对角占优b(k)始终远离零所以没有加保护判断也不会出问题。求出 $M$ 之后在区间 $[x_j, x_{j1}]$ 上代回$$S(x) M_j\frac{(h_j-t)^3}{6h_j} M_{j1}\frac{t^3}{6h_j} \left(y_j - \frac{M_j h_j^2}{6}\right)\frac{h_j-t}{h_j} \left(y_{j1} - \frac{M_{j1}h_j^2}{6}\right)\frac{t}{h_j}$$其中 $t x - x_j$。题目要求打印 $S(i0.5)$在均匀网格 $h\equiv 1$ 下把 $t0.5$ 代进去能化简出一个相当干净的闭式$$S(i0.5) \frac{y_i y_{i1}}{2} - \frac{M_i M_{i1}}{16}$$这个式子给了一个不依赖程序正确性的旁证中点值等于两端函数值的算术平均再减去两端二阶导之和的 1/16。如果 $M$ 整体符号搞反了这个修正项就会变成“加”中点值全部往外鼓画出来肉眼就能看出不对。类似 $S(x_j)y_j$ 恒成立、$S(x_0)0.8$、$S(x_{10})0.2$ 这类不变量写代码时就该顺手加进自检别等交完作业才发现。3. 重积分与外推为什么两个梯形值一组合就跳阶二维积分的题目表面上考“会写循环”实际上考的是外推。设$$I(f) \int_c^d\left(\int_a^b f(x,y),dx\right)dy$$直接上二维复化梯形误差是 $O(h^2)O(k^2)$也就是网格每加密一倍误差降到四分之一。这个速度对 $\varepsilon 0.5\times10^{-5}$ 这种要求来说太慢了要想误差小两个数量级网格得加密十倍二维下计算量涨一百倍。外推干的事情是不加密十倍只用两次相邻网格的结果做一次线性组合把误差阶从二阶抬到四阶。3.1 二维复化梯形张量积结构让两次一维积分等价$a$ 到 $b$ 方向分 $n$ 等份$c$ 到 $d$ 方向分 $m$ 等份节点是 $x_i$ 与 $y_j$ 的笛卡尔积。把二维梯形公式展开会发现它等价于“先对每个 $y_j$ 沿 $x$ 做一次一维复化梯形得到 $G(y_j)\approx\int_a^b f(x,y_j)dx$再对 $G$ 沿 $y$ 做本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
