オイラー法——微分方程式をコードで解く

これまでのepisodeで「x += vx * dt」という更新式を使ってきました。実はこれはオイラー法（Euler’s method）と呼ばれる数値積分の手法です。なぜこの式が成り立つのか、どんな限界があるのかを今回は丁寧に掘り下げます。

数値解法の理解は、物理シミュレーションの精度や安定性に直結します。ゲームエンジンがなぜ特定の積分手法を採用するのかも分かるようになります。

微分方程式とは

物理の世界は「変化の方程式」——微分方程式で記述されます。

dx/dt = vx    （位置の変化率は速度）
dvx/dt = ax   （速度の変化率は加速度）

dx/dt は「時刻tにおけるxの瞬間変化率」です。これを「xはvxの速さで変化している」と読めます。

こういった方程式を解くとは、「時刻tにおけるx, vxの値」を求めることです。解析的に解ける場合もありますが、バネ・衝突・摩擦が絡み合う複雑なシミュレーションでは、コンピューターで数値的に近似解を求めることになります。

オイラー法の直感的説明

オイラー法の考え方は非常にシンプルです。

「今の変化率が少しの間続くと仮定して、次のステップを計算する」

数学的には：

x(t + dt) ≈ x(t) + f(x, t) * dt

ここで f(x, t) は変化率の関数（速度や加速度）です。

グラフで考えると、真の曲線の接線方向に dt だけ進むイメージです。dt が小さいほど接線と曲線の差（誤差）が小さくなります。

// オイラー法の実装（シンプル版）
function eulerStep(x, v, a, dt) {
  var x_new = x + v * dt;       // 今の速度で位置を更新
  var v_new = v + a * dt;       // 今の加速度で速度を更新
  return { x: x_new, v: v_new };
}

dtを小さくすると精度が上がる理由

オイラー法の誤差は O(dt²) のオーダーです。dtを半分にすると、1ステップの誤差は1/4になります。ただし同じ時間をシミュレートするにはステップ数が2倍必要なので、計算量は2倍になります。

この「精度 vs 計算量」のトレードオフがシミュレーション設計の核心です。

// dt = 0.1 のとき（荒い）
var steps = 10;  // 1秒を10ステップで

// dt = 0.001 のとき（精細）
var steps = 1000;  // 1秒を1000ステップで（100倍の計算）

60fpsのアニメーションではdtは約0.016秒です。通常の物理シミュレーションには十分な精度ですが、高速な衝突（例: 弾丸）などには不足することがあります。

オイラー法の誤差の可視化

円運動はオイラー法の誤差を確認するのに適した例題です。真円から外れていく度合いで誤差が分かります。

v = (-y, x) という速度場を使うと中心回りの円運動になります（解析解は完全な円）。

オイラー法の誤差——大きなdtと小さなdtの比較（円運動）

var cx = W/2, cy = H/2, R = 70, omega = 1.0;
var trajs = [
{x: cx+R, y: cy, vx:0, vy:-omega*R, dt:0.08,  color:'#ff6b6b', label:'dt=0.08（粗）', trail:[]},
{x: cx+R, y: cy, vx:0, vy:-omega*R, dt:0.016, color:'#00e5ff', label:'dt=0.016（標準）', trail:[]},
];
var accum = 0;
function loop(ts){
if(!loop.last) loop.last=ts;
var real=Math.min((ts-loop.last)/1000,0.05); loop.last=ts;
accum+=real;
// 各軌道を固有dtで独立に進める
trajs.forEach(function(t){
  while(accum>=t.dt){
    // 速度場: v = omega * (-y+cy, x-cx) で円運動
    t.vx = -omega*(t.y-cy);
    t.vy =  omega*(t.x-cx);
    t.x += t.vx*t.dt;
    t.y += t.vy*t.dt;
    t.trail.push({x:t.x,y:t.y});
    if(t.trail.length>300) t.trail.shift();
  }
});
// accumは最小dtで進める（両方の基準を満たすため）
if(accum>=0.016) accum-=0.016;
ctx.fillStyle='#0d1117'; ctx.fillRect(0,0,W,H);
ctx.beginPath(); ctx.arc(cx,cy,R,0,Math.PI*2);
ctx.strokeStyle='rgba(255,255,255,0.12)'; ctx.lineWidth=1; ctx.stroke();
ctx.beginPath(); ctx.arc(cx,cy,3,0,Math.PI*2); ctx.fillStyle='#fff'; ctx.fill();
trajs.forEach(function(t,idx){
  if(t.trail.length>1){
    ctx.beginPath(); ctx.moveTo(t.trail[0].x,t.trail[0].y);
    for(var i=1;i<t.trail.length;i++) ctx.lineTo(t.trail[i].x,t.trail[i].y);
    ctx.strokeStyle=t.color; ctx.lineWidth=1.5; ctx.stroke();
  }
  ctx.beginPath(); ctx.arc(t.x,t.y,6,0,Math.PI*2); ctx.fillStyle=t.color; ctx.fill();
  ctx.fillStyle=t.color; ctx.font='12px monospace'; ctx.fillText(t.label,8,20+idx*18);
});
requestAnimationFrame(loop);
}
requestAnimationFrame(loop);

赤（dt=0.08）は徐々に円の外側に螺旋を描きます——これが誤差の蓄積です。青（dt=0.016）はほぼ円を保ちます。

改良オイラー法（中点法）の実装

改良オイラー法（別名: 中点法、2次Runge-Kutta法）は、ステップの中間点での傾きを使うことで精度を大幅に向上させます。

k1 = f(x, t)                    （始点の傾き）
k2 = f(x + k1 * dt/2, t + dt/2) （中点の傾き）
x_new = x + k2 * dt              （中点の傾きで更新）

コードで書くと：

function midpointStep(x, v, accelFn, dt) {
  // 始点の加速度
  var a1 = accelFn(x, v);

  // 中点を予測
  var x_mid = x + v * dt / 2;
  var v_mid = v + a1 * dt / 2;

  // 中点の加速度
  var a2 = accelFn(x_mid, v_mid);

  // 中点の速度で位置を更新、中点の加速度で速度を更新
  return {
    x: x + v_mid * dt,
    v: v + a2 * dt
  };
}

オイラー法 vs 中点法の比較（円運動・大きなdt）

var DT = 0.05;
var omega = 1.2;
var R = 70;
var cx = W/2, cy = H/2;

var euler = {x: cx + R, y: cy, vx: 0, vy: -omega * R, trail: []};
var midpt = {x: cx + R, y: cy, vx: 0, vy: -omega * R, trail: []};

var accum = 0;

function accelX(x, y) { return -omega * omega * (x - cx); }
function accelY(x, y) { return -omega * omega * (y - cy); }

function loop(ts) {
if (!loop.last) loop.last = ts;
var realDt = Math.min((ts - loop.last) / 1000, 0.05);
loop.last = ts;
accum += realDt;

while (accum >= DT) {
  accum -= DT;

  var ax1 = accelX(euler.x, euler.y);
  var ay1 = accelY(euler.x, euler.y);
  euler.x += euler.vx * DT;
  euler.y += euler.vy * DT;
  euler.vx += ax1 * DT;
  euler.vy += ay1 * DT;
  euler.trail.push({x: euler.x, y: euler.y});
  if (euler.trail.length > 150) euler.trail.shift();

  var ax1 = accelX(midpt.x, midpt.y);
  var ay1 = accelY(midpt.x, midpt.y);
  var mpx = midpt.x + midpt.vx * DT / 2;
  var mpy = midpt.y + midpt.vy * DT / 2;
  var mvx = midpt.vx + ax1 * DT / 2;
  var mvy = midpt.vy + ay1 * DT / 2;
  var mx2 = accelX(mpx, mpy);
  var my2 = accelY(mpx, mpy);
  midpt.x += mvx * DT;
  midpt.y += mvy * DT;
  midpt.vx += mx2 * DT;
  midpt.vy += my2 * DT;
  midpt.trail.push({x: midpt.x, y: midpt.y});
  if (midpt.trail.length > 150) midpt.trail.shift();
}

ctx.fillStyle = '#0d1117';
ctx.fillRect(0, 0, W, H);

ctx.beginPath();
ctx.arc(cx, cy, R, 0, Math.PI * 2);
ctx.strokeStyle = 'rgba(255,255,255,0.2)';
ctx.lineWidth = 1;
ctx.stroke();

[[euler,'#ff6b6b','オイラー法'], [midpt,'#4ecdc4','中点法']].forEach(function(item, idx) {
  var obj = item[0], col = item[1], name = item[2];
  for (var i = 1; i < obj.trail.length; i++) {
    ctx.beginPath();
    ctx.moveTo(obj.trail[i-1].x, obj.trail[i-1].y);
    ctx.lineTo(obj.trail[i].x, obj.trail[i].y);
    ctx.strokeStyle = col;
    ctx.lineWidth = 2;
    ctx.stroke();
  }
  ctx.beginPath();
  ctx.arc(obj.x, obj.y, 7, 0, Math.PI * 2);
  ctx.fillStyle = col;
  ctx.fill();
  ctx.fillStyle = col;
  ctx.font = '12px monospace';
  ctx.fillText(name, 8, 18 + idx * 18);
});

requestAnimationFrame(loop);
}
requestAnimationFrame(loop);

同じdtで比較すると、中点法（水色）は真円に近い軌道を保ちます。オイラー法（赤）はエネルギーが増大して外に広がります。

シンプレクティックオイラー法（エネルギー保存）

実は、速度と位置の更新順序を変えるだけで「エネルギーが保存される」ようになります。これをシンプレクティックオイラー法（半陰的オイラー法）と呼びます。

// 通常のオイラー法（エネルギーが徐々に増大または減少する）
x_new = x + v * dt
v_new = v + a * dt

// シンプレクティックオイラー法（エネルギーが保存される）
v_new = v + a * dt      // 先に速度を更新
x_new = x + v_new * dt  // 新しい速度で位置を更新（順番が逆！）

ゲームエンジン（Unity等）の Verlet 積分や、多くの物理シミュレーションでこの手法が使われています。計算コストはオイラー法と同じで、精度と安定性が大幅に向上します。

// シンプレクティックオイラー法の実装
function symplecticEulerStep(x, v, a, dt) {
  var v_new = v + a * dt;       // 先に速度を更新
  var x_new = x + v_new * dt;   // 新しい速度で位置を更新
  return { x: x_new, v: v_new };
}

このシリーズでは以降のデモで基本的にシンプレクティックオイラー法を使います。

積分手法の比較表

手法	精度	エネルギー保存	計算コスト
オイラー法	O(dt)	NG（増大）	低
シンプレクティックオイラー	O(dt)	良好	低
中点法（2次RK）	O(dt²)	良好	中
4次Runge-Kutta	O(dt⁴)	良好	高
Verlet積分	O(dt²)	非常に良好	低

ゲームや物理エンジンでは、計算コストの低いシンプレクティックオイラーやVerlet積分がよく使われます。

まとめ

この回でやったこと：

微分方程式の概念と「数値解法で近似する」意味を理解した
オイラー法 x += f(x,t) * dt の仕組みと誤差の原因を学んだ
円運動で大きなdtと小さなdtの誤差を可視化した
中点法（改良オイラー法）でより精度の高い積分ができることを確認した
更新順序を入れ替えるだけのシンプレクティックオイラー法を紹介した

次回は「バネの力」を実装します。フックの法則を使ったバネ振動・減衰・バネチェーンによるロープシミュレーションを作ります。