2
$\begingroup$

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

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$ – Geoff Oxberry Mar 28 '14 at 2:31
  • 1
    $\begingroup$ 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. $\endgroup$ – David Ketcheson Mar 28 '14 at 6:39
3
$\begingroup$

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

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$ – Geoff Oxberry Mar 28 '14 at 3:29
  • $\begingroup$ @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? $\endgroup$ – Ushnik Mukherjee Mar 28 '14 at 4:36
  • $\begingroup$ 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). $\endgroup$ – Geoff Oxberry Mar 28 '14 at 5:49
  • $\begingroup$ @Biswajit, I will give it a try. Thank you for the help. $\endgroup$ – Ushnik Mukherjee Mar 28 '14 at 15:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.