2

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}$.

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
Ushnik Mukherjee
  • 55
  • 4
  • 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 02:31
  • 1
    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 06:39

1 Answers1

3

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.)

Biswajit Banerjee
  • 942
  • 5
  • 16
  • This procedure is typically what is done, and it generally works well. If the Newton iterations don't converge, you could also decrease the time step by a given factor. A good example of how a library adapts to convergence failures can be found in the CVODE user manual. I strongly encourage everyone to call Newton's method from a library to take advantage of more robust tools than one could reasonably implement quickly; they also have the advantage of likely being more efficient and better tested. Implementing Newton yourself is a good exercise, but not especially useful in practice. – Geoff Oxberry Mar 28 '14 at 03:29
  • @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 04: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 05:49
  • @Biswajit, I will give it a try. Thank you for the help. – Ushnik Mukherjee Mar 28 '14 at 15:58