4. Discretization of the differential equation in time

4.1. Introduction of the theta method

A classic way to discretize a first-order differential equation is the \(\mathrm{\Theta }\) -method. Consider the following differential equation:

\(\{\begin{array}{c}\dot{y}(t)=\phi (t,y(t))\\ y(0)={y}_{0}\end{array}\) eq 4.1-1

The \(\Theta\) ‑method consists in discretizing the equation [éq 4.1-1] by a diagram with finite differences

\(\frac{1}{\mathrm{\Delta }t}({y}_{n+1}-{y}_{n})=\theta \mathrm{.}\varphi ({t}_{n+1},{y}_{n+1})+(1-\theta )\mathrm{.}\varphi ({t}_{n},{y}_{n})\) eq 4.1-2

where \({y}_{n+1}\) is an approximation of \(y({t}_{n+1})\), \({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 \ne 0\) the schema is said to be implicit.

4.2. 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)\\ {r}_{\text{vol}}^{+}={r}_{\text{vol}}(r,t+\Delta t)& {r}_{\text{vol}}^{-}={r}_{\text{vol}}(r,t)& {T}^{d+}={T}^{d}(r,t+\Delta t)& {T}^{d-}={T}^{d}(r,t)\\ {g}^{+}=g(r,{T}^{+})& {g}^{-}=g(r,{T}^{-})& & \end{array}\)

where \({T}^{d}(r,t)\) represents the temperature imposed on the domain border, as a function of time and space.

Let’s introduce the following function spaces:

\(\begin{array}{c}{V}_{{t}^{+}}=\left\{v\in {H}^{1}(\Omega ){v}_{\text{/}{\Gamma }_{1}}={T}^{d+}\right\}\\ {V}_{{t}^{-}}=\left\{v\in {H}^{1}(\Omega ){v}_{\text{/}{\Gamma }_{1}}={T}^{d-}\right\}\\ {V}_{0}=\left\{v\in {H}^{1}(\Omega ){v}_{\text{/}{\Gamma }_{1}}=0\right\}\end{array}\)

Since field \({T}^{-}\in {V}_{{t}^{-}}\) is supposed to be known, we are looking for \({T}^{+}\in {V}_{{t}^{+}}\) such as \(\forall v\in {V}_{0}\):

\(\begin{array}{cc}\underset{\Omega }{\mathrm{\int }}\frac{\beta ({T}^{+})\mathrm{-}\beta ({T}^{\mathrm{-}})}{\Delta t}\mathit{v.d}\Omega & +\underset{\Omega }{\mathrm{\int }}(\theta \lambda ({T}^{+})\mathrm{\nabla }{T}^{+}\mathrm{.}\mathrm{\nabla }v+(1\mathrm{-}\theta )\lambda ({T}^{\mathrm{-}})\mathrm{\nabla }{T}^{\mathrm{-}}\mathrm{.}\mathrm{\nabla }v)d\Omega \\ \mathrm{-}\underset{{\Gamma }_{2}}{\mathrm{\int }}(\theta {f}^{+}+(1\mathrm{-}\theta ){f}^{\mathrm{-}})\mathit{v.d}{\Gamma }_{2}& \mathrm{-}\underset{{\Gamma }_{3}}{\mathrm{\int }}(\theta {h}^{+}{T}_{\text{ext}}^{+}+(1\mathrm{-}\theta ){h}^{\mathrm{-}}{T}_{\text{ext}}^{\mathrm{-}}\mathrm{-}\theta {h}^{+}{T}^{+}\mathrm{-}(1\mathrm{-}\theta ){h}^{\mathrm{-}}{T}^{\mathrm{-}})\mathit{v.d}{\Gamma }_{3}\\ \mathrm{-}\underset{{\Gamma }_{4}}{\mathrm{\int }}(\theta {g}^{+}+(1\mathrm{-}\theta ){g}^{\mathrm{-}})\mathit{v.d}{\Gamma }_{4}& \mathrm{=}\\ \underset{\Omega }{\mathrm{\int }}(\theta {r}_{\text{vol}}^{+}+(1\mathrm{-}\theta ){r}_{\text{vol}}^{\mathrm{-}})\mathit{v.d}\Omega & +\underset{\Omega }{\mathrm{\int }}(\theta {r}_{\text{v}}({T}^{+})+(1\mathrm{-}\theta ){r}_{\text{v}}({T}^{\mathrm{-}}))\mathit{v.d}\Omega \end{array}\) eq 4.2-1

In order not to make the writing excessively heavier and insofar as the method is identical to the other terms, the radiation term has not been included in these equations (integral on \({\Gamma }_{5}\)).

By posing:

\(\begin{array}{c}({\text{hT}}_{\text{ext}}{)}^{\theta }\mathrm{=}\theta {h}^{+}{T}_{\text{ext}}^{+}+(1\mathrm{-}\theta ){h}^{\mathrm{-}}{T}_{\text{ext}}^{\mathrm{-}}\\ {f}^{\theta }\mathrm{=}\theta {f}^{+}+(1\mathrm{-}\theta ){f}^{\mathrm{-}}\\ {r}^{\theta }\mathrm{=}\theta {r}_{\text{vol}}^{+}+(1\mathrm{-}\theta ){r}_{\text{vol}}^{\mathrm{-}}\end{array}\) we finally get:

\(\begin{array}{c}\underset{\Omega }{\mathrm{\int }}\frac{\beta ({T}^{+})}{\Delta t}\mathit{v.d}\Omega +\theta \underset{\Omega }{\mathrm{\int }}\lambda ({T}^{+})\mathrm{\nabla }{T}^{+}\mathrm{.}\mathrm{\nabla }\mathit{v.}d\Omega +\theta \underset{{\Gamma }_{3}}{\mathrm{\int }}{h}^{+}{T}^{+}\mathit{v.d}{\Gamma }_{3}\\ \mathrm{-}\theta \underset{{\Gamma }_{4}}{\mathrm{\int }}g({T}^{+})\mathrm{.}\mathit{v.d}{\Gamma }_{4}\mathrm{-}\theta \underset{\Omega }{\mathrm{\int }}{r}_{v}({T}^{+})\mathrm{.}\mathit{v.d}\Omega \mathrm{=}{L}_{1}(v,{T}^{\mathrm{-}})\\ \mathrm{\forall }v\mathrm{\in }{V}_{0}\end{array}\) eq 4.2-2

where we put:

\(\begin{array}{c}{L}_{1}(v,{T}^{\mathrm{-}})\mathrm{=}\underset{\Omega }{\mathrm{\int }}\frac{\beta ({T}^{\mathrm{-}})}{\Delta t}\mathit{v.d}\Omega \mathrm{-}\underset{\Omega }{\mathrm{\int }}(1\mathrm{-}\theta )\lambda ({T}^{\mathrm{-}})\mathrm{\nabla }{T}^{\mathrm{-}}\mathrm{.}\mathrm{\nabla }v\mathrm{.}d\Omega +\underset{{\Gamma }_{2}}{\mathrm{\int }}{f}^{\theta }\mathit{v.d}{\Gamma }_{2}\\ +\underset{{\Gamma }_{3}}{\mathrm{\int }}(({\text{hT}}_{\text{ext}}{)}^{\theta }\mathrm{-}(1\mathrm{-}\theta ){h}^{\mathrm{-}}{T}^{\mathrm{-}})\mathit{v.d}{\Gamma }_{3}+\underset{\Omega }{\mathrm{\int }}{r}^{\theta }\mathit{v.d}\Omega \\ +(1\mathrm{-}\theta )\underset{{\Gamma }_{4}}{\mathrm{\int }}g({T}^{\mathrm{-}})\mathit{v.d}{\Gamma }_{4}+(1\mathrm{-}\theta )\underset{\Omega }{\mathrm{\int }}{r}_{v}({T}^{\mathrm{-}})\mathit{v.d}\Omega \end{array}\) eq 4.2-3

At a given moment of calculation, this term is known. In fact, only the temperature at the previous instant, \({T}^{-}\), as well as the values at the current instant of function known of the time, are involved.

In the case where the temperature distribution in the system at the initial moment is not provided, the stationary problem is solved. The terms of evolution disappear, \(\Theta =1\); the temperature field at the initial instant is given by:

\(\begin{array}{c}\underset{\Omega }{\mathrm{\int }}\lambda ({T}^{t\mathrm{=}0})\mathrm{\nabla }{T}^{t\mathrm{=}0}\mathrm{.}\mathrm{\nabla }\mathit{v.d}\Omega +\underset{{\Gamma }_{3}}{\mathrm{\int }}{h}^{t\mathrm{=}0}{T}^{t\mathrm{=}0}\mathit{v.d}{\Gamma }_{3}\mathrm{-}\underset{{\Gamma }_{4}}{\mathrm{\int }}g({T}^{t\mathrm{=}0})\mathit{v.d}{\Gamma }_{4}\\ \mathrm{=}\underset{{\Gamma }_{2}}{\mathrm{\int }}{f}^{t\mathrm{=}0}\mathit{v.d}{\Gamma }_{2}+\underset{{\Gamma }_{3}}{\mathrm{\int }}{h}^{t\mathrm{=}0}{T}_{\text{ext}}^{t\mathrm{=}0}\mathit{v.d}{\Gamma }_{3}+\underset{\Omega }{\mathrm{\int }}{r}^{t\mathrm{=}0}\mathit{v.d}\Omega \\ \mathrm{\forall }v\mathrm{\in }{V}_{0}\end{array}\) eq 4.2-4

The problem is finally written in a condensed form:

\(\{\begin{array}{c}\mathrm{Soit}{T}^{-}\in {V}_{{t}^{-}}\mathrm{connu},\mathrm{trouver}{T}^{+}\in {V}_{{t}^{+}}\mathrm{tel}\mathrm{que}\\ \forall v\in {V}_{0}\text{:}a(v,{T}^{+})={L}_{1}(v,{T}^{-})\end{array}\) eq 4.2-5