Problem1

In the previous session, C3, we chose $S_c(t) = 1$ and implemented an adaptive solver. In this session, you will take one step further and use the stability factor of the problem.

Create a program that computes the stability factor $S_c$ of a general IVP. This can be acomplished through three steps:

1. Solve the IVP (1).
2. Solve the dual problem (3) for different $t=T$.
3. Compute $S_c(t)$ for each $t=T$.

Instead of solving the dual problem backwards in time, you can solve the problem

$$\begin{array}{rcl} \dot \psi(t)& =& A(T-t)^T \psi (t) \\ \psi(0) &=& \psi^0 \end{array}$$ (8)

Question 1 Check this. Hint: Use the substitution $\psi(t) = \phi(T-t)$.

Compute the stability factor for the ``Akzo-Nobel system of chemical reactions'',

find the concentrations $u(t) = (u_1(t),u_2(t), \dots , u_6(t))$ such that for $0<t<T$,

$$\left\{ \begin{array}{rcl} \dot u_1 &=& -2r_1+r_2-r_3-r_4... ..._5 &=& r_2-r_3+r_5 \\ \dot u_6 &=& -r_5 \end{array} \right.$$ (9)

where $F = 3.3\cdot(0.9/737-u_2)$ and the reaction rates are given by $r_1=18.7\cdot u_1^4\sqrt{u_2}$, $r_2=0.58\cdot u_3u_4$, $r_3=0.58/34.4\cdot u_1u_5$, $r_4=0.09\cdot u_1u_4^2$ and $r_5=0.42\cdot u_6^2\sqrt{u_2}$, with the initial condition $u_0 = (0.437, 0.00123, 0, 0, 0, 0.367)$. Compare your result to figure 1.

Figure 1: The stability factor of the Akzo-Nobel problem. $\phi _0 = (1.0,0,0,0,0)$