---

Скачать книгу

---

initial value problem and is chosen in such a way that it produces an exact solution for a certain model problem. To this end, let us examine the scalar ODE:

      (2.41)

      and let us approximate it using the Theta method:

(2.42)

      where the parameter has not yet been specified. We determine it using the heuristic that this so-called Theta method should be exact for the linear constant-coefficient model problem:

(2.43)

Based on this heuristic and by using the exact solution from (2.43) in scheme (2.42) , we get the value (you should check that this formula is correct; it is a bit of algebra). We get:

      (2.44)

Note: this is a different kind of exponential fitting.

      We need to determine if this scheme is stable (in some sense). To answer this question, we introduce some concepts.

      Definition 2.3 The region of (absolute) stability of a numerical method for an initial value problem is the set of complex values for which all discrete solutions of the model problem (2.43) remain bounded when n approaches infinity.

      Definition 2.4 A numerical method is said to be A-stable if its region of stability is the left-half plane, that is:

      Returning to the exponentially fitted method, we can check that it is A-stable because for all we have , and this condition can be checked using the scheme (2.42) for the model problem (2.43).

      We can generalise the exponential fitting technique to linear and non-linear systems of equations. In the case of a linear system, stiffness is caused by an isolated real negative eigenvalue of the matrix A in the equation:

(2.45)

      where and A is a constant matrix with eigenvalues and eigenvectors

The solution of Equation (2.45) is given by:

      where are arbitrary constants and is a particular integral. In this case we can employ exponential fitting by fitting the dominant eigenvalues which can be computed by the Power method, for example.

      If we assume that the real parts of the eigenvalues are less than zero, we can conclude that the solution tends to the steady-state. Even though this solution is well-behaved the cause of numerical instabilities is the presence of quickly decaying transient components of the solution caused by the dominant eigenvalues of the matrix A in (2.45).

Let us take an example whose matrix has already been given in diagonal form:

      The solution of this system is given by:

Скачать книгу