📘 第 14 章：高等数学应用

这一章在干嘛？

这是本书的「工具箱章节」——把工程师日常会碰到的数学问题（统计、集合、多项式、复数、线性方程组、符号运算、积分微分）一次性铺开。学完你就能用 MATLAB 跑数据拟合、解电路节点方程、画 RC 电路响应、做数值积分，给老板/导师交差。

14.1 统计函数

14.2 集合运算

14.3 多项式、曲线拟合与插值

14.4 复数

14.5 矩阵性质与线性方程组

14.6 符号数学

14.7 数值积分与微分

14.1 统计函数

工程里拿到的数据十有八九是带噪声的，「这一批数大概是什么、波动有多大、有没有异常点」就是统计函数要回答的问题。

14.1.1 一行命令看全家：`max

max 和 min 是最朴素的统计函数。它们也返回极值的索引（多个极值时返回第一个）：

>> x = [9 10 10 9 8 7 3 10 9 8 5 10];
>> [maxval, maxind] = max(x)
maxval = 10
maxind = 2                % 第一个 10 的位置，不是最后一个

矩阵按列操作，要按行请加第二个参数（dim = 2）：

>> mat = [8 9 3; 10 2 3; 6 10 9];
>> mean(mat)              % 默认按列：每列的均值
ans = 8.0000  7.0000  5.0000
>> mean(mat, 2)           % 按行：每行的均值
ans = 6.6667
     5.0000
     8.3333

14.1.2 均值家族

算术平均（arithmetic mean，平常说的「平均」）：

$\bar{x}=\frac{{\sum }_{i=1}^n{x}_i}{n}$

>> mean([9 10 10 9 8 7 3 10 9 8 5 10])
ans = 8.1667

异常点（outlier）会拉偏均值。下面这组数除了中间那个 100 都在 3~10 之间，但均值被拉到 15.23，比任何「正常」值都大：

>> xwithbig = [9 10 10 9 8 100 7 3 10 9 8 5 10];
>> mean(xwithbig)
ans = 15.2308

工程里遇到这种情况，常见做法是先去掉最大最小值再求平均：

>> clean = xwithbig(xwithbig ~= min(xwithbig) & xwithbig ~= max(xwithbig));
>> mean(clean)
ans = 8.5455

调和平均（harmonic mean）和几何平均（geometric mean）用得少但在特定领域很有用 —— 调和平均适合「速率平均」（比如数据下载速率），几何平均适合「增长率平均」（比如投资年化收益）：

$H=\frac{n}{\sum \frac{1}{x_i}},G=\sqrt[n]{x_1{x}_2\cdots x_n}$

>> x = [9 10 10 9 8 7 3 10 9 8 5 10];
>> harmhand = @(x) length(x) / sum(1./x);
>> harmhand(x)
ans = 7.2310
>> geomhand = @(x) nthroot(prod(x), length(x));
>> geomhand(x)
ans = 7.7775

14.1.3 方差与标准差：数据「散」到什么程度

方差（variance）刻画数据围绕均值的离散程度：

$\text{var}=\frac{{\sum }_{i=1}^n(x_i-\bar{x}{)}^2}{n-1}$

>> x = [8 7 5 4 6];
>> var(x)
ans = 2.5000

标准差（standard deviation）就是方差的平方根：

>> std(x)
ans = 1.5811           % = sqrt(2.5)

下面两组数据均值相同（都是 9.5）但离散度天差地别：

>> x1 = [9 10 9.4 9.6]; std(x1)            % 紧凑数据
ans = 0.4163
>> x2 = [2 17 -1.5 20.5]; std(x2)          % 跳跃数据
ans = 10.8704

工程实战：质量控制

如果你是产线工程师，两台设备测得同一直径 16 mm 工件，第一台 std = 0.02 mm（精度高），第二台 std = 0.5 mm（精度差）——只看均值无法判断设备好坏，必须配合标准差看离散度。

14.1.4 众数与中位数

众数（mode）= 出现次数最多的值；多个并列时取最小的那个：

>> mode([9 10 10 9 8 7 3 10 9 8 5 10])
ans = 10                   % 出现了 4 次

中位数（median）= 排序后中间那个值（奇数长）或中间两个的平均（偶数长）：

>> median([1 4 5 9 12])    % 5 个数中间 = 第 3 个
ans = 5
>> median([1 4 5 9 12 33]) % 6 个数中间 = 第 3、4 个的平均 = (5+9)/2
ans = 7

常见坑：均值 vs 中位数在偏态分布上的差异

数据有异常值时，中位数比均值鲁棒得多。例：「大多数员工月薪 1 万，老板月薪 100 万，公司平均月薪 ≈ 9 万」——这种时候报中位数比报均值更反映「真实」水平。

---

14.2 集合运算

把两个向量当作两个集合做运算 —— 工程师常常要回答「两组传感器的告警 ID 里哪些是共同的」「A 班和 B 班的零件编号有哪些不同」这类问题。

>> v1 = 6:-1:2            % v1 = [6 5 4 3 2]
>> v2 = 1:2:7             % v2 = [1 3 5 7]

函数	含义	例子
union(a, b)	并并	[1 2 3 4 5 6 7]
intersect(a, b)	交并	[3 5]
setdiff(a, b)	差并（a 中去掉 b）	[2 4 6]
setxor(a, b)	对称差（只在 a 或只在 b）	[1 2 4 6 7]
unique(v)	去重	[1 2 3 4 5 6]
ismember(a, b)	逻辑向量，a 的每个元素是否在 b 中	[0 1 0 1 0]
issorted(v)	是否升序排列	0 或 1

R2013a 之前返回结果默认升序；之后用 'stable' 保留原始顺序：

>> union(v1, v2, 'stable')
ans = 6  5  4  3  2  1  7   % 先按 v1 的顺序，最后是 v2 中独有的

intersect 还支持返回索引向量，索引回原向量能拿到原始值：

>> [out, i1, i2] = intersect(v1, v2);
>> v1(i1)               % 与 v2 共同的元素在 v1 中的位置
ans = 3  5

用法建议

想做「A 班学生集合有没有完全包含 B 班」这种判定，组合两个 ismember 比手写循环干净得多：

all(ismember(B, A))    % B 中所有元素都在 A 中吗？

---

14.3 多项式、曲线拟合与插值

「数据是采样得到的散点，但我想知道任意时刻的值是多少」—— 这是工程里最常见的数学问题。多项式拟合 + 插值就是答案。

14.3.1 多项式表示法

MATLAB 用系数向量（按降幂排列）表示多项式：

$p(x)=a_n{x}^n+a_{n-1}{x}^{n-1}+\cdots +a_{1}x+a_0⟺[a_n,a_{n-1},\dots ,a_1,a_0]$

>> p = [1 2 -4 3]            % 代表 x^3 + 2x^2 - 4x + 3

两个核心函数：

• roots(p)：求多项式等于 0 的根
• polyval(p, x)：把多项式在 x 处求值

>> roots([4 -2 -8 3])        % 解 4x^3 - 2x^2 - 8x + 3 = 0
ans =
   -1.3660
    1.5000
    0.3660
>> polyval([-2 1 4], 3)      % 求 -2x^2 + x + 4 在 x=3 的值
ans = -11                    % = -2*9 + 3 + 4

常见坑：缺失项要补 0

$2x^4-x^2+5$ 没有 $x^3$ 项和 $x$ 项，要写成 [2 0 -1 0 5]。少一个 0，阶数就错了。

14.3.2 曲线拟合：polyfit + polyval

温度采样案例（工程实例）：一个下午每整点测一次温度，5 个采样点用直线（1阶）来拟合，效果是几乎平的横线 —— 因为这条线代表「下午整体趋势」，没有反映 4 点钟的高温：

>> x = 2:6;                             % 2~6 点钟
>> y = [65 67 72 71 63];                % 对应温度
>> polyfit(x, y, 1)                     % 拟合 1 阶（直线）
ans = 0.0000  67.6000                  % y ≈ 67.6 的水平线

换成 2 阶（二次曲线）效果立刻好了：

>> coeffs = polyfit(x, y, 2)            % 拟合 2 阶（抛物线）
coeffs = -1.8571  14.8571  41.6000
>> curve = polyval(coeffs, x);          % 在每个 x 处求拟合曲线值
>> plot(x, y, 'ro', x, curve); xlabel('Time'); ylabel('Temp');

工程版（画得更密更平滑）：拿一个含 100 点的细网格 linspace(2,6,100) 来画曲线，数据点保持原来的 5 个：

% polytemp.m
x  = 2:6;
y  = [65 67 72 71 63];
coefs = polyfit(x, y, 2);
morex = linspace(min(x), max(x), 100);   % 100 个密集采样点画曲线
plot(x, y, 'ro', morex, polyval(coefs, morex));
xlabel('Time'); ylabel('Temp'); axis([1 7 60 75]);
title('Temperatures one afternoon');

14.3.3 阶数对比：欠拟合、过拟合与「刚刚好」

阶数太低 → 模型太朴素（欠拟合）；阶数太高 → 模型把噪声也学进去了（过拟合）。下面把 1、2、3 阶并列画出来，对比就一目了然：

% polytempsubplot.m
x = 2:6;
y = [65 67 72 71 63];
morex = linspace(min(x), max(x));
for pd = 1:3
    coefs = polyfit(x, y, pd);
    curve = polyval(coefs, morex);
    subplot(1,3,pd);
    plot(x, y, 'ro', morex, curve);
    title(sprintf('Degree %d', pd));
    axis([1 7 60 75]);
end

14.3.4 插值与外推

插值（interpolation）：求采样点之间的值。外推（extrapolation）：求采样点之外的值。

>> polyval(coeffs, 2.5)       % 插值：2:30 时的温度
ans = 67.1357
>> polyval(coeffs, 1.0)       % 外推：1:00 时的温度
ans = 54.6000                 % ⚠️ 外推不可靠，离样本越远越不准

常见坑：外推危险

用一个二次曲线去预测「明天的温度」，相当于让模型在完全没有数据约束的区域发疯。永远把外推视为「仅供参考」，真要做预测用线性回归 + 置信区间或时间序列模型。

14.3.5 工程实例：噪声数据的二次/三次拟合对比

下面两个图是 12 个随机采样点的二阶和三阶拟合曲线对比。点数多时，二阶已经能很好地反映整体趋势；三阶开始把噪声也拟合进去。

---

14.4 复数

电子工程师绕不开复数 —— 交流电路里的阻抗、信号处理里的频域、控制理论里的传递函数，全都是复数。MATLAB 里 i 和 j 都是内置的 $\sqrt{-1}$。

常见坑：循环变量 i 覆盖虚数单位

>> for i = 1:5; end
>> 3*i             % 这时 i=6，结果是 18，不再是复数

永远用 1i 或 1j（永远返回复数 0+1.0000i），不要图省事只用 i/j。

14.4.1 复数创建与基本运算

>> z1 = 4 + 2i              % 标准写法
>> z2 = sqrt(-5)            % MATLAB 自动返回纯虚数
z2 = 0.0000 + 2.2361i
>> z3 = complex(3, -3)      % 用函数构造
>> z4 = 2 + 3j              % 工程习惯写法，结果与 2+3i 相同

实部、虚部、共轭、模长：

>> real(z1), imag(z1)       % 实部、虚部
ans = 4
ans = 2
>> conj(z1)                 % 共轭：a - bi
ans = 4.0000 - 2.0000i
>> abs(z1)                  % 模长：sqrt(a^2 + b^2)
ans = 4.4721

打印陷阱：fprintf('%f\n', z1) 默认只打实部，不会出错但丢一半信息。要打完整得：

>> fprintf('%f + %fi\n', real(z1), imag(z1))
4.000000 + 2.000000i

14.4.2 极坐标形式

复数 $z=a+bi$ 既可以看作复平面上的点 $(a,b)$，也可以看作长度为 $r$、角度 $\theta$ 的向量：

$z=r{e}^{i\theta }=r\cos \theta +ir\sin \theta$

MATLAB 提供直接转换：

>> z = 3 + 4i;
>> r = abs(z);              % 5
>> theta = angle(z);        % 0.9273 rad ≈ 53.13°
>> r * exp(1i*theta)        % 反算回 a+bi
ans = 3.0000 + 4.0000i

下面这个图把一个复数 $z=3+4i$ 画在复平面上（横轴实部，纵轴虚部）：

>> plot(z, '*', 'MarkerSize', 12);
>> xlabel('Real part'); ylabel('Imaginary part');
>> title('Complex number');

14.4.3 复数多项式

多项式系数向量也可以是复数：

>> roots([1 1 -3+2i])       % 解 z^2 + z - 3 + 2i = 0
ans =
  -2.3796 + 0.5320i
   1.3796 - 0.5320i
>> cp = [1 1 -3+2i];
>> polyval(cp, 3)            % 求 z=3 处的值
ans = 9.0000 + 2.0000i

---

14.5 矩阵性质与线性方程组

这是工程数学最硬核的工具 —— 电路节点分析、力学平衡、控制理论状态方程，最底下都是「解一个 $A\mathbf{x}=\mathbf{b}$ 的线性方程组」。

14.5.1 矩阵家族速查

概念	定义	例子
方阵（square）	行数 = 列数	$3\times 3$
对角阵（diagonal）	非主对角元素全 0	diag([4 9 5])
单位阵（identity）	对角线都是 1 的对角阵	eye(3)
三角阵（triangular）	上三角或下三角之外的元素全 0	triu(A) / tril(A)
三对角阵（tridiagonal）	只有主对角 + 上下相邻对角非 0	稀疏存储常用
对称阵（symmetric）	$A=A^T$	力学刚度矩阵
迹（trace）	主对角元素之和	trace(A)

MATLAB 自带判定函数（都是 R2014b 引入的）：

>> A = [1 2 3; 2 5 4; 3 4 6];
>> isdiag(A)        % ans = 0
>> issymmetric(A)   % ans = 1
>> istril(A), istriu(A)

eye 的命名因为发音像字母 i（=identity）—— 故意避开和虚数 i 冲突：

>> eye(4)           % 4x4 单位阵

14.5.2 解线性方程组：两种姿势

把 $A\mathbf{x}=\mathbf{b}$ 这种线性方程组写成矩阵形式后，MATLAB 提供两种等价解法：

姿势 1：直接用反斜杠 \（推荐）

>> A = [4 -2 1; 1 1 5; -2 3 -1];
>> b = [7; 10; 2];
>> x = A\b            % 一行搞定
x =
    3.0244
    2.9512
    0.8049

姿势 2：先求逆再相乘（教学用，生产环境不推荐）

>> x = inv(A) * b     % 结果一样，但慢 + 数值不稳

为什么优先用 A\b 而不是 inv(A)*b？

inv 要算完整个矩阵的逆，O(n³) 还容易在病态矩阵上累积误差；\（backslash）在底层用 LU 分解或 QR 分解，速度快、稳定性高。教材 14.5 节也明确演示了两种做法结果一致。

14.5.3 工程实战：电路节点电压

最常见的工程线性方程是「电路节点电压分析」。下面三方程组求三个节点电压：

$$\{\begin{array}{l}{V}_1=5\\ -6V_1+10V_2-3V_3=0\\ -V_2+51V_3=0\end{array}$$

写成矩阵形式 $A\mathbf{V}=\mathbf{b}$：

>> A = [1 0 0; -6 10 -3; 0 -1 51];
>> b = [5; 0; 0];
>> V = A\b
V =
    5.0000
    3.9286
    0.0770

14.5.4 rref 与增广矩阵

把系数矩阵 $A$ 和常数向量 $\mathbf{b}$ 横向拼起来叫「增广矩阵」[A | b]。用 rref 把它化成「行最简形」后，最后一列就是解：

>> ab = [A b];
>> rref(ab)
ans =
    1     0     0     5.0000
    0     1     0     3.9286
    0     0     1     0.0770

也能反向用 —— 用 rref 把 [A | I] 化成 [I | A^{-1}] 求逆矩阵：

>> A = [1 3 0; 2 1 3; 4 2 3];
>> rref([A eye(size(A))])   % 左侧变 I 时，右侧就是 A^{-1}

---

14.6 符号数学

「让 MATLAB 像 Mathematica 一样把表达式当符号处理」—— 不求具体数值，而是保留 $x$、$y$、$a$、$b$ 这些符号，得到精确的解析解。需要 Symbolic Math Toolbox（额外安装的）。

14.6.1 符号变量与表达式

>> syms a b x y z           % 一次声明多个符号变量
>> myexpr = a*x^2 + b*x     % 仍以符号形式存储
myexpr = a*x^2 + b*x

符号运算自动合并同类项：

>> z = sym('z');
>> z^3 + 2*z^3              % 自动合并
ans = 3*z^3

poly2sym 把系数向量还原成多项式字符串，sym2poly 反过来：

>> poly2sym([1 2 -4 3])
ans = x^3 + 2*x^2 - 4*x + 3
>> sym2poly(ans)
ans = 1 2 -4 3

14.6.2 化简：simplify / collect / expand / factor

函数	用途
simplify(expr)	通用化简（三角恒等式等）
collect(expr)	合并同类项
expand(expr)	展开（乘出来）
factor(expr)	分解因式（不可分解则原样返回）

>> simplify(cos(x)^2 + sin(x)^2)
ans = 1                              % 三角恒等式被识别
>> expand((x+2)*(x-1))
ans = x^2 + x - 2
>> factor(ans)
ans = (x+2)*(x-1)

14.6.3 代入数值：subs

>> myexp = x^3 + 3*x^2 - 2;
>> subs(myexp, 3)                   % 把 x 替换为 3
ans = 52
>> varexp = a*x^2 + b*x;
>> subs(varexp, 'a', 3)             % 指定替换 a
ans = 3*x^2 + b*x

14.6.4 精确分数运算

默认 double 运算下 1/3 + 1/2 = 0.8333，但符号运算下保留精确分数：

>> sym(1/3 + 1/2)
ans = 5/6                           % 精确分数
>> double(ans)
ans = 0.8333                        % 转回 double

用 numden 提取分子分母：

>> [n, d] = numden(sym(1/3 + 1/2))
n = 5
d = 6

14.6.5 解方程：solve

单变量方程：自动把表达式设为 0 求根；多个解全部返回：

>> solve('2*x^2 + x = 6')
ans =
    -2
    3/2                       % 符号解保留为分数
>> double(ans)
ans = -2.0000  1.5000

指定求解变量：多变量表达式默认解 $x$，可用第二个参数指定：

>> solve('a*x^2 + b*x', 'b')   % 求 b
ans = -a*x

多元方程组：返回结构体，.x、.y、.z 字段是各变量解：

>> S = solve('4*x-2*y+z=7', 'x+y+5*z=10', '-2*x+3*y-z=2');
>> S.x, S.y, S.z
x = 124/41
y = 121/41
z = 33/41
>> double([S.x S.y S.z])
ans = 3.0244  2.9512  0.8049

14.6.6 符号画图：ezplot

>> ezplot('x^3 + 3*x^2 - 2')     % 自动选 [-2π, 2π]

14.6.7 细菌分裂：代数求解真实问题

教材练习 36 给出一个对数形式的细菌繁殖公式：

$\log N=\log N_0+\frac{t}{T}\log 2$

已知 $N_0=10^2$、$N=10^8$、$t=8$ h，求倍增时间 $T$：

>> syms T
>> T = solve('log10(10^8) = log10(10^2) + 8/T * log10(2)', T);
>> double(T)
ans = 1.3333                  % ≈ 80 分钟，E. coli 大致这个数量级

---

14.7 数值积分与微分

理论上的微积分有公式，但工程里的 $f(x)$ 往往是采样数据或者无法解析积分的复杂函数 —— 这时候就要数值积分。

14.7.1 梯形法：最朴素的数值积分

${\int }_a^{b}f(x)dx\approx (b-a)\frac{f(a)+f(b)}{2}$

% trapint.m — 一段手写梯形法
function int = trapint(fnh, a, b)
    int = (b-a) * (fnh(a) + fnh(b)) / 2;
end
 
>> f = @(x) 3*x.^2 - 1;          % 被积函数
>> trapint(f, 2, 4)              % 一段梯形近似
ans = 58

内置版：trapz(x, y) 接受 x 和 y 向量：

>> x = 2:4;
>> y = f(x);
>> trapz(x, y)
ans = 58

提高精度：把 $[a,b]$ 切成 $n$ 段，分别用梯形法再求和。x 用等距采样：

>> trapz(2:0.01:4, f(2:0.01:4))   % 100 段
ans = 54.0014                      % 接近真实值 54

真实值 ${\int }_2^4(3x^2-1)dx=[x^3-x{]}_2^4=(64-4)-(8-2)=54$。

现代推荐

R2015b 起更推荐用 integral(f, a, b)（自适应 Simpson/Romberg），速度、精度都比 trapz 强：

>> integral(f, 2, 4)
ans = 54.0000

14.7.2 符号积分：int

如果装了 Symbolic Math Toolbox，直接求原函数：

>> syms x
>> int(3*x^2 - 1)               % 不定积分
ans = x^3 - x
>> int(3*x^2 - 1, 2, 4)         % 定积分
ans = 54

14.7.3 多项式积分与微分：polyint / polyder

对系数向量表示的多项式，积分 = 「每项降一次幂，除以新幂」：

$\int (a_n{x}^n+\dots )dx=\frac{a_n}{n+1}{x}^{n+1}+\dots$

>> origp = [3 4 -4];            % 3x^2 + 4x - 4
>> intp = polyint(origp)
intp = 1 2 -4 0                 % x^3 + 2x^2 - 4x + 0

导数则反过来 ——「每项降一次幂，乘以原幂」：

>> origp = [1 2 -4 3];          % x^3 + 2x^2 - 4x + 3
>> diffp = polyder(origp)
diffp = 3 4 -4                  % 3x^2 + 4x - 4
>> polyval(diffp, 1:3)          % 在 x=1,2,3 处求值
ans = 3 16 35

14.7.4 数值微分：diff

对采样向量，导数 ≈ 相邻差分：

$f^{\prime }(x_i)\approx \frac{y_{i+1}-y_i}{x_{i+1}-x_i}$

>> f = @(x) x.^3 + 2*x.^2 - 4*x + 3;
>> x = 0.5:0.5:3.5;             % 步长 0.5
>> y = f(x);
>> diff(y) ./ diff(x)           % 中点差分
ans = 3.2500  16.2500  35.2500  % 接近 polyder 解析解 [3 16 35]

14.7.5 ODE 求解：ode45

常微分方程（ODE）是工程里最常见的动态系统描述 —— RC 电路、机械振动、化学反应动力学。

RC 电路的一阶 ODE：

$\frac{dV}{dt}=-\frac{V}{RC}+\frac{V_{in}}{RC}$

其中 $V$ 是电容电压，$V_{in}$ 是输入电压（阶跃信号）。

% RC 电路阶跃响应：R=1kΩ, C=1μF, V_in = 5V
R = 1e3;  C = 1e-6;  V_in = 5;
% 状态方程：dV/dt = (V_in - V) / (R*C)
f = @(t, V) (V_in - V) / (R*C);
 
tspan = [0 0.01];              % 仿真 10 ms（约 5 倍 τ=RC）
[V, t] = ode45(f, tspan, 0);   % 初始电压 V(0)=0
 
plot(t*1000, V);
xlabel('时间 (ms)'); ylabel('V_C (V)');
title('RC 电路阶跃响应');
grid on;

运行后会看到一条典型的指数上升曲线：0→5 V，时间常数 $\tau =RC=1$ ms，约 5 ms 达到稳态。

ode 系列函数

• ode45：通用、自适应 Runge-Kutta(4,5)，首选
• ode23：低阶版本，适合刚性问题
• ode15s：刚性方程专用（化学动力学常用）
• 命令窗口敲 odeexamples 看更多官方示例

14.7.6 用 fzero 求方程的根

数值求根（解 $f(x)=0$）：

>> f = @(x) x.^3 - 2*x - 5;
>> r = fzero(f, 2)              % 在 x=2 附近搜索
r = 2.0946

---

实战综合案例：RC 电路放电曲线拟合

假设你测了 RC 电路放电的 10 个采样点，要拟合理论指数曲线 $V(t)=V_0{e}^{-t/RC}$：

% 已知 R = 1kΩ, C = 100μF, τ = 0.1s, 采样 0~0.5s
t_data = 0:0.05:0.5;
V_data = [5.00 4.10 3.30 2.71 2.25 1.85 1.50 1.23 1.02 0.84];
 
% 用对数线性化 → log(V) = log(V0) - t/RC → polyfit 1 阶
coefs = polyfit(t_data, log(V_data), 1);  % y = a*t + b
V0_fit = exp(coefs(2));
RC_fit = -1 / coefs(1);
 
fprintf('拟合 V0 = %.3f V, τ = RC = %.4f s\n', V0_fit, RC_fit);

输出大致是 V0 = 5.00 V, τ = RC = 0.1002 s，误差来自读数与采样频率。这是把理论模型套到实测数据上的标准流程。

---

本章通关标准：

1. 能用 mean / var / std / median / mode 描述一组数据的中心和离散程度，并能识别异常点的影响；
2. 能用 polyfit / polyval 对采样点做多项式拟合，区分插值与外推的安全边界；
3. 能用 A\b 解线性方程组、用 syms + solve 解代数方程、用 integral / trapz 求数值积分、用 ode45 仿真一阶 ODE。

1. 同样一组数据，均值和标准差各描述什么？为什么单看均值不够？

均值描述「整体水平」（数据中心），标准差描述「数据围绕中心的离散程度」。均值相同的两组数据标准差可能差几十倍 —— 例如均值都是 9.5，一组是 [9 10 9.4 9.6]（std ≈ 0.4），另一组是 [2 17 -1.5 20.5]（std ≈ 10.9）。所以质量控制、实验测量、金融风控都必须同时报均值 + 标准差。

2. polyfit(x, y, 2) 用 5 个点拟合二次曲线为什么比拟合三次曲线更稳？

三次多项式有 4 个待定系数 $(a_3,a_2,a_1,a_0)$，5 个数据点提供 5 个方程 —— 1 个自由度留给拟合误差「抖动」；如果采样噪声较大，三次曲线会「追」着噪声走，出现 overfitting 的边界异常扭动（见教材图 14.3 的右侧子图）。二次曲线只有 3 个系数，约束更强、更平滑、更适合这种「采样点少、噪声大」的小数据集。原则：先从低阶试，再考虑加阶。

3. 解 $A\mathbf{x}=\mathbf{b}$，用 A\b 比 inv(A)*b 好在哪？

inv(A) 要先算出整个逆矩阵（O(n³) 计算量 + 数值误差积累），而 \ 在底层根据矩阵性质自动选用 LU 分解、QR 分解等最优算法，对病态矩阵也更稳。教材明确演示两种写法结果一致，但工程实践永远用 A\b。

4. 写一段 5 行代码：用 ode45 仿真 RC 一阶电路的阶跃响应（$R=1k\Omega$，$C=1\mu F$，$V_{in}=5V$）。

f = @(t, V) (5 - V) / (1e3 * 1e-6);   % dV/dt
[V, t] = ode45(f, [0 5e-3], 0);      % 仿真 5 ms 起始 0V
plot(t*1000, V); xlabel('ms'); ylabel('V_C (V)'); grid on;

应该看到从 0V 指数上升到约 4.97 V 的曲线，时间常数约 1 ms。

---

拓展：常用 ODE 模板

下面是一段能直接复用的「任意一阶 ODE 仿真模板」，把状态方程、初始条件、仿真时长抽出来：

% general_ode.m
% 解 dx/dt = f(t, x), 初始 x0, 仿真时长 [t0 tf]
f   = @(t, x) ...;              % 你的状态方程
t0  = 0;   x0 = ...;            % 初始条件
tspan = [t0 tf];                % 仿真区间
 
opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);   % 提高精度
[t, x] = ode45(f, tspan, x0, opts);
 
plot(t, x);
title('ODE 时域响应'); xlabel('t'); ylabel('x'); grid on;

把它从一阶扩展到高阶：把多个状态变量打包成向量 $\mathbf{x}=[x_1;x_2;\dots ]$，状态方程写成向量形式 $\dot{\mathbf{x}}=f(t,\mathbf{x})$，初始条件也用列向量，ode45 完全照搬。

用状态空间模型做控制系统仿真

控制系统传递函数 $G(s)$ 转成状态空间形式 $(\dot{\mathbf{x}}=A\mathbf{x}+Bu,y=C\mathbf{x}+Du)$ 后，就可以直接用 ode45 仿真单位阶跃、单位脉冲、任意输入 $u(t)$ 的响应。这是控制系统工程师最常用的「离线仿真」套路，比 Simulink 灵活，又比手写四阶 Runge-Kutta 省心。