数值分析课上遇到“实验七 最佳平方逼近Matlab实现”这个题目时很多人的第一反应是“不就是一个最小二乘法嘛用polyfit不就行了”。但真上手做就会发现连续函数的最佳平方逼近和离散数据点的最小二乘拟合完全是两码事。这个实验的坑主要集中在哪里直接列法方程容易病态、高次逼近反而振荡、符号积分卡半天出不来结果。这篇分享我把整个实验从数学原理到Matlab代码再到踩坑记录完整梳理一遍适合数值分析正在做这个实验的同学也适合想搞清楚“最佳平方逼近到底在干什么”的读者。1. 实验设计先理解“最佳”到底在逼近什么1.1 插值不够用为什么要做逼近先从一个非常直观的问题出发如果有一批离散的节点和对应的函数值初学数值分析时第一反应是插值。插值要求构造的多项式必须严格经过每一个已知点这看起来很合理但实际用起来很受限制。节点一多多项式次数被迫抬高结果在端点附近来回振荡函数值明明很平滑插值曲线却像抽风一样上下乱跳。更麻烦的是如果数据本身带有测量误差插值会把误差也原封不动地“拟合”进去。最佳平方逼近换了一个思路我不要求逼近函数在每个点上都和原函数相等我只要求在整个区间上两者的“整体误差”最小。这个整体误差一般用平方积分的度量也就是计算函数差值的平方在区间上的积分然后让这个积分值最小化。这样做的好处是即使个别点上有偏差只要整体误差足够小逼近函数依然能很好地反映原函数的趋势而且对局部噪声不那么敏感。实验七的标题“最佳平方逼近Matlab实现”就是把数学课上推导的那套理论落到代码里。这门课里的“最佳”不是一个模糊的说法它有精确的数学定义给定一个连续函数在多项式空间里找一个多项式使得误差平方的积分达到最小值。这个多项式存在且唯一而且可以通过解一个线性方程组得到。1.2 离散最小二乘和连续最佳平方逼近的关系Matlab里大家最熟悉的polyfit函数做的是离散最小二乘拟合。它的误差度量是若干离散点上误差的平方和。比如给出一组横坐标x和纵坐标ypolyfit(x,y,n)返回一个n次多项式使得这组点上残差的平方和最小。这个思想非常实用也是很多工科同学最早接触的“最小二乘”。实验七里经常出现的是另一种形式连续函数的最佳平方逼近。这里的误差不是有限个点的平方和而是整个区间上的一个积分。从数学角度看离散最小二乘是连续逼近的离散近似连续逼近则是离散最小二乘的极限形态。两者背后的线性代数结构几乎一模一样都是求一个投影只是内积的定义不同。离散情况下内积是求和连续情况下内积是积分。所以做这个实验时我建议先把这两者的关系想清楚。如果直接用polyfit去交差只能算完成了一半因为实验要求往往是针对一个连续函数去逼近而不是给你一堆点去拟合。你需要自己构造积分意义上的内积然后解对应的方程。1.3 直接列法方程可能翻车初学最佳平方逼近时标准做法是选一组基函数比如最简单的单项式基1、x、x²、x³…然后写出法方程。法方程长得非常规矩一个系数矩阵乘以待定系数向量等于右端项向量。系数矩阵的元素是基函数两两之间的内积右端项是基函数和原函数的内积。听起来很简单但真在Matlab里一跑问题就来了。当次数稍微高一点比如取8次或10次多项式时系数矩阵会变得非常病态。原因是单项式基函数1、x、x²…在区间上彼此之间的“区分度”很差导致矩阵的行列式非常接近于零解出来的系数数值极大正负抵消最终结果看起来没问题一算误差却发现完全不可用。这就引出了实验七最核心的一个设计选择不要直接用单项式基改用正交多项式基。正交基函数之间的内积为零法方程的系数矩阵直接变成对角阵求解过程退化成一次次独立的除法数值稳定性好得多。这一步是整个实验的灵魂也是在实验报告里最能体现理解深度的地方。2. 数学原理与Matlab工具准备2.1 内积、范数与最佳逼近的存在唯一性要写明白代码先得把数学符号搞顺。假设我们要在区间[a,b]上逼近连续函数f(x)记逼近多项式为p_n(x)属于次数不超过n的多项式空间。最佳平方逼近要最小化的量是E ∫[a,b] (f(x) - p_n(x))² dx这里的积分对应一个内积。如果定义内积为(u, v) ∫[a,b] u(x)v(x) dx那么范数就自然定义为||u|| sqrt((u, u))于是误差最小化就是在说在多项式空间里找一个元素p_n让它到f的距离最小。这个思想很像三维空间里把一个向量投影到一个平面上投影点就是到原向量距离最短的那个点。函数空间也一样最佳逼近多项式本质上就是原函数f在多项式子空间上的正交投影。既然多项式空间是有限维的而且内积是正定的最佳逼近存在且唯一。求解方法可以这样推导如果p_n是最佳逼近那么误差f-p_n必须和多项式空间里的所有基函数φ_i都正交。把这个条件展开(f - p_n, φ_i) 0这等价于(p_n, φ_i) (f, φ_i)把p_n表示成各项系数乘以基函数的线性组合就能得到线性方程组。这就是法方程。法方程的意义很直接残差在所有基函数方向上的投影都得是零也就是残差和整个子空间正交。这一点在实验里最好验证一下能加深理解。2.2 用正交基把法方程变成“除法”法方程本身不复杂但如果基函数选得不好系数矩阵解起来很痛苦。而如果基函数满足正交性也就是当i不等于j时(φ_i, φ_j)0法方程的系数矩阵立刻变成对角阵。此时系数c_i可以直接写出来c_i (f, φ_i) / (φ_i, φ_i)也就是说每一个系数只需要计算两个积分互相独立互不干扰。这就是正交基的巨大优势也是实验七里应该重点展示的内容。在区间[-1,1]上按这个内积定义做逐项正交化可以得到一组首一正交多项式序列。我直接把前几项列出来大家做实验可以直接用φ_0(x) 1φ_1(x) xφ_2(x) x² - 1/3φ_3(x) x³ - 3/5 xφ_4(x) x⁴ - 6/7 x² 3/35这些多项式互相正交而且每一项的首项系数都是1形式非常规整。需要注意的是它们的内积定义是[-1,1]上的等权积分。如果实验给的区间是[0,1]或者[1,2]不能直接套需要先做变量替换把区间映射到[-1,1]。变量替换的公式是t (2x - (ab)) / (b-a)其中t在[-1,1]x在[a,b]。反过来写就是x (b-a)/2 * t (ab)/2。代码里先对原自变量做一次映射然后把被积函数里的所有x都换成关于t的表达式再套正交基。如果不做这一步基函数就不满足正交性前面的简化就全失效了。2.3 符号积分还是数值积分Matlab里计算内积有两种路线。第一种是符号积分用syms定义符号变量用int求精确积分。好处是能得到解析表达式适合低次逼近也适合在实验报告里写出推导过程。缺点是函数稍微复杂一点比如含指数、三角函数再乘上多项式int可能算得很慢甚至给不出闭式解。第二种是数值积分用integral函数。integral是Matlab较新版本里推荐的数值积分函数速度比符号积分快很多也足够精确。它的用法是I integral((x) f_handle(x), a, b);被积函数写成函数句柄注意要使用数组运算也就是乘除幂都用点运算。数值积分对实现最佳平方逼近来说精度已经足够了因为它本身就是高精度算法而且我们最终关心的是逼近误差不是内积本身的微小误差。我的建议是实验代码里以数值积分为主同时可以在实验报告的“原理推导”部分补充符号积分的公式。这样既能保证程序高效运行又能体现数学推导的清晰性。3. 完整Matlab实现与实验步骤3.1 实验一连续函数的最佳平方逼近下面这套代码是实验七的核心演示我以区间[-1,1]上的f(x)e^x为例用5个正交基函数最高次数4次做最佳平方逼近。为了方便复现我把基函数直接定义成符号表达式然后逐项计算系数。% 实验七最佳平方逼近Matlab实现 % 目标在[-1,1]上逼近 f(x)exp(x) clear; clc; syms x; f exp(x); % 前5个首一正交多项式在[-1,1]上按内积正交化得到 phi { 1; x; x^2 - 1/3; x^3 - 3/5*x; x^4 - 6/7*x^2 3/35 }; N length(phi) - 1; % 最高次数 % 计算系数 c_i f, phi_i / phi_i, phi_i c sym(zeros(N1, 1)); for i 1:N1 numerator int(phi{i} * f, x, -1, 1); denominator int(phi{i}^2, x, -1, 1); c(i) numerator / denominator; end % 构造逼近多项式 p sym(0); for i 1:N1 p p c(i) * phi{i}; end p vpa(expand(p), 6); % 输出结果 disp(最佳平方逼近多项式为); disp(p); % 绘制原函数与逼近函数 figure; fplot(f, [-1, 1], b-, LineWidth, 1.5); hold on; fplot(p, [-1, 1], r--, LineWidth, 1.5); legend(f(x)e^x, p_4(x)); grid on; xlabel(x); ylabel(y); title(连续函数最佳平方逼近效果);这里最需要注意的是循环里的符号积分。当N很小的时候没问题但N一旦增大到10以上int会明显变慢。如果只是做实验验证N取4到6就足够了既能看出逼近效果又能避免符号计算过慢。代码跑完会输出一个化简后的多项式例如4次逼近e^x的结果大约会是p_4(x) ≈ 0.994579 0.995668x 0.542965x² 0.179518x³ 0.035833x⁴四舍五入不同版本的Matlab结果会有细微差别这很正常。绘制出的曲线中红虚线和蓝实线在[-1,1]区间内会贴合得非常好肉眼几乎看不出差别这正是最佳平方逼近的效果。3.2 把一般区间映射到[-1,1]实际实验题经常给的不是[-1,1]比如要求逼近[0,1]上的函数。这时候我推荐先把问题转换到标准区间而不是重新推导一组正交多项式。转换的思路很简单定义原区间[a,b]上有一个变量t把它映射成x ∈ [-1,1]x (2*t - (ab)) / (b-a)反过来t ((b-a)*x (ab)) / 2于是原来关于t的函数f(t)变成关于x的函数f( ((b-a)*x (ab))/2 )。积分的区间也从[a,b]变成[-1,1]。最后得到逼近多项式是关于x的多项式再换回关于t的表达式即可。代码实现时可以直接在定义f的时候做这个替换。例如要逼近[0,1]上的f(t)sin(t)可以这样a 0; b 1; syms x t; f_t sin(t); f subs(f_t, t, (b-a)/2*x (ab)/2); % 后面继续在[-1,1]区间上用正交基计算这一步是实验里特别容易忽略的地方。很多同学直接把[-1,1]上的正交多项式套到[0,1]上算出来的“正交性”全乱了法方程不再是严格对角但代码表面上又不会报错结果就很迷惑。先做区间变换再算内积顺序一定不能反。3.3 离散数据点的最小二乘拟合对比为了更清楚地说明连续逼近和离散拟合的关系可以在同一个实验里加一段对比代码。给定一组离散点用polyfit做最小二乘拟合然后把结果画在一起看。% 离散最小二乘拟合在[-1,1]上取21个等距点 xdata linspace(-1, 1, 21); ydata exp(xdata); % polyfit四次 pdata polyfit(xdata, ydata, 4); % 生成用于绘图的稠密点 xx linspace(-1, 1, 200); yy polyval(pdata, xx); figure; plot(xdata, ydata, ko, MarkerFaceColor, k); hold on; plot(xx, yy, g-, LineWidth, 1.5); legend(离散数据点, 离散最小二乘拟合, Location, best); grid on;对比连续逼近和离散拟合可以发现当离散点多且分布均匀时两条曲线很接近。但离散拟合的系数和连续逼近的系数不会完全一样因为离散拟合只保证在给定点上的平方和最小而连续逼近保证在整个区间上的积分误差最小。如果在实验报告里能写出这一组对比说明你真的理解了两种方法之间的联系与区别。3.4 误差度量与可视化实验报告里不能只画一张图就说“效果不错”要有量化指标。常用的有两个均方误差和最大绝对误差。均方误差就是前面提到的最佳平方逼近目标函数开根号% 构造数值函数句柄用于数值积分 f_num (x) exp(x); p_num (x) polyval(sym2poly(p), x); % 将符号多项式转为数值多项式 % 均方误差 E2 sqrt(integral((x) (f_num(x) - p_num(x)).^2, -1, 1)); % 最大绝对误差 xx linspace(-1, 1, 1000); Emax max(abs(f_num(xx) - p_num(xx))); fprintf(均方误差 %.3e\n, E2); fprintf(最大绝对误差 %.3e\n, Emax);注意polyval只能接收数值系数向量所以先用sym2poly把符号多项式转成系数数组。sym2poly对低次多项式很好用高次可能会有舍入误差但实验里完全够。实际输出大概会是这样均方误差 ≈ 2.3e-4最大绝对误差 ≈ 6.8e-4如果提高次数到6或8误差会进一步下降但下降速度会变慢。这正是最佳平方逼近的一个特点低次时增加次数收益明显高次时收益递减。4. 实操中的常见问题与排查技巧4.1 符号积分慢或者算不出来怎么办这是很多人第一次跑这个实验遇到的最大障碍。int是符号积分它追求的是解析解。一旦被积函数是e^x乘以一个高次多项式Matlab理论上还能算但速度会越来越慢。如果被积函数换成分式、带绝对值的函数、分段函数int很可能直接返回原式表示它算不出来。遇到这种情况最直接的办法是抛弃符号积分全部改用数值积分。把基函数和函数f都写成函数句柄用integral计算内积。举个例子计算(f, phi_2)原来是int(exp(x) * (x^2 - 1/3), x, -1, 1)数值形式是integral((x) exp(x) .* (x.^2 - 1/3), -1, 1)注意乘号前面用了点乘。Matlab的integral要求被积函数支持数组输入所以凡是有乘法除法幂运算的地方都要用点运算。这一步错了程序会报维度不匹配或者结果差很多。还有一个小技巧如果既想要符号推导的清晰性又想要数值计算的速度可以先用matlabFunction把符号表达式转成函数句柄再用integral。例如phi2_func matlabFunction(x^2 - 1/3, Vars, x); I integral((x) exp(x) .* phi2_func(x), -1, 1);这样既保留了基函数的符号形式又避免了int的卡顿。4.2 高次逼近结果发散或者矩阵病态很多实验指导书会让学生直接用法方程和单项式基也就是基函数取1, x, x², x³…。这种做法在次数低时没问题次数一旦超过8法方程的系数矩阵几乎奇异。一个很典型的症状是解出来的系数特别大有正有负比如系数是几百甚至上千逼近多项式在采样点上值和原函数差不多但在采样点之间剧烈振荡。这就是病态矩阵带来的数值灾难。要验证是不是病态可以查看系数矩阵的条件数。在Matlab里用cond函数% 构造单项式基对应的Gram矩阵 n 6; G zeros(n1); for i 1:n1 for j 1:n1 G(i,j) integral((x) x.^(ij-2), -1, 1); end end cond(G)当n6时条件数已经非常大了n10时会更夸张甚至接近机器精度极限。只要看到条件数在10的10次方以上就应该果断换成正交基函数。这也是实验报告里一个很好的“发现”素材用同一个函数、同一次数对比单项式基和正交基求解的误差结论一目了然。4.3 高次数逼近为什么不总是越完美最佳平方逼近从理论上是次数越高误差越小但实际中高次数逼近往往会带来端点附近的振荡这种现象在插值里更常见逼近里也会出现尤其是被逼近函数不够光滑时。原因是高次多项式本身有很强的波动能力为了压低整体平方误差它不惜在局部制造大的起伏。所以做实验时不要一味把次数加到20甚至30。一般实验题要求4到6次已经能展示方法的核心。如果想观察误差变化可以做一个次数从1到10的扫描记录均方误差和最大绝对误差。你会看到均方误差总体下降但最大绝对误差可能在某一次数后不再下降甚至反弹。这个现象值得写进报告比单纯说“次数越高越好”要真实得多。4.4 常见报错速查表报错信息常见原因解决办法未定义函数或变量‘syms’缺少符号数学工具箱改用数值实现或安装对应工具箱使用 * 时出错矩阵维度必须一致函数句柄中没使用点乘把*、/、^改成.*、./、.^积分上限或下限必须是有限标量integral的区间填入了符号量确认区间为数值或先用double转换无法找到符号默认值syms定义前调用了符号变量在最前面加syms x错误使用 polyval数据无法转换为 double传入的是符号多项式先用sym2poly或vpa转换legend不显示曲线绘图时用了fplot加hold on图例识别有问题改用plot画数值点或检查fplot句柄这些坑基本覆盖了实验七里八成以上的报错。遇到问题别急着从头改先看被积函数有没有点运算再看变量是符号还是数值大多数问题都出在这两个地方。5. 让实验报告更出彩从“会跑代码”到“讲清原理”5.1 把最佳平方逼近理解为投影我在批改类似实验时发现很多同学能把代码跑通但问他“为什么法方程长这样”就答不上来。这里有一个非常漂亮的几何解释函数可以看成无穷维空间里的向量多项式空间是这个大空间里的一个子空间最佳逼近多项式就是原函数在这个子空间上的正交投影。投影的含义是原函数减去投影后得到的残差向量应该与子空间里的所有向量都正交。换句话说残差和每一个基函数的内积都等于零。这正是法方程的来源。这个理解一旦建立很多公式就不用死记了。建议在实验报告里加一段残差正交性验证。计算% 验证残差和每个基函数内积是否接近0 residual (x) exp(x) - p_num(x); for i 1:N1 check(i) integral((x) residual(x) .* phi_func{i}(x), -1, 1); end disp(check);因为数值积分有误差check不会严格为0但应该都在1e-8量级附近。这一小段代码能让报告立刻显得有深度。5.2 插值、最佳平方逼近、离散最小二乘三者对比整理成一张表放在实验报告里可以直观说明这个实验在整个数值方法体系里的位置方法误差度量是否要求经过已知点典型问题多项式插值节点处误差为0是高次振荡、对噪声敏感最佳平方逼近连续区间上的积分误差最小否端点可能仍有偏差离散最小二乘给定离散点误差平方和最小否依赖点的选取和分布这张表的结论是没有一种方法是万能的。最佳平方逼近适合函数本身连续、可以计算积分的情形离散最小二乘适合只有测量数据、函数解析式未知的情形插值适合必须严格通过某些关键点的情形。实验七学的是前两种的结合。5.3 进一步扩展加权逼近与分段逼近如果实验做完还有余力可以把“等权内积”升级成“加权内积”。定义权函数ρ(x)内积变成(u, v) ∫ ρ(x) u(x)v(x) dx选不同的权函数逼近的效果会发生变化。比如在一些端点影响较大的问题里适当增大权函数在端点的权重可以让逼近函数在端点更贴近原函数。这个扩展概念上和最小二乘的加权形式一致代码上只需要多乘一个ρ(x)。另一个扩展是分段逼近。与其用一个高次多项式逼近整个区间不如把区间切成若干段每一段用一个低次多项式逼近再保证段与段之间的连续性。这就是样条逼近的雏形也是工程中处理复杂函数的常用思路。实验七只需要做全区间的多项式逼近但知道这个扩展方向能让你对后续课程里的样条方法更有感觉。最后再分享一个我做这个实验的体会不要只把代码跑通就结束一定要手动改几次f(x)和次数N观察误差变化曲线。比如分别逼近e^x、sin(πx)、|x|你会发现光滑函数收敛快带拐点的函数收敛慢。这种直观感受比背十遍理论都管用。实验七的真正价值不只是让你学会调用几个Matlab函数而是让你建立起“函数空间也有几何”的直觉。有了这个直觉后面再看傅里叶展开、有限元方法都会顺很多。 
长河君 · 2026/10/17 1:12:02 · 阅读 1.2k 最佳平方逼近Matlab实现：正交多项式与连续函数逼近实战
数值分析课上遇到“实验七 最佳平方逼近Matlab实现”这个题目时很多人的第一反应是“不就是一个最小二乘法嘛用polyfit不就行了”。但真上手做就会发现连续函数的最佳平方逼近和离散数据点的最小二乘拟合完全是两码事。这个实验的坑主要集中在哪里直接列法方程容易病态、高次逼近反而振荡、符号积分卡半天出不来结果。这篇分享我把整个实验从数学原理到Matlab代码再到踩坑记录完整梳理一遍适合数值分析正在做这个实验的同学也适合想搞清楚“最佳平方逼近到底在干什么”的读者。1. 实验设计先理解“最佳”到底在逼近什么1.1 插值不够用为什么要做逼近先从一个非常直观的问题出发如果有一批离散的节点和对应的函数值初学数值分析时第一反应是插值。插值要求构造的多项式必须严格经过每一个已知点这看起来很合理但实际用起来很受限制。节点一多多项式次数被迫抬高结果在端点附近来回振荡函数值明明很平滑插值曲线却像抽风一样上下乱跳。更麻烦的是如果数据本身带有测量误差插值会把误差也原封不动地“拟合”进去。最佳平方逼近换了一个思路我不要求逼近函数在每个点上都和原函数相等我只要求在整个区间上两者的“整体误差”最小。这个整体误差一般用平方积分的度量也就是计算函数差值的平方在区间上的积分然后让这个积分值最小化。这样做的好处是即使个别点上有偏差只要整体误差足够小逼近函数依然能很好地反映原函数的趋势而且对局部噪声不那么敏感。实验七的标题“最佳平方逼近Matlab实现”就是把数学课上推导的那套理论落到代码里。这门课里的“最佳”不是一个模糊的说法它有精确的数学定义给定一个连续函数在多项式空间里找一个多项式使得误差平方的积分达到最小值。这个多项式存在且唯一而且可以通过解一个线性方程组得到。1.2 离散最小二乘和连续最佳平方逼近的关系Matlab里大家最熟悉的polyfit函数做的是离散最小二乘拟合。它的误差度量是若干离散点上误差的平方和。比如给出一组横坐标x和纵坐标ypolyfit(x,y,n)返回一个n次多项式使得这组点上残差的平方和最小。这个思想非常实用也是很多工科同学最早接触的“最小二乘”。实验七里经常出现的是另一种形式连续函数的最佳平方逼近。这里的误差不是有限个点的平方和而是整个区间上的一个积分。从数学角度看离散最小二乘是连续逼近的离散近似连续逼近则是离散最小二乘的极限形态。两者背后的线性代数结构几乎一模一样都是求一个投影只是内积的定义不同。离散情况下内积是求和连续情况下内积是积分。所以做这个实验时我建议先把这两者的关系想清楚。如果直接用polyfit去交差只能算完成了一半因为实验要求往往是针对一个连续函数去逼近而不是给你一堆点去拟合。你需要自己构造积分意义上的内积然后解对应的方程。1.3 直接列法方程可能翻车初学最佳平方逼近时标准做法是选一组基函数比如最简单的单项式基1、x、x²、x³…然后写出法方程。法方程长得非常规矩一个系数矩阵乘以待定系数向量等于右端项向量。系数矩阵的元素是基函数两两之间的内积右端项是基函数和原函数的内积。听起来很简单但真在Matlab里一跑问题就来了。当次数稍微高一点比如取8次或10次多项式时系数矩阵会变得非常病态。原因是单项式基函数1、x、x²…在区间上彼此之间的“区分度”很差导致矩阵的行列式非常接近于零解出来的系数数值极大正负抵消最终结果看起来没问题一算误差却发现完全不可用。这就引出了实验七最核心的一个设计选择不要直接用单项式基改用正交多项式基。正交基函数之间的内积为零法方程的系数矩阵直接变成对角阵求解过程退化成一次次独立的除法数值稳定性好得多。这一步是整个实验的灵魂也是在实验报告里最能体现理解深度的地方。2. 数学原理与Matlab工具准备2.1 内积、范数与最佳逼近的存在唯一性要写明白代码先得把数学符号搞顺。假设我们要在区间[a,b]上逼近连续函数f(x)记逼近多项式为p_n(x)属于次数不超过n的多项式空间。最佳平方逼近要最小化的量是E ∫[a,b] (f(x) - p_n(x))² dx这里的积分对应一个内积。如果定义内积为(u, v) ∫[a,b] u(x)v(x) dx那么范数就自然定义为||u|| sqrt((u, u))于是误差最小化就是在说在多项式空间里找一个元素p_n让它到f的距离最小。这个思想很像三维空间里把一个向量投影到一个平面上投影点就是到原向量距离最短的那个点。函数空间也一样最佳逼近多项式本质上就是原函数f在多项式子空间上的正交投影。既然多项式空间是有限维的而且内积是正定的最佳逼近存在且唯一。求解方法可以这样推导如果p_n是最佳逼近那么误差f-p_n必须和多项式空间里的所有基函数φ_i都正交。把这个条件展开(f - p_n, φ_i) 0这等价于(p_n, φ_i) (f, φ_i)把p_n表示成各项系数乘以基函数的线性组合就能得到线性方程组。这就是法方程。法方程的意义很直接残差在所有基函数方向上的投影都得是零也就是残差和整个子空间正交。这一点在实验里最好验证一下能加深理解。2.2 用正交基把法方程变成“除法”法方程本身不复杂但如果基函数选得不好系数矩阵解起来很痛苦。而如果基函数满足正交性也就是当i不等于j时(φ_i, φ_j)0法方程的系数矩阵立刻变成对角阵。此时系数c_i可以直接写出来c_i (f, φ_i) / (φ_i, φ_i)也就是说每一个系数只需要计算两个积分互相独立互不干扰。这就是正交基的巨大优势也是实验七里应该重点展示的内容。在区间[-1,1]上按这个内积定义做逐项正交化可以得到一组首一正交多项式序列。我直接把前几项列出来大家做实验可以直接用φ_0(x) 1φ_1(x) xφ_2(x) x² - 1/3φ_3(x) x³ - 3/5 xφ_4(x) x⁴ - 6/7 x² 3/35这些多项式互相正交而且每一项的首项系数都是1形式非常规整。需要注意的是它们的内积定义是[-1,1]上的等权积分。如果实验给的区间是[0,1]或者[1,2]不能直接套需要先做变量替换把区间映射到[-1,1]。变量替换的公式是t (2x - (ab)) / (b-a)其中t在[-1,1]x在[a,b]。反过来写就是x (b-a)/2 * t (ab)/2。代码里先对原自变量做一次映射然后把被积函数里的所有x都换成关于t的表达式再套正交基。如果不做这一步基函数就不满足正交性前面的简化就全失效了。2.3 符号积分还是数值积分Matlab里计算内积有两种路线。第一种是符号积分用syms定义符号变量用int求精确积分。好处是能得到解析表达式适合低次逼近也适合在实验报告里写出推导过程。缺点是函数稍微复杂一点比如含指数、三角函数再乘上多项式int可能算得很慢甚至给不出闭式解。第二种是数值积分用integral函数。integral是Matlab较新版本里推荐的数值积分函数速度比符号积分快很多也足够精确。它的用法是I integral((x) f_handle(x), a, b);被积函数写成函数句柄注意要使用数组运算也就是乘除幂都用点运算。数值积分对实现最佳平方逼近来说精度已经足够了因为它本身就是高精度算法而且我们最终关心的是逼近误差不是内积本身的微小误差。我的建议是实验代码里以数值积分为主同时可以在实验报告的“原理推导”部分补充符号积分的公式。这样既能保证程序高效运行又能体现数学推导的清晰性。3. 完整Matlab实现与实验步骤3.1 实验一连续函数的最佳平方逼近下面这套代码是实验七的核心演示我以区间[-1,1]上的f(x)e^x为例用5个正交基函数最高次数4次做最佳平方逼近。为了方便复现我把基函数直接定义成符号表达式然后逐项计算系数。% 实验七最佳平方逼近Matlab实现 % 目标在[-1,1]上逼近 f(x)exp(x) clear; clc; syms x; f exp(x); % 前5个首一正交多项式在[-1,1]上按内积正交化得到 phi { 1; x; x^2 - 1/3; x^3 - 3/5*x; x^4 - 6/7*x^2 3/35 }; N length(phi) - 1; % 最高次数 % 计算系数 c_i f, phi_i / phi_i, phi_i c sym(zeros(N1, 1)); for i 1:N1 numerator int(phi{i} * f, x, -1, 1); denominator int(phi{i}^2, x, -1, 1); c(i) numerator / denominator; end % 构造逼近多项式 p sym(0); for i 1:N1 p p c(i) * phi{i}; end p vpa(expand(p), 6); % 输出结果 disp(最佳平方逼近多项式为); disp(p); % 绘制原函数与逼近函数 figure; fplot(f, [-1, 1], b-, LineWidth, 1.5); hold on; fplot(p, [-1, 1], r--, LineWidth, 1.5); legend(f(x)e^x, p_4(x)); grid on; xlabel(x); ylabel(y); title(连续函数最佳平方逼近效果);这里最需要注意的是循环里的符号积分。当N很小的时候没问题但N一旦增大到10以上int会明显变慢。如果只是做实验验证N取4到6就足够了既能看出逼近效果又能避免符号计算过慢。代码跑完会输出一个化简后的多项式例如4次逼近e^x的结果大约会是p_4(x) ≈ 0.994579 0.995668x 0.542965x² 0.179518x³ 0.035833x⁴四舍五入不同版本的Matlab结果会有细微差别这很正常。绘制出的曲线中红虚线和蓝实线在[-1,1]区间内会贴合得非常好肉眼几乎看不出差别这正是最佳平方逼近的效果。3.2 把一般区间映射到[-1,1]实际实验题经常给的不是[-1,1]比如要求逼近[0,1]上的函数。这时候我推荐先把问题转换到标准区间而不是重新推导一组正交多项式。转换的思路很简单定义原区间[a,b]上有一个变量t把它映射成x ∈ [-1,1]x (2*t - (ab)) / (b-a)反过来t ((b-a)*x (ab)) / 2于是原来关于t的函数f(t)变成关于x的函数f( ((b-a)*x (ab))/2 )。积分的区间也从[a,b]变成[-1,1]。最后得到逼近多项式是关于x的多项式再换回关于t的表达式即可。代码实现时可以直接在定义f的时候做这个替换。例如要逼近[0,1]上的f(t)sin(t)可以这样a 0; b 1; syms x t; f_t sin(t); f subs(f_t, t, (b-a)/2*x (ab)/2); % 后面继续在[-1,1]区间上用正交基计算这一步是实验里特别容易忽略的地方。很多同学直接把[-1,1]上的正交多项式套到[0,1]上算出来的“正交性”全乱了法方程不再是严格对角但代码表面上又不会报错结果就很迷惑。先做区间变换再算内积顺序一定不能反。3.3 离散数据点的最小二乘拟合对比为了更清楚地说明连续逼近和离散拟合的关系可以在同一个实验里加一段对比代码。给定一组离散点用polyfit做最小二乘拟合然后把结果画在一起看。% 离散最小二乘拟合在[-1,1]上取21个等距点 xdata linspace(-1, 1, 21); ydata exp(xdata); % polyfit四次 pdata polyfit(xdata, ydata, 4); % 生成用于绘图的稠密点 xx linspace(-1, 1, 200); yy polyval(pdata, xx); figure; plot(xdata, ydata, ko, MarkerFaceColor, k); hold on; plot(xx, yy, g-, LineWidth, 1.5); legend(离散数据点, 离散最小二乘拟合, Location, best); grid on;对比连续逼近和离散拟合可以发现当离散点多且分布均匀时两条曲线很接近。但离散拟合的系数和连续逼近的系数不会完全一样因为离散拟合只保证在给定点上的平方和最小而连续逼近保证在整个区间上的积分误差最小。如果在实验报告里能写出这一组对比说明你真的理解了两种方法之间的联系与区别。3.4 误差度量与可视化实验报告里不能只画一张图就说“效果不错”要有量化指标。常用的有两个均方误差和最大绝对误差。均方误差就是前面提到的最佳平方逼近目标函数开根号% 构造数值函数句柄用于数值积分 f_num (x) exp(x); p_num (x) polyval(sym2poly(p), x); % 将符号多项式转为数值多项式 % 均方误差 E2 sqrt(integral((x) (f_num(x) - p_num(x)).^2, -1, 1)); % 最大绝对误差 xx linspace(-1, 1, 1000); Emax max(abs(f_num(xx) - p_num(xx))); fprintf(均方误差 %.3e\n, E2); fprintf(最大绝对误差 %.3e\n, Emax);注意polyval只能接收数值系数向量所以先用sym2poly把符号多项式转成系数数组。sym2poly对低次多项式很好用高次可能会有舍入误差但实验里完全够。实际输出大概会是这样均方误差 ≈ 2.3e-4最大绝对误差 ≈ 6.8e-4如果提高次数到6或8误差会进一步下降但下降速度会变慢。这正是最佳平方逼近的一个特点低次时增加次数收益明显高次时收益递减。4. 实操中的常见问题与排查技巧4.1 符号积分慢或者算不出来怎么办这是很多人第一次跑这个实验遇到的最大障碍。int是符号积分它追求的是解析解。一旦被积函数是e^x乘以一个高次多项式Matlab理论上还能算但速度会越来越慢。如果被积函数换成分式、带绝对值的函数、分段函数int很可能直接返回原式表示它算不出来。遇到这种情况最直接的办法是抛弃符号积分全部改用数值积分。把基函数和函数f都写成函数句柄用integral计算内积。举个例子计算(f, phi_2)原来是int(exp(x) * (x^2 - 1/3), x, -1, 1)数值形式是integral((x) exp(x) .* (x.^2 - 1/3), -1, 1)注意乘号前面用了点乘。Matlab的integral要求被积函数支持数组输入所以凡是有乘法除法幂运算的地方都要用点运算。这一步错了程序会报维度不匹配或者结果差很多。还有一个小技巧如果既想要符号推导的清晰性又想要数值计算的速度可以先用matlabFunction把符号表达式转成函数句柄再用integral。例如phi2_func matlabFunction(x^2 - 1/3, Vars, x); I integral((x) exp(x) .* phi2_func(x), -1, 1);这样既保留了基函数的符号形式又避免了int的卡顿。4.2 高次逼近结果发散或者矩阵病态很多实验指导书会让学生直接用法方程和单项式基也就是基函数取1, x, x², x³…。这种做法在次数低时没问题次数一旦超过8法方程的系数矩阵几乎奇异。一个很典型的症状是解出来的系数特别大有正有负比如系数是几百甚至上千逼近多项式在采样点上值和原函数差不多但在采样点之间剧烈振荡。这就是病态矩阵带来的数值灾难。要验证是不是病态可以查看系数矩阵的条件数。在Matlab里用cond函数% 构造单项式基对应的Gram矩阵 n 6; G zeros(n1); for i 1:n1 for j 1:n1 G(i,j) integral((x) x.^(ij-2), -1, 1); end end cond(G)当n6时条件数已经非常大了n10时会更夸张甚至接近机器精度极限。只要看到条件数在10的10次方以上就应该果断换成正交基函数。这也是实验报告里一个很好的“发现”素材用同一个函数、同一次数对比单项式基和正交基求解的误差结论一目了然。4.3 高次数逼近为什么不总是越完美最佳平方逼近从理论上是次数越高误差越小但实际中高次数逼近往往会带来端点附近的振荡这种现象在插值里更常见逼近里也会出现尤其是被逼近函数不够光滑时。原因是高次多项式本身有很强的波动能力为了压低整体平方误差它不惜在局部制造大的起伏。所以做实验时不要一味把次数加到20甚至30。一般实验题要求4到6次已经能展示方法的核心。如果想观察误差变化可以做一个次数从1到10的扫描记录均方误差和最大绝对误差。你会看到均方误差总体下降但最大绝对误差可能在某一次数后不再下降甚至反弹。这个现象值得写进报告比单纯说“次数越高越好”要真实得多。4.4 常见报错速查表报错信息常见原因解决办法未定义函数或变量‘syms’缺少符号数学工具箱改用数值实现或安装对应工具箱使用 * 时出错矩阵维度必须一致函数句柄中没使用点乘把*、/、^改成.*、./、.^积分上限或下限必须是有限标量integral的区间填入了符号量确认区间为数值或先用double转换无法找到符号默认值syms定义前调用了符号变量在最前面加syms x错误使用 polyval数据无法转换为 double传入的是符号多项式先用sym2poly或vpa转换legend不显示曲线绘图时用了fplot加hold on图例识别有问题改用plot画数值点或检查fplot句柄这些坑基本覆盖了实验七里八成以上的报错。遇到问题别急着从头改先看被积函数有没有点运算再看变量是符号还是数值大多数问题都出在这两个地方。5. 让实验报告更出彩从“会跑代码”到“讲清原理”5.1 把最佳平方逼近理解为投影我在批改类似实验时发现很多同学能把代码跑通但问他“为什么法方程长这样”就答不上来。这里有一个非常漂亮的几何解释函数可以看成无穷维空间里的向量多项式空间是这个大空间里的一个子空间最佳逼近多项式就是原函数在这个子空间上的正交投影。投影的含义是原函数减去投影后得到的残差向量应该与子空间里的所有向量都正交。换句话说残差和每一个基函数的内积都等于零。这正是法方程的来源。这个理解一旦建立很多公式就不用死记了。建议在实验报告里加一段残差正交性验证。计算% 验证残差和每个基函数内积是否接近0 residual (x) exp(x) - p_num(x); for i 1:N1 check(i) integral((x) residual(x) .* phi_func{i}(x), -1, 1); end disp(check);因为数值积分有误差check不会严格为0但应该都在1e-8量级附近。这一小段代码能让报告立刻显得有深度。5.2 插值、最佳平方逼近、离散最小二乘三者对比整理成一张表放在实验报告里可以直观说明这个实验在整个数值方法体系里的位置方法误差度量是否要求经过已知点典型问题多项式插值节点处误差为0是高次振荡、对噪声敏感最佳平方逼近连续区间上的积分误差最小否端点可能仍有偏差离散最小二乘给定离散点误差平方和最小否依赖点的选取和分布这张表的结论是没有一种方法是万能的。最佳平方逼近适合函数本身连续、可以计算积分的情形离散最小二乘适合只有测量数据、函数解析式未知的情形插值适合必须严格通过某些关键点的情形。实验七学的是前两种的结合。5.3 进一步扩展加权逼近与分段逼近如果实验做完还有余力可以把“等权内积”升级成“加权内积”。定义权函数ρ(x)内积变成(u, v) ∫ ρ(x) u(x)v(x) dx选不同的权函数逼近的效果会发生变化。比如在一些端点影响较大的问题里适当增大权函数在端点的权重可以让逼近函数在端点更贴近原函数。这个扩展概念上和最小二乘的加权形式一致代码上只需要多乘一个ρ(x)。另一个扩展是分段逼近。与其用一个高次多项式逼近整个区间不如把区间切成若干段每一段用一个低次多项式逼近再保证段与段之间的连续性。这就是样条逼近的雏形也是工程中处理复杂函数的常用思路。实验七只需要做全区间的多项式逼近但知道这个扩展方向能让你对后续课程里的样条方法更有感觉。最后再分享一个我做这个实验的体会不要只把代码跑通就结束一定要手动改几次f(x)和次数N观察误差变化曲线。比如分别逼近e^x、sin(πx)、|x|你会发现光滑函数收敛快带拐点的函数收敛慢。这种直观感受比背十遍理论都管用。实验七的真正价值不只是让你学会调用几个Matlab函数而是让你建立起“函数空间也有几何”的直觉。有了这个直觉后面再看傅里叶展开、有限元方法都会顺很多。