IVP with Systems of First Order ODEs

It was mentioned in the introduction (see Eq. 4) that systems of first order ODEs can arise from a single ODE of order larger than unity. Systems of ODEs also arises naturally from physical modeling. For instance, the following IVP describes the concentrations y_i, $i=1,2,\cdots , n$ of n chemicals in a reactor as a function of time.

$$\frac{dy_1}{dt} = f_1(y_1,y_2,\cdots , y_n, t)$$

$$\frac{dy_2}{dt} = f_2(y_2,y_2,\cdots , y_n, t)$$ ....................................

$$\frac{dy_n}{dt} = f_n(y_2,y_2,\cdots , y_n, t)$$

$$y_1(0) = y_{10}, \:\: y_2(0) =y_{20},\cdots , y_n(0) = y_{n0}.$$ (26)

The methods we discussed for the solution of a single first order ODE can be extended to the system of ODEs. Let's define the n-dimensional vector of solutions ${\bf {y}}(t)={(y_1,y_2,\cdots ,y_n)}^T$ and the n dimensional vector of the right hand side functions ${\bf {f}}({\bf {y}},t)={(f_1, f_2, \cdots , f_n)}^T$. Similarly, the initial conditions can be arranged in a vector ${\bf {y_0}} = {(y_{10},y_{20},\cdots , y_{n0})}^T$. This allows us to write the IVP in the vector notation as

$$\frac{d{\bf {y}}}{dt} = {\bf {f}}({\bf {y}},t), \:\:\: {\bf {y}}(0) = {\bf {y_0}}.$$ (27)

Now, if we apply the forward Euler method for the solution of Eq. 28, with constant time step size h, we get the following explicit equation for the solution at the (k+1)th step given the solution at the k the step:

$${\bf {y}}_{k+1} = {\bf {y}}_{k} + h {\bf {f}}({\bf {y}}_{k}, t_k).$$ (28)

However, the implicit backward Euler method gives

$${\bf {y}}_{k+1} = {\bf {y}}_{k} + h {\bf {f}}({\bf {y}}_{k+1}, t_{k+1}),$$ (29)

which necessitates in general the solution of a system of non-linear algebraic equations at every time step. We can employ Newton's method for a system of non-linear equations to achieve this.