In [1]:
from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read()) 
Out[1]:
In [2]:
import sys #only needed to determine Python version number
print('Python version ' + sys.version)
Python version 3.8.5 (default, Jan 27 2021, 15:41:15) 
[GCC 9.3.0]

M62_CM5 : Convergence et A-stabilité des schémas multipas linéaires

$\newcommand{\R}{\mathbb{R}}$ $\newcommand{\N}{\mathbb{N}}$ $\newcommand{\dt}{\mathrm{d}t}$ $\newcommand{\dx}{\mathrm{d}x}$ $\newcommand{\dtau}{\mathrm{d}\tau}$ $\newcommand{\CC}[1]{\mathscr{C}}$

Convergence d'un schéma multipas linéaire

Considérons le problème de Cauchy

trouver une fonction $y \colon I\subset \mathbb{R} \to \mathbb{R}$ définie sur un intervalle $I=[t_0,T]$ telle que $$\begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=[t_0,T],\\ y(t_0) = y_0, \end{cases}$$ avec $y_0$ une valeur donnée et supposons que l'on ait montré l'existence et l'unicité d'une solution $y$ pour $t\in I$.

Considérons une méthode numérique à $q=p+1$ étapes (=pas) linéaire. Elle s'écrit sous la forme générale $$ u_{n+1} = \sum_{j=0}^p a_ju_{n-j} + h\sum_{j=0}^p b_j\varphi(t_{n-j},u_{n-j}) + hb_{-1}\varphi(t_{n+1},u_{n+1}), \qquad n=p,p+1,\dots,N-1 $$ où les $\{a_k\}$ et $\{b_k\}$ sont des coefficients donnés et $p\ge0$ un entier.

Les données initiales $u_0, u_1, \dots, u_p$ doivent être fournies.
Mis à part $u_0$, qui est égale à $y_0$ donné, les autres valeurs, $u_1, \dots, u_p$ peuvent être obtenues à l’aide de méthodes suffisament précises.
Si $b_{-1}=0$ la méthode est explicite, sinon implicite.

Les méthodes AB$_q$, AM$_q$, N$_q$, MS$_q$ et BFD$_q$ sont toutes des méthodes numériques à $q=p+1$ étapes (=pas) linéaires. En revanche, les méthodes EM et Heun ne rentrent pas dans ce cadre: elles sont des méthodes de la famille RK.

Exemples :

  • Euler explicite
    $$ u_{n+1}=u_n+h\varphi(t_n,u_n) \quad\implies\quad p=0,\ a_0=1,\ b_{-1}=0,\ b_0=1; $$
  • Euler implicite
    $$ u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1}) \quad\implies\quad p=0,\ a_0=1,\ b_{-1}=1,\ b_0=0; $$
  • Crank-Nicolson $$ u_{n+1}=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\right) \quad\implies\quad p=0,\ a_0=1,\ b_{-1}=\frac{1}{2},\ b_0=\frac{1}{2}. $$

Polynômes caractéristiques

Réécrivons le schéma multistep comme suit: $$ u_{n+1} -\sum_{j=0}^p a_ju_{n-j}= h\left(\sum_{j=0}^p b_j\varphi(t_{n-j},u_{n-j})+b_{-1}\varphi(t_{n+1},u_{n+1})\right). $$ Formellement,

  • pour construire le premier polynômes caractéristiques on ne regarde que le LHS et les $u_{n-j}$ deviennent des $r^{p-j}$,
  • pour construire le second polynômes caractéristiques on ne regarde que le RHS divisé par $h$ et les $\varphi(t_{n-j},u_{n-j})$ deviennent des $r^{p-j}$.
    Définition: le polynôme $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j} $$ est appelé premier polynôme caractéristique associé à la méthode numérique.
    Définition: le polynôme $$ \sigma(r)=b_{-1}r^{p+1}+\sum_{j=0}^p b_jr^{p-j} $$ est appelé second polynôme caractéristique associé à la méthode numérique.

Exemples :

  • Euler explicite: $\varrho(r)=r-1$ et $\sigma(r)=1$
  • Euler implicite: $\varrho(r)=r-1$ et $\sigma(r)=r$
  • Crank-Nicolson: $\varrho(r)=r-1$ et $\sigma(r)=\frac{1}{2}r+\frac{1}{2}$
  • AB$_q$ ou AM$_q$: $\varrho(r)=r^{q+1}-r^{q}$
  • N$_q$ ou MS$_q$: $\varrho(r)=r^{q+1}-r^{q-1}$
  • BDF$_2$: $\varrho(r)=r^2-\frac{4}{3}r+\frac{1}{3}$ et $\sigma(r)=\frac{2}{3}$

Consistance

L’erreur de troncature locale de la méthode multi-pas s'écrit $$ \tau_n(h)=\frac{1}{h}\left(y_{n+1}-\sum_{j=0}^p a_jy_{n-j}-h\sum_{j=0}^p b_j\varphi(t_{n-j},y_{n-j})-hb_{-1}\varphi(t_{n+1},y_{n+1})\right). $$ La méthode est dite consistante si $\tau(h)=\max_n |\tau_n(h)|$ tend vers zéro quand $h$ tend vers zéro. Par un développement de Taylor assez fastidieux, on peut montrer que cette condition est équivalente à la condition suivante:

Proposition La méthode est consistante ssi $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=1 \end{cases} $$
Notons que la première condition de consistance équivaut à $\varrho(1)=0$.
Proposition La méthode est d'ordre $\omega\ge1$ ssi elle est consistante et, de plus, $$ \sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j=1 \quad i=2,\dots,\omega. $$

Exemple: schéma AB$_2$

  • Schéma: $u_{n+1}=u_{n}+ \frac{h}{2} \big(3\varphi_{n}-\varphi_{n-1}\big)$
  • Coefficients: $q=2$, $p=1$, $a_0=1$, $a_1=0$, $b_0=\frac{3}{2}$, $b_1=-\frac{1}{2}$, $b_{-1}=0$
  • Consistance: $\begin{cases}\sum_{j=0}^p a_j=1+0=1\\-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-(0a_0+1a_1)+(b_{-1}+b_0+b_1)=1\end{cases}$
  • Ordre $\omega$:
    • pour $i=2$ on a $\sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j=\sum_{j=0}^p j^2 a_j+2\sum_{j=-1}^p j b_j=(0^2a_0+1^2a_1)+2(-1b_{-1}+0b_0+1b_1)=1$ donc $\omega\ge2$,
    • pour $i=3$ on a $\sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j=\sum_{j=0}^p (-j)^3 a_j+3\sum_{j=-1}^p j^2 b_j=(0^3a_0-1^3a_1)+3((-1)^2b_{-1}+0^2b_0+1^2b_1)=3$ donc $\omega<3$
    • donc l'ordre $\omega$ est $2$

Zéro-stabilité

Proposition [condition des racines] Soient $r_0, r_1,\dots, r_p$ les $p+1$ racines du polynôme $\varrho$. On peut montrer que la méthode est **zéro-stable** ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1. \end{cases} $$

En particulier:

  • les méthodes à un pas (EE,EI,CN) vérifient $\varrho(r)=r-1$: le polynôme $\varrho$ n'a qu'une racine, c'est $r=1$. Toutes ces méthodes sont donc zéro-stables. On peut donc conclure que pour les méthodes à un pas, "consistance $=$ convergence";
  • les méthodes de Adam à $q$ pas vérifient $\varrho(r)=r^{p+1}-r^p=r^p(r-1)$: il n'y a que deux racines, la racine $1$ de multiplicité $1$ et la racine $0$ de multiplicité $p$. Toutes ces méthodes sont donc zéro-stables;
  • les méthodes de Nystrom et Milne-Simpson à $q$ pas vérifient $\varrho(r)=r^{p+1}-r^{p-1}=r^{p-1}(r-1)(r+1)$: il n'y a que trois racines, les racines $\pm1$ de multiplicité $1$ et la racine $0$ de multiplicité $p-1$. Toutes ces méthodes sont donc zéro-stables;
  • enfin, on peut montrer que les méthodes BDF sont zéro-stables ssi $p\le5$ (voir les calculs sympy ci-dessous).
Théorème (Dahlquist 1956). Une méthode multipas est convergente si et seulement si elle est consistante et (zero)-stable.

La zero stabilité implique une restriction significative sur l’ordre atteignable pour une méthode multipas.

Première barrière de Dahlquist L’ordre $\omega$ d’une méthode multipas à $q=p+1$ étapes consistante et zéro-stable satisfait les inégalités suivantes: - si $b_{-1}=0$ (la méthode est explicite) alors $\omega\le q=p+1$ - si $b_{-1}\neq0$ (la méthode est implicite) on doit considérér deux cas - si $q=p+1$ est impair alors $\omega\le q+1=p+2$ - si $q=p+1$ est pair alors $\omega\le q+2=p+3$

Exemples :

  • Soit une méthode explicite à un pas, alors $q=1$, $p=0$, $b_{-1}=0$ et son ordre de consistance est au plus $\omega=p+1=1$.
    (Exemple: EE qui a ordre 1)
  • Soit une méthode implicite à un pas, alors $q=1$, $p=0$, $b_{-1}\neq0$ et son ordre de consistance est au plus $\omega=p+2=2$.
    (Exemple: EI qui a ordre 1 et CN qui a ordre 2)
  • Soit une méthode explicite à deux pas, alors $q=2$, $p=1$, $b_{-1}=0$ et son ordre de consistance est au plus $\omega=p+1=2$.
    (Exemple: AB$_2$ qui a ordre 2)
  • Soit une méthode implicite à deux pas, alors $q=2$, $p=1$, $b_{-1}\neq0$ et son ordre de consistance est au plus $\omega=p+3=4$.
    (Exemple: AM$_2$ qui a ordre 3 et MS$_2$ qui a ordre 4)
In [3]:
%reset -f
from matplotlib.pylab import *

from fractions import Fraction

from IPython.display import display, Math
import sympy as sym
sym.init_printing()
r=sym.Symbol('r')


sc = { "EE_AB1" : { "aa":[1], "bb":[ 1], "bm1":0  },
       "EI_AM0_BDF1" : { "aa":[1], "bb":[ 0], "bm1":1  },
       "CN_AM1" : { "aa":[1], "bb":[Fraction(1,2)], "bm1":1/2  }, 
       "AB2": { "aa":[1,0], "bb":[Fraction(3,2),Fraction(-1,2)], "bm1":0  }, 
       "AB3": { "aa":[1,0,0], "bb":[Fraction(23,12),Fraction(-16,12), Fraction(5,12)], "bm1":0  }, 
       "AB4": { "aa":[1,0,0,0], "bb":[Fraction(55,24),Fraction(-59,24), Fraction(37,24), Fraction(-9,24)], "bm1":0  }, 
       "AB5": { "aa":[1,0,0,0,0], "bb":[Fraction(1901,720), Fraction(-2774,720), Fraction(2616,720), Fraction(-1274,720), Fraction(251,720)], "bm1":0  }, 
       "AM2": { "aa":[1,0], "bb":[Fraction(8,12),Fraction(-1,12)], "bm1":Fraction(5,12)  }, 
       "AM3": { "aa":[1,0,0], "bb":[Fraction(19,24),Fraction(-5,24),Fraction(1,24)], "bm1":Fraction(9,24)  }, 
       "AM4": { "aa":[1,0,0,0], "bb":[Fraction(646,720),Fraction(-264,720),Fraction(106,720),Fraction(-19,720)], "bm1":Fraction(251,720) }, 
       "BDF2": { "aa":[Fraction(4,3),Fraction(-1,3)], "bb":[0,0], "bm1":Fraction(2,3) }, 
       "BDF3": { "aa":[Fraction(18,11),Fraction(-9,11),Fraction(2,11)], "bb":[0,0,0], "bm1":Fraction(6,11) }, 
       "BDF4": { "aa":[Fraction(48,25),Fraction(-36,25),Fraction(16,25),Fraction(-3,25)], "bb":[0,0,0,0], "bm1":Fraction(12,25) }, 
       "BDF5": { "aa":[Fraction(300,137),Fraction(-300,137),Fraction(200,137),Fraction(-75,137),Fraction(12,137)], "bb":[0,0,0,0,0], "bm1":Fraction(60,137) }, 
       "BDF6": { "aa":[Fraction(120,49),Fraction(-150,49),Fraction(400,147),Fraction(-75,49),Fraction(24,49),Fraction(-10,147)], "bb":[0,0,0,0,0,0], "bm1":Fraction(20,49) }, 
       "MS2": { "aa":[0,1], "bb":[Fraction(4,3),Fraction(1,3)], "bm1":Fraction(1,3) }
       } 


for schema in sc.keys():
    print("\n",schema)
    q=len(sc[schema]["aa"])
    p=q-1
    print(f"\tq = {q}")
    aa=sc[schema]["aa"]
    bb=sc[schema]["bb"]
    bm1=sc[schema]["bm1"]
    
    if bm1==0:
        print("\tExplicite")
    else:
        print("\tImplicite")
    print("\tConsistance: ",end="")
    
    ordre=[]
    ordre.append( sum(aa)==1 and sum([-j*aa[j] for j in range(len(aa))])+bm1+sum(bb)==1  )
    print(ordre[-1])
    i=1
    while ordre[-1]:
        i+=1
        ordre.append( sum([(-j)**i*aa[j]+i*(-j)**(i-1)*bb[j] for j in range(len(aa))])+i*bm1==1)
    print(f"\tOrdre {len(ordre)-1}")
    
    rho=array([1])
    rho=append(rho,-array(aa))
    display(Math(r'\varrho(r)='+sym.latex(sym.expand(poly1d(rho)(r)))))
#     print(f"\tZéro-stabilité: |r_j|= {abs(roots(rho))}")
#     print("\tRacines: ",end="")
#     display(sym.solve(sym.expand(poly1d(rho)(r))))
    display( Math(r' \{|r_j|\} = '+sym.latex( [abs(rac).simplify().evalf() for rac in sym.solve(sym.expand(poly1d(rho)(r)))]  )))
    
    
        
        
 EE_AB1
	q = 1
	Explicite
	Consistance: True
	Ordre 1
$\displaystyle \varrho(r)=r - 1$
$\displaystyle \{|r_j|\} = \left[ 1.0\right]$
 EI_AM0_BDF1
	q = 1
	Implicite
	Consistance: True
	Ordre 1
$\displaystyle \varrho(r)=r - 1$
$\displaystyle \{|r_j|\} = \left[ 1.0\right]$
 CN_AM1
	q = 1
	Implicite
	Consistance: True
	Ordre 2
$\displaystyle \varrho(r)=r - 1$
$\displaystyle \{|r_j|\} = \left[ 1.0\right]$
 AB2
	q = 2
	Explicite
	Consistance: True
	Ordre 2
$\displaystyle \varrho(r)=r^{2} - r$
$\displaystyle \{|r_j|\} = \left[ 0, \ 1.0\right]$
 AB3
	q = 3
	Explicite
	Consistance: True
	Ordre 3
$\displaystyle \varrho(r)=r^{3} - r^{2}$
$\displaystyle \{|r_j|\} = \left[ 0, \ 1.0\right]$
 AB4
	q = 4
	Explicite
	Consistance: True
	Ordre 4
$\displaystyle \varrho(r)=r^{4} - r^{3}$
$\displaystyle \{|r_j|\} = \left[ 0, \ 1.0\right]$
 AB5
	q = 5
	Explicite
	Consistance: True
	Ordre 5
$\displaystyle \varrho(r)=r^{5} - r^{4}$
$\displaystyle \{|r_j|\} = \left[ 0, \ 1.0\right]$
 AM2
	q = 2
	Implicite
	Consistance: True
	Ordre 3
$\displaystyle \varrho(r)=r^{2} - r$
$\displaystyle \{|r_j|\} = \left[ 0, \ 1.0\right]$
 AM3
	q = 3
	Implicite
	Consistance: True
	Ordre 4
$\displaystyle \varrho(r)=r^{3} - r^{2}$
$\displaystyle \{|r_j|\} = \left[ 0, \ 1.0\right]$
 AM4
	q = 4
	Implicite
	Consistance: True
	Ordre 5
$\displaystyle \varrho(r)=r^{4} - r^{3}$
$\displaystyle \{|r_j|\} = \left[ 0, \ 1.0\right]$
 BDF2
	q = 2
	Implicite
	Consistance: True
	Ordre 2
$\displaystyle \varrho(r)=r^{2} - \frac{4 r}{3} + \frac{1}{3}$
$\displaystyle \{|r_j|\} = \left[ 0.333333333333333, \ 1.0\right]$
 BDF3
	q = 3
	Implicite
	Consistance: True
	Ordre 3
$\displaystyle \varrho(r)=r^{3} - \frac{18 r^{2}}{11} + \frac{9 r}{11} - \frac{2}{11}$
$\displaystyle \{|r_j|\} = \left[ 1.0, \ 0.426401432711221, \ 0.426401432711221\right]$
 BDF4
	q = 4
	Implicite
	Consistance: True
	Ordre 4
$\displaystyle \varrho(r)=r^{4} - \frac{48 r^{3}}{25} + \frac{36 r^{2}}{25} - \frac{16 r}{25} + \frac{3}{25}$
$\displaystyle \{|r_j|\} = \left[ 1.0, \ 0.56086151609339, \ 0.56086151609339, \ 0.38147840911841\right]$
 BDF5
	q = 5
	Implicite
	Consistance: True
	Ordre 5
$\displaystyle \varrho(r)=r^{5} - \frac{300 r^{4}}{137} + \frac{300 r^{3}}{137} - \frac{200 r^{2}}{137} + \frac{75 r}{137} - \frac{12}{137}$
$\displaystyle \{|r_j|\} = \left[ 1.0, \ 0.708710816266406, \ 0.417600758175733, \ 0.708710816266406, \ 0.417600758175733\right]$
 BDF6
	q = 6
	Implicite
	Consistance: True
	Ordre 6
$\displaystyle \varrho(r)=r^{6} - \frac{120 r^{5}}{49} + \frac{150 r^{4}}{49} - \frac{400 r^{3}}{147} + \frac{75 r^{2}}{49} - \frac{24 r}{49} + \frac{10}{147}$
$\displaystyle \{|r_j|\} = \left[ 1.0, \ 0.40612326685391, \ 0.474034857706592, \ 0.474034857706592, \ 0.863380267869827, \ 0.863380267869827\right]$
 MS2
	q = 2
	Implicite
	Consistance: True
	Ordre 4
$\displaystyle \varrho(r)=r^{2} - 1$
$\displaystyle \{|r_j|\} = \left[ 1.0, \ 1.0\right]$

A-stabilité

Dans la section précédente, on a considéré la résolution du problème de Cauchy sur des intervalles bornés. Il existe cependant de nombreuses situations dans lesquelles le problème de Cauchy doit être intégré sur des intervalles en temps très grands ou même infini. Dans ce cas, même pour $h$ fixé, $N$ tend vers l’infini, et un résultat comme la zéro-stabilité n’a plus de sens puisque le membre de droite contient une quantité non bornée ($T-t_0$).
On s’intéresse donc à des méthodes capables d’approcher la solution pour des intervalles en temps arbitrairement grands, même pour des pas de temps $h$ "assez grands".

Soit $\lambda\in\mathbb{C}$ un nombre complexe tel que $\Re(\lambda)=-\beta<0$ et considérons le problème de Cauchy $$\begin{cases} y'(t)=\lambda y(t), &\text{pour }t>0,\\ y(0)=1 \end{cases}$$ Sa solution est $y(t)=e^{\lambda t}$ et comme $\Re(\lambda)<0$ alors $$\lim_{t\to+\infty}|y(t)|=0.$$

Soit $h>0$ un pas de temps donné, $t_n=nh$ pour $n\in\mathbb{N}$ et notons $u_n\approx y(t_n)$ une approximation de la solution $y$ au temps $t_n$.

Si, sous d'éventuelles conditions sur $h$, on a $$ \lim_{n\to+\infty} |u_n| =0, $$ alors on dit que le schéma est A-stable (absolument stable).

Région de stabilité absolue

On appelle région de stabilité absolue $A$ d’une méthode numérique l’ensemble des nombres complexes $z = h\lambda$ pour lesquels la méthode est absolument stable: $$ A=\Big\{z = h\lambda \in \mathbb{C} \ \Big| \ |u_n|\xrightarrow[n\to+\infty]{}0 \Big\} $$ Une méthode numérique est inconditionellement A-stable si $\mathbb{C}^-\subset A$ où $\mathbb{C}^-=\{z\in\mathbb{C}|\Re(z)<0\}$.

Une méthode à $q=p+1$ pas appliquée à ce problème modèle conduit à léquation récurrente à $$ u_{n+1}=\sum_{j=0}^p a_j u_{n-j} + h\lambda\sum_{j=0}^p b_ju_{n-j} + h\lambda b_{-1}u_{n+1} $$ dont le polynôme caracthéristique associé $$ p(r) = \varrho(r) - (h\lambda) \sigma(r) $$ est appelé polynôme de A-stabilité associé à la méthode numérique.

Condition des racines stricte: une méthode multipas linéaire vérifie la condition des racines stricte (éventuellement sous conditions sur $h$) si les racines du polynôme $p$ sont dans le disque unité ouvert de $\mathbb{C}$: $$ \text{il existe $h_0>0$ tel que }\quad |r_j(h\lambda)|<1 \quad\text{pour tout }j=0,\dots,p;\ \forall h\le h_0 $$
Proposition. Une méthode multipas linéaire est absolument stable si et seulement si son polynôme de stabilité vérifie la condition des racines stricte.

Régions de A-stabilité de diverses méthodes d'Adams-Bashforth (à gauche) et d'Adams-Moulton (à droite).
Source: A. Quarteroni, F. Saleri, P. Gervasio Calcul Scientifique: Cours, exercices corrigés et illustrations en Matlab et Octave.

On remarque que les régions de stabilité absolue des méthodes d'Adams-Bashforth et d'Adams-Moulton sont bornées et diminuent quand l'ordre augmente.

A_stab_Adams.png

Régions de A-stabilité de diverses méthodes BDF.
Source: A. Quarteroni, F. Saleri, P. Gervasio Calcul Scientifique: Cours, exercices corrigés et illustrations en Matlab et Octave.

On remarque que les régions de stabilité absolue des méthodes BDF, pour $q\le6$, sont non bornées. En particulier, elles contiennent toujours les réels négatifs.

A_stab_BDF.png

Tracer les régions de A-stabilité.
Posons $z=h\lambda$. Il est en général difficile de déterminer la région de stabilité absolue car on doit décider pour quels $z\in\mathbb{C}$ les racines du polynôme de stabilité vérifient la condition de racines stricte (i.e. $|r| < 1$). Il est plus aisé de déterminer la frontière de cette région car au moins une des racines de $p$ à la frontière est de module $1$. On a donc $|r| = 1$ sur $\partial A$. La frontière de $A$ est alors un sous ensemble des points $z\in\mathbb{C}$ pour lesquels $r=e^{i\vartheta}$, $\vartheta\in\mathbb{R}$. On remplace alors $r =e^{i\vartheta}$ dans $p$ et comme $e^{i\vartheta}$ est racine, $p(e^{i\vartheta}) = 0$ et on résoud l’équation. On obtient $z = z(\vartheta)$ ce qui donne une courbe dans $\mathbb{C}$ en fonction de l’angle $\vartheta$, ce qui décrit donc une courbe polaire. Par conséquent, pour obtenir une représentation approchée de $\partial A$, il suffit d’évaluer le second membre pour diverses valeurs de $r$ sur le cercle unité, diverses valeurs de $\vartheta\in\mathbb{R}$ et $r=e^{i\vartheta}$.

In [33]:
%reset -f
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np

import sympy as sp
sp.init_printing()
def polyCar(schema):
    sp.var('r z')
    p,aa,bb,bm1=schema["p"],schema["aa"],schema["bb"],schema["bm1"]    
    rho=r**(p+1)-sum(aa[j]*r**(p-j)for j in range(p+1))    
    sigma=bm1*r**(p+1)+sum(bb[j]*r**(p-j)for j in range(p+1))
    display(f'rho(r)={rho}')
    display(f'sigma(r)={sigma}')
    return 
def polyCarIMPL(schema):
    sp.var('r z')
    p,aa,bb,bm1=schema["p"],schema["aa"],schema["bb"],schema["bm1"]    
    rho=r**(p+1)-sum(aa[j]*r**(p-j)for j in range(p))    
    sigma=bm1*r**(p+1)+sum(bb[j]*r**(p-j)for j in range(p))
    display(f'rho(r)={rho}')
    display(f'sigma(r)={sigma}')
    return 
    

# Returns > 1 if argument is not in region of absolute stability
def stabilityFunction(z,schema):
    pol=np.array(schema["aa"])+z*np.array(schema["bb"])
    pol=np.append([1-z*schema["bm1"]],-pol)
    return max(abs(np.roots(pol)))


def Astabilite(schema):
    p,aa,bb,bm1=schema["p"],schema["aa"],schema["bb"],schema["bm1"]
    
    # cc liste de points sur le circle unitaire
    rr = np.linspace(0,2*np.pi,201)
    cc = np.cos(rr)+1j*np.sin(rr) 
    ## ff fonction vectorisée qui donne zz
    ff = lambda cc,p,aa,bb,bm1 : (cc**(p+1)-np.polyval(aa,cc))/(bm1*cc**(p+1)+np.polyval(bb,cc)) 
    zz = ff(cc,p,aa,bb,bm1)
    
    plt.figure(figsize=(18,7))
    plt.subplot(1,2,1)
    plt.plot(zz.real,zz.imag,'r')
    plt.grid(True)
#     plt.axis('equal')
    plt.xlabel(r'$\Re(h\beta)$')
    plt.ylabel(r'$\Im(h\beta)$');
    plt.title(schema["Nom"])
    Lx,Ly = plt.xlim(),plt.ylim()
    
    plt.subplot(1,2,2)
    x = np.linspace(Lx[0],Lx[1], 300)
    y = np.linspace(Ly[0],Ly[1], 300)
    [X,Y] = np.meshgrid(x,y)
    Z = np.zeros(X.shape)
    for m in range(X.shape[0]):
        for n in range(X.shape[1]):
            Z[m,n] = stabilityFunction(X[m,n] + 1j * Y[m,n], schema)
    plt.contour(X, Y, Z, [1], colors='k')
    plt.contourf(X, Y, Z, [0,1], colors=[[1, 0.5, 0.8]])
    return None
In [34]:
EE = { "Nom":"EE=AB1", "p":0, "aa":[1], "bb":[1], "bm1":0  }
polyCar(EE)
Astabilite(EE) 
'rho(r)=r - 1'
'sigma(r)=1'
In [35]:
AB2 = { "Nom":"AB2", "p":1, "aa":[1,0], "bb":[3/2,-1/2], "bm1":0  }
polyCar(AB2)
Astabilite(AB2)
'rho(r)=r**2 - r'
'sigma(r)=1.5*r - 0.5'
In [36]:
AB3 = { "Nom":"AB3", "p":2, "aa":[1,0,0], "bb":[23/12,-16/12, 5/12], "bm1":0  }
polyCar(AB3)
Astabilite(AB3)
'rho(r)=r**3 - r**2'
'sigma(r)=1.91666666666667*r**2 - 1.33333333333333*r + 0.416666666666667'
In [37]:
AB4 = { "Nom":"AB4", "p":3, "aa":[1,0,0,0], "bb":[55/24,-59/24, 37/24, -9/24], "bm1":0  }
polyCar(AB4)
Astabilite(AB4)
'rho(r)=r**4 - r**3'
'sigma(r)=2.29166666666667*r**3 - 2.45833333333333*r**2 + 1.54166666666667*r - 0.375'
In [38]:
AB5 = { "Nom":"AB5", "p":4, "aa":[1,0,0,0,0], "bb":[1901/720, -2774/720, 2616/720, -1274/720, 251/720], "bm1":0  }
polyCar(AB5)
Astabilite(AB5)
'rho(r)=r**5 - r**4'
'sigma(r)=2.64027777777778*r**4 - 3.85277777777778*r**3 + 3.63333333333333*r**2 - 1.76944444444444*r + 0.348611111111111'
In [39]:
AM1 = { "Nom":"AM1=CN", "p":1, "aa":[1], "bb":[1/2], "bm1":1/2  }
polyCarIMPL(AM1)
#Astabilite(AM1)
'rho(r)=r**2 - r'
'sigma(r)=0.5*r**2 + 0.5*r'
In [40]:
AM2 = { "Nom":"AM2", "p":1, "aa":[1,0], "bb":[8/12,-1/12], "bm1":5/12  }
polyCar(AM2)
Astabilite(AM2)
'rho(r)=r**2 - r'
'sigma(r)=0.416666666666667*r**2 + 0.666666666666667*r - 0.0833333333333333'
In [41]:
AM3 = { "Nom":"AM3", "p":2, "aa":[1,0,0], "bb":[19/24,-5/24,1/24], "bm1":9/24  }
polyCar(AM3)
Astabilite(AM3)
'rho(r)=r**3 - r**2'
'sigma(r)=0.375*r**3 + 0.791666666666667*r**2 - 0.208333333333333*r + 0.0416666666666667'
In [42]:
AM4 = { "Nom":"AM4", "p":3, "aa":[1,0,0,0], "bb":[646/720,-264/720,106/720,-19/720], "bm1":251/720  }
polyCar(AM4)
Astabilite(AM4)
'rho(r)=r**4 - r**3'
'sigma(r)=0.348611111111111*r**4 + 0.897222222222222*r**3 - 0.366666666666667*r**2 + 0.147222222222222*r - 0.0263888888888889'
In [43]:
MS2 = { "Nom":"MS2", "p":1, "aa":[0,1], "bb":[4/3,1/3], "bm1":1/3  }
polyCar(MS2)
#Astabilite(MS2)
'rho(r)=r**2 - 1'
'sigma(r)=0.333333333333333*r**2 + 1.33333333333333*r + 0.333333333333333'
In [44]:
BDF1 = { "Nom":"BDF1=EI", "p":1, "aa":[1], "bb":[0], "bm1":1  }
polyCarIMPL(BDF1)
Astabilite(BDF1)
'rho(r)=r**2 - r'
'sigma(r)=r**2'
In [45]:
BDF2 = { "Nom":"BDF2", "p":1, "aa":[4/3,-1/3], "bb":[0], "bm1":2/3  }
polyCarIMPL(BDF2)
Astabilite(BDF2)
'rho(r)=r**2 - 1.33333333333333*r'
'sigma(r)=0.666666666666667*r**2'
In [46]:
BDF3 = { "Nom":"BDF3", "p":2, "aa":[18/11,-9/11,2/11], "bb":[0], "bm1":6/11  }
# polyCarIMPL(BDF3)
Astabilite(BDF3)
In [47]:
BDF4 = { "Nom":"BDF4", "p":3, "aa":[48/25,-36/25,16/25,-3/25], "bb":[0], "bm1":12/25  }
# polyCarIMPL(BDF4)
Astabilite(BDF4)