6.1 数值计算

6.1 数值计算

数值计算两把刀：解方程用 DifferentialEquations，看结果用 Plots。本节用同一个物理问题（弹簧振子）把两把刀都开刃，顺带过一遍线性代数标准库。

十分钟解一个 ODE：阻尼弹簧

物理方程：质量块在弹簧与阻尼下振动，位置对时间的二阶导等于负刚度乘位置减阻尼乘速度。转成一阶方程组：

] add DifferentialEquations Plots

using DifferentialEquations, Plots

# 状态 u = [位置, 速度]
function oscillator!(du, u, p, t)
    k, c = p                  # 刚度与阻尼
    du[1] = u[2]
    du[2] = -k * u[1] - c * u[2]
end

u0 = [1.0, 0.0]               # 初位置 1，初速度 0
tspan = (0.0, 10.0)
prob = ODEProblem(oscillator!, u0, tspan, [4.0, 0.5])
sol = solve(prob)

plot(sol, idxs = (1, 2), title = "相图", xlabel = "位置", ylabel = "速度")
​

注意函数名末尾的 !——Julia 惯例：修改参数的函数加感叹号。du 由调用方预分配，这是高性能数值库的通用接口风格。

换个 stiff 问题（方程组快慢尺度悬殊），一行换求解器：

sol2 = solve(prob, Tsit5())        # 显式 Runge-Kutta
sol3 = solve(prob, Rodas5())       # 隐式，刚性系统更稳
​

求解器家族速览

从方程到图像的标准链路

线性代数：标准库就够强

不需要额外安装，LinearAlgebra 是标准库：

using LinearAlgebra

A = [2.0 1.0; 1.0 3.0]
b = [3.0, 5.0]
A \ b                     # 求解 Ax=b，左除是惯用法
eigen(A)                  # 特征值与特征向量
svd(A)                    # 奇异值分解
det(A), norm(b), A' * A   # 行列式、范数、转置乘
​

A \ b 背后是 LU 分解，比手动求逆再乘（inv(A) * b）又快又稳——永远用左除。

💡 关键直觉：解线性方程时不要写 inv(A) * b。数学等价，数值上更慢、舍入误差更大，A \ b 是唯一正确姿势。

Plots 出图的日常三招

x = range(0, 2π, length = 200)
plot(x, sin.(x), label = "sin")                 # 线图
plot!(x, (x -> sin(2x)).(x), label = "sin2x")   # 叠加加感叹号

histogram(randn(10^5), bins = 80)               # 直方图
scatter(randn(200), randn(200), markersize = 3) # 散点
​

案例：调参观察阻尼的三个区间

解出方程只是起点，数值实验的乐趣在扫参数。把阻尼系数从 0.05 拉到 2.5，观察解的形态从"振荡衰减"跨到"过阻尼不振荡"。背景是这条分界线（临界阻尼）在理论上是二倍根号刚度除质量，数值上可以直接抓出来：

using DifferentialEquations

function settle_time(k, c)
    # 返回位置首次进入并保持在 ±0.05 内的时刻
    prob = ODEProblem(oscillator!, [1.0, 0.0], (0.0, 60.0), [k, c])
    sol = solve(prob)
    for t in sol.t
        abs(sol(t)[1]) > 0.05 && return t       # 还在带外，继续
    end
    0.0    # 从未出带（本例不会走到）
end

# 粗扫：欠阻尼、临界、过阻尼三段的稳定时间
[round(settle_time(4.0, c); digits=2) for c in (0.3, 0.4, 0.8, 1.6, 2.5)]
# 临界阻尼理论值 c = 2√k = 4 附近，稳定时间最短
​

操作后把结果排开看：阻尼太小则振荡拖长稳定时间，阻尼太大则像陷进泥潭缓慢爬回，两头都慢、中间最快——这正是工程上"临界阻尼设计"的数值证据。变式：把 settle_time 套上 5.2 的并行映射，对 k–c 平面做二维扫描，一张"稳定时间热力图"就是一份可写进报告的完整数值实验。这个案例同时示范了 sol(t) 可像函数一样插值求值这一便利特性。

线性代数的边界与坑

三个实测过的坑。其一，A \ b 对方阵走 LU、对超定方程组自动转最小二乘，语义随形状变，读旧代码时先确认 A 的形状再解释结果。其二，特殊结构要用特殊求解：对称正定矩阵加 Symmetric(A) 包装后 \ 走乔姆斯基分解，规模大时快一倍且数值更稳。其三，大矩阵乘法别自己写三重循环——A * B 底层是多线程 BLAS，第 5 章的并行知识它内部已经用上了；想控制 BLAS 线程数用 BLAS.set_num_threads(n)。这些细节的共同逻辑是：标准库把数值最优实践包好了，你的任务是选对入口而不是重写内核。

常见报错与读法

ODE 求解的报错集中在三类。dt too small：求解器把步长缩到下限仍不收敛，多半是刚性问题配了显式求解器，换 Rodas5 或 QNDF 通常立刻好转。u modified during integration：右端函数直接改了 u 或 t，违反接口约定，检查函数体里有没有对入参的赋值。收敛慢但不报错：结果看起来"对但特别慢"，多为容差设得过紧，solve(prob; abstol=1e-6, reltol=1e-6) 是科学计算的常用起点，盲目设 1e-12 会把求解时间放大百倍而无精度收益。线性代数侧的高频报错是奇异矩阵的 SingularException，先检查方程是否欠定、数据是否全为常数列，再考虑正则化。排查数值问题的总原则：先怀疑问题形态与求解器不匹配，再怀疑自己代码错，最后才怀疑库。

本节要点回顾

• ODE 三件套：右侧函数（带 !）、ODEProblem、solve，求解器一行可换；
• 刚性系统换隐式求解器，默认 Tsit5 搞不定震荡或龟速时先怀疑刚性；
• A \ b 解方程，永远不用 inv；
• Plots 增量修改用 plot!，检查永远先于出图。