การอินทิเกรตแบบ Verlet ( การออกเสียงภาษาฝรั่งเศส: [vɛʁˈlɛ] ) เป็นวิธีการเชิงตัวเลขที่ใช้ในการอินทิเกรตสมการการเคลื่อนที่ของนิวตัน [ ^{1 ] มัก}ใช้ในการคำนวณ วิถี การเคลื่อนที่ของอนุภาคใน การจำลอง พลศาสตร์ โมเลกุล และกราฟิกคอมพิวเตอร์อัลกอริทึมนี้ถูกใช้ครั้งแรกในปี 1791 โดยJean Baptiste Delambreและได้รับการค้นพบใหม่หลายครั้งนับตั้งแต่นั้นมา โดยล่าสุดโดยLoup Verletในช่วงทศวรรษ 1960 เพื่อใช้ในพลศาสตร์โมเลกุล นอกจากนี้ยังถูกใช้โดยPH CowellและACC Crommelinในปี 1909 เพื่อคำนวณวงโคจรของดาวหางฮัลเลย์และโดยCarl Størmerในปี 1907 เพื่อศึกษาเส้นทางการเคลื่อนที่ของอนุภาคไฟฟ้าในสนามแม่เหล็ก (ดังนั้นจึงเรียกอีกอย่างว่าวิธีของ Størmer ) ^{[ 2 ]} ตัวรวม Verlet ให้ความเสถียรเชิงตัวเลข ที่ดี รวมถึงคุณสมบัติอื่นๆ ที่สำคัญในระบบทางกายภาพเช่นการย้อนกลับของเวลาและการรักษารูปแบบซิมเพล็กติกบนปริภูมิเฟสโดย ไม่มีต้นทุนการคำนวณเพิ่มเติมอย่างมีนัยสำคัญเมื่อเทียบกับวิธี Euler แบบง่าย

สำหรับสมการเชิงอนุพันธ์อันดับสองประเภทที่มีเงื่อนไขเริ่มต้นและสามารถหา คำตอบเชิงตัวเลขโดยประมาณณ เวลา ที่มีขนาดขั้นตอน ได้ โดยใช้วิธีต่อไปนี้:${\displaystyle {\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$${\displaystyle \mathbf {x} (t_{0})=\mathbf {x} _{0}}$${\displaystyle {\dot {\mathbf {x} }}(t_{0})=\mathbf {v} _{0}}$${\displaystyle \mathbf {x} _{n}\approx \mathbf {x} (t_{n})}$${\displaystyle t_{n}=t_{0}+n\,\Delta t}$${\displaystyle \Delta t>0}$

1. ชุด,${\textstyle \mathbf {x} _{1}=\mathbf {x} _{0}+\mathbf {v} _{0}\,\Delta t+{\tfrac {1}{2}}\mathbf {A} (\mathbf {x} _{0})\,\Delta t^{2}}$
2. สำหรับn = 1, 2, ... ทำซ้ำ$${\displaystyle \mathbf {x} _{n+1}=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {A} (\mathbf {x} _{n})\,\Delta t^{2}.}$$

สมการการเคลื่อนที่ของนิวตันสำหรับระบบทางกายภาพแบบอนุรักษ์คือ หรือแยกเป็นรายบุคคล โดยที่ $${\displaystyle {\boldsymbol {M}}{\ddot {\mathbf {x} }}(t)=F{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla V{\bigl (}\mathbf {x} (t){\bigr )},}$$$${\displaystyle m_{k}{\ddot {\mathbf {x} }}_{k}(t)=F_{k}{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla _{\mathbf {x} _{k}}V{\left(\mathbf {x} (t)\right)},}$$

• ${\displaystyle t}$ถึงเวลาแล้ว
• ${\displaystyle \mathbf {x} (t)={\bigl (}\mathbf {x} _{1}(t),\ldots ,\mathbf {x} _{N}(t){\bigr )}}$คือกลุ่มของเวกเตอร์ตำแหน่งของวัตถุ${\displaystyle N}$
• ${\displaystyle V}$คือฟังก์ชันศักย์สเกลาร์
• ${\displaystyle F}$คือค่าความชันเชิงลบของศักยภาพซึ่งให้กลุ่มของแรงที่กระทำต่ออนุภาค
• ${\displaystyle {\boldsymbol {M}}}$คือเมทริกซ์มวลซึ่งโดยทั่วไปจะเป็นเมทริกซ์แนวทแยงที่มีช่องแสดงมวลของอนุภาคแต่ละตัว${\displaystyle m_{k}}$

สมการนี้ เมื่อพิจารณาค่าฟังก์ชันศักย์ที่แตกต่างกันสามารถนำมาใช้อธิบายวิวัฒนาการของระบบทางกายภาพที่หลากหลาย ตั้งแต่การเคลื่อนที่ของโมเลกุล ที่มีปฏิสัมพันธ์กันไปจนถึงวงโคจรของดาวเคราะห์${\displaystyle V}$

หลังจากทำการแปลงเพื่อนำมวลไปไว้ทางด้านขวาและละทิ้งโครงสร้างของอนุภาคหลายตัว สมการอาจลดรูปให้ง่ายขึ้นได้ โดยใช้ฟังก์ชันเวกเตอร์ที่เหมาะสมซึ่งแสดงถึงความเร่งที่ขึ้นอยู่กับตำแหน่ง โดยทั่วไปแล้ว จะมีการกำหนดตำแหน่งเริ่มต้นและความเร็วเริ่มต้นไว้ด้วย $${\displaystyle {\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$$${\displaystyle \mathbf {A} (\mathbf {x} )}$${\displaystyle \mathbf {x} (0)=\mathbf {x} _{0}}$${\displaystyle \mathbf {v} (0)={\dot {\mathbf {x} }}(0)=\mathbf {v} _{0}}$

ในการแปลง ปัญหาค่าเริ่มต้นนี้ให้เป็นแบบไม่ต่อเนื่องและแก้ปัญหาด้วยวิธี เชิงตัวเลข จะต้องเลือก ช่วงเวลาและพิจารณาลำดับจุดสุ่มตัวอย่างงานที่ต้องทำคือการสร้างลำดับจุดที่ติดตามจุดบนเส้นทางของคำตอบที่ถูกต้อง อย่างใกล้ชิด${\displaystyle \Delta t>0}$${\displaystyle t_{n}=n\,\Delta t}$${\displaystyle \mathbf {x} _{n}}$${\displaystyle \mathbf {x} (t_{n})}$

ในขณะที่วิธีของออยเลอร์ใช้ การประมาณ ค่าผลต่างไปข้างหน้าสำหรับอนุพันธ์อันดับแรกในสมการเชิงอนุพันธ์อันดับหนึ่ง การอินทิเกรตแบบเวอร์เลต์สามารถมองได้ว่าเป็นการใช้ การประมาณ ค่าผลต่างส่วนกลางสำหรับอนุพันธ์อันดับสอง: $${\displaystyle {\begin{aligned}{\frac {\Delta ^{2}\mathbf {x} _{n}}{\Delta t^{2}}}&={\frac {{\frac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}-{\frac {\mathbf {x} _{n}-\mathbf {x} _{n-1}}{\Delta t}}}{\Delta t}}\\[6pt]&={\frac {\mathbf {x} _{n+1}-2\mathbf {x} _{n}+\mathbf {x} _{n-1}}{\Delta t^{2}}}\\&=\mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}}$$

การอินทิเกรต Verletในรูปแบบที่ใช้เป็นวิธีการ Størmer ^{[ 3 ]}ใช้สมการนี้เพื่อหาเวกเตอร์ตำแหน่งถัดไปจากสองเวกเตอร์ก่อนหน้าโดยไม่ต้องใช้ความเร็ว $${\displaystyle {\begin{aligned}\mathbf {x} _{n+1}&=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {a} _{n}\,\Delta t^{2},\\[6pt]\mathbf {a} _{n}&=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}}$$

ความสมมาตรเชิงเวลาที่มีอยู่ในวิธีการนี้ช่วยลดระดับของข้อผิดพลาดเฉพาะที่ที่เกิดขึ้นในการอินทิเกรตโดยการแบ่งส่วนย่อยโดยการกำจัดพจน์ที่มีดีกรีคี่ทั้งหมด ในที่นี้คือพจน์ที่มีดีกรีสาม ข้อผิดพลาดเฉพาะที่นี้ถูกวัดโดยการใส่ค่าที่แน่นอนลงในการวนซ้ำและคำนวณการกระจายเทย์เลอร์ณ เวลาของเวกเตอร์ตำแหน่งในทิศทางเวลาต่างๆ โดยที่คือตำแหน่ง คือความเร็วคือความเร่ง และคือค่าความเร่งกระชาก ( อนุพันธ์อันดับสามของตำแหน่งเทียบกับเวลา) ${\displaystyle \Delta t}$${\displaystyle \mathbf {x} (t_{n-1}),\mathbf {x} (t_{n}),\mathbf {x} (t_{n+1})}$${\displaystyle t=t_{n}}$${\displaystyle \mathbf {x} (t\pm \Delta t)}$$${\displaystyle {\begin{aligned}\mathbf {x} (t{+}\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}+{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}{\left(\Delta t^{4}\right)}\\[1ex]\mathbf {x} (t{-}\Delta t)&=\mathbf {x} (t)-\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}-{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}{\left(\Delta t^{4}\right)},\end{aligned}}}$$${\displaystyle \mathbf {x} }$${\displaystyle \mathbf {v} ={\dot {\mathbf {x} }}}$${\displaystyle \mathbf {a} ={\ddot {\mathbf {x} }}}$${\displaystyle \mathbf {b} ={\dot {\mathbf {a} }}={\overset {\dots }{\mathbf {x} }}}$

เมื่อรวมการขยายทั้งสองนี้เข้าด้วยกันจะได้ เราจะเห็นได้ว่าพจน์อันดับที่หนึ่งและอันดับที่สามจากการขยายอนุกรมเทย์เลอร์จะหักล้างกัน ทำให้ตัวอินทิเกรเตอร์ของเวอร์เลต์มีความแม่นยำกว่าการอินทิเกรตโดยใช้การขยายอนุกรมเทย์เลอร์แบบธรรมดาเพียงอย่างเดียวถึงหนึ่งอันดับ $${\displaystyle \mathbf {x} (t{+}\Delta t)=2\mathbf {x} (t)-\mathbf {x} (t{-}\Delta t)+\mathbf {a} (t)\Delta t^{2}+{\mathcal {O}}{\left(\Delta t^{4}\right)}.}$$

ควรระมัดระวังข้อเท็จจริงที่ว่า ค่าความเร่งในที่นี้คำนวณจากคำตอบที่ถูกต้องในขณะที่ในการวนซ้ำ ค่าความเร่งจะคำนวณที่จุดศูนย์กลางของการวนซ้ำในการคำนวณข้อผิดพลาดโดยรวม ซึ่งก็คือระยะห่างระหว่างคำตอบที่ถูกต้องกับลำดับการประมาณค่า สองพจน์นี้จะไม่หักล้างกันอย่างสมบูรณ์ ส่งผลต่อลำดับของข้อผิดพลาดโดยรวม ${\displaystyle \mathbf {a} (t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$${\displaystyle \mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n})}$

เพื่อให้เข้าใจความสัมพันธ์ระหว่างข้อผิดพลาดเฉพาะจุดและข้อผิดพลาดโดยรวม ได้ ดียิ่งขึ้น การพิจารณาตัวอย่างง่ายๆ ที่สามารถแสดงคำตอบที่ถูกต้องและคำตอบโดยประมาณในรูปสูตรที่ชัดเจนได้นั้นจะเป็นประโยชน์ ตัวอย่างมาตรฐานสำหรับงานนี้คือฟังก์ชันเลขชี้กำลัง

Consider the linear differential equation ${\displaystyle {\ddot {x}}(t)=w^{2}x(t)}$ with a constant ${\displaystyle w}$. Its exact basis solutions are ${\displaystyle e^{wt}}$ and ${\displaystyle e^{-wt}}$.

The Størmer method applied to this differential equation leads to a linear recurrence relation$${\displaystyle x_{n+1}-2x_{n}+x_{n-1}=h^{2}w^{2}x_{n},}$$ or $${\displaystyle x_{n+1}-2\left(1+{\tfrac {1}{2}}\left(wh\right)^{2}\right)x_{n}+x_{n-1}=0.}$$ It can be solved by finding the roots of its characteristic polynomial ${\displaystyle q^{2}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)q+1=0}$. These are $${\displaystyle q_{\pm }=1+{\tfrac {1}{2}}\left(wh\right)^{2}\pm wh{\sqrt {1+{\tfrac {1}{4}}\left(wh\right)^{2}}}.}$$ The basis solutions of the linear recurrence are ${\displaystyle x_{n}=q_{+}^{n}}$ and ${\displaystyle x_{n}=q_{-}^{n}}$. To compare them with the exact solutions, Taylor expansions are computed: $${\displaystyle {\begin{aligned}q_{+}&=1+{\tfrac {1}{2}}(wh)^{2}+wh\left(1+{\tfrac {1}{8}}(wh)^{2}-{\tfrac {3}{128}}(wh)^{4}+{\mathcal {O}}{\left(h^{6}\right)}\right)\\&=1+(wh)+{\tfrac {1}{2}}(wh)^{2}+{\tfrac {1}{8}}(wh)^{3}-{\tfrac {3}{128}}(wh)^{5}+{\mathcal {O}}{\left(h^{7}\right)}.\end{aligned}}}$$ The quotient of this series with the exponential ${\displaystyle e^{wh}}$ starts with ${\displaystyle 1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}{\left(h^{5}\right)}}$, so $${\displaystyle {\begin{aligned}q_{+}&=\left(1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}{\left(h^{5}\right)}\right)e^{wh}\\&=e^{-{\frac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)}\,e^{wh}.\end{aligned}}}$$ From there it follows that for the first basis solution the error can be computed as $${\displaystyle {\begin{aligned}x_{n}=q_{+}^{n}&=e^{-{\frac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}\left(h^{4}\right)}\,e^{wt_{n}}\\&=e^{wt_{n}}\left(1-{\tfrac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}{\left(h^{4}\right)}\right)\\&=e^{wt_{n}}+{\mathcal {O}}{\left(h^{2}t_{n}e^{wt_{n}}\right)}.\end{aligned}}}$$ That is, although the local discretization error is of order 4, due to the second order of the differential equation the global error is of order 2, with a constant that grows exponentially in time.

Note that at the start of the Verlet iteration at step ${\displaystyle n=1}$, time ${\displaystyle t=t_{1}=\Delta t}$, computing ${\displaystyle \mathbf {x} _{2}}$, one already needs the position vector ${\displaystyle \mathbf {x} _{1}}$ at time ${\displaystyle t=t_{1}}$. At first sight, this could give problems, because the initial conditions are known only at the initial time ${\displaystyle t_{0}=0}$. However, from these the acceleration ${\displaystyle \mathbf {a} _{0}=\mathbf {A} (\mathbf {x} _{0})}$ is known, and a suitable approximation for the position at the first time step can be obtained using the Taylor polynomial of degree two: $${\displaystyle {\begin{aligned}\mathbf {x} _{1}&=\mathbf {x} _{0}+\mathbf {v} _{0}\Delta t+{\tfrac {1}{2}}\mathbf {a} _{0}\Delta t^{2}\\&\approx \mathbf {x} (\Delta t)+{\mathcal {O}}{\left(\Delta t^{3}\right)}.\end{aligned}}}$$

The error on the first time step then is of order ${\displaystyle {\mathcal {O}}{\left(\Delta t^{3}\right)}}$. This is not considered a problem because on a simulation over a large number of time steps, the error on the first time step is only a negligibly small amount of the total error, which at time ${\displaystyle t_{n}}$ is of the order ${\displaystyle {\mathcal {O}}{\left(e^{Lt_{n}}\Delta t^{2}\right)}}$, both for the distance of the position vectors ${\displaystyle \mathbf {x} _{n}}$ to ${\displaystyle \mathbf {x} (t_{n})}$ as for the distance of the divided differences ${\displaystyle {\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}}$ to ${\displaystyle {\tfrac {\mathbf {x} (t_{n+1})-\mathbf {x} (t_{n})}{\Delta t}}}$. Moreover, to obtain this second-order global error, the initial error needs to be of at least third order.

A disadvantage of the Størmer–Verlet method is that if the time step (${\displaystyle \Delta t}$) changes, the method does not approximate the solution to the differential equation. This can be corrected using the formula^[4]$${\displaystyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+\left(\mathbf {x} _{i}-\mathbf {x} _{i-1}\right){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\Delta t_{i}^{2}.}$$

A more exact derivation uses the Taylor series (to second order) at ${\displaystyle t_{i}}$ for times ${\displaystyle t_{i+1}=t_{i}+\Delta t_{i}}$ and ${\displaystyle t_{i-1}=t_{i}-\Delta t_{i-1}}$ to obtain after elimination of ${\displaystyle \mathbf {v} _{i}}$$${\displaystyle {\frac {\mathbf {x} _{i+1}-\mathbf {x} _{i}}{\Delta t_{i}}}+{\frac {\mathbf {x} _{i-1}-\mathbf {x} _{i}}{\Delta t_{i-1}}}=\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}},}$$ so that the iteration formula becomes $${\displaystyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+(\mathbf {x} _{i}-\mathbf {x} _{i-1}){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}}\,\Delta t_{i}.}$$

ในสมการพื้นฐานของ Størmer นั้นไม่ได้ระบุความเร็วไว้อย่างชัดเจน แต่บ่อยครั้งที่จำเป็นต้องใช้ความเร็วในการคำนวณปริมาณทางกายภาพ บางอย่าง เช่นพลังงานจลน์ซึ่งอาจก่อให้เกิดความท้าทายทางเทคนิคใน การจำลอง พลศาสตร์โมเลกุลเนื่องจากพลังงานจลน์และอุณหภูมิ ณ เวลาใดเวลาหนึ่งจะไม่สามารถคำนวณได้สำหรับระบบจนกว่าจะทราบตำแหน่ง ณ เวลาใดเวลาหนึ่งข้อบกพร่องนี้สามารถแก้ไขได้โดยใช้ อัลกอริทึม ความเร็ว Verletหรือโดยการประมาณความเร็วโดยใช้พารามิเตอร์ตำแหน่งและทฤษฎีบท ค่าเฉลี่ย${\displaystyle t}$${\displaystyle t+\Delta t}$$${\displaystyle \mathbf {v} (t)={\frac {\mathbf {x} (t{+}\Delta t)-\mathbf {x} (t{-}\Delta t)}{2\Delta t}}+{\mathcal {O}}{\left(\Delta t^{2}\right)}.}$$

โปรดสังเกตว่าพจน์ความเร็วนี้ล้าหลังพจน์ตำแหน่งไปหนึ่งขั้น เนื่องจากเป็นความเร็ว ณ เวลาt ไม่ใช่t_p ซึ่งหมายความว่าเป็นค่าประมาณอันดับสองของt_p ด้วยเหตุผลเดียวกัน แต่ลดช่วงเวลาลงครึ่งหนึ่งจะได้ เป็นค่าประมาณอันดับสองของt_p โดยที่ t_p = 0 ${\displaystyle t}$${\displaystyle t+\Delta t}$${\displaystyle \mathbf {v} _{n}={\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n-1}}{2\Delta t}}}$${\displaystyle \mathbf {v} (t_{n})}$${\displaystyle \mathbf {v} _{n+{\frac {1}{2}}}={\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}}$${\displaystyle \mathbf {v} {\left(t_{n+1/2}\right)}}$${\displaystyle t_{n+{\frac {1}{2}}}=t_{n}+{\tfrac {1}{2}}\Delta t}$

เราสามารถลดช่วงเวลาลงเพื่อประมาณความเร็ว ณ เวลาดังกล่าวได้ แต่จะแลกมาด้วยความแม่นยำที่ลดลง: ${\displaystyle t+\Delta t}$$${\displaystyle \mathbf {v} (t{+}\Delta t)={\frac {\mathbf {x} (t{+}\Delta t)-\mathbf {x} (t)}{\Delta t}}+{\mathcal {O}}(\Delta t).}$$

อัลกอริทึมที่เกี่ยวข้องและใช้กันทั่วไปมากกว่าคืออัลกอริทึมความเร็ว Verlet ^{[ 5 ]}ซึ่งคล้ายกับวิธี leapfrogยกเว้นว่าความเร็วและตำแหน่งจะถูกคำนวณที่ค่าเดียวกันของตัวแปรเวลา (leapfrog ไม่ได้ทำเช่นนั้น ดังที่ชื่อบ่งบอก) วิธีนี้ใช้แนวทางที่คล้ายกัน แต่รวมความเร็วไว้อย่างชัดเจน ซึ่งแก้ปัญหาของขั้นตอนเวลาแรกในอัลกอริทึม Verlet พื้นฐาน:

$${\displaystyle {\begin{aligned}\mathbf {x} (t{+}\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\Delta t^{2},\\[6pt]\mathbf {v} (t{+}\Delta t)&=\mathbf {v} (t)+{\frac {\mathbf {a} (t)+\mathbf {a} (t{+}\Delta t)}{2}}\Delta t.\end{aligned}}}$$

สามารถแสดงให้เห็นได้ว่าข้อผิดพลาดในอัลกอริทึม Verlet แบบความเร็วมีขนาดใกล้เคียงกับในอัลกอริทึม Verlet พื้นฐาน โปรดทราบว่าอัลกอริทึมความเร็วไม่จำเป็นต้องใช้หน่วยความจำมากกว่าเสมอไป เพราะใน Verlet พื้นฐาน เราติดตามเวกเตอร์ตำแหน่งสองตัว ในขณะที่ใน Verlet แบบความเร็ว เราติดตามเวกเตอร์ตำแหน่งหนึ่งตัวและเวกเตอร์ความเร็วหนึ่งตัว รูปแบบการใช้งานมาตรฐานของอัลกอริทึมนี้คือ:

1. คำนวณ​${\displaystyle \mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)=\mathbf {v} (t)+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t}$
2. คำนวณ​${\displaystyle \mathbf {x} (t+\Delta t)=\mathbf {x} (t)+\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)\,\Delta t}$
3. ได้มาจากศักยภาพการปฏิสัมพันธ์โดยใช้.${\displaystyle \mathbf {a} (t+\Delta t)}$${\displaystyle \mathbf {x} (t+\Delta t)}$
4. คำนวณ​${\displaystyle \mathbf {v} (t+\Delta t)=\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)+{\tfrac {1}{2}}\,\mathbf {a} (t+\Delta t)\Delta t}$

อัลกอริทึมนี้ยังใช้งานได้กับช่วงเวลาที่แปรผันได้ และเหมือนกับรูปแบบ การ ผสานรวม แบบ 'เตะ-เลื่อน-เตะ' ทุกประการ

หากตัดความเร็วครึ่งก้าวออกไป อัลกอริทึมนี้สามารถย่อให้สั้นลงได้ดังนี้

1. คำนวณ​${\displaystyle \mathbf {x} (t+\Delta t)=\mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t^{2}}$
2. ได้มาจากศักยภาพการปฏิสัมพันธ์โดยใช้.${\displaystyle \mathbf {a} (t+\Delta t)}$${\displaystyle \mathbf {x} (t+\Delta t)}$
3. คำนวณ​${\displaystyle \mathbf {v} (t+\Delta t)=\mathbf {v} (t)+{\tfrac {1}{2}}\,{\bigl (}\mathbf {a} (t)+\mathbf {a} (t+\Delta t){\bigr )}\Delta t}$

อย่างไรก็ตาม โปรดทราบว่าอัลกอริทึมนี้ตั้งอยู่บนสมมติฐานว่าความเร่งขึ้นอยู่กับตำแหน่งเท่านั้นและไม่ขึ้นอยู่กับความเร็ว ${\displaystyle \mathbf {a} (t+\Delta t)}$${\displaystyle \mathbf {x} (t+\Delta t)}$${\displaystyle \mathbf {v} (t+\Delta t)}$

อาจกล่าวได้ว่าผลลัพธ์ระยะยาวของวิธี Velocity Verlet และวิธี Leapfrog นั้นดีกว่าวิธี Semi-Implicit Euler อยู่หนึ่งอันดับ อัลกอริทึมทั้งสองแทบจะเหมือนกันทุกประการ ยกเว้นการเลื่อนความเร็วไปครึ่งช่วงเวลา สามารถพิสูจน์ได้โดยการหมุนลูปข้างต้นให้เริ่มต้นที่ขั้นตอนที่ 3 แล้วสังเกตว่าพจน์ความเร่งในขั้นตอนที่ 1 สามารถกำจัดได้โดยการรวมขั้นตอนที่ 2 และ 4 เข้าด้วยกัน ความแตกต่างเพียงอย่างเดียวคือความเร็ว ณ จุดกึ่งกลางในวิธี Velocity Verlet ถือเป็นความเร็วสุดท้ายในวิธี Semi-Implicit Euler

ข้อผิดพลาดโดยรวมของวิธีการออยเลอร์ทั้งหมดมีลำดับที่หนึ่ง ในขณะที่ข้อผิดพลาดโดยรวมของวิธีการนี้คล้ายกับวิธีการจุดกึ่งกลางมีลำดับที่สอง นอกจากนี้ หากความเร่งเกิดจากแรงในระบบกลไกอนุรักษ์หรือระบบแฮมิลโทเนียนพลังงานของการประมาณจะแกว่งไปมารอบ ๆ พลังงานคงที่ของระบบที่แก้ได้อย่างแม่นยำ โดยมีขอบเขตข้อผิดพลาดโดยรวมอีกครั้งที่มีลำดับที่หนึ่งสำหรับออยเลอร์แบบกึ่งชัดเจนและลำดับที่สองสำหรับเวอร์เลต์-ลีปฟร็อก เช่นเดียวกันนี้ใช้ได้กับปริมาณอนุรักษ์อื่น ๆ ทั้งหมดของระบบ เช่นโมเมนตัมเชิงเส้นหรือเชิงมุมซึ่งจะถูกอนุรักษ์ไว้เสมอหรือเกือบจะถูกอนุรักษ์ไว้ใน ตัวรวมเชิงซิ มเพล็กติก^{[ 6 ]}

วิธี Velocity Verlet เป็นกรณีพิเศษของวิธี Newmark-betaโดยมี และ${\displaystyle \beta =0}$${\textstyle \gamma ={\tfrac {1}{2}}}$

เนื่องจาก อัลกอริทึม Velocity Verletมีประโยชน์โดยทั่วไปในการใช้งาน 3 มิติ ดังนั้นวิธีแก้ปัญหาที่เขียนด้วยภาษา C++ จึงอาจมีลักษณะดังต่อไปนี้ การรวมตำแหน่งแบบนี้จะช่วยเพิ่มความแม่นยำในการจำลองและเกม 3 มิติได้อย่างมาก เมื่อเทียบกับวิธี Euler ทั่วไป

โครงสร้างร่างกาย{Vec3d pos { 0.0 , 0.0 , 0.0 };Vec3d vel { 2.0 , 0.0 , 0.0 }; // 2 เมตร/วินาที ตามแกน xVec3d acc { 0.0 , 0.0 , 0.0 }; // ไม่มีการเร่งความเร็วในตอนแรกมวลสองเท่า= 1.0 ; // 1 กก./** * อัปเดตตำแหน่งและความเร็วโดยใช้การผสานรวม "Velocity Verlet" * @param dt DeltaTime / ขั้นตอนเวลา [เช่น: 0.01] */void update ( double dt ){Vec3d new_pos = pos + vel * dt + acc * ( dt * dt * 0.5 );Vec3d new_acc = apply_forces ();Vec3d new_vel = vel + ( acc + new_acc ) * ( dt * 0.5 );pos = new_pos ;vel = new_vel ;acc = new_acc ;}/** * หากต้องการกำหนดความเร็วให้กับวัตถุ ให้คำนวณเวกเตอร์แรงที่ต้องการแทน * และนำแรงสะสมทั้งหมดมาใช้ที่นี่ */Vec3d apply_forces () const{Vec3d new_acc = Vec3d { 0.0 , 0.0 , -9.81 }; // 9.81 m/s² ลงในแกน z// ใช้แรงอื่นๆ ตามที่ต้องการ...// หมายเหตุ: หลีกเลี่ยงการพึ่งพาตัวแปร `vel` เนื่องจาก Velocity Verlet ถือว่าความเร่งขึ้นอยู่กับตำแหน่งคืนค่าnew_acc ;}};

ข้อผิดพลาดการตัดทอนโดยรวมของวิธี Verlet คือทั้งสำหรับตำแหน่งและความเร็ว ซึ่งแตกต่างจากข้อผิดพลาดเฉพาะที่ในตำแหน่งซึ่งมีค่าเพียงตามที่ได้อธิบายไว้ข้างต้น ความแตกต่างนี้เกิดจากการสะสมของข้อผิดพลาดการตัดทอนเฉพาะที่ตลอดการวนซ้ำทั้งหมด ${\displaystyle {\mathcal {O}}{\left(\Delta t^{2}\right)}}$${\displaystyle {\mathcal {O}}{\left(\Delta t^{4}\right)}}$

ข้อผิดพลาดโดยรวมสามารถหาได้จากข้อสังเกตดังต่อไปนี้:

$${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+\Delta t){\bigr )}={\mathcal {O}}{\left(\Delta t^{4}\right)}}$$

และ

$${\displaystyle x(t_{0}+2\Delta t)=2x(t_{0}+\Delta t)-x(t_{0})+\Delta t^{2}{\ddot {x}}(t_{0}+\Delta t)+{\mathcal {O}}{\left(\Delta t^{4}\right)}.}$$

ดังนั้น

$${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+2\Delta t){\bigr )}=2\cdot \operatorname {error} {\bigl (}x(t_{0}+\Delta t){\bigr )}+{\mathcal {O}}{\left(\Delta t^{4}\right)}=3\,{\mathcal {O}}{\left(\Delta t^{4}\right)}.}$$

ในทำนองเดียวกัน:

$${\displaystyle {\begin{aligned}\operatorname {error} {\bigl (}x(t_{0}+3\Delta t){\bigl )}&=6\,{\mathcal {O}}{\left(\Delta t^{4}\right)},\\[6px]\operatorname {error} {\bigl (}x(t_{0}+4\Delta t){\bigl )}&=10\,{\mathcal {O}}{\left(\Delta t^{4}\right)},\\[6px]\operatorname {error} {\bigl (}x(t_{0}+5\Delta t){\bigl )}&=15\,{\mathcal {O}}{\left(\Delta t^{4}\right)},\end{aligned}}}$$

ซึ่งสามารถสรุปได้โดยทั่วไปดังนี้ (สามารถแสดงได้โดยการอุปมาน แต่ในที่นี้ให้ไว้โดยไม่มีการพิสูจน์):

$${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+n\Delta t){\bigr )}={\frac {n(n+1)}{2}}\,{\mathcal {O}}\left(\Delta t^{4}\right).}$$

หากเราพิจารณาข้อผิดพลาดโดยรวมในตำแหน่งระหว่างและโดยที่ จะเห็นได้ชัดว่า ${\displaystyle x(t)}$${\displaystyle x(t+T)}$${\displaystyle T=n\Delta t}$

$${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+T){\bigr )}=\left({\frac {T^{2}}{2\Delta t^{2}}}+{\frac {T}{2\Delta t}}\right){\mathcal {O}}{\left(\Delta t^{4}\right)},}$$

ดังนั้น ข้อผิดพลาดโดยรวม (สะสม) ในช่วงเวลาคงที่จึงกำหนดโดย

$${\displaystyle \operatorname {error} {\bigr (}x(t_{0}+T){\bigl )}={\mathcal {O}}{\left(\Delta t^{2}\right)}.}$$

เนื่องจากความเร็วถูกกำหนดในลักษณะที่ไม่สะสมจากตำแหน่งในตัวรวม Verlet ดังนั้นข้อผิดพลาดโดยรวมของความเร็วจึงเป็น เช่นกัน${\displaystyle {\mathcal {O}}{\left(\Delta t^{2}\right)}}$

ในการจำลองพลศาสตร์โมเลกุล ข้อผิดพลาดโดยรวมมักมีความสำคัญมากกว่าข้อผิดพลาดเฉพาะจุด ดังนั้นตัวรวมแบบเวอร์เลต์จึงถูกเรียกว่าตัวรวมอันดับสอง

ระบบที่มีอนุภาคหลายตัวและมีข้อจำกัดนั้น การแก้ปัญหาด้วยวิธีการอินทิเกรตแบบเวอร์เลต์จะง่ายกว่าวิธีการของออยเลอร์ ข้อจำกัดระหว่างจุดต่างๆ อาจเป็นศักยภาพที่จำกัดระยะห่างระหว่างจุด หรือแรงดึงดูด ซึ่งสามารถจำลองได้โดย ใช้ สปริงที่เชื่อมต่ออนุภาคเข้าด้วยกัน โดยใช้สปริงที่มีความแข็งอนันต์ จากนั้นจึงสามารถแก้ปัญหาแบบจำลองนี้ได้ด้วยอัลกอริทึมเวอร์เลต์

ในมิติเดียว ความสัมพันธ์ระหว่างตำแหน่งที่ไม่ถูกจำกัดและตำแหน่งจริงของจุดณ เวลา t โดยกำหนดระยะห่างที่ต้องการสามารถหาได้ด้วยอัลกอริทึม ${\displaystyle {\tilde {x}}_{i}^{(t)}}$${\displaystyle x_{i}^{(t)}}$${\displaystyle i}$${\displaystyle t}$${\displaystyle r}$

$${\displaystyle {\begin{aligned}d_{1}&=x_{2}^{(t)}-x_{1}^{(t)},\\[6px]d_{2}&=\|d_{1}\|,\\[6px]d_{3}&={\frac {d_{2}-r}{d_{2}}},\\[6px]x_{1}^{(t+\Delta t)}&={\tilde {x}}_{1}^{(t+\Delta t)}+{\tfrac {1}{2}}d_{1}d_{3},\\[6px]x_{2}^{(t+\Delta t)}&={\tilde {x}}_{2}^{(t+\Delta t)}-{\tfrac {1}{2}}d_{1}d_{3}.\end{aligned}}}$$

การอินทิเกรตแบบเวอร์เลต์มีประโยชน์เพราะมันเชื่อมโยงแรงกับตำแหน่งโดยตรง แทนที่จะแก้ปัญหาโดยใช้ความเร็ว

อย่างไรก็ตาม ปัญหาเกิดขึ้นเมื่อมีแรงจำกัดหลายอย่างกระทำต่ออนุภาคแต่ละตัว วิธีแก้ปัญหาอย่างหนึ่งคือการวนซ้ำผ่านทุกจุดในการจำลอง เพื่อให้ในทุกจุดนั้น การผ่อนคลายข้อจำกัดของจุดก่อนหน้าได้ถูกนำมาใช้แล้วเพื่อเร่งการแพร่กระจายของข้อมูล ในการจำลอง อาจทำได้โดยการใช้ช่วงเวลาการจำลองที่สั้น การใช้จำนวนขั้นตอนการแก้ข้อจำกัดที่คงที่ต่อช่วงเวลา หรือการแก้ข้อจำกัดจนกว่าจะตรงตามเงื่อนไขการเบี่ยงเบนที่กำหนด

เมื่อประมาณค่าข้อจำกัดในระดับท้องถิ่นถึงอันดับแรก วิธีนี้จะเหมือนกับวิธี Gauss–Seidelสำหรับเมทริกซ์ขนาด เล็ก การแยกตัวประกอบ LU นั้นเร็ว กว่าระบบขนาดใหญ่สามารถแบ่งออกเป็นกลุ่ม (ตัวอย่างเช่นตุ๊กตาผ้า แต่ละตัว  = กลุ่ม) ภายในกลุ่มจะใช้วิธี LU ระหว่างกลุ่มจะ ใช้ วิธี Gauss–Seidelโค้ดเมทริกซ์สามารถนำกลับมาใช้ใหม่ได้: ความสัมพันธ์ของแรงกับตำแหน่งสามารถประมาณค่าได้ในระดับท้องถิ่นถึงอันดับแรก และการอินทิเกรตแบบ Verlet สามารถทำให้เป็นแบบแฝงได้มากขึ้น

ซอฟต์แวร์ที่ซับซ้อน เช่น SuperLU ^{[ 7 ]}มีอยู่เพื่อแก้ปัญหาที่ซับซ้อนโดยใช้เมทริกซ์แบบเบาบาง เทคนิคเฉพาะ เช่น การใช้เมทริกซ์ (คลัสเตอร์ของเมทริกซ์) อาจใช้เพื่อแก้ไขปัญหาเฉพาะ เช่น ปัญหาของแรงที่แพร่กระจายผ่านแผ่นผ้าโดยไม่ก่อให้เกิดคลื่นเสียง^{[ 8 ]}

อีกวิธีหนึ่งในการแก้ปัญหาข้อจำกัดแบบโฮโลโนมิกคือการใช้อัลกอริธึมข้อจำกัด

วิธีหนึ่งในการตอบสนองต่อการชนกันคือการใช้ระบบแบบปรับโทษ ซึ่งโดยพื้นฐานแล้วจะใช้แรงที่กำหนดไว้กับจุดที่เกิดการสัมผัส ปัญหาของวิธีนี้คือการเลือกแรงที่ใช้ทำได้ยากมาก หากใช้แรงมากเกินไป วัตถุจะเสียสมดุล หากใช้แรงน้อยเกินไป วัตถุจะทะลุผ่านกัน อีกวิธีหนึ่งคือการใช้ปฏิกิริยาการชนแบบฉายภาพ ซึ่งจะเลือกจุดที่ก่อให้เกิดปัญหาและพยายามเคลื่อนย้ายจุดนั้นไปในระยะทางที่สั้นที่สุดเท่าที่จะเป็นไปได้เพื่อแยกออกจากวัตถุอีกชิ้น

การคำนวณแบบ Verlet จะจัดการกับความเร็วที่เกิดจากการชนในกรณีหลังโดยอัตโนมัติ อย่างไรก็ตาม โปรดทราบว่าไม่รับประกันว่าจะทำในลักษณะที่สอดคล้องกับหลักฟิสิกส์ของการชน (กล่าวคือ การเปลี่ยนแปลงโมเมนตัมไม่รับประกันว่าจะสมจริง) แทนที่จะเปลี่ยนค่าความเร็วโดยปริยาย เราจำเป็นต้องควบคุมความเร็วสุดท้ายของวัตถุที่ชนกันอย่างชัดเจน (โดยการเปลี่ยนตำแหน่งที่บันทึกไว้จากขั้นตอนเวลาที่ผ่านมา)

สองวิธีที่ง่ายที่สุดในการกำหนดความเร็วใหม่คือ การชน แบบยืดหยุ่น สมบูรณ์ และการชนแบบไม่ยืดหยุ่นกลยุทธ์ที่ซับซ้อนกว่าเล็กน้อยซึ่งให้การควบคุมมากกว่าคือการใช้สัมประสิทธิ์การคืนตัว

• สภาพของกูแรนต์–ฟรีดริชส์–ลิวี
• การเคลื่อนตัวของพลังงาน
• ตัวรวมเชิงซิมเพล็กติก
• การบูรณาการแบบ Leapfrog
• อัลกอริทึมของบีแมน

• การสาธิตการผสานรวม Verlet และโค้ดในรูปแบบ Java Applet
• ฟิสิกส์ตัวละครขั้นสูง โดย โทมัส จาคอบเซน
• ทฤษฎีการจำลองพลศาสตร์โมเลกุลถูกเก็บถาวรเมื่อวันที่ 26 เมษายน 2549 ที่Wayback Machine – ด้านล่างของหน้า
• การผสานรวม Verlet ที่พัฒนาด้วย JavaScript สมัยใหม่ – ด้านล่างของหน้า