2b:null 2c:["$","title",null,{"children":"Rotation with different initials | Chuncheng | Observable"}] 31:T420,Plot.plot({ width: width, // height: 400, aspectRatio: 1.0, x: { nice: true, domain: d3.extent(data, (d) => d.x) }, y: { nice: true, domain: d3.extent(data, (d) => d.y) }, color: { legend: true, scheme: "Turbo" }, r: { range: [0.1, 4], legend: true }, marks: [ Plot.dot(data, { x: "x", y: "y", fill: "t", r: "v", fx: "pnt" }), Plot.dot([0], { x: (d) => 0, y: (d) => 0, fill: "black" }), Plot.arrow([0], { x1: (d) => keplerEExtents.ex[0] * -100, y1: (d) => keplerEExtents.ey[0] * -100, x2: (d) => keplerEExtents.ex[0] * 100, y2: (d) => keplerEExtents.ey[0] * 100, opacity: 0.5 }), Plot.text([0], { x: (d) => 0, y: (d) => 0, text: (d) => "center", dx: 4, dy: -4, textAnchor: "start", lineAnchor: "bottom", fontSize: 20 }), Plot.text([0], { x: (d) => data[0].x, y: (d) => data[0].y, text: (d) => "start", dx: 4, dy: -4, textAnchor: "start", lineAnchor: "bottom", fontSize: 20 }) ] })32:Tfb6,**辛积分（Symplectic Integration）** 是一种专门用于模拟**哈密顿系统**（如万有引力下的轨道运动）的数值方法，其核心优势是能够长期保持系统的**能量**、**角动量**等守恒量，避免传统数值方法因误差累积导致的轨迹失真。以下是关键要点和实际应用指南： --- ### **1. 为什么需要辛积分？** - **问题背景**：在牛顿力学中，行星轨道是哈密顿系统（能量守恒）。普通数值方法（如欧拉法、Runge-Kutta）会引入能量漂移，导致轨道“花瓣状”或逐渐发散。 - **辛积分的优势**：通过保持相空间的辛结构（Symplecticity），严格模拟守恒律，即使长时间模拟也能维持轨道稳定性。 --- ### **2. 常用辛积分方法** #### **(1) Verlet 积分（蛙跳法，Leapfrog）** - **公式**： \[ \begin{cases} \mathbf{v}_{n+1/2} = \mathbf{v}_n + \frac{\Delta t}{2} \mathbf{a}(\mathbf{r}_n) \\ \mathbf{r}_{n+1} = \mathbf{r}_n + \Delta t \mathbf{v}_{n+1/2} \\ \mathbf{v}_{n+1} = \mathbf{v}_{n+1/2} + \frac{\Delta t}{2} \mathbf{a}(\mathbf{r}_{n+1}) \end{cases} \] - **特点**：简单高效，适用于保守力场（如万有引力），能量误差有界而非累积。 #### **(2) 辛欧拉法（Symplectic Euler）** - **公式**： \[ \begin{cases} \mathbf{v}_{n+1} = \mathbf{v}_n + \Delta t \mathbf{a}(\mathbf{r}_n) \\ \mathbf{r}_{n+1} = \mathbf{r}_n + \Delta t \mathbf{v}_{n+1} \end{cases} \] - **特点**：比Verlet更简单，但精度略低。 #### **(3) 高阶辛方法** - 如**Yoshida积分**或**Runge-Kutta-Nyström (RKN)** 方法，通过组合多步Verlet实现更高精度。 --- ### **3. 辛积分的实现步骤（以万有引力为例）** #### **步骤1：定义加速度** - 万有引力加速度： \[ \mathbf{a}(\mathbf{r}) = -\frac{GM}{r^3} \mathbf{r}, \quad r = \|\mathbf{r}\| \] #### **步骤2：选择辛算法（以Verlet为例）** ```python def verlet_integration(r0, v0, dt, steps): r = np.zeros((steps, 2)) # 位置轨迹 v = np.zeros((steps, 2)) # 速度轨迹 r[0] = r0 v_half = v0 + 0.5 * dt * acceleration(r0) # 半步速度 for i in range(1, steps): r[i] = r[i-1] + dt * v_half a = acceleration(r[i]) v[i] = v_half + 0.5 * dt * a v_half = v_half + dt * a # 更新半步速度 return r, v ``` #### **步骤3：验证守恒量** - **能量**：计算总能量 $ E = \frac{1}{2}mv^2 - \frac{GMm}{r} $ 是否波动在合理范围内。 - **角动量**：检查 $ \mathbf{L} = \mathbf{r} \times \mathbf{p} $ 是否恒定。 --- ### **4. 辛积分 vs 传统方法** | **特性** | **辛积分** | **显式欧拉/RK4** | |----------------|--------------------------|--------------------------| | **能量守恒** | 长期保持（有界误差） | 误差累积，能量漂移 | | **相位精度** | 轨道形状稳定 | 轨道可能发散或扭曲 | | **计算成本** | 与欧拉法相当 | RK4计算量更高 | | **适用场景** | 哈密顿系统（轨道、分子动力学） | 通用微分方程 | --- ### **5. 常见问题与调试** - **问题1：轨道仍不稳定** → 检查时间步长 $ \Delta t $ 是否过大（建议 $ \Delta t \ll \sqrt{r^3/GM} $）。 - **问题2：能量仍有漂移** → 确保加速度计算无符号错误（如 $ r^3 $ 而非 $ r^2 $）。 - **问题3：初始条件错误** → 行星初速度需满足 $ v_0 \approx \sqrt{GM/r_0} $（圆轨道）或更低（椭圆轨道）。 --- ### **6. 扩展应用** - **多体问题**：辛积分同样适用于三体问题，但需注意引力计算的数值精度。 - **相对论修正**：在强引力场中，需引入后牛顿修正的哈密顿量。 33:[] 2f:{"id":"0f5acfaac78cbb15","slug":"rotation-with-different-initials","title":"Rotation with different initials","thumbnail":"1854e1c6a93db9da4faa0b2c1e6cbad857a2ed749dbebdfb430a27f09324f379","access_level":"public","create_time":"$D2025-08-12T06:19:36.216Z","update_time":"$D2025-08-21T11:44:41.228Z","view_time":null,"first_public_time":"$D2025-08-13T03:53:55.915Z","fork":false,"stars":2,"starred":false,"owner":{"type":"team","id":"bcd65cef1cb4d2ed","login":"listenzcc","name":"Chuncheng","avatar_url":"https://avatars-next.observableusercontent.com/bb3c68e77a708a587ffc2d663bb83421f20d67991dbe6c1623c6e5d2fb73d15e","tier":"starter_2024"},"version":465,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"16328058352682a8","login":"listenzcc","name":"Chuncheng","avatar_url":"https://avatars-next.observableusercontent.com/bb3c68e77a708a587ffc2d663bb83421f20d67991dbe6c1623c6e5d2fb73d15e","tier":null}],"members":[{"type":"individual","id":"16328058352682a8","login":"listenzcc","name":"Chuncheng","avatar_url":"https://avatars-next.observableusercontent.com/bb3c68e77a708a587ffc2d663bb83421f20d67991dbe6c1623c6e5d2fb73d15e","tier":null,"role":"editor"},{"type":"team","id":"bcd65cef1cb4d2ed","login":"listenzcc","name":"Chuncheng","avatar_url":"https://avatars-next.observableusercontent.com/bb3c68e77a708a587ffc2d663bb83421f20d67991dbe6c1623c6e5d2fb73d15e","tier":"starter_2024","role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":0,"collection":null,"type":"notebook","body":{"title":"Rotation with different initials","cells":[{"id":0,"value":"# Rotation with different initials\n\nThe simulation benefits with verlet_integration. (See the doc at last for details.)","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":273,"value":"

Eng extent:

\n Method 1: ${d3.extent(data.filter(d=>d.pnt==\"1\"), (d) => d.e).map(d => d.toFixed(2)).join(' - ')}\n
\n Method 2: ${d3.extent(data.filter(d=>d.pnt==\"2\"), (d) => d.e).map(d => d.toFixed(2)).join(' - ')}\n

\n \n

","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":414,"value":"keplerEExtents = {\n const data1 = data.filter((d) => d.pnt == 2);\n return {\n ex: d3.extent(data1, (d) => d.keplerE.ex),\n ey: d3.extent(data1, (d) => d.keplerE.ey),\n e: d3.extent(data1, (d) => d.keplerE.e)\n };\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":289,"value":"\n

Angular momentum extent:

\n Method 1: ${d3.extent(data.filter(d=>d.pnt==\"1\"), (d) => d.w).map(d => d.toFixed(2)).join(' - ')}\n
\n Method 2: ${d3.extent(data.filter(d=>d.pnt==\"2\"), (d) => d.w).map(d => d.toFixed(2)).join(' - ')}\n

","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":385,"value":"

\n The simulation of rotation\n

\n The left uses the normal simulation method. \n
\n The right uses the verlet_integration method.\n

The color refers the simulation times (in seconds), and the size refers the speed (smaller is slower).

","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":7,"value":"viewof form = Inputs.form({\n angle: Inputs.range([0, 180], { label: \"Angle\", step: 1, value: 41 }),\n v: Inputs.range([0.1, 1], { label: \"V\", step: 0.01, value: 0.48 }),\n dt: Inputs.range([0.001, 0.1], { label: \"dt\", step: 0.001, value: 0.02 }),\n m: Inputs.range([0.1, 10], { label: \"m\", step: 0.1, value: 5.8 })\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":29,"value":"$31","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":404,"value":"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":10,"value":"data = {\n const pnt = {\n x: 10,\n y: 0,\n t: 0,\n vx: Math.cos((form.angle / 180) * Math.PI) * form.v,\n vy: Math.sin((form.angle / 180) * Math.PI) * form.v\n },\n pnt2 = Object.assign({}, pnt),\n data = [];\n\n for (let i = 0; i < 10000; ++i) {\n update(pnt);\n update2(pnt2);\n data.push(Object.assign({ i, pnt: \"1\" }, pnt));\n data.push(Object.assign({ i, pnt: \"2\" }, pnt2));\n }\n return data;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":13,"value":"data","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":25,"value":"function computeR(pnt) {\n return Math.sqrt(pnt.x * pnt.x + pnt.y * pnt.y);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":113,"value":"function computeV(pnt) {\n return Math.sqrt(pnt.vx * pnt.vx + pnt.vy * pnt.vy);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":108,"value":"function rotate(x, y, theta) {\n let cos = Math.cos(theta),\n sin = Math.sin(theta),\n x2 = x * cos - y * sin,\n y2 = x * sin + y * cos;\n\n return { x: x2, y: y2 };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":299,"value":"/*\nverlet_integration method.\n */\nfunction update2(pnt) {\n let { k } = metric(pnt);\n\n let r, a, ax, ay, dt;\n\n dt = form.dt;\n\n r = computeR(pnt);\n a = form.m / r / r;\n ax = (-pnt.x / r) * a;\n ay = (-pnt.y / r) * a;\n pnt.vx += ax * 0.5 * dt;\n pnt.vy += ay * 0.5 * dt;\n\n pnt.x += pnt.vx * dt;\n pnt.y += pnt.vy * dt;\n\n r = computeR(pnt);\n a = form.m / r / r;\n ax = (-pnt.x / r) * a;\n ay = (-pnt.y / r) * a;\n pnt.vx += ax * 0.5 * dt;\n pnt.vy += ay * 0.5 * dt;\n\n pnt.t += dt;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":50,"value":"function update(pnt) {\n let { k } = metric(pnt);\n\n let r, a, ax, ay, dt;\n\n dt = Math.min(form.dt, k / 100);\n\n pnt.x += pnt.vx * dt;\n pnt.y += pnt.vy * dt;\n\n r = computeR(pnt);\n\n a = form.m / r / r;\n ax = (-pnt.x / r) * a;\n ay = (-pnt.y / r) * a;\n pnt.vx += ax * dt;\n pnt.vy += ay * dt;\n pnt.t += dt;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":331,"value":"function metric(pnt) {\n let v, r, e, w, k, keplerE;\n v = computeV(pnt);\n r = computeR(pnt);\n e = v * v * 0.5 - form.m / r;\n w = pnt.vy * pnt.x - pnt.vx * pnt.y;\n k = Math.sqrt((r * r * r) / form.m);\n\n keplerE = computeKeplerE(pnt);\n\n Object.assign(pnt, { v, r, e, w, k, keplerE });\n return pnt;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":410,"value":"function computeKeplerE(pnt) {\n const { x, y, vx, vy } = pnt;\n\n const r = Math.hypot(x, y);\n const v2 = vx * vx + vy * vy;\n const h = x * vy - y * vx; // z分量\n const ex = (vy * h) / form.m - x / r;\n const ey = (-vx * h) / form.m - y / r;\n const e = Math.hypot(ex, ey);\n return { e, ex, ey };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":271,"value":"d3 = require(\"d3\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":381,"value":"$32","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"}],"files":"$W33","stdlib":"1","titling":"auto"}} 30:"$2f:body"