7600 字
38 分钟
计算机图形学笔记(五)：动画与物理仿真
2025-07-24
计算机图形学
C++
/
编程
/
计算机图形学
/
OpenGL

计算机图形学笔记(五)：动画与物理仿真#

这部分用到了物理知识,也是狠狠回忆上了捏.

目录#

  1. 动画基础理论 - 关键帧插值与运动学
  2. 物理仿真数学基础 - 牛顿力学与约束系统
  3. 质点弹簧系统 - 弹性力学与碰撞检测
  4. 数值积分方法 - 欧拉法到Runge-Kutta的演进

动画基础理论#

21.1 动画的数学基础#

21.1.1 插值理论的数学基础#

线性插值的数学定义#

一维线性插值： 给定两个数据点 (t0,f0)(t_0, f_0) 和 (t1,f1)(t_1, f_1)，线性插值函数为： f(t)=f0+t−t0t1−t0(f1−f0)=(1−u)f0+uf1f(t) = f_0 + \frac{t - t_0}{t_1 - t_0}(f_1 - f_0) = (1-u)f_0 + uf_1

其中 u=t−t0t1−t0∈[0,1]u = \frac{t - t_0}{t_1 - t_0} \in [0,1] 是标准化参数。

几何意义： 线性插值在几何上表示连接两点的直线段，满足：

  • 端点插值性：f(t0)=f0f(t_0) = f_0，f(t1)=f1f(t_1) = f_1
  • 线性性：f(αt0+(1−α)t1)=αf0+(1−α)f1f(\alpha t_0 + (1-\alpha)t_1) = \alpha f_0 + (1-\alpha)f_1

多维线性插值#

向量插值： 对于向量 v0,v1∈Rn\mathbf{v}_0, \mathbf{v}_1 \in \mathbb{R}^n： v(t)=(1−t)v0+tv1\mathbf{v}(t) = (1-t)\mathbf{v}_0 + t\mathbf{v}_1

性质：

  • 分量独立性：每个分量独立进行线性插值
  • 仿射不变性：仿射变换与插值可交换

代码实现：

template<typename T>
T lerp(const T& a, const T& b, float t) {
    return a + t * (b - a);
}


// 向量插值
Vector3f lerp_vector(const Vector3f& v1, const Vector3f& v2, float t) {
    return Vector3f(lerp(v1.x(), v2.x(), t),
                   lerp(v1.y(), v2.y(), t),
                   lerp(v1.z(), v2.z(), t));
}


// 颜色插值
Color lerp_color(const Color& c1, const Color& c2, float t) {
    return Color(lerp(c1.r, c2.r, t),
                lerp(c1.g, c2.g, t),
                lerp(c1.b, c2.b, t),
                lerp(c1.a, c2.a, t));
}

球面线性插值（SLERP）的数学理论#

问题背景： 对于单位球面上的点（如单位四元数），线性插值的结果不在球面上，需要特殊的插值方法。

SLERP的数学定义： 对于单位球面上的两点 q1\mathbf{q}_1 和 q2\mathbf{q}_2，球面线性插值为： slerp(q1,q2,t)=sin⁡((1−t)θ)sin⁡θq1+sin⁡(tθ)sin⁡θq2\text{slerp}(\mathbf{q}_1, \mathbf{q}_2, t) = \frac{\sin((1-t)\theta)}{\sin\theta}\mathbf{q}_1 + \frac{\sin(t\theta)}{\sin\theta}\mathbf{q}_2

其中 θ=arccos⁡(q1⋅q2)\theta = \arccos(\mathbf{q}_1 \cdot \mathbf{q}_2) 是两个四元数之间的角度。

几何意义： SLERP沿着连接两点的大圆弧进行插值，保持：

  • 单位长度：插值结果始终在单位球面上
  • 等角速度：角速度恒定
  • 最短路径：沿最短的大圆弧路径

数学推导： 设 q(t)=a(t)q1+b(t)q2\mathbf{q}(t) = a(t)\mathbf{q}_1 + b(t)\mathbf{q}_2，要求：

  1. ∥q(t)∥=1\|\mathbf{q}(t)\| = 1（单位长度）
  2. q(0)=q1\mathbf{q}(0) = \mathbf{q}_1，q(1)=q2\mathbf{q}(1) = \mathbf{q}_2（端点条件）
  3. 等角速度条件

解得：a(t)=sin⁡((1−t)θ)sin⁡θa(t) = \frac{\sin((1-t)\theta)}{\sin\theta}，b(t)=sin⁡(tθ)sin⁡θb(t) = \frac{\sin(t\theta)}{\sin\theta}

Quaternionf slerp(const Quaternionf& q1, const Quaternionf& q2, float t) {
    float dot = q1.dot(q2);


    // 选择最短路径
    Quaternionf q2_corrected = (dot < 0) ? Quaternionf(-q2.coeffs()) : q2;
    dot = std::abs(dot);


    // 如果四元数非常接近，使用线性插值
    if (dot > 0.9995f) {
        Quaternionf result = Quaternionf(q1.coeffs() + t * (q2_corrected.coeffs() - q1.coeffs()));
        result.normalize();
        return result;
    }


    // 球面线性插值
    float theta = std::acos(dot);
    float sin_theta = std::sin(theta);


    float w1 = std::sin((1.0f - t) * theta) / sin_theta;
    float w2 = std::sin(t * theta) / sin_theta;


    return Quaternionf(w1 * q1.coeffs() + w2 * q2_corrected.coeffs());
}

21.1.2 样条插值#

Catmull-Rom样条： 通过控制点的平滑曲线 f(t)=12[2P1+(−P0+P2)t+(2P0−5P1+4P2−P3)t2+(−P0+3P1−3P2+P3)t3]f(t) = \frac{1}{2}[2P_1 + (-P_0 + P_2)t + (2P_0 - 5P_1 + 4P_2 - P_3)t^2 + (-P_0 + 3P_1 - 3P_2 + P_3)t^3]

Vector3f catmull_rom_spline(const Vector3f& p0, const Vector3f& p1,
                           const Vector3f& p2, const Vector3f& p3, float t) {
    float t2 = t * t;
    float t3 = t2 * t;


    Vector3f a = -0.5f * p0 + 1.5f * p1 - 1.5f * p2 + 0.5f * p3;
    Vector3f b = p0 - 2.5f * p1 + 2.0f * p2 - 0.5f * p3;
    Vector3f c = -0.5f * p0 + 0.5f * p2;
    Vector3f d = p1;


    return a * t3 + b * t2 + c * t + d;
}

B样条基函数：

float b_spline_basis(int i, int k, float t, const std::vector<float>& knots) {
    if (k == 0) {
        return (t >= knots[i] && t < knots[i + 1]) ? 1.0f : 0.0f;
    }


    float left_coeff = 0.0f, right_coeff = 0.0f;


    if (knots[i + k] != knots[i]) {
        left_coeff = (t - knots[i]) / (knots[i + k] - knots[i]);
    }


    if (knots[i + k + 1] != knots[i + 1]) {
        right_coeff = (knots[i + k + 1] - t) / (knots[i + k + 1] - knots[i + 1]);
    }


    return left_coeff * b_spline_basis(i, k - 1, t, knots) +
           right_coeff * b_spline_basis(i + 1, k - 1, t, knots);
}

21.2 关键帧动画#

21.2.1 关键帧系统设计#

关键帧数据结构：

template<typename T>
struct Keyframe {
    float time;
    T value;


    // 切线信息（用于Hermite插值）
    T in_tangent;
    T out_tangent;


    // 插值类型
    enum InterpolationType {
        LINEAR,
        CUBIC,
        BEZIER,
        STEP
    } interpolation;
};


template<typename T>
class AnimationCurve {
private:
    std::vector<Keyframe<T>> keyframes;


public:
    void add_keyframe(float time, const T& value,
                     typename Keyframe<T>::InterpolationType interp = Keyframe<T>::LINEAR) {
        Keyframe<T> kf;
        kf.time = time;
        kf.value = value;
        kf.interpolation = interp;


        // 保持时间顺序
        auto it = std::lower_bound(keyframes.begin(), keyframes.end(), kf,
                                  [](const Keyframe<T>& a, const Keyframe<T>& b) {
                                      return a.time < b.time;
                                  });
        keyframes.insert(it, kf);
    }


    T evaluate(float time) const {
        if (keyframes.empty()) return T();
        if (keyframes.size() == 1) return keyframes[0].value;


        // 边界情况
        if (time <= keyframes.front().time) return keyframes.front().value;
        if (time >= keyframes.back().time) return keyframes.back().value;


        // 找到相邻的关键帧
        auto it = std::lower_bound(keyframes.begin(), keyframes.end(), time,
                                  [](const Keyframe<T>& kf, float t) {
                                      return kf.time < t;
                                  });


        const Keyframe<T>& kf1 = *(it - 1);
        const Keyframe<T>& kf2 = *it;


        float t = (time - kf1.time) / (kf2.time - kf1.time);


        switch (kf1.interpolation) {
            case Keyframe<T>::LINEAR:
                return lerp(kf1.value, kf2.value, t);


            case Keyframe<T>::CUBIC:
                return cubic_interpolate(kf1.value, kf1.out_tangent,
                                       kf2.value, kf2.in_tangent, t);


            case Keyframe<T>::STEP:
                return kf1.value;


            default:
                return lerp(kf1.value, kf2.value, t);
        }
    }
};

21.2.2 动画混合#

线性混合：

class AnimationBlender {
private:
    struct AnimationLayer {
        AnimationCurve<Vector3f> position;
        AnimationCurve<Quaternionf> rotation;
        AnimationCurve<Vector3f> scale;
        float weight;
        bool additive;
    };


    std::vector<AnimationLayer> layers;


public:
    void add_layer(const AnimationLayer& layer) {
        layers.push_back(layer);
    }


    Transform evaluate(float time) const {
        Transform result;
        result.position = Vector3f::Zero();
        result.rotation = Quaternionf::Identity();
        result.scale = Vector3f::Ones();


        float total_weight = 0.0f;


        for (const auto& layer : layers) {
            if (layer.weight <= 0.0f) continue;


            Transform layer_transform;
            layer_transform.position = layer.position.evaluate(time);
            layer_transform.rotation = layer.rotation.evaluate(time);
            layer_transform.scale = layer.scale.evaluate(time);


            if (layer.additive) {
                // 加性混合
                result.position += layer_transform.position * layer.weight;
                result.rotation = result.rotation *
                    Quaternionf::Identity().slerp(layer.weight, layer_transform.rotation);
                result.scale += (layer_transform.scale - Vector3f::Ones()) * layer.weight;
            } else {
                // 线性混合
                float normalized_weight = layer.weight;
                if (total_weight > 0) {
                    normalized_weight = layer.weight / (total_weight + layer.weight);
                }


                result.position = lerp(result.position, layer_transform.position, normalized_weight);
                result.rotation = result.rotation.slerp(normalized_weight, layer_transform.rotation);
                result.scale = lerp(result.scale, layer_transform.scale, normalized_weight);


                total_weight += layer.weight;
            }
        }


        return result;
    }
};

21.3 骨骼动画#

21.3.1 骨骼层次结构#

骨骼数据结构：

struct Bone {
    std::string name;
    int parent_index;
    std::vector<int> children_indices;


    // 绑定姿态（T-pose）
    Transform bind_pose;
    Matrix4f inverse_bind_matrix;


    // 当前变换
    Transform local_transform;
    Transform world_transform;
};


class Skeleton {
private:
    std::vector<Bone> bones;
    std::unordered_map<std::string, int> bone_name_to_index;


public:
    void add_bone(const std::string& name, int parent_index, const Transform& bind_pose) {
        Bone bone;
        bone.name = name;
        bone.parent_index = parent_index;
        bone.bind_pose = bind_pose;
        bone.local_transform = bind_pose;
        bone.inverse_bind_matrix = bind_pose.to_matrix().inverse();


        int bone_index = bones.size();
        bones.push_back(bone);
        bone_name_to_index[name] = bone_index;


        // 更新父骨骼的子骨骼列表
        if (parent_index >= 0) {
            bones[parent_index].children_indices.push_back(bone_index);
        }
    }


    void update_world_transforms() {
        for (int i = 0; i < bones.size(); ++i) {
            update_bone_world_transform(i);
        }
    }


private:
    void update_bone_world_transform(int bone_index) {
        Bone& bone = bones[bone_index];


        if (bone.parent_index >= 0) {
            const Bone& parent = bones[bone.parent_index];
            bone.world_transform = parent.world_transform * bone.local_transform;
        } else {
            bone.world_transform = bone.local_transform;
        }
    }


public:
    std::vector<Matrix4f> get_bone_matrices() const {
        std::vector<Matrix4f> matrices;
        matrices.reserve(bones.size());


        for (const auto& bone : bones) {
            Matrix4f bone_matrix = bone.world_transform.to_matrix() * bone.inverse_bind_matrix;
            matrices.push_back(bone_matrix);
        }


        return matrices;
    }
};

21.3.2 蒙皮算法#

线性混合蒙皮（Linear Blend Skinning）：

struct VertexWeight {
    int bone_indices[4];
    float weights[4];


    void normalize() {
        float sum = weights[0] + weights[1] + weights[2] + weights[3];
        if (sum > 0.0f) {
            weights[0] /= sum;
            weights[1] /= sum;
            weights[2] /= sum;
            weights[3] /= sum;
        }
    }
};


class SkinnedMesh {
private:
    std::vector<Vector3f> bind_positions;
    std::vector<Vector3f> bind_normals;
    std::vector<VertexWeight> vertex_weights;


    std::vector<Vector3f> deformed_positions;
    std::vector<Vector3f> deformed_normals;


public:
    void deform(const std::vector<Matrix4f>& bone_matrices) {
        deformed_positions.resize(bind_positions.size());
        deformed_normals.resize(bind_normals.size());


        for (int i = 0; i < bind_positions.size(); ++i) {
            const Vector3f& bind_pos = bind_positions[i];
            const Vector3f& bind_normal = bind_normals[i];
            const VertexWeight& weight = vertex_weights[i];


            Vector3f deformed_pos = Vector3f::Zero();
            Vector3f deformed_normal = Vector3f::Zero();


            for (int j = 0; j < 4; ++j) {
                if (weight.weights[j] > 0.0f) {
                    int bone_index = weight.bone_indices[j];
                    const Matrix4f& bone_matrix = bone_matrices[bone_index];


                    // 变换位置
                    Vector4f pos_homogeneous(bind_pos.x(), bind_pos.y(), bind_pos.z(), 1.0f);
                    Vector4f transformed_pos = bone_matrix * pos_homogeneous;
                    deformed_pos += weight.weights[j] * transformed_pos.head<3>();


                    // 变换法向量
                    Matrix3f normal_matrix = bone_matrix.block<3,3>(0,0);
                    deformed_normal += weight.weights[j] * (normal_matrix * bind_normal);
                }
            }


            deformed_positions[i] = deformed_pos;
            deformed_normals[i] = deformed_normal.normalized();
        }
    }
};

双四元数蒙皮（Dual Quaternion Skinning）： 解决线性混合蒙皮的体积损失问题

struct DualQuaternion {
    Quaternionf real;
    Quaternionf dual;


    DualQuaternion() : real(Quaternionf::Identity()), dual(Quaternionf(0,0,0,0)) {}


    DualQuaternion(const Matrix4f& transform) {
        // 提取旋转和平移
        Matrix3f rotation = transform.block<3,3>(0,0);
        Vector3f translation = transform.block<3,1>(0,3);


        real = Quaternionf(rotation);


        // $dual = 0.5 \times translation \times real$
        Quaternionf trans_quat(0, translation.x(), translation.y(), translation.z());
        dual = Quaternionf(0.5f * (trans_quat * real).coeffs());
    }


    DualQuaternion operator+(const DualQuaternion& other) const {
        DualQuaternion result;
        result.real = Quaternionf(real.coeffs() + other.real.coeffs());
        result.dual = Quaternionf(dual.coeffs() + other.dual.coeffs());
        return result;
    }


    DualQuaternion operator*(float scalar) const {
        DualQuaternion result;
        result.real = Quaternionf(scalar * real.coeffs());
        result.dual = Quaternionf(scalar * dual.coeffs());
        return result;
    }


    void normalize() {
        float norm = real.norm();
        real = Quaternionf(real.coeffs() / norm);
        dual = Quaternionf(dual.coeffs() / norm);
    }


    Vector3f transform_point(const Vector3f& point) const {
        // 从双四元数恢复变换并应用
        Vector3f translation = 2.0f * (dual * real.conjugate()).vec();
        return real * point + translation;
    }
};


void deform_with_dual_quaternions(const std::vector<Matrix4f>& bone_matrices) {
    for (int i = 0; i < bind_positions.size(); ++i) {
        const VertexWeight& weight = vertex_weights[i];


        DualQuaternion blended_dq;


        for (int j = 0; j < 4; ++j) {
            if (weight.weights[j] > 0.0f) {
                int bone_index = weight.bone_indices[j];
                DualQuaternion bone_dq(bone_matrices[bone_index]);


                // 确保四元数在同一半球
                if (j > 0 && blended_dq.real.dot(bone_dq.real) < 0) {
                    bone_dq = bone_dq * (-1.0f);
                }


                blended_dq = blended_dq + bone_dq * weight.weights[j];
            }
        }


        blended_dq.normalize();
        deformed_positions[i] = blended_dq.transform_point(bind_positions[i]);
    }
}

物理仿真数学基础#

22.1 牛顿力学基础#

22.1.1 运动学方程的数学基础#

基本运动量的定义#

位置、速度、加速度的数学关系：

基本物理量定义：

位置： x(t)∈R3\mathbf{x}(t) \in \mathbb{R}^3

速度： v(t)=dxdt\mathbf{v}(t) = \frac{d\mathbf{x}}{dt}

加速度： a(t)=dvdt=d2xdt2\mathbf{a}(t) = \frac{d\mathbf{v}}{dt} = \frac{d^2\mathbf{x}}{dt^2}

匀加速运动的解析解#

常加速度运动方程： 当加速度 a\mathbf{a} 为常数时，运动方程的解析解为：

常加速度运动方程：

速度方程： v(t)=v0+at\mathbf{v}(t) = \mathbf{v}_0 + \mathbf{a}t

位置方程： x(t)=x0+v0t+12at2\mathbf{x}(t) = \mathbf{x}_0 + \mathbf{v}_0 t + \frac{1}{2}\mathbf{a}t^2

推导过程： 从加速度的定义出发： dvdt=a\frac{d\mathbf{v}}{dt} = \mathbf{a}

积分得到： v(t)=∫0ta dτ+v0=at+v0\mathbf{v}(t) = \int_0^t \mathbf{a} \, d\tau + \mathbf{v}_0 = \mathbf{a}t + \mathbf{v}_0

再次积分： x(t)=∫0tv(τ) dτ+x0=∫0t(aτ+v0) dτ+x0=12at2+v0t+x0\mathbf{x}(t) = \int_0^t \mathbf{v}(\tau) \, d\tau + \mathbf{x}_0 = \int_0^t (\mathbf{a}\tau + \mathbf{v}_0) \, d\tau + \mathbf{x}_0 = \frac{1}{2}\mathbf{a}t^2 + \mathbf{v}_0 t + \mathbf{x}_0

代码实现：

struct Particle {
    Vector2D position;
    Vector2D velocity;
    Vector2D acceleration;
    Vector2D forces;
    float mass;
    bool pinned;


    void update_euler(float dt) {
        if (pinned) return;


        acceleration = forces / mass;
        velocity += acceleration * dt;
        position += velocity * dt;


        // 清除力
        forces = Vector2D(0, 0);
    }


    void update_verlet(float dt, Vector2D gravity) {
        if (pinned) return;


        Vector2D new_position = position + (position - last_position) + acceleration * dt * dt;
        last_position = position;
        position = new_position;
    }


private:
    Vector2D last_position;  // Verlet积分需要
};

22.1.2 牛顿第二定律的数学表述#

牛顿第二定律的向量形式#

基本形式： F=ma\mathbf{F} = m\mathbf{a}

其中 F\mathbf{F} 是作用在质量为 mm 的物体上的合外力，a\mathbf{a} 是物体的加速度。

微分方程形式： F=md2xdt2\mathbf{F} = m\frac{d^2\mathbf{x}}{dt^2}

这是一个二阶常微分方程，描述了力与运动的关系。

多力作用的叠加原理#

力的叠加： 当多个力同时作用在物体上时，根据力的叠加原理： ∑i=1nFi=ma\sum_{i=1}^{n} \mathbf{F}_i = m\mathbf{a}

因此加速度为： a=1m∑i=1nFi=F1+F2+⋯+Fnm\mathbf{a} = \frac{1}{m}\sum_{i=1}^{n} \mathbf{F}_i = \frac{\mathbf{F}_1 + \mathbf{F}_2 + \cdots + \mathbf{F}_n}{m}

常见力的数学模型#

重力： Fgravity=mg\mathbf{F}_{gravity} = m\mathbf{g}

弹性力（胡克定律）： Fspring=−k(∥x−x0∥−L0)x−x0∥x−x0∥\mathbf{F}_{spring} = -k(\|\mathbf{x} - \mathbf{x}_0\| - L_0)\frac{\mathbf{x} - \mathbf{x}_0}{\|\mathbf{x} - \mathbf{x}_0\|}

阻尼力： Fdamping=−γv\mathbf{F}_{damping} = -\gamma\mathbf{v}

其中 γ\gamma 是阻尼系数。

在仿真中的应用：

class ForceAccumulator {
private:
    Vector2D total_force;


public:
    void clear() {
        total_force = Vector2D(0, 0);
    }


    void add_force(const Vector2D& force) {
        total_force += force;
    }


    void add_gravity(float mass, const Vector2D& gravity) {
        total_force += mass * gravity;
    }


    void add_spring_force(const Vector2D& position, const Vector2D& other_position,
                         float rest_length, float spring_constant) {
        Vector2D displacement = other_position - position;
        float current_length = displacement.norm();


        if (current_length > 0) {
            Vector2D direction = displacement / current_length;
            float extension = current_length - rest_length;
            Vector2D spring_force = spring_constant * extension * direction;
            total_force += spring_force;
        }
    }


    void add_damping_force(const Vector2D& velocity, float damping_coefficient) {
        total_force -= damping_coefficient * velocity;
    }


    Vector2D get_total_force() const {
        return total_force;
    }
};

22.1.3 能量守恒#

机械能： E=K+U=12mv2+U(x)E = K + U = \frac{1}{2}mv^2 + U(x)

动能： K=12mv2K = \frac{1}{2}mv^2

势能示例：

  • 重力势能：U=mghU = mgh
  • 弹性势能：U=12kx2U = \frac{1}{2}kx^2

能量守恒验证：

class EnergyMonitor {
public:
    float calculate_kinetic_energy(const std::vector<Mass*>& masses) {
        float total_ke = 0.0f;
        for (const auto& mass : masses) {
            if (!mass->pinned) {
                float speed_squared = mass->velocity.norm2();
                total_ke += 0.5f * mass->mass * speed_squared;
            }
        }
        return total_ke;
    }


    float calculate_gravitational_potential(const std::vector<Mass*>& masses,
                                          const Vector2D& gravity) {
        float total_pe = 0.0f;
        for (const auto& mass : masses) {
            // 假设重力向下，y坐标越高势能越大
            total_pe += mass->mass * (-gravity.y) * mass->position.y;
        }
        return total_pe;
    }


    float calculate_elastic_potential(const std::vector<Spring*>& springs) {
        float total_pe = 0.0f;
        for (const auto& spring : springs) {
            Vector2D displacement = spring->m2->position - spring->m1->position;
            float current_length = displacement.norm();
            float extension = current_length - spring->rest_length;
            total_pe += 0.5f * spring->k * extension * extension;
        }
        return total_pe;
    }


    float calculate_total_energy(const std::vector<Mass*>& masses,
                               const std::vector<Spring*>& springs,
                               const Vector2D& gravity) {
        return calculate_kinetic_energy(masses) +
               calculate_gravitational_potential(masses, gravity) +
               calculate_elastic_potential(springs);
    }
};

22.2 弹簧-质点系统#

22.2.1 胡克定律#

基本形式： F=−kxF = -kx 其中k是弹簧常数，x是形变量

向量形式：

Vector2D calculate_spring_force(const Vector2D& pos1, const Vector2D& pos2,
                               float rest_length, float spring_constant) {
    Vector2D displacement = pos2 - pos1;
    float current_length = displacement.norm();


    if (current_length < EPSILON) {
        return Vector2D(0, 0);  // 避免除零
    }


    Vector2D direction = displacement / current_length;
    float extension = current_length - rest_length;


    // 胡克定律：$F = k \times 形变量 \times 方向$
    return spring_constant * extension * direction;
}

22.2.2 阻尼力#

线性阻尼： Fdamping=−γvF_{damping} = -\gamma v

弹簧阻尼：

Vector2D calculate_spring_damping(const Vector2D& vel1, const Vector2D& vel2,
                                 const Vector2D& spring_direction,
                                 float damping_coefficient) {
    // 相对速度
    Vector2D relative_velocity = vel2 - vel1;


    // 沿弹簧方向的速度分量
    float velocity_along_spring = relative_velocity.dot(spring_direction);


    // 阻尼力只作用于弹簧方向
    return -damping_coefficient * velocity_along_spring * spring_direction;
}

22.2.3 完整的弹簧-质点系统#

基于GAMES101 Assignment8的实现：

class Mass {
public:
    Vector2D position;
    Vector2D last_position;  // Verlet积分用
    Vector2D velocity;
    Vector2D forces;
    float mass;
    bool pinned;


    Mass(Vector2D pos, float m, bool pin = false)
        : position(pos), last_position(pos), velocity(0, 0),
          forces(0, 0), mass(m), pinned(pin) {}
};


class Spring {
public:
    Mass* m1;
    Mass* m2;
    float rest_length;
    float k;  // 弹簧常数


    Spring(Mass* mass1, Mass* mass2, float spring_constant)
        : m1(mass1), m2(mass2), k(spring_constant) {
        Vector2D displacement = m2->position - m1->position;
        rest_length = displacement.norm();
    }


    void apply_force() {
        Vector2D displacement = m2->position - m1->position;
        float current_length = displacement.norm();


        if (current_length < EPSILON) return;


        Vector2D direction = displacement / current_length;
        float extension = current_length - rest_length;
        Vector2D spring_force = k * extension * direction;


        // 牛顿第三定律：作用力与反作用力
        m1->forces += spring_force;
        m2->forces -= spring_force;
    }
};


class Rope {
private:
    std::vector<Mass*> masses;
    std::vector<Spring*> springs;


public:
    Rope(Vector2D start, Vector2D end, int num_nodes, float node_mass,
         float spring_k, std::vector<int> pinned_nodes) {


        // 创建质点
        for (int i = 0; i 光线传输的物理意义 num_nodes; ++i) {
            float t = static_cast<float>(i) / (num_nodes - 1);
            Vector2D position = start + t * (end - start);
            bool is_pinned = std::find(pinned_nodes.begin(), pinned_nodes.end(), i)
                           != pinned_nodes.end();
            masses.push_back(new Mass(position, node_mass, is_pinned));
        }


        // 创建弹簧
        for (int i = 0; i < num_nodes - 1; ++i) {
            springs.push_back(new Spring(masses[i], masses[i + 1], spring_k));
        }
    }


    void simulate_euler(float delta_t, Vector2D gravity) {
        // 1. 计算弹簧力
        for (auto& spring : springs) {
            spring->apply_force();
        }


        // 2. 更新质点
        for (auto& mass : masses) {
            if (!mass->pinned) {
                // 添加重力
                mass->forces += mass->mass * gravity;


                // 欧拉积分
                Vector2D acceleration = mass->forces / mass->mass;
                mass->velocity += acceleration * delta_t;
                mass->position += mass->velocity * delta_t;


                // 全局阻尼
                mass->velocity *= 0.99f;
            }


            // 清除力
            mass->forces = Vector2D(0, 0);
        }
    }


    void simulate_verlet(float delta_t, Vector2D gravity) {
        // 1. 计算约束力（弹簧力）
        for (auto& spring : springs) {
            // Verlet积分中，我们直接调整位置来满足约束
            Vector2D displacement = spring->m2->position - spring->m1->position;
            float current_length = displacement.norm();


            if (current_length > EPSILON) {
                Vector2D direction = displacement / current_length;
                float difference = current_length - spring->rest_length;
                Vector2D correction = 0.5f * difference * direction;


                if (!spring->m1->pinned) {
                    spring->m1->position += correction;
                }
                if (!spring->m2->pinned) {
                    spring->m2->position -= correction;
                }
            }
        }


        // 2. Verlet积分
        for (auto& mass : masses) {
            if (!mass->pinned) {
                Vector2D temp_position = mass->position;


                // Verlet公式：$x(t+dt) = 2x(t) - x(t-dt) + a \cdot dt^2$
                Vector2D acceleration = gravity;  // 只考虑重力
                mass->position = 2.0f * mass->position - mass->last_position +
                               acceleration * delta_t * delta_t;


                mass->last_position = temp_position;


                // Verlet阻尼
                mass->position = mass->position * 0.99f + mass->last_position * 0.01f;
            }
        }
    }
};

22.3 数值积分方法#

22.3.1 显式欧拉法#

基本公式：

x(t+h) = x(t) + h·v(t)
v(t+h) = v(t) + h·a(t)

优缺点：

  • ✅ 优点：简单易实现，计算量小
  • ❌ 缺点：数值不稳定，能量不守恒，大时间步长会发散

稳定性分析：

// 弹簧-质点系统的稳定性条件
float calculate_max_stable_timestep(float mass, float spring_constant) {
    // 对于弹簧-质点系统：$dt < \frac{2}{\omega}$，其中$\omega = \sqrt{\frac{k}{m}}$
    float omega = std::sqrt(spring_constant / mass);
    return 2.0f / omega;
}


void adaptive_euler_integration(Rope& rope, float target_dt, Vector2D gravity) {
    float max_safe_dt = calculate_max_stable_timestep(rope.get_min_mass(),
                                                     rope.get_max_spring_constant());


    if (target_dt <= max_safe_dt) {
        rope.simulate_euler(target_dt, gravity);
    } else {
        // 分割时间步
        int substeps = static_cast<int>(std::ceil(target_dt / max_safe_dt));
        float substep_dt = target_dt / substeps;


        for (int i = 0; i < substeps; ++i) {
            rope.simulate_euler(substep_dt, gravity);
        }
    }
}

22.3.2 Verlet积分#

位置Verlet：

$$x(t+h) = 2x(t) - x(t-h) + a(t)h^2$$

速度Verlet： x(t+h)=x(t)+v(t)h+12a(t)h2x(t+h) = x(t) + v(t)h + \frac{1}{2}a(t)h^2 v(t+h)=v(t)+12[a(t)+a(t+h)]hv(t+h) = v(t) + \frac{1}{2}[a(t) + a(t+h)]h

优势：

  • 时间可逆性
  • 更好的能量守恒
  • 对于约束系统更稳定

实现细节：

class VerletIntegrator {
public:
    static void integrate_position(Mass& mass, float dt, const Vector2D& acceleration) {
        if (mass.pinned) return;


        Vector2D new_position = mass.position + (mass.position - mass.last_position) +
                               acceleration * dt * dt;
        mass.last_position = mass.position;
        mass.position = new_position;
    }


    static void integrate_velocity(Mass& mass, float dt, const Vector2D& old_acceleration,
                                  const Vector2D& new_acceleration) {
        if (mass.pinned) return;


        // 速度Verlet：v(t+dt) = v(t) + 0.5*[a(t) + a(t+dt)]*dt
        mass.velocity += 0.5f * (old_acceleration + new_acceleration) * dt;
    }


    static Vector2D calculate_velocity_from_positions(const Mass& mass, float dt) {
        // 从位置差分估算速度
        return (mass.position - mass.last_position) / dt;
    }
};

22.3.3 隐式积分方法#

隐式欧拉：

x(t+h) = x(t) + h·v(t+h)
v(t+h) = v(t) + h·a(t+h)

需要求解线性系统：

class ImplicitEulerSolver {
private:
    // 系统矩阵：$(M - h^2K)\Delta v = h(F + hKv)$
    // M: 质量矩阵, K: 刚度矩阵, F: 外力


public:
    void solve_implicit_step(std::vector<Mass*>& masses,
                           std::vector<Spring*>& springs,
                           float dt, const Vector2D& gravity) {
        int n = masses.size();


        // 构建系统矩阵和右端向量
        Eigen::MatrixXf system_matrix = Eigen::MatrixXf::Zero(2*n, 2*n);
        Eigen::VectorXf rhs = Eigen::VectorXf::Zero(2*n);


        // 质量矩阵部分
        for (int i = 0; i < n; ++i) {
            if (!masses[i]->pinned) {
                float mass = masses[i]->mass;
                system_matrix(2*i, 2*i) = mass;
                system_matrix(2*i+1, 2*i+1) = mass;
            }
        }


        // 刚度矩阵部分
        for (const auto& spring : springs) {
            add_spring_to_system_matrix(system_matrix, spring, dt);
        }


        // 外力
        for (int i = 0; i < n; ++i) {
            if (!masses[i]->pinned) {
                Vector2D force = masses[i]->mass * gravity;
                rhs(2*i) = dt * force.x;
                rhs(2*i+1) = dt * force.y;
            }
        }


        // 求解线性系统
        Eigen::VectorXf delta_v = system_matrix.ldlt().solve(rhs);


        // 更新速度和位置
        for (int i = 0; i < n; ++i) {
            if (!masses[i]->pinned) {
                masses[i]->velocity.x += delta_v(2*i);
                masses[i]->velocity.y += delta_v(2*i+1);
                masses[i]->position += dt * masses[i]->velocity;
            }
        }
    }


private:
    void add_spring_to_system_matrix(Eigen::MatrixXf& matrix, Spring* spring, float dt) {
        // 添加弹簧刚度到系统矩阵
        // 这里需要计算弹簧的雅可比矩阵
        // 实现细节较复杂，涉及非线性弹簧力的线性化
    }
};

22.3.4 Runge-Kutta方法#

四阶Runge-Kutta（RK4）：

class RungeKuttaIntegrator {
public:
    static void rk4_step(Mass& mass, float dt,
                        std::function<Vector2D(const Vector2D&, const Vector2D&)> force_function) {
        if (mass.pinned) return;


        Vector2D x0 = mass.position;
        Vector2D v0 = mass.velocity;


        // k1
        Vector2D k1_v = dt * (force_function(x0, v0) / mass.mass);
        Vector2D k1_x = dt * v0;


        // k2
        Vector2D k2_v = dt * (force_function(x0 + 0.5f * k1_x, v0 + 0.5f * k1_v) / mass.mass);
        Vector2D k2_x = dt * (v0 + 0.5f * k1_v);


        // k3
        Vector2D k3_v = dt * (force_function(x0 + 0.5f * k2_x, v0 + 0.5f * k2_v) / mass.mass);
        Vector2D k3_x = dt * (v0 + 0.5f * k2_v);


        // k4
        Vector2D k4_v = dt * (force_function(x0 + k3_x, v0 + k3_v) / mass.mass);
        Vector2D k4_x = dt * (v0 + k3_v);


        // 更新
        mass.velocity = v0 + (k1_v + 2.0f * k2_v + 2.0f * k3_v + k4_v) / 6.0f;
        mass.position = x0 + (k1_x + 2.0f * k2_x + 2.0f * k3_x + k4_x) / 6.0f;
    }
};

22.4 约束求解#

22.4.1 位置约束#

距离约束： C(x1,x2)=∥x2−x1∥−L=0C(\mathbf{x}_1, \mathbf{x}_2) = \|\mathbf{x}_2 - \mathbf{x}_1\| - L = 0

约束投影方法：

void satisfy_distance_constraint(Mass* m1, Mass* m2, float rest_length) {
    Vector2D displacement = m2->position - m1->position;
    float current_length = displacement.norm();


    if (current_length < EPSILON) return;


    Vector2D direction = displacement / current_length;
    float difference = current_length - rest_length;
    Vector2D correction = 0.5f * difference * direction;


    // 根据质量分配修正量
    float total_mass = m1->mass + m2->mass;
    float w1 = m2->mass / total_mass;  // 质量越大，修正越小
    float w2 = m1->mass / total_mass;


    if (!m1->pinned) {
        m1->position += w1 * correction;
    }
    if (!m2->pinned) {
        m2->position -= w2 * correction;
    }
}

22.4.2 Position Based Dynamics (PBD)#

基本算法流程：

class PBDSolver {
public:
    void simulate_step(std::vector<Mass*>& masses,
                      std::vector<Constraint*>& constraints,
                      float dt, const Vector2D& gravity) {
        // 1. 预测位置
        for (auto& mass : masses) {
            if (!mass->pinned) {
                mass->velocity += dt * gravity;
                mass->predicted_position = mass->position + dt * mass->velocity;
            }
        }


        // 2. 迭代求解约束
        for (int iter = 0; iter < solver_iterations; ++iter) {
            for (auto& constraint : constraints) {
                constraint->project(masses);
            }
        }


        // 3. 更新速度和位置
        for (auto& mass : masses) {
            if (!mass->pinned) {
                mass->velocity = (mass->predicted_position - mass->position) / dt;
                mass->position = mass->predicted_position;
            }
        }
    }


private:
    int solver_iterations = 5;
};


class DistanceConstraint : public Constraint {
private:
    int index1, index2;
    float rest_length;
    float stiffness;


public:
    void project(std::vector<Mass*>& masses) override {
        Mass* m1 = masses[index1];
        Mass* m2 = masses[index2];


        Vector2D displacement = m2->predicted_position - m1->predicted_position;
        float current_length = displacement.norm();


        if (current_length < EPSILON) return;


        Vector2D direction = displacement / current_length;
        float constraint_value = current_length - rest_length;


        // 计算修正量
        float w1 = m1->pinned ? 0.0f : 1.0f / m1->mass;
        float w2 = m2->pinned ? 0.0f : 1.0f / m2->mass;
        float lambda = -constraint_value / (w1 + w2);


        Vector2D correction1 = stiffness * lambda * w1 * direction;
        Vector2D correction2 = -stiffness * lambda * w2 * direction;


        if (!m1->pinned) {
            m1->predicted_position += correction1;
        }
        if (!m2->pinned) {
            m2->predicted_position += correction2;
        }
    }
};

质点弹簧系统#

23.1 系统建模#

23.1.1 质点弹簧系统的数学建模#

系统的数学描述#

状态向量： 质点弹簧系统的状态可以用状态向量描述： s=(x1v1x2v2⋮xnvn)∈R6n\mathbf{s} = \begin{pmatrix} \mathbf{x}_1 \\ \mathbf{v}_1 \\ \mathbf{x}_2 \\ \mathbf{v}_2 \\ \vdots \\ \mathbf{x}_n \\ \mathbf{v}_n \end{pmatrix} \in \mathbb{R}^{6n}

其中 xi∈R3\mathbf{x}_i \in \mathbb{R}^3 是第 ii 个质点的位置，vi∈R3\mathbf{v}_i \in \mathbb{R}^3 是速度。

系统动力学方程： dsdt=(v1a1v2a2⋮vnan)\frac{d\mathbf{s}}{dt} = \begin{pmatrix} \mathbf{v}_1 \\ \mathbf{a}_1 \\ \mathbf{v}_2 \\ \mathbf{a}_2 \\ \vdots \\ \mathbf{v}_n \\ \mathbf{a}_n \end{pmatrix}

其中加速度 ai=Fimi\mathbf{a}_i = \frac{\mathbf{F}_i}{m_i} 由牛顿第二定律确定。

拓扑结构的数学表示#

邻接矩阵： 系统的连接关系可以用邻接矩阵 A∈{0,1}n×n\mathbf{A} \in \{0,1\}^{n \times n} 表示：

邻接矩阵元素定义：

当质点 ii 和 jj 之间有弹簧连接时：Aij=1A_{ij} = 1

其他情况：Aij=0A_{ij} = 0

链式结构（绳子）的数学模型：

class RopeTopology {
public:
    static std::vector<Spring*> create_chain(std::vector<Mass*>& masses, float k) {
        std::vector<Spring*> springs;


        for (int i = 0; i < masses.size() - 1; ++i) {
            springs.push_back(new Spring(masses[i], masses[i + 1], k));
        }


        return springs;
    }
};

网格结构（布料）：

class ClothTopology {
public:
    static std::vector<Spring*> create_cloth_springs(
        std::vector<std::vector<Mass*>>& mass_grid, float k) {


        std::vector<Spring*> springs;
        int rows = mass_grid.size();
        int cols = mass_grid[0].size();


        // 结构弹簧（相邻质点）
        for (int i = 0; i < rows; ++i) {
            for (int j = 0; j < cols; ++j) {
                // 水平弹簧
                if (j < cols - 1) {
                    springs.push_back(new Spring(mass_grid[i][j], mass_grid[i][j + 1], k));
                }
                // 垂直弹簧
                if (i < rows - 1) {
                    springs.push_back(new Spring(mass_grid[i][j], mass_grid[i + 1][j], k));
                }
            }
        }


        // 剪切弹簧（对角线）
        for (int i = 0; i < rows - 1; ++i) {
            for (int j = 0; j < cols - 1; ++j) {
                springs.push_back(new Spring(mass_grid[i][j], mass_grid[i + 1][j + 1], k * 0.5f));
                springs.push_back(new Spring(mass_grid[i + 1][j], mass_grid[i][j + 1], k * 0.5f));
            }
        }


        // 弯曲弹簧（跨越一个质点）
        for (int i = 0; i < rows; ++i) {
            for (int j = 0; j < cols - 2; ++j) {
                springs.push_back(new Spring(mass_grid[i][j], mass_grid[i][j + 2], k * 0.25f));
            }
        }
        for (int i = 0; i < rows - 2; ++i) {
            for (int j = 0; j < cols; ++j) {
                springs.push_back(new Spring(mass_grid[i][j], mass_grid[i + 2][j], k * 0.25f));
            }
        }


        return springs;
    }
};

23.1.2 材料属性#

弹性模量与弹簧常数的关系： k=E×ALk = \frac{E \times A}{L} 其中E是杨氏模量，A是截面积，L是原长

不同材料的参数：

struct MaterialProperties {
    float youngs_modulus;     // 杨氏模量 (Pa)
    float density;            // 密度 $(kg/m^3)$
    float poisson_ratio;      // 泊松比
    float damping_coefficient; // 阻尼系数


    static MaterialProperties steel() {
        return {200e9f, 7850.0f, 0.3f, 0.01f};
    }


    static MaterialProperties rubber() {
        return {0.01e9f, 1000.0f, 0.5f, 0.1f};
    }


    static MaterialProperties cotton() {
        return {5e9f, 1500.0f, 0.4f, 0.05f};
    }
};


float calculate_spring_constant(const MaterialProperties& material,
                               float cross_section_area, float length) {
    return (material.youngs_modulus * cross_section_area) / length;
}

23.2 高级弹簧模型#

23.2.1 非线性弹簧#

指数弹簧模型：

class NonlinearSpring : public Spring {
private:
    float exponential_factor;


public:
    Vector2D calculate_force() override {
        Vector2D displacement = m2->position - m1->position;
        float current_length = displacement.norm();


        if (current_length < EPSILON) return Vector2D(0, 0);


        Vector2D direction = displacement / current_length;
        float strain = (current_length - rest_length) / rest_length;


        // 非线性力：$F = k \times strain \times \exp(\alpha \times |strain|)$
        float force_magnitude = k * strain * std::exp(exponential_factor * std::abs(strain));


        return force_magnitude * direction;
    }
};

分段线性弹簧：

class PiecewiseLinearSpring : public Spring {
private:
    struct ForceSegment {
        float strain_start;
        float strain_end;
        float stiffness;
    };


    std::vector<ForceSegment> segments;


public:
    Vector2D calculate_force() override {
        Vector2D displacement = m2->position - m1->position;
        float current_length = displacement.norm();


        if (current_length < EPSILON) return Vector2D(0, 0);


        Vector2D direction = displacement / current_length;
        float strain = (current_length - rest_length) / rest_length;


        // 找到对应的刚度段
        float stiffness = k;  // 默认刚度
        for (const auto& segment : segments) {
            if (strain >= segment.strain_start && strain < segment.strain_end) {
                stiffness = segment.stiffness;
                break;
            }
        }


        float force_magnitude = stiffness * strain * rest_length;
        return force_magnitude * direction;
    }
};

23.2.2 弹塑性模型#

塑性变形：

class PlasticSpring : public Spring {
private:
    float yield_strain;      // 屈服应变
    float plastic_modulus;   // 塑性模量
    float accumulated_plastic_strain;


public:
    Vector2D calculate_force() override {
        Vector2D displacement = m2->position - m1->position;
        float current_length = displacement.norm();


        if (current_length < EPSILON) return Vector2D(0, 0);


        Vector2D direction = displacement / current_length;
        float total_strain = (current_length - rest_length) / rest_length;
        float elastic_strain = total_strain - accumulated_plastic_strain;


        // 检查是否超过屈服点
        if (std::abs(elastic_strain) > yield_strain) {
            float excess_strain = std::abs(elastic_strain) - yield_strain;
            float strain_sign = (elastic_strain > 0) ? 1.0f : -1.0f;


            // 更新塑性应变
            accumulated_plastic_strain += strain_sign * excess_strain;
            elastic_strain = strain_sign * yield_strain;
        }


        float force_magnitude = k * elastic_strain * rest_length;
        return force_magnitude * direction;
    }
};

23.3 碰撞检测与响应#

23.3.1 质点-平面碰撞#

碰撞检测：

struct Plane {
    Vector2D point;
    Vector2D normal;


    float distance_to_point(const Vector2D& p) const {
        return (p - point).dot(normal);
    }
};


bool check_particle_plane_collision(const Mass& mass, const Plane& plane,
                                   float& penetration_depth) {
    float distance = plane.distance_to_point(mass.position);


    if (distance < 0) {
        penetration_depth = -distance;
        return true;
    }


    return false;
}

碰撞响应：

void resolve_particle_plane_collision(Mass& mass, const Plane& plane,
                                     float restitution_coefficient) {
    float penetration;
    if (!check_particle_plane_collision(mass, plane, penetration)) return;


    // 位置修正
    mass.position += penetration * plane.normal;


    // 速度修正
    float velocity_normal = mass.velocity.dot(plane.normal);
    if (velocity_normal < 0) {  // 只处理朝向平面的速度
        Vector2D velocity_normal_component = velocity_normal * plane.normal;
        Vector2D velocity_tangential = mass.velocity - velocity_normal_component;


        // 法向速度反弹
        Vector2D new_velocity_normal = -restitution_coefficient * velocity_normal_component;


        // 切向摩擦
        float friction_coefficient = 0.3f;
        Vector2D friction_force = -friction_coefficient * velocity_tangential;


        mass.velocity = new_velocity_normal + velocity_tangential + friction_force;
    }
}

23.3.2 自碰撞检测#

空间哈希：

class SpatialHash {
private:
    float cell_size;
    std::unordered_map<int, std::vector<Mass*>> hash_table;


    int hash_position(const Vector2D& pos) {
        int x = static_cast<int>(pos.x / cell_size);
        int y = static_cast<int>(pos.y / cell_size);
        return x * 73856093 ^ y * 19349663;  // 大质数哈希
    }


public:
    void clear() {
        hash_table.clear();
    }


    void insert(Mass* mass) {
        int hash = hash_position(mass->position);
        hash_table[hash].push_back(mass);
    }


    std::vector<Mass*> query_nearby(const Vector2D& position, float radius) {
        std::vector<Mass*> nearby_masses;


        // 检查周围的9个格子
        for (int dx = -1; dx <= 1; ++dx) {
            for (int dy = -1; dy <= 1; ++dy) {
                Vector2D offset_pos = position + Vector2D(dx * cell_size, dy * cell_size);
                int hash = hash_position(offset_pos);


                auto it = hash_table.find(hash);
                if (it != hash_table.end()) {
                    for (Mass* mass : it->second) {
                        float distance = (mass->position - position).norm();
                        if (distance <= radius) {
                            nearby_masses.push_back(mass);
                        }
                    }
                }
            }
        }


        return nearby_masses;
    }
};

连续碰撞检测：

bool continuous_collision_detection(const Mass& mass1, const Mass& mass2,
                                   float dt, float& collision_time) {
    Vector2D relative_position = mass2.position - mass1.position;
    Vector2D relative_velocity = mass2.velocity - mass1.velocity;
    float min_distance = 0.1f;  // 最小安全距离


    // 求解二次方程：$|p + v \cdot t|^2 = min\_distance^2$
    float a = relative_velocity.norm2();
    float b = 2.0f * relative_position.dot(relative_velocity);
    float c = relative_position.norm2() - min_distance * min_distance;


    float discriminant = b * b - 4 * a * c;


    if (discriminant >= 0 && a > EPSILON) {
        float t1 = (-b - std::sqrt(discriminant)) / (2 * a);
        float t2 = (-b + std::sqrt(discriminant)) / (2 * a);


        // 选择在时间步内的最早碰撞时间
        if (t1 >= 0 && t1 <= dt) {
            collision_time = t1;
            return true;
        } else if (t2 >= 0 && t2 <= dt) {
            collision_time = t2;
            return true;
        }
    }


    return false;
}

数值积分方法#

24.1 稳定性分析#

24.1.1 线性稳定性理论#

线性系统的稳定性分析#

线性常微分方程系统： 考虑线性系统： dxdt=Ax\frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x}

其中 A∈Rn×n\mathbf{A} \in \mathbb{R}^{n \times n} 是系统矩阵。

解析解的稳定性： 解析解为 x(t)=eAtx0\mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0，系统稳定当且仅当 A\mathbf{A} 的所有特征值 λi\lambda_i 满足 Re(λi)≤0\text{Re}(\lambda_i) \leq 0。

数值方法的稳定性#

显式欧拉法的放大因子： 显式欧拉法：xn+1=xn+hAxn=(I+hA)xn\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{A}\mathbf{x}_n = (\mathbf{I} + h\mathbf{A})\mathbf{x}_n

放大因子为：G=1+hλG = 1 + h\lambda

稳定性条件： 数值解稳定当且仅当： ∣1+hλ∣≤1|1 + h\lambda| \leq 1

对于所有特征值 λ\lambda。

稳定域分析：

  • 实特征值：λ∈R\lambda \in \mathbb{R}，稳定条件为 −2≤hλ≤0-2 \leq h\lambda \leq 0
  • 复特征值：λ=α+iβ\lambda = \alpha + i\beta，稳定域为复平面上以 (−1,0)(-1, 0) 为圆心，半径为 1 的圆盘

时间步长的选择#

最大稳定时间步长： hmax=2∣λmax∣h_{max} = \frac{2}{|\lambda_{max}|}

其中 λmax\lambda_{max} 是系统矩阵的最大特征值（按模长）。

稳定性测试：

class StabilityAnalyzer {
public:
    static bool is_euler_stable(float eigenvalue_real, float eigenvalue_imag, float dt) {
        std::complex<float> lambda(eigenvalue_real, eigenvalue_imag);
        std::complex<float> amplification_factor = 1.0f + lambda * dt;
        return std::abs(amplification_factor) <= 1.0f;
    }


    static float max_stable_timestep_euler(float max_eigenvalue_magnitude) {
        // 对于实特征值：$dt < \frac{2}{|\lambda_{max}|}$
        return 2.0f / max_eigenvalue_magnitude;
    }


    static void analyze_spring_system_stability(float mass, float spring_constant) {
        float omega = std::sqrt(spring_constant / mass);
        float max_dt = 2.0f / omega;


        std::cout << "弹簧系统稳定性分析:\n";
        std::cout << "  固有频率: " << omega << " rad/s\n";
        std::cout << "  最大稳定时间步: " << max_dt << " s\n";
        std::cout << "  建议时间步: " << max_dt * 0.5f << " s\n";
    }
};

24.1.2 能量守恒性#

能量漂移监控：

class EnergyConservationMonitor {
private:
    float initial_energy;
    std::vector<float> energy_history;


public:
    void initialize(const std::vector<Mass*>& masses,
                   const std::vector<Spring*>& springs,
                   const Vector2D& gravity) {
        EnergyMonitor monitor;
        initial_energy = monitor.calculate_total_energy(masses, springs, gravity);
        energy_history.clear();
    }


    void record_energy(const std::vector<Mass*>& masses,
                      const std::vector<Spring*>& springs,
                      const Vector2D& gravity) {
        EnergyMonitor monitor;
        float current_energy = monitor.calculate_total_energy(masses, springs, gravity);
        energy_history.push_back(current_energy);
    }


    float calculate_energy_drift() const {
        if (energy_history.empty()) return 0.0f;


        float current_energy = energy_history.back();
        return std::abs(current_energy - initial_energy) / initial_energy;
    }


    void print_energy_statistics() const {
        if (energy_history.size() < 2) return;


        float min_energy = *std::min_element(energy_history.begin(), energy_history.end());
        float max_energy = *std::max_element(energy_history.begin(), energy_history.end());
        float energy_range = max_energy - min_energy;


        std::cout << "能量守恒统计:\n";
        std::cout << "  初始能量: " << initial_energy << "\n";
        std::cout << "  能量范围: [" << min_energy << ", " << max_energy << "]\n";
        std::cout << "  能量漂移: " << calculate_energy_drift() * 100 << "%\n";
        std::cout << "  能量振荡: " << energy_range / initial_energy * 100 << "%\n";
    }
};

24.2 自适应时间步长#

24.2.1 误差估计#

Richardson外推法：

class AdaptiveTimestepper {
private:
    float tolerance;
    float min_dt;
    float max_dt;


public:
    float estimate_optimal_timestep(Rope& rope, float current_dt, const Vector2D& gravity) {
        // 保存当前状态
        auto saved_state = rope.save_state();


        // 用当前时间步积分一步
        rope.simulate_euler(current_dt, gravity);
        auto state_full_step = rope.save_state();


        // 恢复状态，用两个半步积分
        rope.restore_state(saved_state);
        rope.simulate_euler(current_dt * 0.5f, gravity);
        rope.simulate_euler(current_dt * 0.5f, gravity);
        auto state_half_steps = rope.save_state();


        // 计算误差
        float error = calculate_state_difference(state_full_step, state_half_steps);


        // 恢复原始状态
        rope.restore_state(saved_state);


        // 根据误差调整时间步
        float safety_factor = 0.8f;
        float new_dt = current_dt * safety_factor * std::pow(tolerance / error, 0.2f);


        return std::clamp(new_dt, min_dt, max_dt);
    }


private:
    float calculate_state_difference(const RopeState& state1, const RopeState& state2) {
        float max_position_diff = 0.0f;


        for (int i = 0; i < state1.positions.size(); ++i) {
            float diff = (state1.positions[i] - state2.positions[i]).norm();
            max_position_diff = std::max(max_position_diff, diff);
        }


        return max_position_diff;
    }
};

24.2.2 CFL条件#

Courant-Friedrichs-Lewy条件：

class CFLCondition {
public:
    static float calculate_cfl_timestep(const std::vector<Spring*>& springs,
                                       const std::vector<Mass*>& masses) {
        float min_dt = std::numeric_limits<float>::max();


        for (const auto& spring : springs) {
            // 计算弹簧的波速
            float reduced_mass = (spring->m1->mass * spring->m2->mass) /
                               (spring->m1->mass + spring->m2->mass);
            float wave_speed = std::sqrt(spring->k / reduced_mass);


            // CFL条件：dt < dx / c
            float element_length = spring->rest_length;
            float cfl_dt = element_length / wave_speed;


            min_dt = std::min(min_dt, cfl_dt);
        }


        return 0.5f * min_dt;  // 安全系数
    }


    static void print_cfl_analysis(const std::vector<Spring*>& springs,
                                  const std::vector<Mass*>& masses) {
        float cfl_dt = calculate_cfl_timestep(springs, masses);


        std::cout << "CFL稳定性分析:\n";
        std::cout << "  CFL时间步限制: " << cfl_dt << " s\n";
        std::cout << "  对应频率: " << 1.0f / cfl_dt << " Hz\n";
    }
};

24.3 高阶积分方法#

24.3.1 多步法#

Adams-Bashforth方法：

class AdamsBashforthIntegrator {
private:
    std::deque<Vector2D> acceleration_history;
    int order;


public:
    AdamsBashforthIntegrator(int method_order) : order(method_order) {}


    void integrate_step(Mass& mass, float dt, const Vector2D& current_acceleration) {
        if (mass.pinned) return;


        acceleration_history.push_back(current_acceleration);


        if (acceleration_history.size() > order) {
            acceleration_history.pop_front();
        }


        Vector2D velocity_increment(0, 0);


        if (acceleration_history.size() == 1) {
            // 一阶（显式欧拉）
            velocity_increment = dt * acceleration_history[0];
        } else if (acceleration_history.size() == 2) {
            // 二阶Adams-Bashforth
            velocity_increment = dt * (1.5f * acceleration_history[1] -
                                     0.5f * acceleration_history[0]);
        } else if (acceleration_history.size() >= 3) {
            // 三阶Adams-Bashforth
            velocity_increment = dt * (23.0f/12.0f * acceleration_history[2] -
                                     16.0f/12.0f * acceleration_history[1] +
                                     5.0f/12.0f * acceleration_history[0]);
        }


        mass.velocity += velocity_increment;
        mass.position += dt * mass.velocity;
    }
};

24.3.2 预测-校正方法#

Adams-Bashforth-Moulton方法：

class PredictorCorrectorIntegrator {
public:
    void integrate_step(Mass& mass, float dt,
                       std::function<Vector2D(const Mass&)> force_function) {
        if (mass.pinned) return;


        // 预测步（Adams-Bashforth）
        Vector2D current_acceleration = force_function(mass) / mass.mass;
        Vector2D predicted_velocity = mass.velocity + dt * current_acceleration;
        Vector2D predicted_position = mass.position + dt * predicted_velocity;


        // 创建预测状态
        Mass predicted_mass = mass;
        predicted_mass.position = predicted_position;
        predicted_mass.velocity = predicted_velocity;


        // 校正步（Adams-Moulton）
        Vector2D predicted_acceleration = force_function(predicted_mass) / mass.mass;
        Vector2D corrected_velocity = mass.velocity +
            0.5f * dt * (current_acceleration + predicted_acceleration);
        Vector2D corrected_position = mass.position +
            0.5f * dt * (mass.velocity + corrected_velocity);


        mass.velocity = corrected_velocity;
        mass.position = corrected_position;
    }
};

24.4 专用物理积分器#

24.4.1 辛积分器#

Störmer-Verlet方法：

class SymplecticIntegrator {
public:
    // 保持哈密顿系统的辛结构
    static void stormer_verlet_step(Mass& mass, float dt,
                                   std::function<Vector2D(const Vector2D&)> force_function) {
        if (mass.pinned) return;


        // 位置半步更新
        mass.position += 0.5f * dt * mass.velocity;


        // 计算新位置的力
        Vector2D force = force_function(mass.position);
        Vector2D acceleration = force / mass.mass;


        // 速度全步更新
        mass.velocity += dt * acceleration;


        // 位置半步更新
        mass.position += 0.5f * dt * mass.velocity;
    }


    // Leapfrog方法（等价于Störmer-Verlet）
    static void leapfrog_step(Mass& mass, float dt,
                             std::function<Vector2D(const Vector2D&)> force_function) {
        if (mass.pinned) return;


        // 速度在半时间步更新
        Vector2D force = force_function(mass.position);
        Vector2D acceleration = force / mass.mass;
        mass.velocity += 0.5f * dt * acceleration;


        // 位置全时间步更新
        mass.position += dt * mass.velocity;


        // 计算新位置的力并完成速度更新
        force = force_function(mass.position);
        acceleration = force / mass.mass;
        mass.velocity += 0.5f * dt * acceleration;
    }
};

24.4.2 约束保持积分器#

SHAKE算法：

class SHAKEIntegrator {
private:
    float tolerance;
    int max_iterations;


public:
    void integrate_with_constraints(std::vector<Mass*>& masses,
                                   std::vector<DistanceConstraint*>& constraints,
                                   float dt, const Vector2D& gravity) {
        // 1. 无约束的Verlet步
        for (auto& mass : masses) {
            if (!mass->pinned) {
                Vector2D acceleration = gravity;
                Vector2D new_position = 2.0f * mass->position - mass->last_position +
                                       acceleration * dt * dt;
                mass->last_position = mass->position;
                mass->position = new_position;
            }
        }


        // 2. 迭代满足约束
        for (int iter = 0; iter < max_iterations; ++iter) {
            bool converged = true;


            for (auto& constraint : constraints) {
                float constraint_error = constraint->evaluate(masses);


                if (std::abs(constraint_error) > tolerance) {
                    converged = false;
                    constraint->apply_correction(masses, dt);
                }
            }


            if (converged) break;
        }


        // 3. 更新速度
        for (auto& mass : masses) {
            if (!mass->pinned) {
                mass->velocity = (mass->position - mass->last_position) / dt;
            }
        }
    }
};


class DistanceConstraint {
private:
    int index1, index2;
    float target_distance;


public:
    float evaluate(const std::vector<Mass*>& masses) {
        Vector2D displacement = masses[index2]->position - masses[index1]->position;
        float current_distance = displacement.norm();
        return current_distance - target_distance;
    }


    void apply_correction(std::vector<Mass*>& masses, float dt) {
        Mass* m1 = masses[index1];
        Mass* m2 = masses[index2];


        Vector2D displacement = m2->position - m1->position;
        float current_distance = displacement.norm();


        if (current_distance < EPSILON) return;


        Vector2D direction = displacement / current_distance;
        float constraint_error = current_distance - target_distance;


        // 计算拉格朗日乘数
        float w1 = m1->pinned ? 0.0f : 1.0f / m1->mass;
        float w2 = m2->pinned ? 0.0f : 1.0f / m2->mass;
        float lambda = -constraint_error / (w1 + w2);


        // 应用位置修正
        Vector2D correction1 = lambda * w1 * direction;
        Vector2D correction2 = -lambda * w2 * direction;


        if (!m1->pinned) m1->position += correction1;
        if (!m2->pinned) m2->position += correction2;
    }
};

这个完整的数值积分框架为物理仿真提供了坚实的数学基础，确保了仿真的稳定性和精度。

计算机图形学笔记(五)：动画与物理仿真
https://kyc001.github.io/posts/计算机图形学笔记五/
作者
kyc001
发布于
2025-07-24
许可协议
CC BY-NC-SA 4.0
计算机图形学笔记(六)：现代图形学前沿
计算机图形学笔记(四)：光线追踪与全局光照
© 2025 kyc001. All Rights Reserved. / RSS / Sitemap
Powered by Astro & Fuwari