Back to Subject

Unit 4: Numerical Solution of Ordinary Differential Equations

ECE183 — Mathematics For Robotics 8 min read

I. Orientation: Differential Equations as Robot Motion Models

An ordinary differential equation (ODE) relates an unknown function of one independent variable to its derivatives. In robotics, the independent variable is usually time (t), while the dependent variables describe position, velocity, joint angle, or actuator state. Numerical methods approximate the solution when an exact analytical solution is difficult or unavailable.

  • General first-order form:
    [
    \frac{dy}{dt}=f(t,y), \qquad y(t_0)=y_0
    ]
    Here, (t) is the independent variable, (y(t)) is the unknown state, (f) is a known function, and (y_0) is the initial value.

  • Initial-value assumption: A unique solution generally requires (f(t,y)) and its relevant derivative with respect to (y) to be continuous near ((t_0,y_0)).

  • Discrete grid: Numerical computation uses
    [
    t_n=t_0+nh
    ]
    where (h) is the step size and (y_n) approximates (y(t_n)).

  • Error convention: Local truncation error is introduced in one step; global error accumulates over many steps. Smaller (h) usually improves accuracy but increases computation.

  • Robot-state interpretation: For a state vector (\mathbf{x}),
    [
    \dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x})
    ]
    may contain joint positions, velocities, and sensor or actuator states.

II. Numerical Methods for Initial Value Problem

A. Numerical methods for initial value problem

Numerical methods for an initial value problem generate successive approximations from a known starting state.

  • Problem statement: Given
    [
    y'=f(t,y), \qquad y(t_0)=y_0,
    ]
    calculate (y_1,y_2,\ldots) at (t_1,t_2,\ldots).

  • Integral foundation: Integrating between (tn) and (t{n+1}=tn+h) gives
    [
    y    {n+1}=yn+\int{tn}^{t{n+1}}f(t,y(t))\,dt.
    ]
    Numerical methods approximate this integral.

  • One-step methods: Euler and Runge-Kutta methods calculate (y_{n+1}) primarily from (t_n,y_n).

  • Multi-step methods: These use several previous values such as (yn,y{n-1},y_{n-2}) to predict the next value.

  • Algorithmic pattern:

    TEXT
      Set t = t0 and y = y0
  For n = 0, 1, ..., N-1:
      calculate an approximation to y at t + h
      replace t by t + h


    The step size (h) must be selected according to the required accuracy and stability.

  • Worked example: For (y'=y,\ y(0)=1), with (h=0.1), the exact value at (t=0.1) is (e^{0.1}\approx1.10517). Numerical methods estimate this value from the initial slope (f(0,1)=1).

III. Euler and Modified Euler Method

A. Euler and Modified Euler method

Euler’s method uses the slope at the beginning of an interval, whereas Modified Euler uses a more accurate average of beginning and ending slopes.

  • Euler principle: Approximate the curve by its tangent at ((t_n,yn)):
    [
    y    {n+1}=y_n+h f(t_n,y_n).
    ]
    The method has first-order accuracy, with global error (O(h)).

  • Geometric meaning: The derivative (f(t_n,y_n)) is held constant across the interval ([t_n,t_n+h]).

  • Modified Euler predictor-corrector: First predict:
    [
    y_{n+1}^{(p)}=y_n+h f(t_n,yn).
    ]
    Then correct using the average of two slopes:
    [
    y    {n+1}=y_n+\frac{h}{2}
    \left[f(t_n,yn)+f(t{n+1},y_{n+1}^{(p)})\right].
    ]

  • Accuracy comparison:

    1. Euler: Uses one slope and is second-order locally, first-order globally.
    2. Modified Euler: Uses two slopes and is second-order globally, (O(h^2)).

  • Worked example: For (y'=y,\ y_0=1,\ h=0.1):
    [
    y_1^{(p)}=1+0.1(1)=1.1,
    ]
    [
    y_1=1+\frac{0.1}{2}(1+1.1)=1.105.
    ]
    This is closer to (e^{0.1}\approx1.10517) than Euler’s (1.1).

  • Practical limitation: Both methods may become unstable for rapidly changing or stiff equations. Step-size reduction improves accuracy but may not resolve severe stability restrictions.

IV. Runge-Kutta Methods

A. Runge-Kutta methods

Runge-Kutta methods obtain high accuracy by evaluating several slopes within each step, without explicitly calculating higher derivatives.

  • General purpose: They approximate the integral of (f(t,y)) by weighted slope samples inside the step.

  • Fourth-order Runge-Kutta method: The common RK4 formula is
    [
    \begin{aligned}
    k_1&=f(t_n,y_n),\
    k_2&=f\left(t_n+\frac h2,y_n+\frac h2k_1\right),\
    k_3&=f\left(t_n+\frac h2,y_n+\frac h2k_2\right),\
    k_4&=f(t_n+h,y_n+hk3),
    \end{aligned}
    ]
    [
    y    {n+1}=y_n+\frac h6(k_1+2k_2+2k_3+k_4).
    ]

  • Symbol definitions: (k_1,k_2,k_3,k_4) are slope estimates; (h) is the step size; (t_n,y_n) are the current grid point and approximation.

  • Accuracy: RK4 has local truncation error (O(h^5)) and global error (O(h^4)), assuming sufficient smoothness.

  • Computational trade-off: Four evaluations of (f) are required per step, but no prior derivative values are stored.

  • Worked example: For (y'=y,\ y(0)=1,\ h=0.1), the RK4 estimate is approximately
    [
    y(0.1)\approx1.1051708,
    ]
    which closely matches (e^{0.1}=1.1051702).

  • Robot application: RK4 can integrate a robot’s state equation
    [
    \dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x},\mathbf{u})
    ]
    when the control input (\mathbf{u}) is known over the step.

V. Multi-step Methods

A. Multi-step methods

Multi-step methods use previously computed solution values and derivative values to advance the solution efficiently.

  • Basic principle: Instead of evaluating (f) repeatedly within one interval, approximate the integral
    [
    y_{n+1}=yn+\int{tn}^{t{n+1}}f(t,y(t))\,dt
    ]
    using values at several earlier grid points.

  • Adams-Bashforth example: The explicit four-step formula is
    [
    y_{n+1}=y_n+\frac h{24}
    \left(55fn-59f{n-1}+37f{n-2}-9f{n-3}\right),
    ]
    where (f_j=f(t_j,y_j)).

  • Adams-Moulton correction: An implicit four-step formula is
    [
    y_{n+1}=yn+\frac h{24}
    \left(9f    {n+1}+19fn-5f{n-1}+f{n-2}\right).
    ]
    Since (f    {n+1}) depends on the unknown (y_{n+1}), prediction or iteration is required.

  • Starting values: A four-step method needs (y_0,y_1,y_2,y_3). These are commonly generated by RK4.

  • Advantages: Multi-step methods can achieve high order with only one new function evaluation per step after initialization.

  • Limitations: They require stored history and are less convenient when the step size changes. A restart is often needed after discontinuities or failed convergence.

VI. Milne Method

A. Milne method

Milne’s method is a predictor-corrector multi-step method that uses four equally spaced previous points to estimate the next solution.

  • Milne predictor: For (y'=f(t,y)),
    [
    y{n+1}^{(p)}
    =y    {n-3}+\frac{4h}{3}
    \left(2fn-f{n-1}+2f_{n-2}\right).
    ]

  • Milne corrector:
    [
    y{n+1}^{(c)}
    =y    {n-1}+\frac h3
    \left(f_{n+1}^{(p)}+4fn+f{n-1}\right).
    ]
    Here,
    [
    f{n+1}^{(p)}=f(t{n+1},y_{n+1}^{(p)}).
    ]

  • Iteration: The corrected value can be substituted again into (f_{n+1}) if greater accuracy is required.

  • Required data: Values (y{n-3},y{n-2},y_{n-1},y_n), usually at equal spacing (h), are necessary. Initial values may be generated using RK4.

  • Procedure:

    TEXT
      Calculate f_n, f_(n-1), f_(n-2)
  Compute the Milne predictor y_(n+1)^p
  Evaluate f_(n+1)^p
  Compute the Milne corrector y_(n+1)^c
  Repeat correction if required
  Shift the four-point history

  • Accuracy and limitation: The predictor-corrector pair is high-order for smooth problems, but accumulated errors in stored history can affect later values. It is also sensitive to irregular step sizes.

VII. Lagrange's Equation

A. Lagrange's equation

Lagrange’s equation provides the differential equations of a robot by expressing dynamics in terms of energy and generalized coordinates.

  • General equation:
    [
    \frac{d}{dt}\left(\frac{\partial L}{\partial \dot q_i}\right)
    -\frac{\partial L}{\partial q_i}
    =Q_i,
    ]
    where (q_i) is a generalized coordinate, (\dot q_i) its velocity, (L=T-V) is the Lagrangian, (T) is kinetic energy, (V) is potential energy, and (Q_i) is a generalized non-conservative force.

  • Robot dynamics form: For an (n)-joint manipulator,
    [
    \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}
    +\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}
    +\mathbf{g}(\mathbf{q})
    =\boldsymbol{\tau},
    ]
    where (\mathbf{q}) is the joint-position vector, (\mathbf{M}) the inertia matrix, (\mathbf{C}\dot{\mathbf{q}}) Coriolis and centrifugal terms, (\mathbf{g}) gravity, and (\boldsymbol{\tau}) actuator torque.

  • Conversion to an ODE: Define the state
    [
    \mathbf{x}=
    \begin{bmatrix}
    \mathbf{q}\
    \dot{\mathbf{q}}
    \end{bmatrix}.
    ]
    Then
    [
    \dot{\mathbf{x}}=
    \begin{bmatrix}
    \dot{\mathbf{q}}\
    \mathbf{M}^{-1}
    \left(\boldsymbol{\tau}
    -\mathbf{C}\dot{\mathbf{q}}
    -\mathbf{g}\right)
    \end{bmatrix}.
    ]
    This first-order system can be integrated by Euler, Modified Euler, RK4, or a multi-step method.

  • Simple pendulum example: For angle (q), mass (m), length (l), and gravity (g),
    [
    ml^2\ddot q+mgl\sin q=\tau.
    ]
    With (x_1=q,\ x_2=\dot q),
    [
    \dot x_1=x_2,\qquad
    \dot x_2=\frac{\tau}{ml^2}-\frac{g}{l}\sin x_1.
    ]

  • Numerical significance: Lagrange’s equation supplies the physical model; numerical ODE methods supply its time evolution. Accuracy depends on the model parameters, force inputs, integration step size, and numerical stability.

Previous
Unit 3
Next
Unit 5