EE107b | Numerical Integration Schemes

Explanations and examples related to the two integration methods below, as well as many other methods for numerically solving PDEs can be found in the excellent and practical book Applied Numerical Modelling for Engineers .

Euler-Cauchy Method

The first method I used is an adaptation of the Euler-Cuachy method for first-order initial value problems. The method is iterative and is derived from the second-order Taylor approximation. I have included a short outline of the adapted Euler-Cauchy prodecure.

Given a second order equation of form y'' = f(t,y,y') and initial values y(0) and y'(0), iterate the following procedure for t = n*h [0 < n < n_max] where Fs = 1/h is the simulation sampling rate:
  1. y''(t) = f(t,y,y')
  2. y(t+h) = y(t) + h*y'(t) + (h^2/2!)*y''(t)
  3. y'(t+h) = y'(t) + hy''(t)

After implementing this scheme, I quickly discovered that the convergence of this routine was much too poor to generate realistic convergent solutions. Luckily there is a much better way of doing this available...

Runge-Kutta-Nystrom

Because RKN is derived from the fourth-order Taylor approximation, it has much better convergence properties than Euler-Cauchy. For the purposes of modeling synthesis this means that higher eigenmodes of the system will stably converge. RKN requires the computation of several 'auxiliary' quantities, used to estimate higher-order terms of the series expansion.

Given a second-order problem of form y'' = f(t,y,y'') with initial values y(0) and y'(0), iterate the following procedure for t = n*h [0 < n < n_max] where Fs = 1/h is the simulation sampling rate:
  1. k1 = (h/2)*f(t,y,y')
  2. K = (h/2)*(y'(t) + k1/2)
  3. k2 = (h/2)*f(t+h/2, y(t)+K, y'(t)+k1)
  4. k3 = (h/2)*f(t+h/2, y(t)+K, y'(t)+k2)
  5. L = h*(y'(t) + h3)
  6. k4 = (h/2)*f(t+h/2, y(t)+L, y'(t)+2*k3)
  7. y(t+h) = y(t) + h*(y'(t) + (1/3)*(k1 + k2 + k3))
  8. y'(t+h) = y'(t) + (1/3)*(k1 + 2*k2 + 2*k3 + k4)
BACK to "Physical Modeling Synthesis"