Chapter 8 Exercise

8.5. Chapter 8 Exercise#

You should try the following exercise questions first, then check with the answers.

Exercise 8.1

A certain chemical reaction is described by a stiff system

\[\begin{split} \begin{aligned} y'_1 ~&=~ -1000y_1 + y_2, \quad y_1(0) = 1 \\ y'_2 ~&=~ 1000y_1 - 2y_2, \quad y_2(0) = 0. \end{aligned} \end{split}\]

Show that \(\,y_1\,\) decays rapidly while \(\,y_2\,\) decays slowly.
Find the eigenvalues of the associated Jacobian matrix and the stiffness ratio.

Solution to

Matlab live script solution for Exercise 8.1

Exercise 8.2

Consider the problem

\[ \begin{aligned} && y' ~&=~ -100y + \cos x, & y(0) &= 1. && \end{aligned} \]

The general solution is \(y(x) = -A\cos x - B\sin x + Ce^{-100x}\). In what sense could this single equation be said to exhibit stiffness? If Euler’s rule were used to obtain a numerical solution, what would be the maximum allowable steplength.

Solution to

From both the equation and the analytical solution it can be seen that the solution has a fast decaying element as well as a slowly varying element, hence the problem can be considered as stiff.

The stability region for Euler method is \([-2, 0]\), this means

\[h\lambda \in [-2, 0]\]

as \(\lambda =-100\), so

\[ h \in (0, 0.02)\]

Exercise 8.3

Detertermine whether the following problem is stiff

\[ \begin{aligned} && y'' ~&=~ -20y' -19y, & y(0) &= 2,\quad y'(0) = -20. && \end{aligned} \]

[Hint: reduce the equation to a system of two first order ODEs i.e. let \(y = y_1\), \(y' = y_2\), then differentiation gives \(y' = y_1'\), and \(y'' = y'_2\,\)].

Solution to

Method 1: Find the analytical solution

Rearrange the equation as

\[ y'' + 20 y' + 19 y=0 \]

The characteristic equation is

\[ \lambda^2 + 20 \lambda + 19=0 \]

so the roots of the characteristic equation are

\[ \lambda_1 = -1, \quad \lambda_2 = -19. \]

so the general solution is

\[ y= c_1 e^{-x} + c_2 e^{-19x} \]

Substituting the initial conditions, we can determine the coefficient \(c_1=1\) and \(c_2=1\). Therefore the solution is

\[ y=e^{-x}+e^{-19x} \]

The solution has two components, \(e^{-19x}\) changing much faster than \(e^{-x}\), so the problem can be considered stiff.

Method 2: Reduce the equation to a first-order system

\[\begin{split} \begin{aligned} y'_1 = & ~ y_2 \\ y'_2 = & -19 y_1 - 20 y_2 \end{aligned} \end{split}\]

The Jacobian matrix is

\[\begin{split} J = \begin{bmatrix} 0 & 1 \\ -19 & -20 \end{bmatrix} \end{split}\]

its eigenvalues are

\[ \lambda_1 = -1, \quad \lambda_2 = -19 \]

so the stiffness ratio of the system is

\[ S = 19 \]

Exercise 8.4

Consider the non-linear chemical system in Example 8.3. Find time-dependent solution of the species \(y_1\), \(y_2\), and \(y_3\) for a duration of 60 seconds using the stiff equation solver routine ode15s from the MATLAB ODE solver routines. Try to run the model using ode113 solver and comment on your results. Change both tolerances RelTol and AbsTol and comment on your results.

Solution to

Matlab code using ode15s solver

Matlab code using ode113 solver

Exercise 8.5

Consider the Robertson Problem in the MATLAB Program 4. Find time-dependent solution of the species \(y_1\), \(y_2\), and \(y_3\) for a duration of \(t = 100\) seconds using the stiff equation solver routine ode15s from MATLAB ODE solver routines. Use tic toc commands to record the run time. Try to run the model using ode113 solver and comment on your results. Change the tolerances RelTol and AbsTol and comment on your results. Using the results from the ode15s, find the eigenvalues of the the Jacobian Matrix at \(t = 0\), \(t = 1.0e-5\) and \(t = 60\) seconds, and by finding the corresponding stiffness ratio, comment on the stability of the system. Produce a table to record the run times and your computational stats.

Solution to

You will find the running time for ode15s is much shorter than that for ode113. Check your results and comment.

Exercise 8.6

Find the interval of absolute stability for the backward differentiation formulae (BDF) for \(k = 1,2,\dots,4\) methods using root locus method, and show whether they are convergent. For the same BDF methods plot the region of absolute stability using boundary root locus method.

Solution to

Matlab code for Exercise 8.6

Solution are shown in the lecture notes Fig. 8.1.

Exercise 8.7

Consider the following system of ODEs which represents the mathematical formulation of an electric circuit:

\[\begin{split} \begin{aligned} \frac{d}{dt} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} ~=~ \begin{pmatrix} -1 & -1\phantom{000}\\ \phantom{-}1 & -5000 \end{pmatrix} ~ \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}, \qquad\text{where}~ \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \end{pmatrix} \end{aligned} \end{split}\]

The above system can be re-arranged as a system of two ODEs:

\[\begin{split} \begin{aligned} && y'_1 ~&=~ -y_1 - y_2 & y_1(0) &= 1 && \\ && y'_2 ~&=~ \phantom{-}y_1 - 5000\,\,y_2 & y_2(0) &= 1 && \end{aligned} \end{split}\]

Find time-dependent solution of \(y_1\) and \(y_2\) for a duration of 5 seconds using the stiff equation solver routine ode15s from the MATLAB ODE solver routines. Try to run the model using ode113 solver and comment on your results. Change both tolerances RelTol and AbsTol and comment on your results. Set the option ‘Stats’ ‘on’ in the options=odeset command, e.g.

options=odeset('Stats','on','RelTol',1e-6,'AbsTol',1.e-8)

and prepare a table for comparison of the stats using the two different solvers. Comment on your results.

Solution to

Matlab code for Exercise 8.7

The plot of the solutions is shown below: