5. Discretization of the differential equation in time#
A classic way to discretize a first-order differential equation is to use a \(\theta\) -method. Consider the following differential equation:
\(\{\begin{array}{c}\frac{\partial }{\partial t}y(t)=\varphi (t,y(t))\\ y(0)={y}_{0}\end{array}\)
The \(\theta -\) method consists in discretizing the equation by a scheme with finite differences
\(\frac{1}{\Delta t}({y}_{n+1}-{y}_{n})=\theta \varphi ({t}_{n+1},{y}_{n+1})+(1-\theta )\varphi ({t}_{n},{y}_{n})\)
where \({y}_{n+1}\) is an approximation of \(y({t}_{n+1})\), with \({y}_{n}\) being assumed to be known
and \(\theta\) is the method parameter, \(\theta \in \left[\mathrm{0,1}\right]\).
Note:
if \(\theta =0\) the schema is said to be explicit,
if \(\theta \ge 0\) the schema is said to be implicit.
5.1. Accuracy of the method#
Assume \(y\) that is sufficiently regular (at least 3 times differentiable), by a Taylor expansion at point \({t}_{n}\) we obtain:
\(y({t}_{n+1})-y({t}_{n})=\Delta ty\text{'}({t}_{n})+\frac{\Delta {t}^{2}}{2}y\text{''}({t}_{n})+O(\Delta {t}^{2})\)
and
\(\begin{array}{cc}\theta \varphi ({t}_{n+1},y({t}_{n+1}))+(1-\theta )\varphi ({t}_{n},y({t}_{n}))& =\theta y\text{'}({t}_{n+1})+(1-\theta )y\text{'}({t}_{n})\\ & =y\text{'}({t}_{n+1})+\theta (y\text{'}({t}_{n+1})-y\text{'}({t}_{n}))\\ & =y\text{'}({t}_{n})+\theta \Delta ty\text{''}({t}_{n})+O(\Delta {t}^{2})\end{array}\)
The solution therefore verifies approximately:
\(\frac{1}{\Delta t}(y({t}_{n+1})-y({t}_{n}))=\theta \varphi ({t}_{n+1},y({t}_{n+1}))+(1-\theta )\varphi ({t}_{n},y({t}_{n}))+(\frac{1}{2}-\theta )\Delta ty\text{''}({t}_{n})+O(\Delta {t}^{2})\)
The diagram is of order 1 in time if \(\theta \ne \frac{1}{2}\), and of order 2 if \(\theta =\frac{1}{2}\) (Crank-Nicolson diagram).
5.2. Stability of the method#
Consider the following differential equation:
\(\{\begin{array}{ccc}y\text{'}=-\lambda y& t\ge 0& \lambda \in R\\ y(0)={y}_{0}& & \end{array}\)
Using the \(\theta\) -method in this differential equation we get:
\({y}_{n+1}=\frac{1-(1-\theta )\lambda \Delta t}{1+\theta \lambda \Delta t}{y}_{n}\text{}0\le n\le N-1\)
Or again:
\({y}_{n+1}={r}^{n}(\lambda \Delta t){y}_{0}\text{avec}r(x)=\frac{1-(1-\theta )x}{1+\theta x}\)
The approximate solution \({y}_{n}\) must be limited (the exact solution of the initial problem being so), which imposes the following condition:
\(\mid r(\lambda \Delta t)\mid \le 1\)
By studying the variations of function \(r(x)\), it is easy to see that:
if \(\theta \ge\) the condition is true regardless of \(\Delta t\), the schema is unconditionally stable;
if \(\theta <\) the condition is only true if \(\Delta t\le \frac{2}{\lambda (1-2\theta )}\), the schema is conditionally stable.
In the THER_LINEAIRE [U4.33.01] command, the parameter \(\theta\) is data that can be provided by the user, the default value is set to 0.57. This value has the reputation of being preferable to the Crank-Nicolson value (0.5) and « optimal » for quadratic interpolations, but we have not found a trace of the justifications.
5.3. Application to the heat equation#
Let’s use the \(\theta\) -method in the variational formulation of the heat equation, where we asked:
\(\begin{array}{cccc}{T}^{+}=T(r,t+\Delta t)& {T}^{-}=T(r,t)& {h}^{+}=h(r,t+\Delta t)& {h}^{-}=h(r,t)\\ {f}^{+}=f(r,t+\Delta t)& {f}^{-}=f(r,t)& {T}_{\text{ext}}^{+}={T}_{\text{ext}}(r,t+\Delta t)& {T}_{\text{ext}}^{-}={T}_{\text{ext}}(r,t)\\ {s}^{+}=s(r,t+\Delta t)& {s}^{-}=s(r,t)& {T}_{1}^{+}={T}_{1}(r,t+\Delta t)& {T}_{1}^{-}={T}_{1}(r,t)\end{array}\)
Let’s introduce the following function spaces:
\(\begin{array}{c}{V}_{{t}^{+}}=\left\{v\in {H}^{1}(\Omega )\text{}{v}_{/{\Gamma }_{1}}={T}_{1}(r,{t}^{+})\right\}\\ {V}_{{t}^{-}}=\left\{v\in {H}^{1}(\Omega )\text{}{v}_{/{\Gamma }_{1}}={T}_{1}(r,{t}^{-})\right\}\\ {V}_{0}=\left\{v\in {H}^{1}(\Omega )\text{}{v}_{/{\Gamma }_{1}}=0\right\}\end{array}\)
Since field \({T}^{-}\in {V}_{{t}^{-}}\) is supposed to be known, we’re looking for \({T}^{+}\in {V}_{{t}^{+}}\):
\(\begin{array}{}\underset{\Omega }{\int }\rho {C}_{p}\frac{{T}^{+}-{T}^{-}}{\Delta t}vd\Omega +\underset{\Omega }{\int }(\theta \lambda \nabla {T}^{+}\text{.}\nabla v+(1-\theta )\lambda \nabla {T}^{-}\text{.}\nabla v)d\Omega \\ -\underset{{\Gamma }_{2}}{\int }(\theta {f}^{+}+(1-\theta ){f}^{-})vd{\Gamma }_{2}-\underset{{\Gamma }_{3}}{\int }(\theta {h}^{+}{T}_{\text{ext}}^{+}+(1-\theta ){h}^{-}{T}_{\text{ext}}^{-}-\theta {h}^{+}{T}^{+}-(1-\theta ){h}^{-}{T}^{-})vd{\Gamma }_{3}\\ =\underset{\Omega }{\int }(\theta {s}^{+}+(1-\theta ){s}^{-})vd\Omega \\ \forall v\in {V}_{0}\end{array}\)
By posing:
\(\begin{array}{}({\text{hT}}_{\text{ext}}{)}^{\theta }=\theta {h}^{+}{T}_{\text{ext}}^{+}+(1-\theta ){h}^{-}{T}_{\text{ext}}^{-}\\ {f}^{\theta }=\theta {f}^{+}+(1-\theta ){f}^{-}\end{array}\)
we finally get:
\(\begin{array}{}\underset{\Omega }{\int }\frac{\rho {C}_{p}}{\Delta t}{T}^{+}vd\Omega +\underset{\Omega }{\int }\theta \lambda \nabla {T}^{+}\text{.}\nabla vd\Omega +\underset{{\Gamma }_{3}}{\int }\theta {h}^{+}{T}^{+}vd{\Gamma }_{3}\\ =\underset{\Omega }{\int }\frac{\rho {C}_{p}}{\Delta t}{T}^{-}vd\Omega -\underset{\Omega }{\int }(1-\theta )\lambda \nabla {T}^{-}\text{.}\nabla vd\Omega +\underset{{\Gamma }_{2}}{\int }{f}^{\theta }vd{\Gamma }_{2}\\ +\underset{{\Gamma }_{3}}{\int }(({\text{hT}}_{\text{ext}}{)}^{\theta }-(1-\theta ){h}^{-}{T}^{-})vd{\Gamma }_{3}+\underset{\Omega }{\int }(\theta {s}^{+}+(1-\theta ){s}^{-})vd\Omega \\ \forall v\in {V}_{0}\end{array}\)