Problem

Create a program for solving general IVP with backward Euler (a frame of the program can be found in BackwardEulerNewton.m). You should use Newton iterations for solving the implicit difference equation that arises when computing $u_{n+1}$ from $u_n$.

You should test your solver on three different IVP:s. First

$$\begin{array}{rcl} \dot u (t) &=& -u(t) \quad 0<t<T \\ u(0) &=& 1 \end{array}$$ (3)

Compare your result to the exact solution.

The second example is a stiff IVP, ``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.$$ (4)

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 computed solution to the Akzo-Nobel problem.

The third and last example is the ``Volterra-Lotka'' equations, also known as the prey-predator equations,

$$\begin{array}{rcl} \dot u_1 &=& u_1(a-bu_2) \\ \dot u_2 &=& -u_2(c-du_1) \end{array}$$ (5)

where $u_1$ is the number of some prey (for example rabbits) and $u_2$ is the number of its predator (for example foxes). $a,b,c,d$ are parameters representing the interaction of the two species.

Compute solutions for different parameters. The solution should have a periodic behaviour.

Question 2 Can you notice any damping of the solution? What may be the cause of this?