trapazoidal rule

ECE1254H Modeling of Multiphysics Systems. Lecture 15: Nonlinear differential equations. Taught by Prof. Piero Triverio

November 12, 2014 ece1254 , , , , , , , ,

[Click here for a PDF of this post with nicer formatting]


Peeter’s lecture notes from class. These may be incoherent and rough.

Nonlinear differential equations

Assume that the relationships between the zeroth and first order derivatives has the form

F\lr{ x(t), \dot{x}(t) } = 0
x(0) = x_0

The backward Euler method where the derivative approximation is

\dot{x}(t_n) \approx \frac{x_n – x_{n-1}}{\Delta t},

can be used to solve this numerically, reducing the problem to

F\lr{ x_n, \frac{x_n – x_{n-1}}{\Delta t} } = 0.

This can be solved with Newton’s method. How do we find the initial guess for Newton’s? Consider a possible system in fig. 1.


fig. 1. Possible solution points


One strategy for starting each iteration of Newton’s method is to base the initial guess for \( x_1 \) on the value \( x_0 \), and do so iteratively for each subsequent point. One can imagine that this may work up to some sample point \( x_n \), but then break down (i.e. Newton’s diverges when the previous value \( x_{n-1} \) is used to attempt to solve for \( x_n \)). At that point other possible strategies may work. One such strategy is to use an approximation of the derivative from the previous steps to attempt to get a better estimate of the next value. Another possibility is to reduce the time step, so the difference between successive points is reduced.

Analysis, accuracy and stability (\(\Delta t \rightarrow 0\))

Consider a differential equation

\dot{x}(t) = f(x(t), t)
x(t_0) = x_0

A few methods of solution have been considered

  • (FE) \( x_{n+1} – x_n = \Delta t f(x_n, t_n) \)
  • (BE) \( x_{n+1} – x_n = \Delta t f(x_{n+1}, t_{n+1}) \)
  • (TR) \( x_{n+1} – x_n = \frac{\Delta t}{2} f(x_{n+1}, t_{n+1}) + \frac{\Delta t}{2} f(x_{n}, t_{n}) \)

A common pattern can be observed, the generalization of which are called
\textit{linear multistep methods}
(LMS), which have the form

\sum_{j=-1}^{k-1} \alpha_j x_{n-j} = \Delta t \sum_{j=-1}^{k-1} \beta_j f( x_{n-j}, t_{n-j} )

The FE (explicit), BE (implicit), and TR methods are now special cases with

  • (FE) \( \alpha_{-1} = 1, \alpha_0 = -1, \beta_{-1} = 0, \beta_0 = 1 \)
  • (BE) \( \alpha_{-1} = 1, \alpha_0 = -1, \beta_{-1} = 1, \beta_0 = 0 \)
  • (TR) \( \alpha_{-1} = 1, \alpha_0 = -1, \beta_{-1} = 1/2, \beta_0 = 1/2 \)

Here \( k \) is the number of timesteps used. The method is explicit if \( \beta_{-1} = 0 \).

Definition: Convergence


  • \(x(t)\) : exact solution
  • \(x_n\) : computed solution
  • \(e_n\) : where \( e_n = x_n – x(t_n) \), is the global error

The LMS method is convergent if

\max_{n, \Delta t \rightarrow 0} \Abs{ x_n – t(t_n) } \rightarrow 0 %\xrightarrow[t \rightarrow 0 ]{} 0

Convergence: zero-stability and consistency (small local errors made at each iteration),

where zero-stability is “small sensitivity to changes in initial condition”.

Definition: Consistency

A local error \( R_{n+1} \) can be defined as

R_{n+1} = \sum_{j = -1}^{k-1} \alpha_j x(t_{n-j}) – \Delta t \sum_{j=-1}^{k-1} \beta_j f(x(t_{n-j}), t_{n-j}).

The method is consistent if

\lim_{\Delta t} \lr{
\max_n \Abs{ \inv{\Delta t} R_{n+1} } = 0

or \( R_{n+1} \sim O({\Delta t}^2) \)

ECE1254H Modeling of Multiphysics Systems. Lecture 14: Backward Euler method and trapezoidal methods. Taught by Prof. Piero Triverio

November 10, 2014 ece1254 , , ,

[Click here for a PDF of this post with nicer formatting]


Peeter’s lecture notes from class. These may be incoherent and rough.

Backward Euler method

Discretized time dependent partial differential equations were seen to have the form

G \Bx(t) + C \dot{\Bx}(t) = B \Bu(t),

where \( G, C, B \) are matrices, and \( \Bu(t) \) is a vector of sources.

The backward Euler method augments \ref{eqn:multiphysicsL14:20} with an initial condition. For a one dimensional system such an initial condition could a zero time specification

G x(t) + C \dot{x}(t) = B u(t),
x(0) = x_0

Discretizing time as in fig. 1.


fig. 1. Discretized time

The discrete derivative, using a backward difference, is

\dot{x}(t = t_n) \approx \frac{ x_n – x_{n-1} }{\Delta t}

Evaluating \ref{eqn:multiphysicsL14:40} at \( t = t_n \) is

G x_n + C \dot{x}(t = t_n) = B u(t_n),

or approximately

G x_n + C \frac{x_n – x_{n-1}}{\Delta t} = B u(t_n).


\lr{ G + \frac{C}{\Delta t} } x_n = \frac{C}{\Delta t} x_{n-1}
B u(t_n).

Assuming that matrices \( G, C \) are constant, and \( \Delta t \) is fixed, a matrix inversion can be avoided, and a single LU decomposition can be used. For \( N \) sampling points (not counting \( t_0 = 0 \)), \( N \) sets of backward and forward substitutions will be required to compute \( x_1 \) from \( x_0 \), and so forth.

Backwards Euler is an implicit method.

Recall that the forward Euler method gave

x_{n+1} =
x_n \lr{ I – C^{-1} \Delta t G }
+ C^{-1} \Delta t B u(t_n)

This required

  • \( C \) must be invertible.
  • \( C \) must be cheap to invert, perhaps \( C = I \), so that
    x_{n+1} =
    \lr{ I – \Delta t G } x_n
    + \Delta t B u(t_n)
  • This is an explicit method
  • This can be cheap but unstable.

Trapezoidal rule (TR)

The derivative can be approximated using an average of the pair of derivatives as illustrated in fig. 2.


fig. 2. Trapezoidal derivative approximation

\frac{x_n – x_{n-1}}{\Delta t} \approx \frac{

Application to \ref{eqn:multiphysicsL14:40} for \( t_{n-1}, t_n \) respectively gives

G x_{n-1} + C \dot{x}(t_{n-1}) &= B u(t_{n-1}) \\
G x_{n} + C \dot{x}(t_{n}) &= B u(t_{n}) \\

Averaging these

G \frac{ x_{n-1} + x_n }{2} + C
= B
u(t_{n}) }{2},

and inserting the trapezoidal approximation

G \frac{ x_{n-1} + x_n }{2}
x_{n} –
}{\Delta t}
= B
u(t_{n}) }{2},

and a final rearrangement yields

\lr{ G + \frac{2}{\Delta t} C } x_n

\lr{ G – \frac{2}{\Delta t} C } x_{n_1}
+ B
u(t_{n}) }.

This is

  • also an implicit method.
  • requires LU of \( G – 2 C /\Delta t \).
  • more accurate than BE, for the same computational cost.

In all of these methods, accumulation of error is something to be very careful of, and in some cases such error accumulation can even be exponential.

This is effectively a way to introduce central differences. On the slides this is seen to be more effective at avoiding either artificial damping and error accumulation that can be seen in backwards and forwards Euler method respectively.