# How can I implement the implicit Euler method for a small nonlinear system of ODEs?

I am trying to solve a system of coupled ODEs:

\begin{align} \frac{dn_A}{dt} & = e\left[j(t) - f\, θ_H\sinh\left(\frac{g\,n_A}{T}\right)\right] \\ \frac{dθ_H}{dt} & = a\left[bP\,(1 - θ_H )^2 - c~(θ_H)^2 - f\,θ_H\,\sinh\left(\frac{g\,n_A}{T}\right)\right] \\ \frac{dP}{dt} & =h(i - P) + d\,T\,(P-P_a)\left[bP~(1 - θ_H )^2 - c~(θ_H)^2\right]. \end{align}

Here $e, f, g, T, a, b, c,h,i,d$ are constants, $j$ is a time dependent variable, and $P_a$ = 101325. The initial values for the three parameters $n_A, θ_H, P$ are $0, 0.5$, and $8106$ respectively.

I was able to solve the equations using the explicit Euler method, but I am having trouble applying the implicit Euler method to do the same. For the first equation I can use the explicit Euler method to estimate the initial value for $n_A$ at time point $t_{n+1}$, but I am not able to figure out what should I take as a value for the $θ_H$ term that also gets involved in the Newton-Raphson iteration scheme for calculation at $t_{n+1}$.

• Again, please don't add "thank you"-type comments at the end of the question; these statements tend to be "filler". If people answer your question to your satisfaction, thank them by voting up their answers, and leaving comments thanking them. – Geoff Oxberry Mar 28 '14 at 2:31
• If you're taking really large time steps with implicit Euler, then using explicit Euler as a predictor might be significantly worse than just taking the last solution value as your initial guess. – David Ketcheson Mar 28 '14 at 6:39

Define \begin{align} A^{n+1} & = e\left[j^{n+1} - f\, θ_H^{n+1}\sinh\left(\frac{g\,n_A^{n+1}}{T}\right)\right] \\ B^{n+1} & = a\left[bP^{n+1}\,(1 - θ_H^{n+1} )^2 - c~(θ_H^{n+1})^2 - f\,θ_H^{n+1}\,\sinh\left(\frac{g\,n_A^{n+1}}{T}\right)\right] \\ C^{n+1} & = h(i - P^{n+1}) + d\,T\,(P^{n+1}-P_a)\left[bP^{n+1}~(1 - θ_H^{n+1} )^2 - c~(θ_H^{n+1})^2\right] \end{align} Then backward Euler gives you \begin{align} n_A^{n+1} & = n_A^n + \Delta t\,A^{n+1} \\ \theta_H^{n+1} & = \theta_H^n + \Delta t\,B^{n+1} \\ P^{n+1} & = P^n + \Delta t\,C^{n+1} \,. \end{align} Reorganize into the form \begin{align} F^{n+1} &:= n_A^{n+1} - \Delta t\,A^{n+1} - n_A^n = 0 \\ G^{n+1} &:= \theta_H^{n+1} - \Delta t\,B^{n+1} - \theta_H^n = 0 \\ H^{n+1} &:= P^{n+1} - \Delta t\,C^{n+1} - P^n = 0 \end{align} Let us assume that you have an equation for $j^{n+1}$ that can be solved explicitly. Then the above system of equations has the form $$\mathbf{r}(\mathbf{u}^{n+1}) = 0$$ where $\mathbf{u}^{n+1} = (n_A^{n+1}, \theta_H^{n+1}, P^{n+1})$. Newton's method can then be expressed as $$\mathbf{u}^{n+1}_{k+1} = \mathbf{u}^{n+1}_k - \mathbf{T}^{-1}(\mathbf{u}^{n+1}_k)\,\mathbf{r}(\mathbf{u}^{n+1}_k)$$ where $$\mathbf{T}(\mathbf{u}^{n+1}_k) = \frac{\partial \mathbf{r}(\mathbf{u}^{n+1}_k)}{\partial \mathbf{u}^{n+1}_k} \,.$$ The components of the matrix $\mathbf{T}$ are found using the usual relations $$T_{ij} = \frac{\partial r_i}{\partial u_j} \,.$$ You already have initial values for $\mathbf{u}$. I'm not sure why you can't use those for your Newton iterations. If Newton iterations don't converge you may have to use an alternative algorithm (such as bisection/line search etc.)
• @GeoffOxberry Please correct me if I am wrong, the initial guess for $(\mathbf{u_k}^{n+1})$ is obtained using the forward/explicit euler method, right? – Ushnik Mukherjee Mar 28 '14 at 4:36
• A decent initial guess is the solution from the current time step (that is, $\mathbf{u}^{n}$). According to Ascher and Petzold, better guesses are available (so I assume they're probably given in a text like Hairer and Wanner, Lambert, or Butcher). – Geoff Oxberry Mar 28 '14 at 5:49