第1章:计算模拟的数学基础

说实话,很多人一听到「数学基础」四个字就想打退堂鼓。我当年也是这样。但做了十几年模拟项目后,我慢慢发现——真正在代码里用到的数学,其实就那么几块。今天咱们就把它们捋清楚。

核心观点:计算模拟的本质,就是用数学语言描述物理世界。你不需要成为数学家,但得知道工具箱里有什么工具。

1.1 微积分基础:变化率的语言

微积分这东西,说白了就是研究「变化」的学问。模拟里最常见的就是跟踪某个量随时间怎么变。

导数:瞬间的变化率

导数描述的是「这一刻的变化速度」。比如你开车,速度表显示的就是位移对时间的导数。

# 数值求导的简单实现
def derivative(f, x, h=1e-6):
    return (f(x + h) - f(x - h)) / (2 * h)

# 例子:求 x^2 在 x=3 处的导数
f = lambda x: x**2
print(derivative(f, 3))  # 输出接近 6.0

我的经验:实际项目中,我很少手算导数。大部分时候用自动微分库(比如PyTorch的autograd)。但理解导数的几何意义——切线的斜率——能帮你调试很多奇怪的问题。

积分:累积的总量

积分是导数的逆运算。如果说导数看「速度」,积分就算「总路程」。模拟里经常用积分来累加某个量,比如总热量、总流量。

# 数值积分:梯形法则
def integrate(f, a, b, n=1000):
    h = (b - a) / n
    total = 0.5 * (f(a) + f(b))
    for i in range(1, n):
        total += f(a + i * h)
    return total * h

# 计算 sin(x) 从 0 到 pi 的积分
import math
print(integrate(math.sin, 0, math.pi))  # 接近 2.0

避坑指南:我曾经在模拟流体时,积分步长选得太大,结果总质量不守恒。后来才意识到——数值积分不是精确的,步长越小越准,但计算量也越大。这是个经典的权衡问题。

1.2 线性代数基础:高维数据的骨架

你想想看,模拟里处理的数据很少是单个数字。温度场有几千个点,每个点有坐标和温度值——这就是向量和矩阵的用武之地。

向量:数据的自然表示

向量就是一串有序的数字。在模拟里,一个粒子的位置是三维向量,一个网格点的物理量也是向量。

# 向量运算示例
import numpy as np

v1 = np.array([1.0, 2.0, 3.0])
v2 = np.array([4.0, 5.0, 6.0])

# 点积:衡量两个向量的相似度
dot_product = np.dot(v1, v2)
print(f"点积: {dot_product}")

# 叉积:得到垂直于两个向量的新向量
cross_product = np.cross(v1, v2)
print(f"叉积: {cross_product}")

矩阵:线性变换的载体

矩阵可以看作是对向量做「操作」的工具。旋转、缩放、投影——这些在图形学和物理模拟里天天用。

# 矩阵乘法:旋转二维向量
import numpy as np

theta = np.pi / 4  # 45度
rotation_matrix = np.array([
    [np.cos(theta), -np.sin(theta)],
    [np.sin(theta),  np.cos(theta)]
])

vector = np.array([1.0, 0.0])
rotated = rotation_matrix @ vector
print(f"旋转后的向量: {rotated}")

特征值:系统的「固有频率」

特征值和特征向量听起来高大上,其实它们描述的是:一个矩阵作用在某个方向上,只是拉伸或压缩,不改变方向。这在结构力学里特别重要——建筑物的固有频率就是特征值问题。

关键理解:特征值告诉你系统「最敏感」的方向。我在做地震模拟时,就是靠特征值分析找出建筑最容易损坏的位置。

1.3 概率论与数理统计:不确定性的量化

模拟不是算命。真实世界充满随机性——材料有缺陷、测量有误差、边界条件不确定。概率统计就是处理这些「不确定性」的工具。

基本概念:随机变量与分布

随机变量描述一个不确定的量。分布告诉我们这个量取不同值的可能性有多大。

# 生成正态分布随机数
import numpy as np
import matplotlib.pyplot as plt

# 均值0,标准差1的正态分布
samples = np.random.normal(0, 1, 10000)

# 计算统计量
mean = np.mean(samples)
std = np.std(samples)
print(f"均值: {mean:.3f}, 标准差: {std:.3f}")

蒙特卡洛方法:用随机数解决确定性问题

这个方法我特别喜欢。它的思路很简单:如果一个问题太难直接求解,那就随机采样,用统计结果逼近真实值。

# 用蒙特卡洛方法估算圆周率
import random

def estimate_pi(num_points=100000):
    inside = 0
    for _ in range(num_points):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        if x*x + y*y <= 1:
            inside += 1
    return 4 * inside / num_points

print(f"估算的圆周率: {estimate_pi()}")

我的习惯:做蒙特卡洛模拟时,我总会先跑一个小样本看看收敛情况。如果结果波动太大,就增加采样数。别一上来就跑百万级样本——那是浪费计算资源。

1.4 数值计算方法简介

理论公式很美,但计算机只能做加减乘除。数值计算方法就是「怎么把连续的问题离散化,让计算机能算」。

离散化:从连续到离散

现实世界是连续的——温度场处处有值。但计算机只能存有限个点。离散化就是把连续空间切成网格,每个网格存一个值。

# 一维热传导方程的离散化
# 理论: dT/dt = alpha * d²T/dx²
# 离散: T_new[i] = T[i] + dt * alpha * (T[i+1] - 2*T[i] + T[i-1]) / dx²

def heat_1d(T, alpha, dt, dx):
    n = len(T)
    T_new = T.copy()
    for i in range(1, n-1):
        T_new[i] = T[i] + dt * alpha * (T[i+1] - 2*T[i] + T[i-1]) / (dx*dx)
    return T_new

迭代法:逐步逼近真解

很多方程没有解析解,只能一步步逼近。迭代法的核心思想是:猜一个初始值,然后不断修正,直到结果足够精确。

# 牛顿法求方程的根
def newton(f, df, x0, tol=1e-6, max_iter=100):
    x = x0
    for i in range(max_iter):
        fx = f(x)
        if abs(fx) < tol:
            return x
        x = x - fx / df(x)
    return x

# 求 x² - 2 = 0 的根(即 sqrt(2))
f = lambda x: x**2 - 2
df = lambda x: 2*x
root = newton(f, df, 1.0)
print(f"sqrt(2) ≈ {root:.6f}")

避坑指南:我曾经用牛顿法求解一个非线性方程组,结果迭代发散。后来发现初始值选得太离谱。记住——迭代法不保证收敛,初始猜测很关键。如果发散,换个初始值或者用更稳健的二分法。

知识体系总览

下面这张图把本章的核心内容串起来了。你可以把它当作一个「数学工具箱」的索引——遇到什么问题,就去翻对应的工具。

计算模拟数学基础 微积分基础 导数 → 变化率 积分 → 累积量 数值微分/积分 线性代数 向量 → 数据表示 矩阵 → 线性变换 特征值 → 固有特性 概率与统计 随机变量与分布 蒙特卡洛方法 数值计算方法 离散化技术 迭代逼近法

嗯,以上就是计算模拟需要的数学基础。别被这些名词吓到——实际用起来,每个工具都有固定的套路。后面几章我们会把这些工具一个个用起来,到时候你就知道它们有多实用了。


公众号:蓝海资料掘金营,微信deep3321