工科数学完全征服第4篇：数值分析

当工程问题无法用解析方法求解，或者需要对复杂系统进行数值仿真时，就需要用到数值分析(Numerical Analysis)。晶体管非线性特性分析、复杂电路网络的暂态响应、有限元分析等，全都建立在数值分析之上。本文将从原理出发，把核心数值方法一路讲到完整的Python实现。

1. 误差分析 (Error Analysis)

1.1 误差的种类

数值计算中不存在"完美解"。理解误差的种类，是可靠数值计算的起点。

舍入误差(Round-off Error): 这是由于计算机用有限位数表示实数而产生的。

IEEE 754 双精度(double)：约15-16位有效数字

$1/3 = 0.333333\ldots \to 0.3333333333333333 \text{ (有限表示)}$

截断误差(Truncation Error): 这是在用有限项近似无穷级数时产生的。

$e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots \approx 1 + x + \frac{x^2}{2!} \text{ (截断高次项)}$

绝对误差(Absolute Error): $E_a = |x_{\text{true}} - x_{\text{approx}}|$

相对误差(Relative Error): $E_r = \frac{|x_{\text{true}} - x_{\text{approx}}|}{|x_{\text{true}}|}$

百分比误差： $\epsilon_r = E_r \times 100\%$

有效数字(Significant Figures): 若 $\epsilon_r < 0.5\%$ ，则可保证有效数字达到 $d$ 位： $d = 2 - \log_{10}(E_r \times 100)$

1.2 数值稳定性 (Numerical Stability)

误差不应在计算过程中被放大。

良态(Well-conditioned): 输入误差不会在输出中被显著放大

病态(Ill-conditioned): 输入的微小变化会导致输出发生巨大变化

条件数(Condition Number): $\kappa(A) = \|A\| \cdot \|A^{-1}\|$

在线性系统 $Ax = b$ 中，若 $\kappa(A)$ 很大，数值解就会不稳定。

import numpy as np

# 希尔伯特矩阵(声名狼藉的病态矩阵)
n = 6
H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
kappa = np.linalg.cond(H)
print(f"希尔伯特矩阵 ({n}x{n}) 条件数: {kappa:.2e}")

# 简单矩阵
A = np.array([[2, 1], [1, 3]], dtype=float)
print(f"简单矩阵条件数: {np.linalg.cond(A):.2f}")

# 误差传播演示
b_true = np.array([3, 4], dtype=float)
x_true = np.linalg.solve(A, b_true)

# 加入微小误差
b_perturbed = b_true + np.array([0.01, 0.0])
x_perturbed = np.linalg.solve(A, b_perturbed)

print(f"\n输入误差: {np.linalg.norm(b_perturbed - b_true):.4f}")
print(f"输出误差: {np.linalg.norm(x_perturbed - x_true):.4f}")

1.3 大O记号与收敛阶数

收敛阶数(Order of Convergence):

$E_{n+1} \leq C \cdot E_n^p$

$p = 1$ : 一阶(线性)收敛
$p = 2$ : 二阶收敛(快得多)
$p = 1.618$ : 超线性收敛(割线法)

大O记号(Big-O Notation):

$f(h) = O(h^p) \Leftrightarrow |f(h)| \leq C h^p \text{ as } h \to 0$

例: 前向差分误差为 $O(h)$ ，中心差分误差为 $O(h^2)$

2. 求解非线性方程

2.1 二分法 (Bisection Method)

原理: 波尔查诺定理 — 若 $f(a) \cdot f(b) < 0$ ，则 $[a, b]$ 之间存在根。

算法:

确认 $f(a) \cdot f(b) < 0$
计算中点 $c = (a + b)/2$
若 $f(c) = 0$ 则结束；若 $f(a) \cdot f(c) < 0$ 则令 $b = c$ ；否则令 $a = c$
若 $|b - a| < \epsilon$ 则结束

收敛性: 每次迭代区间减半，为一阶收敛。 $n$ 次迭代后的误差： $E_n \leq (b-a)/2^n$

所需迭代次数： $n \geq \log_2\!\left(\frac{b-a}{\epsilon}\right)$

def bisection(f, a, b, tol=1e-10, max_iter=100):
    """用二分法求 f(x) = 0 的根"""
    if f(a) * f(b) > 0:
        raise ValueError("f(a) 和 f(b) 必须符号相反")

    history = []
    for i in range(max_iter):
        c = (a + b) / 2
        history.append({'iter': i+1, 'a': a, 'b': b, 'c': c,
                        'f(c)': f(c), 'error': abs(b-a)/2})

        if abs(f(c)) < tol or (b - a) / 2 < tol:
            print(f"收敛: 第 {i+1} 次迭代")
            return c, history

        if f(a) * f(c) < 0:
            b = c
        else:
            a = c

    return (a + b) / 2, history

# 示例: x^3 - 2x - 5 = 0
def equation(x):
    return x**3 - 2*x - 5

root, history = bisection(equation, 1, 3)
print(f"根: {root:.10f}")
print(f"验证: f({root:.6f}) = {equation(root):.2e}")

# 输出收敛历史
print("\n迭代 | 近似值 | 误差")
for h in history[:8]:
    print(f"  {h['iter']:2d} | {h['c']:.8f} | {h['error']:.2e}")

2.2 牛顿-拉夫森法 (Newton-Raphson Method)

原理: 在 $x_n$ 处做一阶泰勒近似

$f(x) \approx f(x_n) + f'(x_n)(x - x_n) = 0$

$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

几何解释: $x_n$ 处切线与x轴的交点即为 $x_{n+1}$

收敛性: 理想情况下为二次收敛 — 每次迭代有效数字位数翻倍！

$E_{n+1} \approx \frac{f''(x^*)}{2f'(x^*)}E_n^2$

局限与注意事项:

若 $f'(x_n) = 0$ 则发散(切线与x轴平行)
不保证收敛(初始点的选择很重要)
在重根处会退化为一阶收敛

def newton_raphson(f, df, x0, tol=1e-10, max_iter=50):
    """牛顿-拉夫森法"""
    x = x0
    history = [{'iter': 0, 'x': x, 'f(x)': f(x), 'error': float('inf')}]

    for i in range(max_iter):
        fx = f(x)
        dfx = df(x)

        if abs(dfx) < 1e-15:
            print("警告: 导数接近0，存在发散风险")
            break

        x_new = x - fx / dfx
        error = abs(x_new - x)
        history.append({'iter': i+1, 'x': x_new,
                        'f(x)': f(x_new), 'error': error})

        if error < tol:
            print(f"收敛: 第 {i+1} 次迭代")
            return x_new, history

        x = x_new

    return x, history

# f(x) = x^3 - 2x - 5, f'(x) = 3x^2 - 2
def f(x):
    return x**3 - 2*x - 5

def df(x):
    return 3*x**2 - 2

root, history = newton_raphson(f, df, x0=2.0)
print(f"牛顿-拉夫森根: {root:.10f}")

print("\n迭代 | 近似值 | 误差")
for h in history:
    err_str = f"{h['error']:.2e}" if h['error'] != float('inf') else "   -"
    print(f"  {h['iter']:2d} | {h['x']:.10f} | {err_str}")

收敛速度比较:

import numpy as np
import matplotlib.pyplot as plt

# 二分法 vs 牛顿-拉夫森法 收敛速度比较
true_root = 2.0945514815423265  # x^3 - 2x - 5 = 0 的真实根

_, hist_bisect = bisection(f, 1, 3, tol=1e-12, max_iter=50)
_, hist_newton = newton_raphson(f, df, x0=2.0, tol=1e-12, max_iter=20)

errors_b = [abs(h['c'] - true_root) for h in hist_bisect]
errors_n = [abs(h['x'] - true_root) for h in hist_newton if h['error'] != float('inf')]

plt.figure(figsize=(10, 5))
plt.semilogy(range(1, len(errors_b)+1), errors_b, 'b-o', label='二分法 (一阶收敛)')
plt.semilogy(range(1, len(errors_n)+1), errors_n, 'r-s', label='牛顿-拉夫森法 (二阶收敛)')
plt.xlabel('迭代次数')
plt.ylabel('绝对误差 (对数坐标)')
plt.title('二分法 vs 牛顿-拉夫森法 收敛速度')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('convergence_comparison.png', dpi=150)
plt.show()

2.3 割线法 (Secant Method)

当不知道 $f'(x)$ 时，用有限差分来近似:

$x_{n+1} = x_n - f(x_n)\cdot\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$

收敛阶数： $p = (1+\sqrt{5})/2 \approx 1.618$ (黄金比例，超线性)

def secant_method(f, x0, x1, tol=1e-10, max_iter=50):
    """割线法"""
    for i in range(max_iter):
        f0, f1 = f(x0), f(x1)
        if abs(f1 - f0) < 1e-15:
            break
        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
        if abs(x2 - x1) < tol:
            print(f"收敛: 第 {i+1} 次迭代")
            return x2
        x0, x1 = x1, x2
    return x1

root_s = secant_method(f, 1.0, 3.0)
print(f"割线法根: {root_s:.10f}")

2.4 使用 scipy.optimize

from scipy.optimize import fsolve, brentq, minimize_scalar
import numpy as np

# 求解单一方程
def equation(x):
    return x**3 - 2*x - 5

# 布伦特法 (二分法 + 割线法 + 反二次插值的混合方法)
root = brentq(equation, 1, 3, xtol=1e-12)
print(f"brentq 根: {root:.12f}")

# fsolve (基于牛顿-拉夫森法)
root_fs = fsolve(equation, 2.0)[0]
print(f"fsolve 根: {root_fs:.12f}")

# 非线性方程组
from scipy.optimize import fsolve

def system(vars):
    x, y = vars
    eq1 = x**2 + y**2 - 4   # 圆
    eq2 = x*y - 1             # 双曲线
    return [eq1, eq2]

solutions = []
for x0, y0 in [(1, 1), (-1, 1), (1, -1), (-1, -1)]:
    sol = fsolve(system, [x0, y0], full_output=True)
    if sol[2] == 1:  # 收敛成功
        solutions.append(sol[0])

print("\n方程组的解:")
for sol in solutions:
    print(f"  (x, y) = ({sol[0]:.6f}, {sol[1]:.6f})")

# 电路方程: 二极管特性 (非线性)
def diode_equation(V, V_s=5, R=1000, I_s=1e-12, n=1, V_T=0.02585):
    """V_s = V + R * I_s * (exp(V/(n*V_T)) - 1)"""
    I = I_s * (np.exp(V / (n * V_T)) - 1)
    return V_s - V - R * I

V_diode = brentq(diode_equation, 0, 1)
I_diode = 1e-12 * (np.exp(V_diode / 0.02585) - 1)
print(f"\n二极管工作点: V = {V_diode:.4f} V, I = {I_diode*1000:.4f} mA")

3. 数值微分与积分

3.1 有限差分 (Finite Difference)

从泰勒级数推导导数的数值近似。

前向差分(Forward Difference):

$f(x+h) = f(x) + hf'(x) + \frac{h^2}{2}f''(x) + \cdots$

$f'(x) \approx \frac{f(x+h) - f(x)}{h} + O(h)$

一阶精度(误差与 $h$ 成正比)。

后向差分(Backward Difference):

$f'(x) \approx \frac{f(x) - f(x-h)}{h} + O(h)$

中心差分(Central Difference):

$f(x+h) = f(x) + hf'(x) + \frac{h^2}{2}f''(x) + \frac{h^3}{6}f'''(x) + \cdots$

$f(x-h) = f(x) - hf'(x) + \frac{h^2}{2}f''(x) - \frac{h^3}{6}f'''(x) + \cdots$

两式相减： $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h} + O(h^2)$

二阶精度 — 在相同的 $h$ 下，比前向差分精确得多。

二阶导数:

$f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} + O(h^2)$

这被用于热方程、波动方程的有限差分解法中。

import numpy as np
import matplotlib.pyplot as plt

def numerical_derivative(f, x, h=1e-5, method='central'):
    """数值微分"""
    if method == 'forward':
        return (f(x + h) - f(x)) / h
    elif method == 'backward':
        return (f(x) - f(x - h)) / h
    elif method == 'central':
        return (f(x + h) - f(x - h)) / (2 * h)

def f(x):
    return np.sin(x)

def df_true(x):
    return np.cos(x)

x0 = np.pi / 4
true_deriv = df_true(x0)

# 按步长 h 分析误差
h_values = np.logspace(-16, 0, 100)
errors_fwd = []
errors_cen = []

for h in h_values:
    e_fwd = abs(numerical_derivative(f, x0, h, 'forward') - true_deriv)
    e_cen = abs(numerical_derivative(f, x0, h, 'central') - true_deriv)
    errors_fwd.append(e_fwd if e_fwd > 0 else 1e-17)
    errors_cen.append(e_cen if e_cen > 0 else 1e-17)

plt.figure(figsize=(10, 6))
plt.loglog(h_values, errors_fwd, 'b-', linewidth=2, label='前向差分 O(h)')
plt.loglog(h_values, errors_cen, 'r-', linewidth=2, label='中心差分 O(h^2)')
plt.loglog(h_values, h_values, 'b--', alpha=0.5, label='O(h) 参考线')
plt.loglog(h_values, h_values**2, 'r--', alpha=0.5, label='O(h^2) 参考线')
plt.xlabel('h (步长)')
plt.ylabel('绝对误差')
plt.title('数值微分误差分析 (误差随h的变化)')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('numerical_derivative_error.png', dpi=150)
plt.show()

3.2 数值积分 (Numerical Integration)

梯形法则(Trapezoidal Rule):

$\int_a^b f(x)\,dx \approx h\left[\frac{f(a)}{2} + f(x_1) + f(x_2) + \cdots + f(x_{n-1}) + \frac{f(b)}{2}\right]$

其中 $h = (b-a)/n$ 。误差： $O(h^2)$

辛普森1/3法则(Simpson's 1/3 Rule):

将每两个子区间用二次多项式近似：

$\int_a^b f(x)\,dx \approx \frac{h}{3}\left[f(a) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \cdots + 4f(x_{n-1}) + f(b)\right]$

$n$ 必须为偶数。误差： $O(h^4)$ — 比梯形法则精确得多！

辛普森3/8法则：三次多项式，误差 $O(h^4)$ ， $n$ 须为3的倍数。

高斯求积法(Gaussian Quadrature):

通过选取最优的积分点和权重，使精度最大化。

$n$ 点高斯-勒让德求积法可精确积分次数不超过 $(2n-1)$ 的多项式。

$\int_{-1}^{1}f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)$

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def trapezoidal(f, a, b, n):
    """梯形法则"""
    x = np.linspace(a, b, n+1)
    h = (b - a) / n
    y = f(x)
    return h * (y[0]/2 + np.sum(y[1:-1]) + y[-1]/2)

def simpsons(f, a, b, n):
    """辛普森1/3法则 (n须为偶数)"""
    if n % 2 != 0:
        n += 1
    x = np.linspace(a, b, n+1)
    h = (b - a) / n
    y = f(x)
    return h/3 * (y[0] + 4*np.sum(y[1::2]) + 2*np.sum(y[2:-1:2]) + y[-1])

# 示例: sin(x) 从 0 到 pi 的积分 = 2
f_test = np.sin
a, b = 0, np.pi
true_val = 2.0

print("数值积分精度比较")
print(f"精确值: {true_val}")
print("-" * 60)
print(f"{'方法':<20} {'n':>6} {'近似值':>15} {'相对误差':>12}")
print("-" * 60)

for n in [4, 8, 16, 32]:
    trap = trapezoidal(f_test, a, b, n)
    simp = simpsons(f_test, a, b, n)
    print(f"{'梯形':<20} {n:>6} {trap:>15.10f} {abs(trap-true_val)/true_val:>12.2e}")
    print(f"{'辛普森':<20} {n:>6} {simp:>15.10f} {abs(simp-true_val)/true_val:>12.2e}")

# scipy.integrate.quad - 自适应求积法
result, error = integrate.quad(f_test, a, b)
print(f"{'scipy.quad':<20} {'自适应':>6} {result:>15.10f} {'<'+str(error):>12}")

# 较难的积分: int_0^inf exp(-x^2) dx = sqrt(pi)/2
def gaussian(x):
    return np.exp(-x**2)

result_gauss, _ = integrate.quad(gaussian, 0, np.inf)
print(f"\n高斯积分: {result_gauss:.10f}")
print(f"精确值 (sqrt(pi)/2): {np.sqrt(np.pi)/2:.10f}")

# 二重积分: int_0^1 int_0^1 x*y dx dy = 0.25
def f2d(y, x):
    return x * y

result_2d, _ = integrate.dblquad(f2d, 0, 1, 0, 1)
print(f"\n二重积分: {result_2d:.10f} (精确值: 0.25)")

4. ODE数值解法

4.1 欧拉法 (Euler's Method)

思路: 用当前斜率预测下一个点

$y_{n+1} = y_n + hf(x_n, y_n)$

局部截断误差： $O(h^2)$ ，整体误差： $O(h)$ (一阶方法)

几何解释: 在每一点沿切线方向移动 $h$

import numpy as np
import matplotlib.pyplot as plt

def euler_method(f, x0, y0, h, x_end):
    """欧拉法"""
    x_vals = [x0]
    y_vals = [y0]
    x = x0
    y = y0

    while x < x_end - 1e-10:
        h_actual = min(h, x_end - x)
        y = y + h_actual * f(x, y)
        x = x + h_actual
        x_vals.append(x)
        y_vals.append(y)

    return np.array(x_vals), np.array(y_vals)

4.2 改进欧拉法 / 休恩法 (Heun's Method)

思路: 先预测(predict)后校正(correct)

$k_1 = f(x_n, y_n)$ $\tilde{y}_{n+1} = y_n + hk_1 \quad \text{(预测)}$ $k_2 = f(x_n + h, \tilde{y}_{n+1})$ $y_{n+1} = y_n + \frac{h}{2}(k_1 + k_2) \quad \text{(校正)}$

整体误差： $O(h^2)$ — 比欧拉法精确得多。

4.3 四阶龙格-库塔法 (RK4)

这是工程领域中应用最广泛的ODE数值解法。

$k_1 = hf(x_n, y_n)$ $k_2 = hf\!\left(x_n + \frac{h}{2},\, y_n + \frac{k_1}{2}\right)$ $k_3 = hf\!\left(x_n + \frac{h}{2},\, y_n + \frac{k_2}{2}\right)$ $k_4 = hf(x_n + h,\, y_n + k_3)$ $y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$

直观理解: 对区间内四个斜率做加权平均。

$k_1$ ：起点斜率
$k_2$ 、 $k_3$ ：中点斜率(两个估计值，分别用 $k_1$ 和 $k_2$ )
$k_4$ ：终点斜率

整体误差： $O(h^4)$ — 误差与 $h^4$ 成正比

比较: 若将 $h$ 减半

欧拉法：误差减为1/2
RK4：误差减为1/16

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def rk4_method(f, x0, y0, h, x_end):
    """四阶龙格-库塔法"""
    x_vals = [x0]
    y_vals = [y0]
    x = x0
    y = y0

    while x < x_end - 1e-10:
        h_actual = min(h, x_end - x)
        k1 = h_actual * f(x, y)
        k2 = h_actual * f(x + h_actual/2, y + k1/2)
        k3 = h_actual * f(x + h_actual/2, y + k2/2)
        k4 = h_actual * f(x + h_actual, y + k3)
        y = y + (k1 + 2*k2 + 2*k3 + k4) / 6
        x = x + h_actual
        x_vals.append(x)
        y_vals.append(y)

    return np.array(x_vals), np.array(y_vals)

# 完整的RLC电路仿真 - 比较多种方法
R, L, C = 50, 0.2, 100e-6
Vs = 100

def rlc_ode(t, state):
    q, dq = state
    d2q = (Vs - R * dq - q/C) / L
    return [dq, d2q]

def rlc_scalar(t, y):
    """将耦合ODE合并为单一函数"""
    return rlc_ode(t, y)

# 初始条件
y0 = [0.0, 0.0]
t_end = 0.02
h_values = [0.001, 0.0005, 0.0001]

# scipy solve_ivp (基准解)
t_ref = np.linspace(0, t_end, 5000)
sol_ref = solve_ivp(rlc_ode, [0, t_end], y0,
                    t_eval=t_ref, method='RK45',
                    rtol=1e-10, atol=1e-12)

v_ref = sol_ref.y[0] / C  # 电容电压

# 用RK4计算
def rlc_2nd_order(t, y):
    """耦合一阶系统"""
    q, dq = y[0], y[1]
    return [dq, (Vs - R*dq - q/C) / L]

fig, axes = plt.subplots(2, 1, figsize=(12, 10))

axes[0].plot(sol_ref.t * 1000, v_ref, 'k-', linewidth=2.5, label='基准 (RK45高精度)')

for h in h_values:
    n_steps = int(t_end / h)
    t_rk4 = np.linspace(0, t_end, n_steps + 1)

    # 手动实现RK4
    y = np.array(y0, dtype=float)
    y_hist = [y.copy()]
    for i in range(n_steps):
        t_i = t_rk4[i]
        k1 = h * np.array(rlc_2nd_order(t_i, y))
        k2 = h * np.array(rlc_2nd_order(t_i + h/2, y + k1/2))
        k3 = h * np.array(rlc_2nd_order(t_i + h/2, y + k2/2))
        k4 = h * np.array(rlc_2nd_order(t_i + h, y + k3))
        y = y + (k1 + 2*k2 + 2*k3 + k4) / 6
        y_hist.append(y.copy())

    y_arr = np.array(y_hist)
    v_rk4 = y_arr[:, 0] / C
    axes[0].plot(t_rk4 * 1000, v_rk4, '--',
                 label=f'RK4 (h={h*1000:.1f}ms)', alpha=0.8)

axes[0].set_xlabel('时间 (ms)')
axes[0].set_ylabel('电容电压 (V)')
axes[0].set_title('RLC串联电路 - RK4数值解')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].axhline(y=Vs, color='gray', linestyle=':', alpha=0.7, label='稳态')

# 电流曲线
i_ref = sol_ref.y[1]
axes[1].plot(sol_ref.t * 1000, i_ref * 1000, 'k-', linewidth=2, label='电流 i(t)')
axes[1].set_xlabel('时间 (ms)')
axes[1].set_ylabel('电流 (mA)')
axes[1].set_title('RLC电路电流')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rlc_numerical.png', dpi=150)
plt.show()

# 检查阻尼特性
omega_0 = 1 / np.sqrt(L * C)
alpha = R / (2 * L)
omega_d = np.sqrt(max(omega_0**2 - alpha**2, 0))
print(f"谐振频率: {omega_0/(2*np.pi):.1f} Hz")
print(f"阻尼系数: {alpha:.1f} s^-1")
if omega_0 > alpha:
    print(f"欠阻尼，阻尼振荡频率: {omega_d/(2*np.pi):.1f} Hz")
elif omega_0 == alpha:
    print("临界阻尼")
else:
    print("过阻尼")

4.4 刚性ODE与自适应步长

刚性ODE(Stiff ODE): 指时间常数相差极大的多个模态同时存在的情形。

例如：快速RC瞬态响应与缓慢温度变化并存的系统

虽然必须使用很小的 $h$ ，但实际上起决定作用的是缓慢变化的成分。

from scipy.integrate import solve_ivp
import numpy as np

# 刚性方程示例: Robertson化学反应
def robertson(t, y):
    """
    dy0/dt = -0.04*y0 + 1e4*y1*y2
    dy1/dt = 0.04*y0 - 1e4*y1*y2 - 3e7*y1^2
    dy2/dt = 3e7*y1^2
    """
    dy0 = -0.04*y[0] + 1e4*y[1]*y[2]
    dy1 = 0.04*y[0] - 1e4*y[1]*y[2] - 3e7*y[1]**2
    dy2 = 3e7*y[1]**2
    return [dy0, dy1, dy2]

y0 = [1, 0, 0]
t_span = (0, 1e11)
t_eval = np.logspace(-6, 11, 1000)

# 非刚性方法 (RK45): 太慢
# 刚性方法 (Radau, BDF): 高效

sol = solve_ivp(robertson, t_span, y0,
                method='Radau',  # 刚性方法
                t_eval=t_eval,
                rtol=1e-6, atol=1e-9)

print(f"函数求值次数: {sol.nfev}")
print(f"雅可比矩阵求值次数: {sol.njev}")

5. 线性系统数值解法

5.1 高斯消元法 (Gaussian Elimination)

前向消元: 通过行运算构造上三角矩阵

回代: 求解上三角系统

运算次数： $O(n^3)$

部分主元法(Partial Pivoting): 在每一列中选取绝对值最大的元素作为主元，以提升数值稳定性。

import numpy as np

def gaussian_elimination_pivoting(A, b):
    """带部分主元法的高斯消元法"""
    n = len(b)
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)

    # 前向消元
    for i in range(n):
        # 部分主元法: 在第i列中找出最大元素
        max_row = i + np.argmax(abs(A[i:, i]))
        A[[i, max_row]] = A[[max_row, i]]
        b[[i, max_row]] = b[[max_row, i]]

        if abs(A[i, i]) < 1e-15:
            raise ValueError("奇异矩阵！")

        for j in range(i+1, n):
            factor = A[j, i] / A[i, i]
            A[j, i:] -= factor * A[i, i:]
            b[j] -= factor * b[i]

    # 回代
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]

    return x

# 示例: 三回路电路网络 (KVL)
A = np.array([
    [10, -3, -1],
    [-3, 12, -2],
    [-1, -2,  8]
], dtype=float)
b = np.array([100, 50, -30], dtype=float)

x_gauss = gaussian_elimination_pivoting(A, b)
x_numpy = np.linalg.solve(A, b)

print("电路回路电流 (高斯消元法):")
for i, v in enumerate(x_gauss):
    print(f"  I{i+1} = {v:.6f} A")

print(f"\n与numpy结果的误差: {np.max(abs(x_gauss - x_numpy)):.2e}")

5.2 LU分解 (LU Decomposition)

$A = LU$

$L$ ：下三角矩阵， $U$ ：上三角矩阵

优点: 对同一个 $A$ ，可以高效求解多个不同的 $b$ 。

$Ax = b \to Ly = b$ (前代)， $Ux = y$ (回代)

from scipy.linalg import lu, lu_factor, lu_solve
import numpy as np

A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]], dtype=float)
b = np.array([8, -11, -3], dtype=float)

# LU分解
P, L, U = lu(A)
print("L 矩阵:")
print(np.round(L, 4))
print("U 矩阵:")
print(np.round(U, 4))

# 用LU分解求解
lu_and_piv = lu_factor(A)
x = lu_solve(lu_and_piv, b)
print(f"\n解: {x}")

# 对多个右端向量高效求解
for b_new in [b, b+1, b*2]:
    x_new = lu_solve(lu_and_piv, b_new)  # 复用 L, U

5.3 迭代法 (Iterative Methods)

在大型稀疏矩阵上比直接法更高效。

雅可比法(Jacobi Method): $x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j \neq i}a_{ij}x_j^{(k)}\right)$

同时更新所有变量。

高斯-赛德尔法(Gauss-Seidel Method): $x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j < i}a_{ij}x_j^{(k+1)} - \sum_{j > i}a_{ij}x_j^{(k)}\right)$

立即使用已更新的值 → 收敛速度快于雅可比法。

收敛条件: 对角占优矩阵(diagonal dominant)

$|a_{ii}| > \sum_{j \neq i}|a_{ij}| \quad \forall i$

import numpy as np

def gauss_seidel(A, b, x0=None, tol=1e-10, max_iter=1000):
    """高斯-赛德尔迭代法"""
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()

    for iteration in range(max_iter):
        x_old = x.copy()
        for i in range(n):
            sum1 = np.dot(A[i, :i], x[:i])
            sum2 = np.dot(A[i, i+1:], x_old[i+1:])
            x[i] = (b[i] - sum1 - sum2) / A[i, i]

        residual = np.linalg.norm(x - x_old)
        if residual < tol:
            print(f"收敛: 第 {iteration+1} 次迭代, 残差: {residual:.2e}")
            return x

    print(f"在最大迭代次数 {max_iter} 处终止, 残差: {residual:.2e}")
    return x

# 对角占优矩阵示例 (节点电压分析)
A = np.array([
    [ 4, -1, -1,  0],
    [-1,  4,  0, -1],
    [-1,  0,  4, -1],
    [ 0, -1, -1,  4]
], dtype=float)
b = np.array([5, 10, 15, 20], dtype=float)

x_gs = gauss_seidel(A, b)
x_exact = np.linalg.solve(A, b)
print(f"高斯-赛德尔: {x_gs}")
print(f"精确解: {x_exact}")
print(f"最大误差: {np.max(abs(x_gs - x_exact)):.2e}")

5.4 使用 numpy.linalg

import numpy as np

A = np.array([[3, -0.1, -0.2],
              [0.1, 7, -0.3],
              [0.3, -0.2, 10]], dtype=float)
b = np.array([7.85, -19.3, 71.4], dtype=float)

# 直接求解
x = np.linalg.solve(A, b)
print(f"解: {x}")

# 特征值
eigenvalues, eigenvectors = np.linalg.eig(A)
print(f"\n特征值: {eigenvalues}")

# 逆矩阵 (不推荐 - 数值误差大)
A_inv = np.linalg.inv(A)
print(f"逆矩阵验证 (A*A_inv - I 的最大值): {np.max(abs(A @ A_inv - np.eye(3))):.2e}")

# 条件数与秩
print(f"条件数: {np.linalg.cond(A):.2f}")
print(f"秩: {np.linalg.matrix_rank(A)}")

# 最小二乘解 (超定系统)
A_over = np.vstack([A, [1, 1, 1]])
b_over = np.append(b, 10)
x_ls, residuals, rank, sv = np.linalg.lstsq(A_over, b_over, rcond=None)
print(f"\n最小二乘解: {x_ls}")

6. 插值与近似

6.1 拉格朗日插值 (Lagrange Interpolation)

经过 $n+1$ 个点 $(x_0, y_0), \ldots, (x_n, y_n)$ 的 $n$ 次多项式:

$P_n(x) = \sum_{k=0}^n y_k L_k(x)$

基函数(Lagrange basis): $L_k(x) = \prod_{j=0, j\neq k}^n \frac{x - x_j}{x_k - x_j}$

龙格现象(Runge's Phenomenon): 用等间距点做高次插值时，在两端会出现振荡。

解决办法：使用切比雪夫点(Chebyshev nodes)

$x_k = \cos\!\left(\frac{(2k+1)\pi}{2(n+1)}\right), \quad k = 0, 1, \ldots, n$

6.2 样条插值 (Spline Interpolation)

在每个子区间使用低次多项式。三次样条最为常用。

相邻区间上函数值、一阶导数、二阶导数均连续 → 曲线光滑。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import lagrange, CubicSpline, interp1d

# 测量数据 (温度 vs 电阻)
x_data = np.array([0, 20, 40, 60, 80, 100], dtype=float)  # 摄氏度
y_data = np.array([100, 108, 116, 125, 134, 143], dtype=float)  # 欧姆

x_fine = np.linspace(0, 100, 300)

# 拉格朗日插值
poly = lagrange(x_data, y_data)
y_lagrange = poly(x_fine)

# 三次样条
cs = CubicSpline(x_data, y_data)
y_spline = cs(x_fine)

# 线性插值
f_linear = interp1d(x_data, y_data)
y_linear = f_linear(x_fine)

plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, s=100, zorder=5, label='测量数据', color='black')
plt.plot(x_fine, y_linear, 'g--', linewidth=1.5, label='线性插值')
plt.plot(x_fine, y_lagrange, 'r-', linewidth=2, label='拉格朗日插值')
plt.plot(x_fine, y_spline, 'b-', linewidth=2, label='三次样条')
plt.xlabel('温度 (C)')
plt.ylabel('电阻 (Ohm)')
plt.title('温度-电阻插值比较')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('interpolation.png', dpi=150)
plt.show()

# 预测50度时的电阻
print(f"50度时的电阻预测:")
print(f"  线性: {f_linear(50):.2f} Ohm")
print(f"  样条: {cs(50):.2f} Ohm")

6.3 最小二乘法 (Least Squares Method)

用 $n$ 次多项式( $m > n+1$ )最佳近似 $m$ 个数据点。

$\min \sum_{i=1}^m [y_i - p(x_i)]^2$

正规方程(Normal Equations): $A^T A \mathbf{c} = A^T \mathbf{y}$

其中 $A$ 是范德蒙德矩阵。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 带噪声的指数衰减数据
np.random.seed(42)
t_data = np.linspace(0, 5, 30)
y_true = 3 * np.exp(-0.7 * t_data)
y_noisy = y_true + 0.2 * np.random.randn(len(t_data))

# 多项式最小二乘拟合
degrees = [1, 2, 3, 5]
t_fine = np.linspace(0, 5, 200)

plt.figure(figsize=(12, 8))
plt.scatter(t_data, y_noisy, s=50, color='black', zorder=5, label='测量数据')
plt.plot(t_fine, 3*np.exp(-0.7*t_fine), 'k-', linewidth=2, label='真实信号')

for deg in degrees:
    coeffs = np.polyfit(t_data, y_noisy, deg)
    y_fit = np.polyval(coeffs, t_fine)
    plt.plot(t_fine, y_fit, '--', linewidth=1.5, label=f'{deg}次多项式')

plt.ylim([-1, 4])
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.title('最小二乘多项式拟合')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('polynomial_fit.png', dpi=150)
plt.show()

# 非线性最小二乘: 指数函数拟合
def exp_model(t, A, k):
    return A * np.exp(-k * t)

popt, pcov = curve_fit(exp_model, t_data, y_noisy, p0=[3, 0.5])
A_fit, k_fit = popt
A_std, k_std = np.sqrt(np.diag(pcov))

print(f"指数拟合结果:")
print(f"  A = {A_fit:.4f} +/- {A_std:.4f} (真值: 3)")
print(f"  k = {k_fit:.4f} +/- {k_std:.4f} (真值: 0.7)")

# 残差与 R^2
y_pred = exp_model(t_data, *popt)
ss_res = np.sum((y_noisy - y_pred)**2)
ss_tot = np.sum((y_noisy - np.mean(y_noisy))**2)
r_squared = 1 - ss_res / ss_tot
print(f"  R^2 = {r_squared:.6f}")

7. 特征值问题数值解法

7.1 幂法 (Power Method)

求矩阵 $A$ 的最大绝对值特征值 $\lambda_1$ 及其对应的特征向量 $\mathbf{v}_1$ 。

算法:

选取初始向量 $\mathbf{q}_0$
$\mathbf{z}_{k+1} = A\mathbf{q}_k$
$\mathbf{q}_{k+1} = \mathbf{z}_{k+1}/\|\mathbf{z}_{k+1}\|$
$\lambda_{k+1} = \mathbf{q}_{k+1}^T A \mathbf{q}_{k+1}$ (瑞利商)
重复直至收敛

收敛速度： $|\lambda_2/\lambda_1|^k$ (特征值比率的幂)

import numpy as np

def power_method(A, tol=1e-10, max_iter=1000):
    """幂法 (最大特征值)"""
    n = A.shape[0]
    q = np.ones(n) / np.sqrt(n)  # 初始向量

    for k in range(max_iter):
        z = A @ q
        eigenvalue = np.dot(q, z)  # 瑞利商
        q_new = z / np.linalg.norm(z)

        if np.linalg.norm(q_new - q) < tol or np.linalg.norm(q_new + q) < tol:
            print(f"收敛: 第 {k+1} 次迭代")
            return eigenvalue, q_new

        q = q_new

    return eigenvalue, q

# 振动系统的固有频率分析
# 双振子质量-弹簧系统的刚度矩阵
k = 1000  # N/m
m = 1     # kg

K = np.array([[2*k, -k], [-k, k]], dtype=float)
M = np.array([[m, 0], [0, m]], dtype=float)

# 广义特征值问题: K*v = lambda*M*v
# 变换: M^(-1/2) * K * M^(-1/2) * u = lambda * u
M_inv_sqrt = np.diag(1/np.sqrt(np.diag(M)))
A_sym = M_inv_sqrt @ K @ M_inv_sqrt

lam_max, v_max = power_method(A_sym)
omega_max = np.sqrt(lam_max)

# 全部特征值 (numpy)
eigenvalues_all = np.linalg.eigvalsh(A_sym)
omega_all = np.sqrt(eigenvalues_all)

print("振动系统固有频率:")
for i, omega in enumerate(omega_all):
    print(f"  模态 {i+1}: {omega/(2*np.pi):.3f} Hz")

print(f"\n幂法求得的最大固有频率: {omega_max/(2*np.pi):.3f} Hz")

7.2 QR算法

可同时求出所有特征值，是现代数值线性代数的核心算法。

思路: 反复做QR分解，矩阵最终收敛为上三角形式，其对角元素即为特征值。

import numpy as np

def qr_algorithm(A, max_iter=1000, tol=1e-10):
    """基本QR算法"""
    A_k = A.copy().astype(float)
    n = A.shape[0]

    for k in range(max_iter):
        Q, R = np.linalg.qr(A_k)
        A_k_new = R @ Q

        # 检查下三角元素的大小
        off_diag_norm = np.max(np.abs(np.tril(A_k_new, -1)))
        if off_diag_norm < tol:
            print(f"QR收敛: 第 {k+1} 次迭代")
            break

        A_k = A_k_new

    return np.diag(A_k)

# 对称矩阵 (实特征值)
A = np.array([
    [4, 1, 0],
    [1, 3, 1],
    [0, 1, 2]
], dtype=float)

eig_qr = qr_algorithm(A)
eig_numpy = np.linalg.eigvalsh(A)

print("QR算法特征值:", np.sort(eig_qr))
print("numpy特征值: ", np.sort(eig_numpy))

7.3 完整的特征值分析：结构共振分析

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg

def vibration_analysis(K, M, n_modes=None):
    """
    广义特征值问题: K*v = omega^2 * M * v
    K: 刚度矩阵, M: 质量矩阵
    """
    # 用scipy的linalg.eigh求解广义特征值问题
    eigenvalues, eigenvectors = linalg.eigh(K, M)

    omega = np.sqrt(np.abs(eigenvalues))  # 角频率 (rad/s)
    freq = omega / (2 * np.pi)           # 频率 (Hz)

    if n_modes is None:
        n_modes = len(eigenvalues)

    print("固有频率分析结果:")
    print("-" * 40)
    for i in range(n_modes):
        print(f"模态 {i+1}: {freq[i]:.3f} Hz ({omega[i]:.2f} rad/s)")

    return omega[:n_modes], eigenvectors[:, :n_modes]

# 三自由度振动系统
k1, k2, k3 = 1000, 1500, 800   # N/m
m1, m2, m3 = 1.0, 1.5, 0.8    # kg

K = np.array([
    [k1+k2,   -k2,    0],
    [-k2,   k2+k3,  -k3],
    [0,       -k3,   k3]
], dtype=float)

M = np.diag([m1, m2, m3])

omegas, modes = vibration_analysis(K, M)

# 模态振型可视化
positions = np.arange(1, 4)
fig, axes = plt.subplots(1, 3, figsize=(14, 5))

for i, ax in enumerate(axes):
    mode_shape = modes[:, i] / np.max(abs(modes[:, i]))
    ax.plot([0] + list(positions), [0] + list(mode_shape),
            'bo-', linewidth=2, markersize=10)
    ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
    ax.set_title(f'模态 {i+1}\n{omegas[i]/(2*np.pi):.2f} Hz')
    ax.set_xlabel('自由度')
    ax.set_ylabel('归一化模态振型')
    ax.grid(True, alpha=0.3)
    ax.set_xticks([0, 1, 2, 3])
    ax.set_xticklabels(['固定端', '质量1', '质量2', '质量3'])

plt.suptitle('三自由度振动系统模态振型')
plt.tight_layout()
plt.savefig('vibration_modes.png', dpi=150)
plt.show()

8. 综合实战案例：MOSFET偏置电路数值分析

import numpy as np
from scipy.optimize import fsolve
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# MOSFET简化模型参数
Vth = 1.5   # 阈值电压 (V)
k = 2e-3    # 跨导参数 (A/V^2)
VDD = 15    # 电源电压 (V)
RD = 3000   # 漏极电阻 (Ohm)
RS = 1000   # 源极电阻 (Ohm)
R1 = 100e3  # 分压电阻1 (Ohm)
R2 = 47e3   # 分压电阻2 (Ohm)

def mosfet_bias_equations(vars):
    """MOSFET偏置非线性方程组"""
    VGS, ID = vars

    # 由分压电路得到栅极电压
    VG = VDD * R2 / (R1 + R2)

    # 漏极电流 (饱和区)
    if VGS > Vth:
        ID_calc = k / 2 * (VGS - Vth)**2
    else:
        ID_calc = 0

    # KVL: VGS = VG - VS = VG - ID*RS
    VGS_calc = VG - ID * RS

    return [VGS - VGS_calc, ID - ID_calc]

VG = VDD * R2 / (R1 + R2)
initial_guess = [VG - 1, 1e-3]

solution = fsolve(mosfet_bias_equations, initial_guess, full_output=True)
VGS_dc, ID_dc = solution[0]
VDS_dc = VDD - ID_dc * (RD + RS)

print("MOSFET DC偏置工作点:")
print(f"  VGS = {VGS_dc:.3f} V")
print(f"  ID  = {ID_dc*1000:.3f} mA")
print(f"  VDS = {VDS_dc:.3f} V")
print(f"  饱和区检验: VDS > VGS - Vth = {VGS_dc - Vth:.3f} V -> {'OK' if VDS_dc > VGS_dc - Vth else '线性区'}")

# 传输特性曲线
VGS_range = np.linspace(0, 5, 200)
ID_range = np.where(VGS_range > Vth,
                    k/2 * (VGS_range - Vth)**2,
                    0)

plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.plot(VGS_range, ID_range * 1000, 'b-', linewidth=2)
plt.axvline(x=VGS_dc, color='r', linestyle='--', label=f'工作点 VGS={VGS_dc:.2f}V')
plt.axhline(y=ID_dc*1000, color='g', linestyle='--', label=f'工作点 ID={ID_dc*1000:.2f}mA')
plt.xlabel('VGS (V)')
plt.ylabel('ID (mA)')
plt.title('MOSFET传输特性')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(1, 2, 2)
VDS_range = np.linspace(0, VDD, 200)

for VGS_val in [2.0, 2.5, 3.0, 3.5]:
    if VGS_val <= Vth:
        continue
    ID_sat = k/2 * (VGS_val - Vth)**2
    ID_curve = np.where(
        VDS_range < VGS_val - Vth,
        k * ((VGS_val - Vth)*VDS_range - VDS_range**2/2),
        ID_sat
    )
    plt.plot(VDS_range, ID_curve * 1000,
             label=f'VGS={VGS_val}V')

# 负载线
ID_load = (VDD - VDS_range) / (RD + RS)
plt.plot(VDS_range, ID_load * 1000, 'k--', linewidth=2, label='负载线')
plt.axvline(x=VDS_dc, color='r', linestyle=':', alpha=0.7)
plt.xlabel('VDS (V)')
plt.ylabel('ID (mA)')
plt.title('MOSFET输出特性与负载线')
plt.legend(fontsize=8)
plt.ylim([0, None])
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('mosfet_analysis.png', dpi=150)
plt.show()

总结与后续

本篇(第4篇)涵盖的内容:

误差分析: 舍入/截断误差、绝对/相对误差、数值稳定性、条件数
非线性方程: 二分法(一阶)、牛顿-拉夫森法(二阶)、割线法(1.618阶)、scipy.optimize
数值微分: 前向/后向/中心差分、误差分析、最优步长
数值积分: 梯形法则、辛普森法则、高斯求积法、scipy.integrate
ODE数值解法: 欧拉法、休恩法、完整的RK4实现、刚性ODE
线性系统: 高斯消元法、LU分解、迭代法、numpy.linalg
插值: 拉格朗日插值、样条插值、龙格现象、最小二乘法、非线性曲线拟合
特征值问题: 幂法、QR算法、振动分析

通过工科数学4篇系列，我们从ODE、傅里叶分析、PDE、复变分析、Z变换一路走到数值分析，完整征服了电子/电气/机械工程的核心数学工具。

参考资料

Burden, R. & Faires, J. "Numerical Analysis", 10th Edition, Cengage
Heath, M. "Scientific Computing: An Introductory Survey", 2nd Edition
Press, W. et al. "Numerical Recipes: The Art of Scientific Computing", 3rd Edition
Kreyszig, E. "Advanced Engineering Mathematics", 10th Edition, Wiley
NumPy 线性代数文档
SciPy 优化文档
SciPy 积分文档