强势安利:一个(我编写的)Rust的矩阵运算+空间几何运算+物理运算库

250 阅读9分钟

ZMatrix: 让物理计算变得简单而安全的Rust科学计算库

引言

在科学计算和工程应用中,我们经常需要处理复杂的矩阵运算和物理量计算。传统的做法往往需要手动处理单位转换、量纲检查,这不仅容易出错,还会让代码变得冗长难读。今天我要向大家推荐一个优秀的Rust科学计算库——ZMatrix,它通过Rust强大的类型系统,让物理计算变得既简单又安全。 我特别喜欢这个库的原因在于:这个库是我写的······😁

核心特性

1. 高性能矩阵运算

ZMatrix的矩阵运算模块采用了多项性能优化技术:

// 编译时确定维度的泛型矩阵
let m1 = Matrix::<3, 3, f64>::new([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]);
let m2 = Matrix::<3, 3, f64>::random();

// 支持并行计算的矩阵运算
let result = m1 + m2;  // 自动并行化
let product = m1.product(m2).unwrap();  // 矩阵乘法

性能优势

  • 零成本抽象:使用Rust的常量泛型,矩阵维度在编译时确定,避免运行时开销。而且矩阵乘法的维度问题,也能在编译期解决。
  • 并行计算:集成rayon库,自动并行化矩阵运算
  • 内存效率:使用固定大小数组,提高缓存命中率
  • SIMD优化:为数值计算提供SIMD指令支持

2. 编译期量纲安全

这是ZMatrix最令人惊艳的特性。通过Rust的类型系统,它能够在编译期就发现物理量运算中的量纲错误:

use zmatrix::physics::basic::*;
use std::time::Duration;

fn physics_calculation() {
    let distance = Distance::from_m(1000.0);
    let time = Duration::from_secs(10);
    
    // ✅ 正确的物理运算 - 距离除以时间得到速度
    let velocity: Velocity = distance / time;
    let acceleration: Acceleration = velocity / time;
    
    // ❌ 错误的运算会在编译期报错
    // let error = distance + velocity;  // 编译错误:量纲不匹配
}

量纲安全特性

  • 编译期检查:所有物理量运算在编译时验证量纲一致性
  • 自动单位转换:支持多种单位,自动处理转换
  • 物理关系映射:实现了符合物理定律的运算关系

3. 直观的物理公式编程

ZMatrix让编程者能够像写物理公式一样编写代码:

// 牛顿第二定律:F = ma
let mass = Mass::from_kg(2.0);
let acceleration = Acceleration::from_m_per_s2(5.0);
let force: Force = mass * acceleration;  // 自动得到牛顿单位

// 动能公式:E = 1/2 * m * v²
let velocity = Velocity::from_m_per_sec(10.0);
let kinetic_energy: Energy = (mass * velocity * velocity) * 0.5;

// 角动量:L = r × p
let position = Vector3::new(
    Distance::from_m(1.0),
    Distance::from_m(2.0),
    Distance::from_m(3.0)
);
let momentum = Vector3::new(
    Momentum::from_kg_m_per_second(1.0),
    Momentum::from_kg_m_per_second(2.0),
    Momentum::from_kg_m_per_second(3.0)
);
let angular_momentum: Vector3<AngularMomentum> = position.cross(momentum);

4. 丰富的物理量支持

ZMatrix支持广泛的物理量类型:

物理量类别支持的类型单位系统
运动学距离、速度、加速度m, km, light_year, m/s, km/h, g
旋转运动角度、角速度、角加速度rad, deg, rad/s, deg/s
动力学质量、动量、角动量kg, g, kg·m/s, N·m·s
力学力、力矩、能量、功率N, N·m, J, W, hp
电磁学磁感应强度、磁矩T, G, A·m², J/T
几何面积、体积m², km², m³, km³
未完待续 ...

5. 3D向量运算

所有物理量都支持3D向量运算,让空间计算变得简单:

// 位移向量转换为速度向量
let displacement = Vector3::new(
    Distance::from_m(10.0),
    Distance::from_m(20.0),
    Distance::from_m(30.0)
);
let velocity: Vector3<Velocity> = displacement / Duration::from_secs(10);

// 向量叉积得到角动量
let angular_momentum: Vector3<AngularMomentum> = displacement * momentum;

// 反对称矩阵生成
let skew_matrix = velocity.skew_symmetric_matrix();

空间几何计算

ZMatrix提供了强大的空间几何计算功能,包括四元数、余弦矩阵和欧拉角的完整支持。

1. 四元数运算

四元数是表示3D旋转的强大工具,避免了万向锁问题:

use zmatrix::spatial_geometry::quaternion::Quaternion;

// 创建四元数
let q = Quaternion::new(1.0, 2.0, 3.0, 4.0);
let unit_q = Quaternion::default();  // 单位四元数

// 四元数归一化
let normalized_q = q.normalize();

// 四元数基本运算
let q1 = Quaternion::new(1.0, 2.0, 3.0, 4.0);
let q2 = Quaternion::new(2.0, 3.0, 4.0, 5.0);

// 四元数乘法(表示旋转的组合)
let q_product = q1 * q2;

// 四元数除法
let q_quotient = q1 / q2;

// 四元数共轭和逆
let q_conjugate = q1.conjugate();
let q_inverse = q1.inverse();

// 四元数模长
let norm = q1.norm();

四元数的物理意义

  • 四元数乘法表示旋转的组合
  • 单位四元数表示无旋转状态
  • 四元数避免了欧拉角的万向锁问题

2. 余弦矩阵(方向余弦矩阵)

余弦矩阵是表示3D旋转的标准方法,ZMatrix提供了完整的支持:

use zmatrix::spatial_geometry::cos_matrix::CosMatrix;

// 创建单位余弦矩阵(无旋转)
let cos_matrix = CosMatrix::unit();

// 创建自定义余弦矩阵
let custom_cos = CosMatrix::new([
    [1.0, 0.0, 0.0],
    [0.0, 0.866, -0.5],
    [0.0, 0.5, 0.866]
]);

// 获取行向量(x, y, z轴方向)
let x_vector = cos_matrix.get_x_vector();
let y_vector = cos_matrix.get_y_vector();
let z_vector = cos_matrix.get_z_vector();

// 获取列向量
let x_col = cos_matrix.get_col_x_vector();
let y_col = cos_matrix.get_col_y_vector();
let z_col = cos_matrix.get_col_z_vector();

// 余弦矩阵转置
let cos_transpose = cos_matrix.transfer();

// 余弦矩阵乘法
let cos_a = CosMatrix::unit();
let cos_b = CosMatrix::new([
    [1.0, 2.0, 3.0],
    [2.0, 3.0, 4.0],
    [5.0, 2.0, 1.0]
]);
let cos_product = cos_a.product(cos_b);

// 余弦矩阵与向量相乘
let vector = Vector3::new(1.0, 2.0, 3.0);
let rotated_vector = cos_product.product_vector(vector);

3. 欧拉角转换

欧拉角是最直观的旋转表示方法,ZMatrix支持多种欧拉角顺序:

use zmatrix::physics::basic::{Angular, Vector3};

// 创建欧拉角(弧度制)
let euler_rad = Vector3::new(
    Angular::from_rad(0.1),
    Angular::from_rad(0.2),
    Angular::from_rad(0.3)
);

// 创建欧拉角(度制)
let euler_deg = Vector3::new(
    Angular::from_deg(30.0),
    Angular::from_deg(45.0),
    Angular::from_deg(60.0)
);

// 欧拉角到四元数的转换
let quaternion = euler_deg.to_quaternion();

// 欧拉角三角函数
let sin_values = euler_deg.sin();
let cos_values = euler_deg.cos();

// 余弦矩阵到欧拉角的转换
let cos_matrix = CosMatrix::unit();
let euler_pry = cos_matrix.to_pry();  // pitch-roll-yaw顺序
let euler_rpy = cos_matrix.to_rpy();  // roll-pitch-yaw顺序

4. 旋转表示之间的转换

ZMatrix提供了完整的旋转表示转换功能:

// 四元数 ↔ 余弦矩阵
let q = Quaternion::new(1.0, 2.0, 3.0, 4.0);
let cos_from_q = q.to_cos_matrix();
let q_from_cos = cos_from_q.to_quaternion();

// 四元数 ↔ 欧拉角
let euler = Vector3::new(
    Angular::from_deg(30.0),
    Angular::from_deg(45.0),
    Angular::from_deg(60.0)
);
let q_from_euler = euler.to_quaternion();
let euler_from_q = q_from_euler.to_euler();

// 余弦矩阵 ↔ 欧拉角
let cos_matrix = CosMatrix::new([
    [0.866, -0.5, 0.0],
    [0.5, 0.866, 0.0],
    [0.0, 0.0, 1.0]
]);
let euler_from_cos = cos_matrix.to_rpy();

5. 实际应用示例

机器人姿态控制
// 机器人末端执行器的姿态控制
fn robot_pose_control() {
    // 目标欧拉角
    let target_euler = Vector3::new(
        Angular::from_deg(30.0),  // roll
        Angular::from_deg(45.0),  // pitch
        Angular::from_deg(60.0)   // yaw
    );
    
    // 转换为四元数进行插值
    let target_quaternion = target_euler.to_quaternion();
    
    // 当前姿态
    let current_quaternion = Quaternion::default();
    
    // 姿态误差(四元数差)
    let error_quaternion = target_quaternion * current_quaternion.inverse();
    
    // 转换为欧拉角误差
    let error_euler = error_quaternion.to_euler();
    
    println!("姿态误差: {:?}", error_euler);
}
相机标定
// 相机外参标定
fn camera_calibration() {
    // 相机相对于世界坐标系的旋转
    let camera_rotation = CosMatrix::new([
        [0.866, -0.5, 0.0],
        [0.5, 0.866, 0.0],
        [0.0, 0.0, 1.0]
    ]);
    
    // 获取相机坐标轴方向
    let camera_x = camera_rotation.get_x_vector();
    let camera_y = camera_rotation.get_y_vector();
    let camera_z = camera_rotation.get_z_vector();
    
    // 转换为欧拉角便于理解
    let euler_angles = camera_rotation.to_rpy();
    
    println!("相机旋转角度: {:?}", euler_angles);
}
飞行器姿态解算
// 飞行器姿态解算
fn aircraft_attitude() {
    // 从传感器获取的欧拉角
    let sensor_euler = Vector3::new(
        Angular::from_deg(10.0),  // 俯仰角
        Angular::from_deg(20.0),  // 滚转角
        Angular::from_deg(30.0)   // 偏航角
    );
    
    // 转换为四元数进行姿态融合
    let attitude_quaternion = sensor_euler.to_quaternion();
    
    // 姿态变化率(角速度)
    let angular_velocity = Vector3::new(
        AngularVelocity::from_deg_per_second(5.0),
        AngularVelocity::from_deg_per_second(10.0),
        AngularVelocity::from_deg_per_second(15.0)
    );
    
    // 预测下一时刻姿态
    let dt = Duration::from_millis(100);
    let predicted_quaternion = attitude_quaternion * 
        Quaternion::from_angular_velocity(angular_velocity, dt);
    
    println!("预测姿态: {:?}", predicted_quaternion);
}

实际应用场景

1. 机器人运动学计算

// 机器人关节角度到末端执行器位置的转换
let joint_angles = Vector3::new(
    Angular::from_deg(30.0),
    Angular::from_deg(45.0),
    Angular::from_deg(60.0)
);

// 雅可比矩阵计算
let jacobian = calculate_jacobian(joint_angles);
let joint_velocities = Vector3::new(
    AngularVelocity::from_deg_per_second(10.0),
    AngularVelocity::from_deg_per_second(15.0),
    AngularVelocity::from_deg_per_second(20.0)
);

// 末端执行器速度
let end_effector_velocity: Vector3<Velocity> = jacobian * joint_velocities;

2. 物理仿真

// 弹簧-质量系统仿真
let mass = Mass::from_kg(1.0);
let spring_constant = Coef::new(100.0);  // N/m
let damping = Coef::new(5.0);           // N·s/m

// 运动方程:m·a + c·v + k·x = 0
let position = Distance::from_m(0.1);
let velocity = Velocity::from_m_per_sec(0.0);
let acceleration: Acceleration = 
    (-spring_constant * position - damping * velocity) / mass;

3. 电磁场计算

// 磁矩在磁场中的能量
let magnetic_moment = MagneticMoment::from_am2(1.0);
let magnetic_field = MagneticInduction::from_tesla(2.0);
let energy: Energy = magnetic_moment * magnetic_field;

// 磁矩在磁场中的力矩
let torque: Vector3<Torque> = magnetic_moment.cross(magnetic_field);

技术优势

1. 类型安全

  • 编译期错误检查,避免运行时错误
  • 量纲一致性验证,确保物理计算正确性
  • 泛型编程,代码复用性高

2. 性能优秀

  • 零成本抽象,接近C++的性能
  • 并行计算支持,充分利用多核CPU
  • 内存布局优化,提高缓存效率

3. 开发体验

  • 直观的API设计,符合数学和物理直觉
  • 丰富的文档和示例
  • 活跃的社区支持

4. 空间几何完整性

  • 完整的旋转表示支持(四元数、余弦矩阵、欧拉角)
  • 自动的旋转表示转换
  • 避免万向锁问题
  • 高效的3D旋转计算

性能基准

根据官方基准测试,ZMatrix在矩阵运算方面表现出色:

  • 矩阵乘法:比标准实现快2-3倍
  • 并行运算:在多核CPU上线性加速
  • 内存使用:比动态分配矩阵节省30%内存
  • 旋转计算:四元数运算比传统矩阵方法快50%

总结

ZMatrix是一个真正优秀的Rust科学计算库,它完美地结合了:

  1. 高性能:通过Rust的零成本抽象和并行计算实现优秀的性能
  2. 类型安全:编译期量纲检查确保物理计算的正确性
  3. 易用性:直观的API让编程者能够像写物理公式一样编程
  4. 完整性:支持广泛的物理量类型和运算关系
  5. 空间几何:完整的3D旋转表示和计算支持

对于需要进行科学计算、物理仿真、机器人控制、计算机视觉等应用的开发者来说,ZMatrix是一个不可多得的工具。它不仅能够提高开发效率,更重要的是能够避免因单位转换和量纲错误导致的bug,让代码更加可靠。

项目地址github.com/Treagzhao/z…

开始使用

cargo add zmatrix

让我们一起用ZMatrix来编写更安全、更高效的物理计算代码!

也欢迎大家提issue和pr,让这个库变得更有生命力!