None .. An html document created by ipypublish

outline: ipypublish.templates.outline_schemas/rst_outline.rst.j2 with segments: - nbsphinx-ipypublish-content: ipypublish sphinx content

4. Résolution d’EDO en temps¶

Basile Marchand (basile.marchand@mines-paristech.fr)

4.1. Introduction¶

Intégration numérique ?

Considérons la forme générale d’une EDO du premier ordre

\[ \begin{align}\begin{aligned} \dot{\mathbf{q}} = \mathbf{B}( \mathbf{q}, t), \;\; \mathbf{q}\in \mathbb{R}^n\\On cherche alors à déterminer\end{aligned}\end{align} \]

$\mathbf{q}(t)\in\mathbb{R}^n,\; \forall t \in [0,T]$.

La question :

Comment faire travailler l’ordinateur à notre place ?

Pour information :

Il est possible également de traiter le cas d’EDO du second order par des méthodes dédiées (Newmark par exemple)

Mais ça sert à quoi ?

De manière générale, il existe de nombreux cas de figures où l’on ne sait pas trouver de solution analytique à une EDO.

En mécanique de nombreuses lois de comportement s’expriment sous la forme d’EDO en temps. Traduisant de fait l’impact de l’histoire sur l’état actuel.

Il existe toute un zoologie de méthodes d’intégration temporelle. Mais dans le cadre de ce cours on se limite aux méthodes les plus utilisées en mécanique à savoir :

La $\theta$-méthode
La méthode de Runge-Kutta

4.1.1. La notion de discrétisation en temps¶

Pour tous les schémas d’intégration il est nécessaire de se donner une discrétisation en temps. On découpe donc l’intervalle de temps $[t\_{init}, t\_{end}]$ en $n$ intervalles de même taille ou non

\[\Delta t_{i} = t_{i+1} - t_{i}\]

4.2. La $\theta$-méthode¶

L’idée de base

La $\theta$-méthode repose sur une évaluation de la dérivée temporelle par différences finies complétée par une évaluation de l’EDO en $t_\theta$

\[t_\theta = (1-\theta) t_i + \theta t_{i+1} \; ; \; \mathbf{q}_\theta = (1-\theta) \mathbf{q}_i + \theta \mathbf{q}_{i+1}\]

Le schéma numérique ainsi construit s’exprime :

\[\dfrac{ \mathbf{q}_{i+1} - \mathbf{q}_i }{\Delta t} = \mathbf{B}( \mathbf{q}_\theta , t_\theta )\]

La $\theta$-méthode regroupe entre autres les trois méthodes suivantes : - Euler explicite $ (:raw-latex:`\theta`=0) $ - Cranck-Nicholson $ (:raw-latex:`theta`=0.5 )$ - Euler implicite $ (:raw-latex:`theta`=1) $

La méthode d’Euler explicite :math:`(theta=0)`

On considère la dérivée au début de l’intervalle de temps

\[ \begin{align}\begin{aligned} \dfrac{ \mathbf{q}_{i+1} - \mathbf{q}_i }{\Delta t} = \mathbf{B}( \mathbf{q}_i , t_i )\\Ce qui permet d'écrire la suite :\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \mathbf{q}_{i+1} = \mathbf{q}_{i} + \Delta t \mathbf{B}( \mathbf{q}_i , t_i )\\La résolution est immédiate ! Connaissant :math:`\mathbf{q}_{i}` on\end{aligned}\end{align} \]

peut immédiatement déterminer $\mathbf{q}_{i+1}$.

En revanche attention au pas de temps. Car la méthode est conditionellemnt stable.

La méthode de Cranck-Nicholson :math:`(theta = 0.5)`

La méthode de Crank-Nicholson est similaire à une formule de quadrature par la méthode des trapèzes pour les EDO linéaires.

\[ \begin{align}\begin{aligned} \mathbf{q}_{i+1} = \mathbf{q}_i + \Delta t \mathbf{B}\left( \frac{\mathbf{q}_i+\mathbf{q}_{i+1}}{2}, \frac{t_i + t_{i+1}}{2} \right)\\Cette méthode est d'ordre 2\end{aligned}\end{align} \]

La méthode d’Euler implicite :math:`(theta = 1)`

On considère la dérivée à la fin de l’intervalle de temps

\[ \begin{align}\begin{aligned} \dfrac{ \mathbf{q}_{i+1} - \mathbf{q}_i }{\Delta t} = \mathbf{B}( \mathbf{q}_{i+1} , t_i )\\Ce qui permet d'écrire la suite :\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \mathbf{q}_{i+1} = \mathbf{q}_{i} + \Delta t \mathbf{B}( \mathbf{q}_{i+1} , t_{i+1} )\\Problème **non-linéaire** !!\end{aligned}\end{align} \]

Nécessite la résolution par une méthode de Newton et donc le calcul d’une matrice jacobienne

Synthèse dans le cas d’une EDO linéaire

Dans le cas où l’opérateur $\mathbf{B}$ est linéaire on peut nécessairement écrire l’EDO sous la forme :

\[\mathbf{C} \mathbf{\dot{q}} = \mathbf{K}\cdot \mathbf{q} + \mathbf{F}(t)\]

Dans ce cas la $\theta$-méthode se ramène à la résolution de systèmes linéaire :

\[\left( \frac{1}{\Delta t} \mathbf{C} + \theta \mathbf{K} \right) \cdot \mathbf{q}_{i+1} = \left( \frac{1}{\Delta t} \mathbf{C} + (1-\theta) \mathbf{K} \right) \cdot \mathbf{q}_{i} + (1-\theta) \mathbf{F}_{i} + \theta \mathbf{F}_{i+1}\]

Dans le cas de pas de temps constant il est fortement recommandé de calculer une factorisation de $\left( \frac{1}{\Delta t} \mathbf{C} + \theta \mathbf{K} \right)$ en phase préliminaire

Retour sur le Euler implicite

\[ \begin{align}\begin{aligned} \mathbf{q}_{i+1} = \mathbf{q}_{i} + \Delta t \mathbf{B}( \mathbf{q}_{i+1} , t_{i+1} )\\Le problème à résoudre est non-linéaire !\end{aligned}\end{align} \]

A chaque incrément on doit résoudre

\[\mathbf{G}(\mathbf{q}_{i+1} ) = \mathbf{q}_{i+1} - \mathbf{q}_{i} - \Delta t \mathbf{B}( \mathbf{q}_{i+1} , t\_{i+1} ) = 0\]

On applique alors une méthode de Newton :

\[\mathbf{q}_{i+1}^{(k+1)} = \mathbf{q}_{i+1}^{(k)} + \mathcal{J}^{-1}( \mathbf{G}( \mathbf{q}_{i+1}^{(k)} ) ) \cdot \mathbf{G}(\mathbf{q}_{i+1}^{(k)})\]

Avec $\mathcal{J}( \mathbf{F}( \mathbf{q}_{i+1}^{(k)} ) )$ la matrice Jacobienne définie par :

\[\mathcal{J}( \mathbf{G}( \mathbf{q}_{i+1}^{(k)} )) = \left. \frac{ \partial G }{\partial \mathbf{q}}\right\vert_{ \mathbf{q}_{i+1}^{(k)}}\]

Dans la pratique deux approches pour calculer la matrice jacobienne : - Calcul analytique - Par perturbation

Pour plus de détails revoir le cours précédent sur la résolution de systèmes non-linéaires.

4.2.1. Quelques propriétés mathématiques¶

Existence, unicité et convergence

Existence et unicité : en associant une condition initiale $\mathbf{q}_0$ ont obtient un problème de Cauchy

\[\begin{cases} \dot{\mathbf{q}} &= \mathbf{B}( \mathbf{q}, t) \newline \mathbf{q}(t_0) &= \mathbf{q}_0 \end{cases}\]

Le thérorème de Cauchy-Lipchitz stipule alors que :

Si on connait :math:`mathbf{q}_0` et si :math:`mathbf{B}(mathbf{q}, t)` est localement lipchitzienne par rapport à la première variable alors la solution existe et elle est unique

\[\exists k \in \mathbb{R}^{+},\, \forall t > t\_0,\, \forall (\mathbf{q}^1, \mathbf{q}^2)\]

\[\Vert \mathbf{B}( \mathbf{q}^2, t ) - \mathbf{B}(\mathbf{q}^1, t) \Vert \leq k \Vert \mathbf{q}^2 - \mathbf{q}^1 \Vert\]

Théorème de Lax :

Un schéma numérique aux différences finis est convergent s’il est **consistant* et stable.*

Consistance d’un schéma numérique

Une méthode d’intégration temporelle est consistante si et seulement si elle engendre une suite $\lbrace \mathbf{q}_{i} \rbrace_{i=1,...,n}$ telle que

\[\lim_{\Delta t \rightarrow 0} \frac{ \mathbf{q}(t_{i+1}) - \mathbf{q}(t_{i}) }{\Delta t} = \mathbf{\dot{q}}(t_i)\]

La consistance d’une méthode d’intégration est une condition nécessaire de convergence

Stabilité d’un schéma numérique

Une méthode d’intégration est stable si une petite variation du vecteur d’état à l’instant $t_i$ n’entraîne que des variations consécutives qui ne sont pas croissantes au cours du temps

Stabilité de la :math:`theta`-méthode : cas de la diffusion

Soit un barreau soumis en $x=0$ à une température imposée et en $x=L$ à un flux imposé.

La discrétisation d’un tel problème permet d’écrire le système

\[ \begin{align}\begin{aligned} \mathbf{C}\cdot \mathbf{\dot{q}} + \mathbf{K}\cdot \mathbf{q} = \mathbf{F}\\Où les différents opérateurs sont définis :\end{aligned}\end{align} \]

\[C_{ij} = \rho C_p \frac{h}{4} ( 2 \delta_{ij} + \delta_{(i-1)j} + \delta_{(i+1)j} )\; ; \; C_{1j} = \rho C_p \frac{h}{4} (\delta_{1j} + \delta_{2j} )\]

\[K_{ij} = \frac{k}{h} ( 2 \delta_{ij} - \delta_{(i-1)j} - \delta_{(i+1)j} )\; ; \; K_{1j} = \frac{k}{h} (\delta_{1j} - \delta_{2j} ) \; ; \; F_{i} = T^d \delta_{(N)i}\]

Avec $\delta_{ij}$ le symbole de Kronecker

Illustration

*TODO*

Astuce : changement de variable par projection dans une base de vecteurs propres:

\[\exists \boldsymbol{\phi}_k \in \mathbb{R}^N, \: \lambda_k \in \mathbb{R}^+, \; - \lambda_k \: \mathbf{C} \: \boldsymbol{\phi}_k + \mathbf{K} \boldsymbol{\phi}_k = 0\]

\[ \begin{align}\begin{aligned} \exists \mathbf{y} \in \mathbb{R}^N, \; \mathbf{q} = \sum_{k=1}^N \phi_k \: y_k = \boldsymbol{\phi} \mathbf{y}, \; \boldsymbol{\phi} = [\phi_1, ...,\phi_N] \in \mathbb{R}^{N\times N}\\On en déduit:\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \dot{\mathbf{y}} = \boldsymbol{\phi}^T \mathbf{C} \mathbf{B}( \boldsymbol{\phi} \mathbf{y} , t)\\Donc:\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \dot{\mathbf{y}} = - \boldsymbol{\phi}^T \mathbf{K} \boldsymbol{\phi} \mathbf{y} + \boldsymbol{\phi}^T \mathbf{F}\\avec\end{aligned}\end{align} \]

$(\boldsymbol{\phi}^T \mathbf{K} \boldsymbol{\phi})_{ij} = \lambda_{i} \delta_{ij}$. Il s’agit donc d’un système d’équations découplées.

Considérons l’effet d’une perturbation à la ligne $k$ :

\[\dot{y}^{(k)} = - \lambda^{(k)} y^{(k)} + F^{(k)}\]

Appliquons la $\theta$-méthode à la ligne $k$:

\[ \begin{align}\begin{aligned} \dot{y}^{(k)} = - \lambda^{(k)} y^{(k)} + F^{(k)}\\devient\end{aligned}\end{align} \]

\[\frac{y_{i+1}^{(k)} - y_{i}^{(k)}}{\Delta t} = -\lambda^{(k)} ( (1-\theta) y_{i}^{(k)} + \theta y_{i+1}^{(k)}) + (1-\theta) F_{i}^{(k)} + \theta F_{i+1}^{(k)}\]

\[ \begin{align}\begin{aligned} \Rightarrow ( 1 + \Delta t \theta \lambda^{(k)} ) y_{i+1}^{(k)} = (1 - \Delta t \: (1-\theta) \: \lambda^{(k)}) \: y_{i}^{(k)} + \Delta t \: ( (1-\theta) \: F_{i}^{(k)} + \theta \: F_{i+1}^{(k)})\\Considérons l'effet d'un petite perturbation :math:`\delta y_{i}^{(k)}`\end{aligned}\end{align} \]

en $t_i$ pour la ligne $k$. On obtient l’équation d’évolution de la perturbation:

\[ \begin{align}\begin{aligned} ( 1 + \Delta t \: \theta \: \lambda^{(k)} ) \: \delta y_{i+1}^{(k)} = (1 - \Delta t \: (1-\theta) \: \lambda^{(k)}) \: \delta y_{i}^{(k)}\\Le coefficient d'amplification est donc:\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \mathcal{A}^{(k)} = \frac{1 + \Delta t \: (\theta - 1) \: \lambda^{(k)}}{1 + \Delta t \: \theta \: \lambda^{(k)}}, \quad \delta y_{i+1}^{(k)} = \mathcal{A}^{(k)} \: \delta y_{i}^{(k)}\\Si\end{aligned}\end{align} \]

$\max_{k} | \frac{1 + \Delta t \: (\theta - 1) \: \lambda^{(k)}}{1 + \Delta t \: \theta \: \lambda^{(k)}} | < 1$ alors il n’y a pas d’amplification de perturbation et donc le schéma est stable.

Cas linéaire

\[\frac{1 + \Delta t \: (\theta - 1) \: \lambda^{(k)}}{1 + \Delta t \: \theta \: \lambda^{(k)}} = 1 - \frac{\lambda^{(k)}}{1 + \Delta t \: \theta \: \lambda^{(k)}} < 1\]

Existe-t-il une valeur de $\Delta t > 0$, notée $\Delta t_{c}^{(k)}$ telle que:

\[\frac{1 + \Delta t_{c}^{(k)} \: (\theta - 1) \: \lambda^{(k)}}{1 + \Delta t_{c}^{(k)} \: \theta \: \lambda^{(k)}} = -1\]

\[\Rightarrow \: 1 + \Delta t_{c}^{(k)} \: (\theta - 1) \: \lambda^{(k)} + 1 + \Delta t_{c}^{(k)} \: \theta \: \lambda^{(k)} = 0\]

\[ \begin{align}\begin{aligned} \Rightarrow \: 2 + \Delta t_{c}^{(k)} \: ( 2 \: \theta - 1) \: \lambda^{(k)} = 0\\Si $ :raw-latex:`\theta < \frac{1}{2}$, il existe $\Delta t_{c} >` 0 $\end{aligned}\end{align} \]

tel que:

\[\Delta t_{c} = \min_k \frac{1}{(\frac{1}{2} - \theta) \: \lambda^{(k)} } = \frac{1}{(\frac{1}{2} - \theta) \: \max_k \lambda^{(k)} }\]

Si $\theta \ge \frac{1}{2}$, on a $\max_k | \mathcal{A}_k | < 1$. Le schéma est stable.

Si $\theta < \frac{1}{2}$, comment choisir $\Delta t$?

Nous souhaitons avoir $\mathcal{A}^{(k)} > -1$:

\[1 + \Delta t \: (\theta - 1) \: \lambda^{(k)} > -1 - \Delta t \: \theta \: \lambda^{(k)}\]

\[\Rightarrow \: 2 + \Delta t \: (2 \: \theta - 1) \: \lambda^{(k)} > 0\]

\[\Rightarrow \: \Delta t \: \lambda^{(k)} < \frac{-2}{(2 \: \theta - 1)}\]

\[\Rightarrow \: \Delta t < \frac{-2}{(2 \: \theta - 1) \: \lambda^{(k)}} \: \forall k=1, ...,N\]

Il faut donc \(\Delta t < \Delta t_{c}\).
La \(\theta\)-méthode est stable pour
\(\theta \ge \frac{1}{2}\). Donc Euler implicite et
Crank-Nicholson sont des schémas stables.
La \(\theta\)-méthode est conditionnellement stable pour
\(\theta < \frac{1}{2}\). Il faut déterminer un pas de temps
critique \(\Delta t\_c\) et choisir \(\Delta t < \Delta t\
_c\).

Retour sur le problème de diffusion

Dans le cas du problème parabolique de diffusion thermique, les valeurs propres $\lambda^{(k)}$ dépendent de $h$, la taille des éléments. On peut alors introduire le coefficient CFL (Courant–Friedrichs–Lewy):

\[ \begin{align}\begin{aligned} \beta = \frac{k}{\rho \: C_p} \: \frac{\Delta t}{h^2}\\Une estimation de :math:`\max_k \lambda^{(k)}` peut être obtenue avec\end{aligned}\end{align} \]

la plus grande valeur propre du plus petit élément pris de façon isolée (théorème de Irons et Threhane).

On obtient ici:

\[\max_k \lambda^{(k)} \approx \frac{8 \: k}{\rho \: C_p \: h^2}\]

On en déduit:

\[ \begin{align}\begin{aligned} \Delta t_c = \frac{ \rho \: C_p \: h^2 }{\left(\frac{1}{2} - \theta\right) \: 8 \: k }\\Plus le plus petit élément est petit, plus il faut réduire le pas de\end{aligned}\end{align} \]

temps pour $\theta < \frac{1}{2}$.

Dans le cas étudié, on a:

\[ \begin{align}\begin{aligned} \frac{\Delta t}{\Delta t_c} = \left(\frac{1}{2} - \theta \right) \: 8 \: \beta\\La condition CFL est donc:\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \left(\frac{1}{2} - \theta\right) \: 8 \: \frac{k}{\rho \: C_p} \: \frac{\Delta t}{h^2} < 1\\pour avoir un schéma stable. Si :math:`\theta < \frac{1}{2}` plus\end{aligned}\end{align} \]

$h$ est petit pour capter de forts gradients, plus il faut $\Delta t$ petit.

4.3. Méthode de Runge-Kutta¶

L’idée de base

Parfois l’emploi de la méthode d’Euler implicite est impossible et la résolution avec un Euler explicite couteraît trop cher.

Principe de la méthode : - Prédiction - Correction

*Runge-Kutta d’ordre 2 (RK2)*

Le principe est d’estimer la dérivée au milieu du pas d’intégration

La prévision est construite à partir de deux évaluations de la fonction $\mathbf{B}$

\[\begin{cases} \delta_1 \mathbf{q} &=& \Delta t \: \mathbf{B}( \mathbf{q}_i, \: t_i) \newline \delta_2 \mathbf{q} &=& \Delta t \: \mathbf{B}( \mathbf{q}_i + \frac{\delta_1 \mathbf{q}}{2} , \: t_i+\frac{ \Delta t}{2}) \newline \mathbf{q}_{i+1} &=& \mathbf{q}_i + \delta_2 \mathbf{q} + o(\Delta t^2) \end{cases}\]

$\delta_1 \mathbf{q}$ est l’incrément calculé par la méthode d’Euler explicite.

Ce schéma annule le premier ordre de l’erreur.

*Runge-Kutta d’ordre 4 (RK4)*

La prévision est construite avec quatre évaluations de la fonction $\mathbf{B}$:

\[\begin{cases} \delta_1 \mathbf{q} &=& \Delta t \: \mathbf{B}( \mathbf{q}_i, \: t_i) \newline \delta_2 \mathbf{q} &=& \Delta t \: \mathbf{B}( \mathbf{q}_i + \frac{\delta_1 \mathbf{q}}{2} , \: t_i+\frac{ \Delta t }{2}) \newline \delta_3 \mathbf{q} &=& \Delta t \: \mathbf{B}( \mathbf{q}_i + \frac{\delta_2 \mathbf{q}}{2} , \: t_i+\frac{ \Delta t }{2}) \newline \delta_4 \mathbf{q} &=& \Delta t \: \mathbf{B}( \mathbf{q}_i + \delta_3 \mathbf{q}, \: t_i+\Delta t) \newline \mathbf{q}_{i+1} &=& \mathbf{q}_i + \frac{ \delta_1 \mathbf{q}}{6} + \frac{ \delta_2 \mathbf{q}}{3} + \frac{ \delta_3 \mathbf{q}}{3} + \frac{ \delta_4 \mathbf{q}}{6} + o(\Delta t^4) \end{cases}\]

Remarque sur la complexité numérique:

Les pas de temps utilisés pour ce schéma ne sont pas nécessairement deux fois plus grands que ceux utilisés pour le schéma du second ordre. Mais s’ils le sont la complexité numérique du schéma est moins grande que celle de RK2.

Contrôle du pas de temps

Approche empirique pour évaluer la convergence de l’approximation: 1. diviser la taille du pas de temps par 2 2. comparer les prévisions 3. recommencer tant que l’écart entre les deux dernières prévisions est jugé trop important.

Mais en général l’erreur d’approximation ne converge pas de façon uniforme sur tout l’intervalle de temps.

Il est parfois utile de concentrer les petits pas de temps en certains endroits de l’intervalle de temps.

Une méthode adaptative de pas de temps permet de mieux répartir sur l’intervalle de temps les erreurs d’approximation en ayant le souci de réduire les temps de calcul. L’adaptation des pas de temps permet aussi de faciliter la linéarisation de certains problèmes d’évolution non linéaires.

*Approche par extrapolation locale [Press94].*

En doublant le pas de temps, on estime la solution exacte $\overline{\mathbf{q}}$ de deux façons par RK4:

\[\overline{\mathbf{q}}(t_i + (2\Delta t)) = \mathbf{q}_{i+1} + o( (2 \Delta t)^4)\]

\[ \begin{align}\begin{aligned} \overline{\mathbf{q}}(t_i + \Delta t + \Delta t) = \mathbf{q}_{i+1} + o( \Delta t^4)\\On suppose qu'il existe $ c :raw-latex:`\in `:raw-latex:`\mathbb{R}` $\end{aligned}\end{align} \]

tel que:

\[\overline{\mathbf{q}}(t_i + (2\Delta t)) = \mathbf{q}_{i+1} + c \: (2 \Delta t)^5 + o( (2 \Delta t)^5)\]

\[ \begin{align}\begin{aligned} \overline{\mathbf{q}}(t_i + \Delta t + \Delta t) = \widetilde{\mathbf{q}}_{i+1} + 2 \: c \: \Delta t^5 + o( \Delta t^5)\\En supposant négligeable l'ordre 5, on estime l'erreur selon:\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned} \Delta \mathbf{q} = \widetilde{\mathbf{q}}_{i+1} - \mathbf{q}_{i+1}, \: \eta = \max_i \vert \Delta \mathbf{q}_i \vert, \quad \eta = 30 \: c \: \Delta t^5\\On souhaite connaitre le pas de temps :math:`\Delta t^\star` qui permet\end{aligned}\end{align} \]

d’avoir une erreur $\eta^\star = \epsilon\_{tol}$ si $\eta > \epsilon\_{tol}$, tel que:

\[ \begin{align}\begin{aligned} \eta^\star = 30 \: c \: \Delta t^{\star 5}\\On en déduit le pas de temps adapté :\end{aligned}\end{align} \]

$\Delta t^\star = \Delta t \: (\frac{\eta}{\epsilon\_{tol}})^{1/5}$

Le contrôle du pas de temps passe par une estimation de l’erreur

*Il y a donc un coût additionnel !!*

C’est néanmoins très utile voir impératif

*Cela permet d’optimiser les coûts de calcul*

4.4. Modèle de Maxwell¶

\[ \begin{align}\begin{aligned} \dot{\varepsilon} = \frac{1}{E} \dot{\sigma} + \frac{1}{\eta} \sigma\\avec :math:`\varepsilon(t)=10^{-3} t` imposé et :math:`\sigma(t=0) = 0`\end{aligned}\end{align} \]

Calculer la solution analytique
Déterminer la solution par une approche Euler explicite
Déterminer la solution par une approche Euler implicite
Faire de même avec la méthode RK4
Utiliser scipy.integrate.ode pour calculer la solution de ce problème et comparer

Données : $\eta = 1000\,MPa.s$, $E=20000\,MPa$, $T=10\,s$