摘要

物理 AI 是人工智能发展的下一个阶段：继感知 AI、生成式 AI 与代理 AI 之后，机器开始被要求理解摩擦、惯性、因果关系等物理概念，并在真实世界中完成感知、决策与行动。机器人是物理 AI 最重要的载体，但机器人学习长期受制于数据稀缺——真实数据的采集成本高、速度慢、风险大，且难以覆盖多样的长尾场景。仿真的价值正在于此：它能在低成本、零风险的虚拟环境中生成近乎无限的数据，并让策略先在仿真中训练、再迁移到真实机器人上，因而成为支撑机器人学习与 Sim-to-Real 迁移的关键基础设施。为此，NVIDIA、Google DeepMind 与 Disney Research 联合开发了开源 GPU 物理仿真框架 Newton。它底层依托 Warp，并采用 USD（通用场景描述）统一描述场景与物体，内置多种物理求解器，覆盖刚体、多体、可变形体、布料和颗粒等对象；除物理仿真外，还提供 IK 逆运动学求解器，将笛卡尔空间目标转换为关节空间配置。针对多物理场景，它支持代理耦合与 ADMM 耦合两种多求解器协同模式。借助 Warp 的 GPU 并行能力与 Worlds 多世界机制，Newton 能够高效生成大规模合成数据，服务于控制策略学习、算法验证与 Sim-to-Real 迁移，目前也可通过 MXMACA 软件栈在沐曦 GPU 上运行。

图 1 Newton 的技术与能力总览：Newton 负责仿真与数据生成，下游训练可由其他框架完成

1. 用于机器人模拟的物理仿真引擎

在 GTC 2025 主题演讲中，人工智能的发展被概括为几个阶段[1]：始于感知 AI，涵盖计算机视觉与语音识别；随后是生成式 AI，专注于在不同模态之间转换；近年来出现的代理 AI 具备自主性，能够理解上下文、推理并规划行动；而正在到来的下一波是物理 AI——它理解摩擦、惯性、因果关系与物体恒常性等概念，让机器人真正能够在物理世界中行动。具身智能，尤其是机器人，正是物理 AI 的重要载体，但机器人学习面临一个突出的瓶颈——数据稀缺：真实世界中的数据采集成本高、速度慢、风险大，且难以覆盖多样化的长尾场景。仿真的价值正在于此：它能在低成本、零风险的虚拟环境中生成近乎无限的数据，并让策略先在仿真中训练、再迁移到真实机器人上，因而成为支撑机器人学习与 Sim-to-Real 迁移的关键基础设施。

在机器人仿真中，物理仿真引擎处于核心位置：数据合成、策略训练、算法验证与 Sim-to-Real 迁移都依赖它。具体而言，仿真可用于训练机器人控制策略，通过强化学习和模仿学习获得可迁移的控制能力；仿真生成的数据还可用于训练更高层的世界模型和基础模型，例如视觉-语言-行动模型；它也能在算法部署到昂贵的真实机器人之前，先在仿真中进行大规模、可重复的测试，验证性能与鲁棒性；此外，Sim-to-Real 迁移研究通过教师-学生蒸馏等技术，可以训练出只依赖真实传感器可用信息（如摄像头、关节角度）的学生策略，从而弥合仿真与现实的差距。

Newton 物理引擎[2,3]正是为此设计的开源 GPU 物理仿真框架，由 NVIDIA、Google DeepMind 与 Disney Research 联合开发，面向机器人学习与 Physical AI。它底层依托 Warp，并采用 USD 作为场景描述标准，目标是在同一体系中容纳大规模 GPU 仿真、多种物理求解器、可扩展的碰撞与接触模型，并与机器人学习框架衔接。基于 Newton 的机器人仿真，希望为机器人控制策略学习、算法验证和 Sim-to-Real 迁移提供统一、可扩展且 GPU 加速的基础设施。

2. 基于USD格式的场景和物体描述

USD（Universal Scene Description，通用场景描述）是由皮克斯（Pixar）提出、现由 OpenUSD 联盟（AOUSD）负责标准化的第三方场景描述标准[4]。Newton 用它来定义仿真场景与其中的物体，同一份 USD 场景也可以被其他工具与引擎使用。

在 Newton 中，USD 是连接资产世界与物理仿真的统一场景描述层。它以层级化、可组合的方式描述复杂的三维世界，能够在一个场景里组织机器人、物体、环境、材质与语义信息，也能描述刚体、柔体、关节、碰撞体、质量与惯性、摩擦和阻尼等物理属性，使求解器能够识别场景中的对象并参与仿真。无论资产来自第三方资产库、由 CAD 或 URDF 等格式转换而来，还是通过 AI 生成、真实扫描重建或手工创建，最终都可以统一成 USD 场景供 Newton 使用。

这样做的好处是"一次描述、多处复用"。机器人、物体、环境以及它们的物理属性只需描述一次，物理仿真求解器与逆运动学求解器就能在同一份场景上工作：前者负责模拟，后者负责规划。USD 还通过引用、实例和分层等机制支持资产的组合与复用，使构建大规模、复杂场景不必从零开始。由此，Newton 把场景构建与物理仿真解耦——上游可以灵活地生产、转换和组合资产，下游则在统一的场景描述上完成仿真、规划、控制策略学习与验证。

图 2 在 Newton 中添加一个 USD 格式的机器人代码片段

3. 物理仿真求解器

物理学在机器人仿真中发挥着至关重要的作用，为在现实环境中准确虚拟呈现机器人行为和交互奠定了基础。机器人仿真的可信度取决于是否遵循基本物理规律——质量与动量守恒、接触与摩擦、执行器与传感器建模。Newton 内置 8 个求解器，各自擅长的领域如下。

1. MuJoCo：面向带接触与约束的关节系统，如工业机械臂、灵巧手，是 Newton 中适配 MuJoCo 模型最直接的选择。它通过 mujoco_warp 提供 MuJoCo 物理引擎的 GPU 加速[5]，并支持 MuJoCo 与 mujoco_warp 两种后端；默认沿用 MuJoCo 自带的碰撞管线。它用广义坐标描述关节结构，只做刚体仿真。

2. Kamino：面向结构复杂的机器人，尤其是带闭环连杆的机构，如并联腿、四连杆传动。它用 PADMM 迭代优化方法求解受约束的刚体运动与受力[6]，关节的硬连接、关节限位、带摩擦的接触与碰撞都能一并处理。它把每个刚体的位置与朝向独立记录（最大坐标表示），因此闭环结构与普通关节同等对待，无需特殊处理。只支持刚体，不支持粒子与布料。

3. Featherstone：面向常规的关节机器人，是刚体仿真的基础选项。它用广义坐标描述关节结构，因此关节约束在数学上被精确满足[7]，不会出现关节漂移；代价是只支持由关节连接的刚体，无法直接模拟自由漂浮的物体。

4. Style3D：面向布料，如衣物、旗帜、软性薄片。它通过迭代方式求解布料运动[8]，可选碰撞处理，并可通过 newton.solvers.style3d 提供的辅助函数把布料网格加进模型。只处理布料，不支持刚体与关节；其中的粒子仅作为布料网格的顶点使用。

5. VBD：一个统一的求解器，粒子和刚体都能处理，适合布料、软体与刚体相互作用的场景。粒子部分使用块下降迭代方法（VBD）[9]，刚体部分使用其增强版本（AVBD）[10]，因此刚体关节、接触以及粒子与刚体的耦合可以在同一次求解中一并完成。它通过软约束处理关节，因此约束精度不如 Featherstone 或 MuJoCo。

6. XPBD：面向刚体与粒子的通用场景，尤其适合需要快速、稳定结果的接触与堆叠问题。它以位置约束的方式处理关节和接触[11]，支持刚体、粒子以及实验性的软体，但布料不支持自碰撞。

7. SemiImplicit：面向粒子、布料与软体的通用场景，并支持基础的梯度求解，可用于可微仿真。它采用辛欧拉积分[12]，容易实现、计算开销小，但并非无条件稳定，需要足够小的时间步长来支撑所需的刚度和阻尼，因此更适合刚度不高的软性物体。支持刚体、关节、粒子和布料（布料不支持自碰撞）。

8. ImplicitMPM：面向颗粒与弹塑性材料，如沙土、粮食、橡皮类材料。它采用隐式物质点法（MPM）[13]，支持屈服、粘性、硬化等常见的颗粒与弹塑性行为；对时间步长无条件稳定，因此很适合模拟很硬或接近完全非弹性的材料。不支持刚体、关节和布料。

这些求解器共同覆盖刚体、关节、可变形体、布料与颗粒材料，使 Newton 能在统一框架下满足多样化的机器人仿真需求。

4. IK逆运动学规划

除物理仿真求解器外，Newton 还提供 IKSolver，即逆运动学（Inverse Kinematics, IK）求解器。它负责运动学层面的问题：用户给定末端执行器（如指尖、手掌）的目标位置和/或朝向，并可附带关节限位等约束，IKSolver 求解出一组关节角度，使末端尽可能接近目标且不违反限位。返回的关节角度可写入机器人模型状态，作为下一步仿真的初始状态或控制参考，也可用于控制策略学习、轨迹规划与验证。需要注意的是，IK 只决定“关节该转到什么角度”，不计算力或力矩；若要在仿真中真正执行这些动作，还需控制器或逆动力学模块把关节目标转换为驱动力或力矩。

在实现上，IKSolver 通过一组加权目标来构造优化问题，让各项目标的总误差最小。目标主要有三类：位置目标让机器人的某个点（如指尖）移动到指定位置，姿态目标让某个连杆（如手掌）转到指定朝向，关节限位则以软约束的方式防止结果超出关节的物理范围。用户可以同时添加多个目标，例如同时要求位置和姿态达标，求解器会在所有可行解中保留整体误差最小的那组关节角度。它从模型中自动读取关节结构与限位范围，并可选两种优化后端：LM 算法适合大多数场景，L-BFGS 更适合关节数量很多的机器人。

由此，IKSolver 把“末端要去哪里”这类直观的机器人任务，转换为可求解的关节空间目标，为后续的仿真、规划或控制提供输入。

5. 多物理仿真的耦合模式

在真实场景中，机器人运动常会与场景中的其他物体发生接触和交互。此时，单一求解器往往不足以完成整个仿真，需要多种不同类型的求解器协同模拟。为此，Newton 框架在其求解器中提供了两类多物理仿真耦合模式：代理耦合（SolverCoupledProxy）和 ADMM（Alternating Direction Method of Multipliers，交替方向乘子法）耦合（SolverCoupledADMM）。需要注意的是，该耦合框架目前仍处于实验（experimental）阶段，其公共 API 与行为可能发生变化，且当前的通用代理与 ADMM 耦合路径最多支持两个求解器入口（Entry）。

代理耦合采用主从式代理策略：将源求解器中的真实端点复制为目标求解器内部的代理体（Proxy Body / Proxy Particle），并为其赋予缩放后的虚拟质量与惯性，使目标求解器能够以代理对象的形式参与耦合计算。在正向同步阶段，源端点的位姿和速度被推送给这些代理体，使目标求解器能够像处理自身元素一样与它们进行碰撞或约束交互，而无需修改目标求解器的内部逻辑。目标求解器完成子步计算后，系统通过“收获”（Harvest）钩子收集代理体上积累的所有接触力与力矩，将其汇总并经松弛因子（Relaxation）滤波后反馈给源求解器。该模式支持“延迟”（Lagged，使用上一步状态以保证稳定）和“交错”（Staggered，使用更新后状态以加速收敛）两种同步方式。代理耦合适合交互以一方为主导的场景，例如将少数刚体作为障碍物投射到密集的变形体（布料、MPM）上；其本质上是求解器之间的松弛不动点迭代，因此虚拟质量缩放、接触刚度与求解器迭代次数都会影响阻尼与收敛性，需要仔细调参以防止数值振荡。

ADMM 耦合则是一种对称的、基于数学约束的框架，它在两个不同求解器拥有的真实端点之间显式建立“界面行”（Interface Rows），而不创建任何虚拟克隆体。在每个耦合子步中，两个子求解器各自独立向前推进，随后 ADMM 耦合器统一求解这些界面约束（如关节、附着或接触），计算出修正力并更新对偶变量（拉格朗日乘子）。最终，耦合器将等大反向的力分别施加回两个真实端点，完成一次迭代；整个过程在固定迭代次数内重复，形成“预测-校正”闭环。该方式需要针对每种交互类型显式实现对应的界面行求解逻辑：当前通用 ADMM 路径支持跨求解器模型关节（仅球铰 BALL、固定 FIXED 和转动 REVOLUTE 三类）、粒子-刚体附着，以及刚-刚、刚-粒、粒-粒接触；其中球铰生成锚点重合的平移行，固定关节额外增加角度行，转动关节保持铰链轴并可附加干摩擦行。其调参（如惩罚权重 rho 与近端权重 gamma）较为敏感，适合追求精确约束且允许迭代收敛的刚-柔或刚-刚耦合场景。

这些耦合机制在概念和数值实现上都较为复杂，但 Newton 已将其封装在 SolverCoupled 系列实验性接口中，使用者无需从零实现子求解器之间的状态同步、力收集与约束迭代。实际使用时，通常只需声明耦合关系并设置相应参数：代理耦合需要配置入口（SolverCoupled.Entry）与代理映射（SolverCoupledProxy.Proxy），ADMM 耦合需要配置入口与 ADMM 配置（SolverCoupledADMM.Config）。框架会自动完成正向同步、子步推进、代理力收集或界面约束求解，并在迭代中完成反馈与校正，从而把多物理耦合的复杂度留在底层。这样，用户既能保持各专用求解器的独立性，又能以较低门槛构建跨求解器的多物理仿真；对于需要精细调节的场景，仍可通过松弛因子、虚拟质量缩放、惩罚权重 rho 与近端权重 gamma 等参数调整稳定性与收敛性。总体而言，Newton 的多物理耦合机制在扩展性与易用性之间取得了平衡，使机器人仿真能够覆盖刚体、软体、布料、颗粒材料等多种物理对象及其相互作用。

图 3 Proxy 耦合示例：入口定义与代理参数配置（mass_scale 为虚拟质量缩放，mode 可选 lagged 或 staggered，proxy_relaxation 为力滤波松弛）

图 4 ADMM 耦合示例：入口定义与配置参数（iterations 为耦合迭代次数，rho 为惩罚刚度，gamma 为近端牵引权重，另含接触对与接触匹配等选项）

6. 沐曦GPU和MXMACA软件栈支持基于Warp的实现

Warp 是 Newton 实现 GPU 加速与可微仿真的底层计算框架。开发者可以用接近 Python 的写法编写并行计算程序，这些程序在运行时被即时编译（Just-In-Time, JIT）为 GPU 代码，从而高效运行仿真。Warp 具备自动微分能力，能够承担物理求解、碰撞检测、粒子计算等任务，让原本需要数天的仿真缩短到几分钟。它还能与 PyTorch 等机器学习框架衔接，使仿真可以纳入端到端训练流程，支持基于梯度的优化、系统辨识与策略学习。因此，Warp 既是 Newton 的 GPU 加速引擎，也是连接物理仿真与机器学习的重要桥梁。

Warp 的基础编程模型是"一线程一元素"的 SIMT 形式：内核在每个线程上独立执行，内核启动时的网格与块维度直接决定并行的粒度。这一模型足以覆盖碰撞检测、粒子积分、接触求解等大多数计算，Newton 的求解器也主要建立在其上。在此之上，Warp 还扩展出了块级编程模型（tile）：把数据划分成块，由块内所有线程协作完成加载、归约与矩阵运算，从而更充分地利用共享内存与张量核心等硬件单元，并减少手工索引与共享内存管理；tile 程序同样支持自动微分，可一并纳入训练回路。在 Newton 中，这一模型通常以选项的形式暴露——例如 Featherstone 求解器的 use_tile_gemm 参数打开后，矩阵运算就会切换到 tile 路径。Featherstone 求解器求解关节加速度时，需要对每个机器人实例的广义惯性矩阵做 Cholesky 分解，这一步在 tile 路径下由如下内核完成（源码见 newton/_src/solvers/featherstone/kernels.py）：

wp.tile_load 把整个矩阵搬入共享内存，wp.tile_diag_add 加上关节阻尼与惯量项，wp.tile_cholesky 完成分解，wp.tile_transpose 调整存储顺序后写回。每一步都由块内线程协作完成，而 articulation 取自 wp.tid()，即每个线程块负责一个机器人实例的关节系统。

支撑这些 tile 运算的是 Warp 自带的块级协作实现，不依赖任何外部数学库即可运行。在此基础上，Warp 还能接入 NVIDIA 的 MathDx 设备端数学库以进一步提升性能：接入后，tile 的矩阵乘法与 Cholesky 分解会改由其中的 cuBLASDx、cuSolverDx 等接口执行，从而借助经过深度优化的硬件实现获得更好的运行时表现。这一接入是可选的，代价在于 MathDx 的编译过程较慢，因此实际使用中需要在编译时间与运行效率之间做取舍。

目前，MXMACA 软件栈已成功支持 Warp，在沐曦 GPU 上可以使用 Newton 物理引擎进行机器人仿真。

 

图 5 展示了一个方块摆放任务在沐曦 GPU 上的执行效果：场景用 MuJoCo 求解器仿真，5 个方块的初始位置在平面上随机给定；机械臂的每个动作——接近、抓取、抬起、移动到放置点、松开、退回——都由 IK 逆运动学求解器规划，位置与姿态目标再加上关节限位构成其优化目标，最终把方块逐个水平放置、排成 METAX 字样。由于初始摆位随机给定，每一次运行都是一条不同的轨迹，而不是同一段动作的重复播放。

 

图 6 给出的是上面摆放方块算例在沐曦 GPU 与 国际主流 GPU 产品上的执行结果对比：第 2000 帧时各方块的位置(x, y, z)。以欧氏距离衡量，此时前三个方块已完成放置，与目标位置的偏差分别为 0.0042、0.0035 与 0.0023；其余两个方块仍在搬运途中，距目标尚有一段距离。逐分量比较，同一方块在两平台上的坐标差最大为 2×10⁻⁶ m，属于浮点舍入的量级，说明 Newton 在沐曦 GPU 上的数值表现与国际主流GPU平台接近，可以为国内团队开展机器人仿真与强化学习训练提供一个可用的平台选择。

 

不过，单个机器人的仿真算例通常包含的物体并不多，计算量有限，很难把 GPU 用满。为此 Newton 引入了"Worlds（世界）"机制：可以在一个模型里同时运行多个互不干扰的仿真，每个世界都有自己的机器人与场景物体，世界之间不会发生碰撞。这一机制主要面向大规模并行仿真，尤其服务于强化学习（Reinforcement Learning, RL）——成千上万个机器人实例可以在同一块 GPU 上同时训练，数据采集效率随之大幅提升。它也为域随机化提供了天然载体：不同世界可以设置不同的物理参数与初始条件，包括方块的摆放位置，一次仿真就能生成大量多样化的轨迹与交互数据。实际使用时，通过 world_count 参数设置并行世界的数量即可：把 world_count 从默认的十几个调高到数十甚至数百，一次仿真即可产出海量轨迹数据，图 7 即为多世界并行时的场景。与此同时，单世界下闲置的计算资源也被填满：提高并行世界数量可以显著提升任务并行度，减少小规模仿真造成的资源空闲。

图 5 用 MuJoCo 求解器仿真机械臂、由逆运动学规划抓取动作，在沐曦 GPU 上将 5 个方块依次放置成 METAX 字样（方块初始位置在平面上随机给定）

图 6 把 world_count 调高后，同一场景在单块 GPU 上并行推进大量世界，多个机器人实例同时执行抓放任务，各世界的方块初始位置随机不同

7. 总结与展望

本文介绍了一款面向机器人仿真领域、功能强大且支持 GPU 加速的开源 Newton 物理仿真引擎，并阐述其核心能力：基于 USD 格式统一描述场景与物体，多种物理仿真求解器与逆运动学（IK）求解器，多求解器协同的多物理耦合机制，以及基于 Warp 的 GPU 并行计算与 Worlds 多世界机制。其中，Warp 将并行计算程序即时编译到 GPU，Worlds 使成千上万个相互独立的世界可在一次仿真中并行推进，从而为强化学习等场景大规模生成合成数据。目前，Newton 已可通过 MXMACA 软件栈在沐曦 GPU 上运行，可用于机器人控制策略学习与验证等下游任务。从更长的视角看，Newton 这类 GPU 物理引擎还有相当大的成长空间，三个方向值得关注。

其一是求解精度与可信度的提升。当前机器人仿真求解器大致处在游戏物理与工程仿真之间：它追求毫秒级实时性，数值精度以"策略能训出来"为准。随着机器人从实验室走向产线与户外，接触、摩擦、磨损、装配公差等问题的可信度会越来越重要，物理引擎的精度目标也会随之向工程仿真靠拢。

其二是多物理场耦合的深化。Newton 目前的代理耦合与 ADMM 耦合本质上是把不同求解器拼接起来，交换的物理量主要还局限在力这一项上；而高可信仿真常常需要多个物理场同时参与，例如力与热、流与固、电磁与力学。要把机器人的电机发热、线缆形变、气流扰动与控制器的电学响应一并算进去，耦合机制需要从"传力"扩展到"传多个物理量"——这正是成熟 CAE 工具多年积累的地方。

其三是与决策优化求解器的结合。物理仿真回答"这样做会发生什么"，决策优化回答"应该怎么做"，两者天然互补，又都是 GPU 并行擅长的工作负载。按任务的时空尺度，这种结合可分三层：最上层是任务与路径规划，在仓库、产线等大空间内统筹任务分配、多机调度与路径，决策变量动辄成千上万，解的优劣直接决定吞吐；中间层是轨迹优化，在避障与动力学约束之间寻找最优轨迹，其解空间连续可导，适合梯度类方法；最下层是运动规划与精细操作，如装配中的插入、拧紧与抓取姿态选择，目标更局部，但对接触、摩擦与力的描述要求最严，恰由 Newton 的求解器兜住。文献中把这一分层称作任务与运动规划（Task and Motion Planning, TAMP）。三层能够衔接，靠的是仿真同时充当优化问题的模型来源与验证环境：几何、惯性、雅可比等信息来自仿真，可微仿真还提供梯度，优化解出的动作再回注仿真中验证，形成"上层规划—下层验证"的闭环。

这三个方向的共同前提，是算力仍有余量。多数机器人仿真算例规模有限，很难独自用满 GPU，目前主要靠 Worlds 多世界并行把算力填满；反过来说，这也意味着在不增加硬件的前提下仍有算力可用于更昂贵的工作：更密的接触迭代、更小的积分步长、更多物理场的联立求解，以及同一决策问题上大量候选方案的并行搜索。GPU 利用率从"被并行喂饱"转向"被复杂度喂饱"，可能是物理引擎下一个阶段的主要特征。

把上述方向落到工程实践上，沐曦在 GPU 加速 CAE 与智能决策两方面都有长期的研究积累与工程经验，恰好覆盖求解精度、多物理场耦合与决策优化这几条路径，可为物理仿真引擎的下一阶段演进提供直接借鉴。

参考资料

[1] NVIDIA. "GTC 2025 Keynote with Jensen Huang." NVIDIA GTC, Mar. 2025. https://www.nvidia.com/gtc/

[2] S. Huang, M. Macklin, and G. State. "Announcing Newton, an Open-Source Physics Engine for Robotics Simulation." NVIDIA Technical Blog, Mar. 2025.

[3] The Newton Developers. "Newton Physics Documentation." https://newton-physics.github.io/newton/stable/

[4] OpenUSD. "USD Frequently Asked Questions." Universal Scene Description Documentation. https://openusd.org/release/usdfaq.html

[5] E. Todorov, T. Erez, and Y. Tassa. "MuJoCo: A physics engine for model-based control." 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5026–5033, 2012.

[6] V. Tsounis, G. Maloisel, C. Schumacher, R. Grandia, A. Serifi, D. Müller, C. Amevor, T. Widmer, and M. Bächer. "Kamino: GPU-based Massively Parallel Simulation of Multi-Body Systems with Challenging Topologies." arXiv preprint arXiv:2603.16536, 2026.

[7] R. Featherstone. Rigid Body Dynamics Algorithms. Springer, 2008.

[8] L. Lan, et al. "Efficient GPU Cloth Simulation with Non-distance Barriers and Subspace Reuse." ACM Transactions on Graphics, vol. 43, no. 6, Article 228, 2024.

[9] A. H. Chen, Z. Liu, Y. Yang, and C. Yuksel. "Vertex Block Descent." ACM Transactions on Graphics, vol. 43, no. 4, Article 116, 2024.

[10] C. Giles, E. Diaz, and C. Yuksel. "Augmented Vertex Block Descent." ACM Transactions on Graphics, vol. 44, no. 4, Article 90, 2025.

[11] M. Macklin, M. Müller, and N. Chentanez. "XPBD: Position-Based Simulation of Compliant Constrained Dynamics." Proceedings of the 9th International Conference on Motion in Games (MIG), pp. 49–54, 2016.

[12] R. D. Skeel and J. L. Cieśliński. "On the famous unpublished preprint 'Methods of integration which preserve the contact transformation property of the Hamilton equations' by René De Vogelaere." arXiv:2003.12268, 2020.

[13] D. Sulsky, Z. Chen, and H. L. Schreyer. "A particle method for history-dependent materials." Computer Methods in Applied Mechanics and Engineering, vol. 118, pp. 179–196, 1994.