2.5 แคลคูลัสและระเบียบวิธีเชิงตัวเลข

ทุกเฟรม เกมต้องขยับสิ่งของ: ตัวละครตกลงมา ลูกบอลเด้ง เชือกแกว่งไปมา กระสุนปืนพุ่งเป็นเส้นโค้งในอากาศ เบื้องหลังทั้งหมดนี้คือคำถามเดิมที่ถามซ้ำ ๆ ว่า ถ้ารู้ว่าตอนนี้สิ่งของอยู่ตรงไหนและเคลื่อนที่ยังไง มันควรจะอยู่ตรงไหนในอีกเสี้ยววินาทีข้างหน้า? การตอบคำถามนี้ซ้ำ ๆ หกสิบครั้งขึ้นไปต่อวินาที คือสิ่งที่บทนี้ว่าด้วย คุณรู้เรื่อง struct, pointer, และ loop มาจากบทก่อนหน้าแล้ว บทนี้เพิ่มคณิตศาสตร์เข้ามานิดหน่อย (แตะ calculus แบบเบา ๆ ซึ่งเป็นสาขาคณิตศาสตร์ที่ว่าด้วยการเปลี่ยนแปลง) และแสดงให้เห็นว่าเกมเอนจินแปลงคณิตศาสตร์นั้นให้กลายเป็น loop ที่รันทุกเฟรมได้ยังไง

ไม่ต้องกังวลถ้า calculus ฟังดูน่ากลัว คุณไม่ต้องแก้สมการด้วยมือ สิ่งที่ต้องรู้มีแค่สองไอเดีย — derivative (อนุพันธ์) กับ integral (อินทิกรัล) — และทักษะหนึ่งอย่าง คือแปลงมันให้เป็น loop โค้ดสั้น ๆ ทุกอย่างในบทนี้สร้างต่อจากสองไอเดียนี้

1. Derivative: ตอนนี้มันเปลี่ยนแปลงเร็วแค่ไหน

Derivative (อนุพันธ์) ตอบคำถาม "ปริมาณนี้กำลังเปลี่ยนแปลงเร็วแค่ไหน ณ ขณะนี้เป๊ะ ๆ" มันคืออัตราการเปลี่ยนแปลง ถ้า x คือตำแหน่งของวัตถุ derivative ของมัน dx/dt ("dx ต่อ dt" คือการเปลี่ยนแปลงของ x ต่อการเปลี่ยนแปลงเวลา t ที่เล็กมาก ๆ) คือ velocity (ความเร็ว) ของมัน ถ้า v คือความเร็ว derivative ของมัน dv/dt คือ acceleration (ความเร่ง) — ความเร็วเองกำลังเปลี่ยนเร็วแค่ไหน

ไม่ต้องใช้สัญลักษณ์ calculus ก็รู้สึกได้ ถ้าไมล์วัดระยะทางของรถอ่านได้ 100 กม. ตอนบ่ายโมง และ 160 กม. ตอนบ่ายสอง ความเร็วเฉลี่ยตลอดชั่วโมงนั้นคือ 60 กม./ชม. — นั่นคือ (ระยะทางที่เปลี่ยน) / (เวลาที่เปลี่ยน) derivative ก็ไอเดียเดียวกัน แต่ย่อให้เล็กลง: แทนที่จะเป็นทั้งชั่วโมง คุณย่อช่วงเวลาให้เล็กลงเรื่อย ๆ จนมันกลายเป็นคำถามเกี่ยวกับขณะเดียว มาดูการย่อนี้เกิดขึ้นในโค้ดกัน โดยใช้ฟังก์ชันตำแหน่ง x(t) = t^2 (วัตถุที่เร่งความเร็วขึ้นเรื่อย ๆ อย่างสม่ำเสมอ เหมือนอะไรบางอย่างที่เร่งขึ้นภายใต้แรงขับคงที่):

#include <cstdio>

double position(double t) { return t * t; }   // example: x(t) = t^2

int main() {
    double t = 2.0;                 // check the rate of change AT t = 2
    double steps[] = {1.0, 0.1, 0.01, 0.0001};

    for (double dt : steps) {
        double avgV = (position(t + dt) - position(t)) / dt;  // slope of a nearby secant
        printf("dt = %8.4f  average velocity = %.6f\n", dt, avgV);
    }
}

Output:

dt =   1.0000  average velocity = 5.000000
dt =   0.1000  average velocity = 4.100000
dt =   0.0100  average velocity = 4.010000
dt =   0.0001  average velocity = 4.000100

เมื่อ dt (ช่องว่างเวลาเล็ก ๆ) หดเข้าใกล้ศูนย์ ความเร็วเฉลี่ยจะเข้าใกล้ 4.0 ค่าที่มันกำลังมุ่งเข้าหานั้น — คือ derivative ที่แท้จริง ณ t = 2 สำหรับ x(t) = t^2 calculus ให้สูตรลัดสำหรับ derivative ไว้คือ dx/dt = 2t แทนค่า t = 2 จะได้ 4 พอดี ตรงกับที่ dt ที่หดลงเรื่อย ๆ ลู่เข้าไปหา ในเกมคุณแทบไม่ต้องใช้สูตรลัดนี้เลย — physics ให้ force และ acceleration มาตรง ๆ อยู่แล้ว — แต่นี่คือไอเดียที่ซ่อนอยู่หลังประโยค "velocity คือ derivative ของ position" ทุกประโยคที่คุณจะเจอ

2. Integral: velocity สะสมกลายเป็น position

Integral (อินทิกรัล) คือการย้อนกลับ: แทนที่จะถามว่า "สิ่งนี้เปลี่ยนแปลงเร็วแค่ไหน" มันถามว่า "ถ้ารู้อัตราการเปลี่ยนแปลงในทุกขณะ ปริมาณนี้สะสมรวมได้เท่าไหร่" มันคือ accumulation (การสะสม) ถ้ารู้ velocity ของวัตถุในทุกขณะ การอินทิเกรต velocity นั้นตลอดเวลาจะบอกว่ามันเดินทางไปได้ไกลแค่ไหน — คือ position ของมัน ถ้ารู้ acceleration ในทุกขณะ การอินทิเกรตมันจะบอก velocity

ลองนึกภาพกราฟ velocity เทียบกับเวลา ระยะทางที่เดินทางระหว่างสองช่วงเวลาคือ พื้นที่ใต้กราฟนั้น คุณประมาณพื้นที่นี้ได้โดยหั่นเวลาเป็นแถบบาง ๆ แล้วสำหรับแต่ละแถบวาดสี่เหลี่ยมผืนผ้าบาง ๆ สูง v กว้าง dt — พื้นที่ของมัน v * dt คือระยะทางโดยประมาณที่เดินทางในเสี้ยวเวลานั้น รวมพื้นที่สี่เหลี่ยมทั้งหมดเข้าด้วยกัน คุณก็ได้ประมาณค่า integral แล้ว วิธีนี้เรียกว่า Riemann sum และมันคือสิ่งที่เกมทำทุกเฟรมเป๊ะ ๆ เพียงแต่ทำทีละสี่เหลี่ยม ทีละเฟรม

นี่คือไอเดียนั้นในรูปโค้ด: วัตถุที่มี acceleration คงที่ a เริ่มจากหยุดนิ่ง สะสมทีละช่วงเวลาเล็ก ๆ dt เราเทียบผลลัพธ์ที่สะสมได้กับสูตรปิด (closed-form) จาก physics คือ x = 0.5 * a * t^2:

#include <cstdio>

int main() {
    double a = 2.0;      // constant acceleration
    double T = 3.0;      // total time
    double dt = 0.001;   // tiny time slices

    double v = 0.0, x = 0.0;
    for (double t = 0.0; t < T; t += dt) {
        x += v * dt;      // accumulate position from velocity (one thin slice)
        v += a * dt;      // accumulate velocity from acceleration
    }

    double exactX = 0.5 * a * T * T;   // closed-form: x = 1/2 a t^2
    printf("numeric x  = %.4f\n", x);
    printf("exact x    = %.4f\n", exactX);
}

Output:

numeric x  = 9.0030
exact x    = 9.0000

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

3. จากคณิตศาสตร์ต่อเนื่อง สู่ game loop

ตัวอย่างข้างบนรู้ acceleration ล่วงหน้าเป็นสูตรที่ชัดเจน เกมจริงไม่มีความหรูหราแบบนั้น acceleration ของตัวละครขึ้นอยู่กับ input ของผู้เล่น, gravity, สิ่งที่มันกำลังชนอยู่, และ spring หรือ force ใดก็ตามที่กำลังดึงมันอยู่ตอนนี้ — มันเปลี่ยนไปจากเฟรมหนึ่งไปอีกเฟรมได้ในแบบที่ไม่มีสูตรไหนจับได้ ดังนั้นแทนที่จะแก้สมการที่แท้จริง เกมเอนจินจะทำสูตรเล็ก ๆ ซ้ำ ๆ ทุกเฟรม สูตรที่ทำซ้ำนี้เรียกว่า numerical integration (การอินทิเกรตเชิงตัวเลข — numerical แปลว่าใช้เลขคณิตแบบก้าวทีละขั้นโดยประมาณ แทนที่จะใช้สูตรที่แท้จริง)

วัตถุทุกตัวที่เกมจำลองจะเก็บ state (สถานะ) เล็ก ๆ ไว้: อย่างน้อยก็ position x และ velocity v ทุกเฟรม เอนจินจะคำนวณ acceleration a ปัจจุบันจาก force ที่มีผลอยู่ แล้วใช้ stepping rule (กฎการก้าว) เพื่อแปลง (x, v, a) ที่ "ตอนนี้" ให้กลายเป็น (x, v) ที่ "ตอนนี้ + dt" ที่เหลือของบทนี้ทั้งหมดว่าด้วย stepping rule ตัวเดียวนี้ — ปรากฏว่ามีสูตรหลายแบบ และมันไม่ได้ทำงานเหมือนกันหมด แม้ว่าในโค้ดจะดูเกือบเหมือนกันก็ตาม

เพื่อเปรียบเทียบสูตรต่าง ๆ อย่างยุติธรรม ตัวอย่างทั้งหมดตั้งแต่นี้ไปจะใช้เคสทดสอบเดียวกัน: มวลติดกับสปริง ไม่มีแรงเสียดทาน สปริงดึงกลับเข้าหาจุดกึ่งกลางด้วยแรงที่แปรผันตรงกับระยะที่มันถูกยืดออก — กฎของฮุก (Hooke's law), F = -k * x โดย k คือความแข็งของสปริง Acceleration คือ a = F / m = -(k / m) * x สปริงที่ไม่มีแรงเสียดทานควรแกว่งตลอดไปที่ความสูงเดิม ไม่ได้พลังงานเพิ่มหรือเสียพลังงานไปเลย — ซึ่งทำให้มันเป็นวิธีที่สมบูรณ์แบบในการจับ stepping rule ที่แอบรั่วพลังงานเข้าหรือออก

4. Explicit (Forward) Euler: กฎการก้าวที่ง่ายที่สุด

วิธีตรงที่สุดในการแปลง "position กับ velocity ตอนนี้" ให้กลายเป็น "position กับ velocity อีกขณะข้างหน้า" เรียกว่า explicit Euler (เรียกอีกอย่างว่า forward Euler ตั้งชื่อตามนักคณิตศาสตร์ Leonhard Euler) แต่ละก้าว คำนวณ acceleration จาก state ปัจจุบัน แล้วขยับทั้ง velocity และ position ไปข้างหน้าโดยใช้ค่าที่มีอยู่ ก่อน ก้าวนั้น:

นี่คือมันตอนก้าวสปริงไร้แรงเสียดทานของเรา (k = 4, m = 1, dt = 0.1 เริ่มจากยืดไปที่ x = 1 อยู่นิ่ง) นอกจาก position กับ velocity เราพิมพ์ energy (พลังงาน) รวมของสปริงด้วย (0.5*k*x^2 + 0.5*m*v^2 คือผลรวมของพลังงานสะสมในสปริงกับพลังงานจลน์) — physics บอกว่าตัวเลขนี้ควรคงที่เป๊ะ ๆ เพราะไม่มีอะไรเอาพลังงานออกหรือเติมเข้าไป:

#include <cstdio>

int main() {
    const double k = 4.0, m = 1.0, dt = 0.1;
    double x = 1.0, v = 0.0;

    printf("step     x        v      energy\n");
    for (int i = 0; i <= 6; i++) {
        double energy = 0.5*k*x*x + 0.5*m*v*v;
        printf("%3d   %7.4f  %7.4f  %6.4f\n", i, x, v, energy);
        double a = -k/m * x;
        double xNew = x + v*dt;   // explicit Euler: uses OLD v
        double vNew = v + a*dt;
        x = xNew;
        v = vNew;
    }
}

Output:

step     x        v      energy
  0    1.0000   0.0000  2.0000
  1    1.0000  -0.4000  2.0800
  2    0.9600  -0.8000  2.1632
  3    0.8800  -1.1840  2.2497
  4    0.7616  -1.5360  2.3397
  5    0.6080  -1.8406  2.4333
  6    0.4239  -2.0838  2.5306

ตัวเลข position กับ velocity ดูสมเหตุสมผลดี — สปริงกำลังแกว่งกลับ ตรงตามที่คาดไว้เป๊ะ แต่ดูที่คอลัมน์ energy สิ: มันเริ่มที่ 2.0000 แล้วไต่ขึ้นทุกก้าว โดยไม่มีอะไรป้อนเข้าไปเลย นี่ไม่ควรเกิดขึ้น และมันคือหัวข้อของส่วนถัดไป

5. ทำไม forward Euler ถึงระเบิด

Explicit Euler ใช้ position หรือ velocity จากก่อนก้าวนั้นเสมอในการขยับสิ่งของไปข้างหน้า รายละเอียดข้อนี้แปลว่ามันมองย้อนหลังนิดหน่อยอยู่เสมอ ใช้ข้อมูลที่เก่าไปแล้ว สำหรับระบบที่แกว่งไปมาแบบสปริง ความล่าช้านี้มีผลลัพธ์ที่เฉพาะเจาะจงและน่าเกลียดมาก: แต่ละก้าวจะ overshoot (เลยจุดที่ควรจะเป็น) เล็กน้อย และทุก overshoot จะเติมพลังงานเล็ก ๆ ที่ไม่เคยมีอยู่จริงเข้าไป ถ้าก้าวต่อไปเรื่อย ๆ พลังงานก็จะสะสมเพิ่มขึ้นเรื่อย ๆ เรื่องนี้มักถูกเรียกว่าซิมูเลชัน "ระเบิด" หรือ "ได้พลังงานเพิ่ม" และถ้าก้าวมากพอ มันจะทำให้ position พุ่งไปเป็นค่าที่ไร้สาระ

วิธีที่สะอาดที่สุดในการเห็นเรื่องนี้คือพล็อต velocity เทียบกับ position — เรียกว่า phase space สำหรับสปริงไร้แรงเสียดทานจริง กราฟนั้นจะเป็นวงปิด วนซ้ำตลอดไป เพราะพลังงานไม่เคยเปลี่ยน Explicit Euler กลับวาดเกลียวที่ค่อย ๆ กว้างขึ้นแทน:

6. Semi-implicit (symplectic) Euler: ทางแก้บรรทัดเดียว

Semi-implicit Euler (เรียกอีกอย่างว่า symplectic Euler) เปลี่ยนแค่อย่างเดียว: มันอัปเดต velocity ก่อน แล้วใช้ velocity ใหม่ ในการอัปเดต position แทนที่จะใช้ตัวเก่า

นั่นคือการเปลี่ยนแปลงทั้งหมด — สลับลำดับสองบรรทัด นี่คือสปริงตัวเดิมเป๊ะ ตัวเลขเดิมเป๊ะ มีแค่การสลับนั้น:

#include <cstdio>

int main() {
    const double k = 4.0, m = 1.0, dt = 0.1;
    double x = 1.0, v = 0.0;

    printf("step     x        v      energy\n");
    for (int i = 0; i <= 6; i++) {
        double energy = 0.5*k*x*x + 0.5*m*v*v;
        printf("%3d   %7.4f  %7.4f  %6.4f\n", i, x, v, energy);
        double a = -k/m * x;
        v = v + a*dt;      // update velocity FIRST
        x = x + v*dt;      // then use the NEW velocity
    }
}

Output:

step     x        v      energy
  0    1.0000   0.0000  2.0000
  1    0.9600  -0.4000  1.9232
  2    0.8816  -0.7840  1.8618
  3    0.7679  -1.1366  1.8254
  4    0.6236  -1.4438  1.8199
  5    0.4542  -1.6932  1.8462
  6    0.2667  -1.8749  1.9000

ตอนนี้ energy สั่นเล็กน้อย (2.00 -> 1.92 -> 1.86 -> 1.83 -> 1.82 -> 1.85 -> 1.90) แต่ไม่วิ่งหนีไปไหน — มันอยู่ในช่วงแคบ ๆ รอบค่าจริง และจะแกว่งอยู่รอบ ๆ นั้นตลอดไปแทนที่จะไต่ขึ้นเรื่อย ๆ ใน phase space มันอยู่บนวงเดิม (เกือบ) แทนที่จะเป็นเกลียวที่กว้างขึ้น:

วิธีที่มีคุณสมบัติ "พลังงานอยู่ในขอบเขต ไม่วิ่งหนี" แบบนี้เรียกว่า symplectic การเปลี่ยนบรรทัดเดียวนี้คือเหตุผลว่าทำไม semi-implicit Euler ไม่ใช่ explicit Euler ถึงเป็นสิ่งที่เอนจินจริงใช้สำหรับ physics: Rigidbody ของ Unity, Box2D, และ physics engine ของเกมส่วนใหญ่อัปเดต velocity ก่อนแล้วค่อย position ทีหลัง เพราะความเสถียรนี้โดยเฉพาะ มันใช้เลขคณิตต่อก้าวเท่ากับ explicit Euler เป๊ะ — คุณได้ซิมูเลชันที่เสถียรกว่ามากแบบฟรี ๆ แค่เลือกลำดับสองบรรทัดให้ถูก

7. Verlet integration: ข้ามการเก็บ velocity ไปเลย

Verlet integration (จาก Loup Verlet นักฟิสิกส์ที่ใช้มันจำลองโมเลกุล) ใช้วิธีที่ต่างออกไป: แทนที่จะเก็บ velocity มันจำ position ก่อนหน้า กับ position ปัจจุบัน ของวัตถุไว้ แล้วก้าวไปข้างหน้าโดยใช้ทั้งสองค่านั้นบวกกับ acceleration ปัจจุบัน:

#include <cstdio>

int main() {
    const double k = 4.0, m = 1.0, dt = 0.1;
    double x = 1.0, v = 0.0;              // start the same as the earlier examples
    double xPrev = x - v*dt;              // reconstruct a starting "previous position"

    printf("step     x      energy(approx)\n");
    for (int i = 0; i <= 6; i++) {
        double vApprox = (x - xPrev) / dt;  // velocity is IMPLIED, not stored
        double energy = 0.5*k*x*x + 0.5*m*vApprox*vApprox;
        printf("%3d   %7.4f   %6.4f\n", i, x, energy);
        double a = -k/m * x;
        double xNext = 2*x - xPrev + a*dt*dt;
        xPrev = x;
        x = xNext;
    }
}

Output:

step     x      energy(approx)
  0    1.0000   2.0000
  1    0.9600   1.9232
  2    0.8816   1.8618
  3    0.7679   1.8254
  4    0.6236   1.8199
  5    0.4542   1.8462
  6    0.2667   1.9000

สังเกตว่าตัวเลขพวกนี้เหมือนกับ semi-implicit Euler เป๊ะ นั่นไม่ใช่เรื่องบังเอิญ — ถ้าคุณขยายสูตร Verlet โดยใช้ x_prev = x - v*dt มันจะย่อลงเหลือสองบรรทัดเดียวกันกับ semi-implicit Euler เป๊ะ พวกมันเป็นญาติใกล้ชิดกัน: ทั้งคู่อยู่ในตระกูล "symplectic" เดียวกัน ทั้งคู่รักษาพลังงานให้อยู่ในขอบเขตแทนที่จะรั่วออกไป ข้อได้เปรียบจริง ๆ ของ Verlet ปรากฏที่อื่น: constraints (ข้อจำกัด)

ทำไม Verlet ถึงเป็นตัวเลือกหลักสำหรับผ้าและ ragdoll

ผ้า เชือก และ ragdoll มักถูกจำลองเป็นตารางหรือโซ่ของจุดที่เชื่อมกันด้วย "ไม้" แข็งที่มีความยาวพักคงที่ เพราะ Verlet แตะแค่ position เท่านั้น (ไม่มี velocity ที่ต้องคอย sync) คุณขยับทุกจุดไปข้างหน้าด้วยสูตรข้างบนได้ แล้วค่อยดันจุดตรง ๆ ให้ตรงตาม distance constraint ("สองจุดนี้ต้องห่างกันพอดี 1 หน่วย") โดยไม่ต้องกังวลเรื่องอัปเดต velocity แยกต่างหาก — velocity ที่แฝงอยู่ของก้าวถัดไปจะออกมาเองจากตำแหน่งที่จุดนั้นลงเอย วิธีนี้บางทีเรียกว่า position-based dynamics

#include <cstdio>
#include <cmath>

struct Point { float x, y; };

void satisfyDistance(Point& a, Point& b, float restLength) {
    float dx = b.x - a.x, dy = b.y - a.y;
    float dist = std::sqrt(dx*dx + dy*dy);
    float diff = (dist - restLength) / dist;
    a.x += dx * 0.5f * diff;   // pull a toward (or push away from) b
    a.y += dy * 0.5f * diff;
    b.x -= dx * 0.5f * diff;   // and the same for b, in the opposite direction
    b.y -= dy * 0.5f * diff;
}

int main() {
    Point a{0.0f, 0.0f};
    Point b{3.0f, 0.0f};    // stretched too far apart
    float rest = 1.0f;      // the "stick" between them should be length 1

    float before = b.x - a.x;
    satisfyDistance(a, b, rest);
    float after = b.x - a.x;

    printf("distance before = %.2f\n", before);
    printf("a = (%.2f, %.2f)   b = (%.2f, %.2f)\n", a.x, a.y, b.x, b.y);
    printf("distance after  = %.2f\n", after);
}

Output:

distance before = 3.00
a = (1.00, 0.00)   b = (2.00, 0.00)
distance after  = 1.00

สองจุดเริ่มห่างกัน 3 หน่วย แล้วถูกดึงเข้าหากันคนละครึ่ง ลงเอยที่ห่างกันพอดี 1 หน่วย — คือความยาวพัก ซิมูเลชันผ้าจริง ๆ จะรัน Verlet position step แล้ววนลูปทุก stick constraint แล้ว apply การดันนี้สักหลายรอบ (รอบเดียวมักไม่พอเมื่อไม้หลายอันใช้จุดร่วมกัน) ซึ่งเป็นเหตุผลที่ผ้าในเกมมักดูค่อนข้าง "หลวม" หรือใช้เวลาสักเฟรมสองเฟรมกว่าจะนิ่ง — มันกำลังแก้การดันเล็ก ๆ พวกนี้แบบวนซ้ำ (iterative) อยู่

8. RK4: แม่นยำกว่ามาก แต่มีต้นทุนจริง

RK4 (Runge-Kutta อันดับสี่ ตั้งชื่อตามนักคณิตศาสตร์สองคน) เป็น stepping rule ที่แม่นยำกว่า แทนที่จะสุ่มตัวอย่าง slope (คือ acceleration กับ velocity) แค่ครั้งเดียวต่อก้าว มันสุ่มตัวอย่าง สี่ ครั้ง — ตอนเริ่ม สองครั้งใกล้จุดกึ่งกลาง และตอนจบ — แล้วผสมตัวอย่างทั้งสี่เข้าด้วยกันด้วยน้ำหนักคงที่:

#include <cstdio>

int main() {
    const double k = 4.0, m = 1.0, dt = 0.1;
    double x = 1.0, v = 0.0;

    auto accel = [&](double xx) { return -k/m * xx; };

    printf("step     x        v      energy\n");
    for (int i = 0; i <= 6; i++) {
        double energy = 0.5*k*x*x + 0.5*m*v*v;
        printf("%3d   %7.4f  %7.4f  %6.4f\n", i, x, v, energy);

        double x1 = x,            v1 = v;
        double a1 = accel(x1);
        double x2 = x + v1*dt/2,  v2 = v + a1*dt/2;
        double a2 = accel(x2);
        double x3 = x + v2*dt/2,  v3 = v + a2*dt/2;
        double a3 = accel(x3);
        double x4 = x + v3*dt,    v4 = v + a3*dt;
        double a4 = accel(x4);

        x = x + (dt/6.0)*(v1 + 2*v2 + 2*v3 + v4);
        v = v + (dt/6.0)*(a1 + 2*a2 + 2*a3 + a4);
    }
}

Output:

step     x        v      energy
  0    1.0000   0.0000  2.0000
  1    0.9801  -0.3973  2.0000
  2    0.9211  -0.7788  2.0000
  3    0.8253  -1.1293  2.0000
  4    0.6967  -1.4347  2.0000
  5    0.5403  -1.6829  2.0000
  6    0.3624  -1.8641  2.0000

คอลัมน์ energy อ่านได้ 2.0000 ตลอดทาง — ที่ความละเอียดนี้ RK4 ไม่รั่วหรือสั่นให้เห็นเลย มันแม่นยำกว่า Euler ทั้งสองแบบมากต่อก้าว แล้วทำไมแทบไม่มีใครใช้ RK4 สำหรับ physics เกมแบบเรียลไทม์? ดูที่โค้ดสิ: ทุกก้าวเรียก accel สี่ครั้ง ไม่ใช่ครั้งเดียว สำหรับฉากที่มีวัตถุไม่กี่ตัวนั่นไม่มีความหมายอะไร แต่สำหรับ physics body หลายพันตัว งาน acceleration สี่เท่าต่อตัวต่อเฟรมคือต้นทุนจริง และเกมต้องจ่ายต้นทุนนั้นทุกเฟรมตลอดไป ไม่ใช่แค่ครั้งเดียว การสั่นของ energy เล็ก ๆ ที่สม่ำเสมอของ semi-implicit Euler (หัวข้อ 6) คือราคาที่เกมเกือบทุกเกมยินดีจ่ายแลกกับงานแค่หนึ่งในสี่

9. ทั้งสี่วิธีเทียบกัน บนสปริงตัวเดิม

ตัวอย่างข้างบนทั้งหมดใช้สปริงตัวเดิมเป๊ะ (k = 4, m = 1, dt = 0.1, เริ่มที่ x = 1, v = 0) เอาคอลัมน์ energy จากแต่ละตัวมาเรียงกันจะเล่าเรื่องทั้งหมดในตารางเดียว จำไว้ว่า: energy จริงควรอยู่ที่ 2.0000 เป๊ะตลอดไป เพราะไม่มีอะไรเอาพลังงานออกหรือเติมเข้าไปในสปริงไร้แรงเสียดทาน

สามบทเรียนซ่อนอยู่ในตารางเดียวนี้ อย่างแรก "โค้ดที่ดูแม่นยำกว่า" ไม่เหมือนกับ "โค้ดที่เสถียร" — explicit Euler กับ semi-implicit Euler ทำเลขคณิตต่อก้าวเท่ากันเป๊ะ แต่ตัวหนึ่งระเบิดอีกตัวไม่ระเบิด ลำดับสองบรรทัดคือความต่างทั้งหมด อย่างที่สอง Verlet กับ semi-implicit Euler ทำงานเหมือนกันสำหรับเคสง่าย ๆ นี้ นั่นคือเหตุผลว่าทำไมทั้งสองตัวเป็นค่าเริ่มต้นที่สมเหตุสมผล ขึ้นอยู่กับว่าคุณคิดในรูป velocity หรือ position ได้ง่ายกว่ากัน อย่างที่สาม RK4 แม่นยำกว่าจริง ๆ แต่ความแม่นยำนั้นไม่ฟรี — คุณกำลังแลกเวลา CPU ที่ต้องจ่ายทุกเฟรมตลอดไป กับความแม่นยำที่ gameplay ส่วนใหญ่ไม่มีทางสังเกตเห็น

10. Fixed timestep + accumulator: physics ต้องไม่ขึ้นกับ frame rate

Stepping rule ทุกตัวข้างบนรับ dt เป็น input สิ่งที่ดูเป็นเรื่องธรรมดาคือวัดว่าเฟรมล่าสุดใช้เวลานานแค่ไหน แล้วส่งตัวเลขนั้นตรง ๆ เป็น dt นี่คือกับดัก เวลาเฟรมไม่เคยนิ่งสมบูรณ์แบบ — มันสั่นไปตามความซับซ้อนของฉาก งาน OS เบื้องหลัง และ lag spike ตรง ๆ ถ้า dt คือเวลาเฟรมดิบ ๆ gameplay เดียวกันเป๊ะจะให้ผลลัพธ์ physics ต่างกันบนเครื่องเร็วเทียบกับเครื่องช้า และผลลัพธ์ยังต่างกันได้จากการรันครั้งหนึ่งไปอีกครั้งบนเครื่องเดียวกัน แย่กว่านั้น หัวข้อ 5 แสดงให้เห็นว่า dt ที่ใหญ่ขึ้นทำให้ error ของการก้าวแบบ explicit โตเร็วขึ้น — ดังนั้น lag spike ครั้งเดียวอาจทำให้ physics ทำงานต่างไปมากในหนึ่งเฟรม หรือในกรณีร้ายอาจทำให้วัตถุที่เคลื่อนที่เร็วทะลุกำแพงบาง ๆ ไปเลย เพราะมันกระโดดครั้งใหญ่ครั้งเดียวแทนที่จะเป็นก้าวเล็ก ๆ หลายก้าว

ทางแก้คือแยกเวลา physics ออกจากเวลา render ให้ขาดจากกันเลย เลือก dt คงที่ตัวเดียวสำหรับ physics (ทั่วไปคือ 1/60 วินาที) แล้วเก็บ accumulator ไว้ — ยอดรวมที่วิ่งอยู่ของ "เวลาจริงที่ค้างซิมูเลชันอยู่" ทุกเฟรม บวกเวลาเฟรมที่วัดได้เข้า accumulator แล้วรัน physics step เป็นลูป กินทีละ fixedDt จนกว่าจะเหลือไม่พอสำหรับอีกก้าวเต็ม ๆ เฟรมที่ช้าแค่ทำให้มี physics step รันเพิ่มมากขึ้นในเฟรมนั้น เพื่อให้ซิมูเลชันตามทัน — แต่ physics step แต่ละตัวใช้ dt เดิมเป๊ะเสมอ ดังนั้นตัวซิมูเลชันเองจึงกำหนดผลได้แน่นอน (deterministic) และไม่ขึ้นกับ frame rate แม้ว่าการ render จะไม่เป็นแบบนั้นก็ตาม

#include <cstdio>

void stepPhysics() {
    // moves every body forward by exactly fixedDt -- always the same dt
}

int main() {
    const double fixedDt = 1.0 / 60.0;     // ~0.016667s, ALWAYS this value
    double accumulator = 0.0;

    double frameTimes[] = {0.021, 0.018, 0.052, 0.017}; // "measured" real frame time

    for (double frameTime : frameTimes) {
        if (frameTime > 0.25) frameTime = 0.25;   // clamp a huge spike (see below)
        accumulator += frameTime;
        printf("frame took %.3fs -> accumulator = %.5f\n", frameTime, accumulator);

        int steps = 0;
        while (accumulator >= fixedDt) {
            stepPhysics();
            accumulator -= fixedDt;
            steps++;
        }
        double alpha = accumulator / fixedDt;   // 0..1, for render interpolation
        printf("  ran %d physics step(s), leftover = %.5f, alpha = %.2f\n",
               steps, accumulator, alpha);
    }
}

Output:

frame took 0.021s -> accumulator = 0.02100
  ran 1 physics step(s), leftover = 0.00433, alpha = 0.26
frame took 0.018s -> accumulator = 0.02233
  ran 1 physics step(s), leftover = 0.00567, alpha = 0.34
frame took 0.052s -> accumulator = 0.05767
  ran 3 physics step(s), leftover = 0.00767, alpha = 0.46
frame took 0.017s -> accumulator = 0.02467
  ran 1 physics step(s), leftover = 0.00800, alpha = 0.48

ดูเฟรมที่สาม: มันใช้เวลา 0.052 วินาทีซึ่งนานผิดปกติ (lag spike) ลูปเลยรัน physics step สามครั้งติดกันเพื่อตามให้ทัน แทนที่จะเป็นก้าวใหญ่ก้าวเดียวด้วย dt ที่บวมขึ้น physics ยังคงเสถียรเท่ากับตอน frame rate นิ่ง ๆ เพียงแค่รันบ่อยขึ้นในเฟรมนั้นเฟรมเดียว ค่า accumulator ที่เหลืออยู่หลังลูป (ไม่เคยเต็ม fixedDt) ถูกใช้เป็น alpha เศษส่วนสำหรับ interpolate ตำแหน่งตอน render เพื่อให้ภาพลื่นระหว่าง physics step สองตัว ภาพบนจอเลยไม่ดูกระตุก แม้ว่า physics เองจะอัปเดตเป็นก้อนคงที่เท่านั้นก็ตาม

11. กับดักของ floating-point ในโค้ด physics

ตัวเลขทุกตัวในทุกตัวอย่างข้างบนเป็น float หรือ double — การประมาณค่าจำนวนจริงของคอมพิวเตอร์ ไม่ใช่ค่าที่แน่นอน โค้ด physics รันการประมาณค่าพวกนี้หลายล้านครั้งตลอดเซสชันการเล่น ดังนั้นสองนิสัยนี้สำคัญกว่าโค้ดทั่วไปมาก

อย่าเทียบ float ด้วย ==

คณิตศาสตร์ค่าจริงแบบเดียวกัน ที่ทำผ่านลำดับ floating-point operation ต่างกัน มักจะลงเอยที่ค่าที่คลาดเคลื่อนในหลักไบนารีท้าย ๆ นั่นแปลว่าการเช็คความเท่ากันแบบ t == 0.3 อาจล้มเหลวแบบเงียบ ๆ ได้ แม้ว่า t จะเท่ากับ 0.3 เป๊ะในทางปฏิบัติก็ตาม:

#include <cstdio>

int main() {
    double t = 0.0;
    const double dt = 0.1;
    for (int i = 0; i < 3; i++) t += dt;

    printf("t after 3 steps of 0.1  = %.17f\n", t);
    printf("t == 0.3 ?               %s\n", (t == 0.3) ? "true" : "false");

    const double eps = 1e-9;
    printf("fabs(t - 0.3) < eps ?    %s\n", (t > 0.3 - eps && t < 0.3 + eps) ? "true" : "false");
}

Output:

t after 3 steps of 0.1  = 0.30000000000000004
t == 0.3 ?               false
fabs(t - 0.3) < eps ?    true

การบวก 0.1 สามครั้งได้ 0.30000000000000004 ไม่ใช่ 0.3 — 0.1 เองเก็บให้แม่นยำเป๊ะในระบบ binary floating point ไม่ได้ ความคลาดเคลื่อนเล็ก ๆ เลยฝังอยู่ตั้งแต่การบวกครั้งแรก การเช็ค == ตรง ๆ รายงานว่า false แม้ว่าค่าทั้งสองจะเท่ากันในทางปฏิบัติทุกแบบก็ตาม ทางแก้คือ epsilon comparison: แทนที่จะถามว่า "ค่าสองตัวนี้เท่ากันเป๊ะไหม" ให้ถามว่า "ค่าสองตัวนี้ห่างกันไม่เกินค่าคลาดเคลื่อนเล็ก ๆ (eps ย่อจาก "epsilon") ไหม" ใช้แบบนี้กับการเทียบค่า physics ใด ๆ — เช็คว่าลูกบอล "หยุดแล้ว" หรือยัง (speed < eps) แทนที่จะเป็น speed == 0, หรือสองตำแหน่ง "ตรงกัน" แทนที่จะต้องเหมือนกันทุกบิต

ความแม่นยำแย่ลงยิ่งอยู่ไกลจากศูนย์

ตัวเลข floating-point ไม่ได้กระจายความแม่นยำให้เท่ากันทุกค่า — มันอัดความแม่นยำไว้ใกล้ศูนย์เยอะกว่า และมีน้อยลงเมื่อไกลออกไป นั่นแปลว่าการบวกเลขเล็ก ๆ ตัวเดิมเข้ากับ position มีผลต่างกันมาก ขึ้นอยู่กับว่า position นั้นใหญ่แค่ไหนอยู่แล้ว:

#include <cstdio>

int main() {
    float smallStep = 0.01f;
    float positions[] = {0.0f, 1000.0f, 100000.0f, 10000000.0f};
    for (float start : positions) {
        float p = start;
        for (int i = 0; i < 100; i++) p += smallStep;   // should move by 1.0 total
        printf("start = %10.1f  moved by %.6f (expected 1.000000)\n", start, p - start);
    }
}

Output:

start =        0.0  moved by 0.999999 (expected 1.000000)
start =     1000.0  moved by 1.000977 (expected 1.000000)
start =   100000.0  moved by 0.781250 (expected 1.000000)
start = 10000000.0  moved by 0.000000 (expected 1.000000)

ใกล้จุดกำเนิด การบวก 0.01f ร้อยครั้งขยับวัตถุเกือบเท่ากับ 1.0 เป๊ะ ตามที่คาดไว้ ห่างจากจุดกำเนิดสิบล้านหน่วย มันไม่ขยับเลย — 0.01 เล็กกว่าช่องว่างระหว่างค่า float ที่แทนได้สองค่าตรงนั้น ทุกการบวกเลยถูกปัดทิ้งเงียบ ๆ จนเหลือศูนย์ นี่คือบั๊กเบื้องหลังรายงานแบบ "ตัวละครโลกเปิดของฉันเริ่มสั่น หรือไม่ตอบสนองต่อการขยับเล็ก ๆ เมื่อเดินออกไปไกลจากจุดกึ่งกลางแผนที่พอ" เกมโลกเปิดขนาดใหญ่จัดการเรื่องนี้ด้วยเทคนิค floating-origin: ย้ายจุดศูนย์กลางของโลกทั้งหมด (หรืออย่างน้อยทุกอย่างใกล้กล้อง) กลับไปที่ (0, 0, 0) เป็นระยะ ๆ หรือเก็บตำแหน่งโลกเป็น double แล้วแปลงเป็น float เทียบกับกล้องแค่ตอนก่อน render เท่านั้น

12. คำศัพท์

• Derivative (อนุพันธ์) — อัตราการเปลี่ยนแปลงของปริมาณ ณ ขณะหนึ่ง; velocity คือ derivative ของ position
• Integral (อินทิกรัล) — การสะสมของอัตราการเปลี่ยนแปลงตลอดเวลา; position คือ integral ของ velocity
• Riemann sum — การประมาณ integral โดยบวกสี่เหลี่ยมผืนผ้าบาง ๆ เข้าด้วยกัน (อัตรา x ช่วงเวลาเล็ก ๆ)
• Numerical integration — การประมาณสมการการเคลื่อนที่แบบต่อเนื่องด้วยก้าวเล็ก ๆ แบบไม่ต่อเนื่องที่ทำซ้ำ
• Timestep (dt) — ช่วงเวลาเล็ก ๆ ที่ simulation step หนึ่งก้าวไปข้างหน้า
• Explicit (forward) Euler — ก้าว position โดยใช้ velocity เก่า; ง่าย แต่รั่วพลังงานและระเบิดได้
• Semi-implicit (symplectic) Euler — ก้าว velocity ก่อน แล้วก้าว position ด้วย velocity ใหม่; เสถียร ต้นทุนต่ำ ค่าเริ่มต้นปกติในเกม
• Symplectic — คุณสมบัติของ integrator ที่พลังงานอยู่ในขอบเขต (สั่นไปมา) แทนที่จะดริฟท์หนีไปตามเวลา
• Phase space — กราฟ velocity เทียบกับ position; วิธีที่ดีในการดูว่า integrator เสถียรหรือกำลังระเบิด
• Verlet integration — ก้าวโดยใช้ position ก่อนหน้ากับปัจจุบัน แทนที่จะเก็บ velocity; เหมาะกับผ้าและ ragdoll ที่ใช้ constraint
• Position-based dynamics — ซิมูเลชันโดยขยับ position ตรง ๆ ให้ตรงตาม constraint (เช่นระยะทางคงที่) แทนที่จะผ่าน force กับ velocity
• RK4 (Runge-Kutta อันดับสี่) — integrator ที่สุ่มตัวอย่าง slope สี่ครั้งต่อก้าวแล้วผสมกัน; แม่นยำมาก ต้นทุนประมาณสี่เท่าต่อก้าว
• Fixed timestep — ก้าว physics simulation ด้วย dt เดิมเสมอ ไม่ว่าเฟรมที่ render จะใช้เวลานานแค่ไหนจริง ๆ
• Accumulator — ยอดรวมที่วิ่งอยู่ของ "เวลาจริงที่ค้างซิมูเลชันอยู่" ถูกใช้ทีละ fixed timestep
• Spiral of death (เกลียวมรณะ) — ลูป feedback ที่วิ่งหนีตัวเอง เฟรมช้าทำให้มี step ตามให้ทันเพิ่ม ซึ่งกินเวลาเพิ่ม ทำให้เฟรมถัดไปมี step มากขึ้นอีก
• Interpolation (alpha) — ผสม state physics ก่อนหน้ากับปัจจุบันโดยใช้เศษที่เหลือของ accumulator เพื่อให้ render ดูลื่นระหว่าง fixed physics step
• Epsilon comparison — เช็คว่าค่า floating-point สองตัวห่างกันไม่เกินค่าคลาดเคลื่อนเล็ก ๆ แทนที่จะเช็คความเท่ากันเป๊ะ
• Floating-origin — ย้ายจุดศูนย์กลางของโลกขนาดใหญ่กลับมาใกล้จุดกำเนิดเป็นระยะ ๆ เพื่อเลี่ยงการเสียความแม่นยำของ floating-point เมื่อไกลจาก (0,0,0)

13. แบบฝึกหัด

step 0: x=1.0000  v=0.0000  a=-2.0000  E=1.0000
step 1: x=1.0000  v=-1.0000            E=1.5000
step 2: x=0.5000  v=-2.0000            E=2.2500

step 0: x=1.0000  v=0.0000   E=1.0000
step 1: x=0.5000  v=-1.0000  E=0.7500
step 2: x=-0.2500 v=-1.5000  E=1.1875

นั่นคือชุดเครื่องมือ numerical methods ที่โปรแกรมเมอร์เกมใช้จริง: derivative กับ integral ให้คำศัพท์ (อัตราการเปลี่ยนแปลง, การสะสม), explicit Euler แสดงวิธีตรงไปตรงมาและทำไมมันถึงล้มเหลว, semi-implicit Euler กับ Verlet คือค่าเริ่มต้นที่เสถียรและต้นทุนต่ำที่ใช้อยู่ในเอนจินจริง, RK4 มีไว้สำหรับเคสหายากที่ต้องการความแม่นยำจริง ๆ และ fixed-timestep accumulator ทำให้ทั้งหมดนี้ทำงานเหมือนกันไม่ว่าเฟรมจะ render เร็วหรือช้าแค่ไหน จำสองนิสัยจากหัวข้อ 11 ไว้ — อย่าเทียบ float ด้วย == และระวังความแม่นยำเมื่อไกลจากจุดกำเนิด — แล้วคุณก็มีทุกอย่างที่ต้องใช้ในการก้าวการเคลื่อนที่ไปข้างหน้าอย่างปลอดภัย เฟรมแล้วเฟรมเล่า