Vérification de CFL de Maxwell
===

Il s'agit ici de vérifier un calcul de CFL lié à Maxwell dans RK(3,3) suite à un mail de Nicolas du 4 novembre 2020 à 16h21.


```python
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
```


```python
k = sp.symbols("k",real=True)
dt = sp.symbols(r"\Delta\ t",real=True,positive=True)
```

On cherche à résoudre uniquement la partie Maxwell :
$$
  \partial_t U = A\cdot U
$$
avec $U$ le vecteur de variable (uniquement $B$ et $E$, seules variables intervenant dans Maxwell) et la matrice $A$ suivante :  


```python
A = sp.Matrix([
        [ 0      ,  0      ,  0     , sp.I*k ],
        [ 0      ,  0      , -sp.I*k, 0      ],
        [ 0      , -sp.I*k ,  0     , 0      ],
        [ sp.I*k ,  0      ,  0     , 0      ],
    ])
U = sp.MatrixSymbol("U^n",4,1)

A
```

On souhaite résoudre ce problème avec une méthode RK(3,3) ou RK(4,4) pour déterminer la CFL dans ce problème pour tout mode de Fourier $k$.


```python
def rk33(un,L,dt):
  """
    méthode RK(3,3) de Shu-Osher
  """
  u1 = un + dt*L*un
  u2 = sp.Rational(3,4)*un + sp.Rational(1,4)*u1 + sp.Rational(1,4)*dt*L*u1
  return sp.Rational(1,3)*un + sp.Rational(2,3)*u2 + sp.Rational(2,3)*dt*L*u2
```


```python
def rk44(un,L,dt):
  """
    méthode RK(4,4)
  """
  u1 = un + sp.Rational(1,2)*dt*L*un
  u2 = un + sp.Rational(1,2)*dt*L*u1
  u3 = un + dt*L*u2
  return sp.Rational(1,3)*( -un + u1 + 2*u2 + u3 ) + sp.Rational(1,6)*dt*L*u3
```

On vérifie nos méthodes de Runge-Kutta, pour cela on vérifie que l'on obtient bien sur une itération le développement de Taylor évalué en $\Delta t \lambda$ ($\lambda$ étant une linéarisation de notre fonction dans notre problème).


```python
# vérification de notre méthode RK
L = sp.symbols("L")
un = sp.symbols("u_n")

rk33(un,L,dt).expand().collect(un).subs(L,sp.symbols(r"\lambda"))
```


```python
rk44(un,L,dt).expand().collect(un).subs(L,sp.symbols(r"\lambda"))
```

Maintenant on test avec la matrice $A$ à la place du symbole $L$, et le vecteur $U$ en lieu et place de $u_n$.


```python
sp.Matrix(rk33(U,A,dt).expand().simplify())
```


```python
sp.Matrix(rk44(U,A,dt).expand().simplify())
```

Il n'est pas possible avec `sympy` d'écrire cela sous la forme $B_{A,\Delta t}\cdot U$. Donc pour RK(3,3) on trouve à la main la matrice :

$$
  \begin{pmatrix}
      1 - \frac{\Delta t^2k^2}{2} & 0 & 0 & i\Delta tk - \frac{i\Delta t^3k^3}{6} \\
      0 & 1 - \frac{\Delta t^2k^2}{2} & -i\Delta tk + \frac{i\Delta t^3k^3}{6} & 0 \\
      0 & -i\Delta tk + \frac{i\Delta t^3k^3}{6} & 1 - \frac{\Delta t^2k^2}{2} & 0 \\
      i\Delta tk - \frac{i\Delta t^3k^3}{6} & 0 & 0 & 1 - \frac{\Delta t^2k^2}{2} \\
  \end{pmatrix}
$$

Mais on sait que l'on doit obtenir le développement de Taylor de la fonction exponentielle :


```python
l = sp.symbols(r"\lambda")
sp.exp(dt*l).series(dt*l,n=5)
```


```python
B_rk33 = sp.eye(4) + dt*A + sp.Rational(1,2)*(A*dt)**2 + sp.Rational(1,6)*(dt*A)**3
B_rk33
```


```python
B_rk44 = sp.eye(4) + dt*A + sp.Rational(1,2)*(A*dt)**2 + sp.Rational(1,6)*(dt*A)**3 + sp.Rational(1,24)*(dt*A)**4
B_rk44
```


```python
# on vérifie bien le résultat
sp.Matrix(B_rk33*U).expand()
```


```python
sp.Matrix(B_rk44*U).expand()
```

Maintenant calculons les valeurs propres de nos matrices $B_{rk(n,n)}$


```python
ev_rk33 = B_rk33.eigenvals()

for v,m in ev_rk33.items():
    display(v.collect(sp.I))
    print("multiplicity of ",m)
```


```python
ev_rk44 = B_rk44.eigenvals()

for v,m in ev_rk44.items():
    display(v.collect(sp.I))
    print("multiplicity of ",m)
```

Maintenant on trace le module de ces valeurs propres pour différentes valeurs de $\Delta tk$.


```python
nK = np.linspace(-2,2,100)

for i,v in enumerate(ev_rk33):
    evk = sp.lambdify(k,v.subs(dt,1))(nK)
    kmin = np.amin(nK[np.where(np.abs(evk)<1.0)])
    kmax = np.amax(nK[np.where(np.abs(evk)<1.0)])
    print(kmin,kmax)
    plt.plot(nK,np.abs(evk),linewidth=2-i)
    plt.vlines(kmin,np.amin(np.abs(evk))*0.99,np.max(np.abs(evk))*1.01,colors="gray",linestyles="dotted")
    plt.vlines(kmax,np.amin(np.abs(evk))*0.99,np.max(np.abs(evk))*1.01,colors="gray",linestyles="dotted")
plt.hlines(1.0,-2,2,colors="gray",linestyles="dotted")
```

J'en conclus qu'avec RK(3,3), il faut que $\Delta tk \leq 1.71$ pour assurer la stabilité. Sachant que $k$ sont les fréquences de Fourier, donc $k = \frac{2\pi n}{L}$ avec $L = \pi$ et $n\in[\![-\frac{N_z}{2},\frac{N_z}{2}]\!]$. On a donc $k_\text{max} = N_z$, il faut donc $\Delta t \leq \frac{1.71}{N_z}$, avec typiquement $N_z = 15$ pour les simus grossières : $\Delta t\leq0.114$ ou $N_z = 27$ pour les simus plus fines : $\Delta t \leq 0.063$.


```python
nK = np.linspace(-3,3,100)

for i,v in enumerate(ev_rk44):
    evk = sp.lambdify(k,v.subs(dt,1))(nK)
    kmin = np.amin(nK[np.where(np.abs(evk)<1.0)])
    kmax = np.amax(nK[np.where(np.abs(evk)<1.0)])
    print(kmin,kmax)
    plt.plot(nK,np.abs(evk),linewidth=2-i)
    plt.vlines(kmin,np.amin(np.abs(evk))*0.99,np.max(np.abs(evk))*1.01,colors="gray",linestyles="dotted")
    plt.vlines(kmax,np.amin(np.abs(evk))*0.99,np.max(np.abs(evk))*1.01,colors="gray",linestyles="dotted")
plt.hlines(1.0,-3,3,colors="gray",linestyles="dotted")
```

J'en conclus qu'avec RK(4,4), il faut que $\Delta tk \leq 2.81$ pour assurer la stabilité. Sachant que $k$ sont les fréquences de Fourier, donc $k = \frac{2\pi n}{L}$ avec $L = \pi$ et $n\in[\![-\frac{N_z}{2},\frac{N_z}{2}]\!]$. On a donc $k_\text{max} = N_z$, il faut donc $\Delta t \leq \frac{2.81}{N_z}$, avec typiquement $N_z = 15$ pour les simus grossières : $\Delta t\leq0.187$ ou $N_z = 27$ pour les simus plus fines : $\Delta t \leq 0.104$.

Maintenant étudions le résultat en fixant la valeur de $N_z$ donc de $k$ à 27 (pour se rapprocher de l'approche de Nicolas).


```python
B27 = B_rk33.subs(k,27)
B27
```

On cherche maintenant à trouver le $\Delta t$ tel que la valeur propre vaut 1.


```python
for v in B27.eigenvals():
  display(sp.solve(sp.Eq(sp.Abs(v)**2,1),dt))
```


```python
np.sqrt(3.)/27.
```

Résultat très similaire à l'étude numérique.

Pour $N_z=15$


```python
B15 = B_rk33.subs(k,15)
B15
```


```python
for v in B15.eigenvals():
  display(sp.solve(sp.Eq(sp.Abs(v)**2,1),dt))
```


```python
np.sqrt(3.)/15.
```

On a visiblement quelque chose en $\frac{\sqrt{3}}{N_z}$, on vérifie ça directement sur $B_{RK(3,3)}$.


```python
for v in B_rk33.eigenvals():
  display(sp.solve(sp.Eq(sp.Abs(v)**2,1),dt))
```

Maintenant avec RK(4,4)


```python
B27 = B_rk44.subs(k,27)
B27
```


```python
for v in B27.eigenvals():
  x = sp.solve(sp.Eq(sp.Abs(v)**2,1),dt)[0]
  display(x,x.evalf())
```

Maintenant avec $N_z=15$ :


```python
B15 = B_rk44.subs(k,15)

for v in B15.eigenvals():
  x = sp.solve(sp.Eq(sp.Abs(v)**2,1),dt)[0]
  display(x,x.evalf())
```


```python
for v in B_rk44.eigenvals():
  display(sp.solve(sp.Eq(sp.Abs(v)**2,1),dt))
```

Donc pour RK(3,3) on a une CFL de Maxwell en $\frac{\sqrt{3}}{N_z}$ et pour RK(4,4) une CFL en $\frac{2\sqrt{2}}{N_z}$.

En fait il s'agit de chercher la CFL d'un transport pur non diffusif, donc trouver ce qu'on avait nommé $y_\text{max}$ pour la stabilité de CD2. Donc les valeurs ont déjà été calculé sur [Ponio](http://jmassot.perso.math.cnrs.fr/ponio.html), j'ai un ami webdev qui s'est aussi amusé à faire une visualisation, un peu plus optimisé [ici](https://runge-kutta-method-viewer.netlify.app/runge-kutta-methods/compare?left=10&right=15).


```python

```

209b0d701f614eceaa94f0507bcb4effca16a857
ipynb
Jupyter Notebook
notes/notebooks/Maxwell_cfl.ipynb
kivvix/draft
33b605be27e556df061f856be8e84e5b3f49a219
MIT
notes/notebooks/Maxwell_cfl.ipynb
kivvix/draft
33b605be27e556df061f856be8e84e5b3f49a219
MIT
notes/notebooks/Maxwell_cfl.ipynb
kivvix/draft
33b605be27e556df061f856be8e84e5b3f49a219
MIT
Qwen/Qwen-72B
1. YES
2. YES